Two issues from trying an electromagnetism problem (complex numbers and vector curl)
Hello,
I'm been aware of SfePy for a very long time, but only recently have actually been trying it. I've been surprised at how *simple* it is to get things going, it's a great package. Thanks a lot for making it.
What I'm trying to do is solve a scattering problem from electromagnetism. I'm working with a 3D vector electric field (of complex numbers). So far as I can tell, I've been able to get the weak form of the problem, and it's boundary conditions, implemented correctly.
There are two things I'd like that are not working at the moment. I'd could really use some guidance, so would appreciate any help. These are:
I'd like to evaluate a surface integral of the electric field intensity I = |E|^2 = E dot conj(E), where E is a vector, dot is the dot product, and conj is the complex conjugate. I'm implementing my own term for this, which I think should be nearly the same as the existing IntegrateSurfaceTerm. However, it appears that the complex numbers are complicating things. The function method of my new term is apparently called twice in 'eval' mode, once for the real part and once for the imag part. Looking at the source code of Term, this seems to be the implemented behavior. Is it possible to get the function method called with the complex numbers themselves? If not, what is the recommended technique for doing operations in which the real parts and the imag parts do not separate cleanly (like for the magnitude of a complex vector)? If it helps, I only need things like these in eval mode, not as part of the solved equations.
I'd like to evaluate the magnetic field H from the electric field E. This is simply H = -1j * c * curl(E), where c is a constant (material parameter) and curl is the vector curl operation. The vector field E is solved for as part of my problem. What is the best way to add the equation for H into the problem? I don't see a curl operation anywhere. Do I need a custom term or can I do this with existing terms? It is fully determined by E, as above.
Thanks for your help,
Irwin
Hello Irwin,
On 08/27/2018 07:17 PM, Irwin Zaid wrote:
Hello,
I'm been aware of SfePy for a very long time, but only recently have actually been trying it. I've been surprised at how *simple* it is to get things going, it's a great package. Thanks a lot for making it.
Thanks for the kind words!
What I'm trying to do is solve a scattering problem from electromagnetism. I'm working with a 3D vector electric field (of complex numbers). So far as I can tell, I've been able to get the weak form of the problem, and it's boundary conditions, implemented correctly.
Before answering your question, let me point out this. When doing electromagnetic simulations, AFAIK, vector finite elements from the H(curl) space are needed (e.g. Nedelec elements). I have no experience with such elements, and those are not implemented in SfePy. So unless you are sure that H1 elements are OK for your problem, SfePy might not be the right tool to use.
There are two things I'd like that are not working at the moment. I'd could really use some guidance, so would appreciate any help. These are:
- I'd like to evaluate a surface integral of the electric field intensity I = |E|^2 = E dot conj(E), where E is a vector, dot is the dot product, and conj is the complex conjugate. I'm implementing my own term for this, which I think should be nearly the same as the existing IntegrateSurfaceTerm. However, it appears that the complex numbers are complicating things. The function method of my new term is apparently called twice in 'eval' mode, once for the real part and once for the imag part. Looking at the source code of Term, this seems to be the implemented behavior. Is it possible to get the function method called with the complex numbers themselves? If not, what is the recommended technique for doing operations in which the real parts and the imag parts do not separate cleanly (like for the magnitude of a complex vector)? If it helps, I only need things like these in eval mode, not as part of the solved equations.
It is not currently possible, but it could be rather easily implemented, assuming your term evaluation function would be pure Python (or Cython). The reason of the current implementation is more or less historical - a number of term evaluation functions are implemented in C, so they expect just simple float64 arrays.
- I'd like to evaluate the magnetic field H from the electric field E. This is simply H = -1j * c * curl(E), where c is a constant (material parameter) and curl is the vector curl operation. The vector field E is solved for as part of my problem. What is the best way to add the equation for H into the problem? I don't see a curl operation anywhere. Do I need a custom term or can I do this with existing terms? It is fully determined by E, as above.
If this is meant just for post-processing, then the curl operator could be added. However, for solving the full coupled Maxwell equations, H(curl) elements would be required, as mentioned above. This is currently out of scope for me to add.
Best regards, r.
Before answering your question, let me point out this. When doing electromagnetic simulations, AFAIK, vector finite elements from the H(curl) space are needed (e.g. Nedelec elements). I have no experience with such elements, and those are not implemented in SfePy. So unless you are sure that H1 elements are OK for your problem, SfePy might not be the right tool to use.
This is very clear, thanks. To be clear on my end, because this is a scattering problem, this actually reduces to solving the vector Helmholtz equation for a single field, either E or H. It's not the full Maxwell's equations, this is a standard approach for light scattering. The other field is easily determined from the first, as I mentioned above for H.
My understanding right now is that I can solve the vector Helmholtz equation within SfePy, and I believe I've implemented it correctly, but could be wrong.
It is not currently possible, but it could be rather easily implemented, assuming your term evaluation function would be pure Python (or Cython). The reason of the current implementation is more or less historical - a number of term evaluation functions are implemented in C, so they expect just simple float64 arrays.
I'm happy to do pure Python or Cython. How would you suggest I proceed to get the complex numbers? Is there any way I can make this work with SfePy as it is now? I'd be willing to do PRs for SfePy in the near future if I am sure I can solve my problem within it -- which is why I'd like to work with it as it is now.
If this is meant just for post-processing, then the curl operator could be added. However, for solving the full coupled Maxwell equations, H(curl) elements would be required, as mentioned above. This is currently out of scope for me to add.
Yes, it is meant for just post-processing. Basically I solve for E, and get from H = const * curl(E). I'd also be interested in using the curl operator on parameter fields in the problem, but definitely not for test or trial fields. Again, is there a way to do this now with custom terms, or does it require changes to the SfePy source?
I think we're on the same page -- I *believe* I can use SfePy for this because I'm solving the equations for light scattering (vector Helmholtz equation), rather than the full Maxwell equations. But please correct me if I'm wrong. :)
Irwin
On 08/28/2018 12:43 AM, Irwin Zaid wrote:
Before answering your question, let me point out this. When doing electromagnetic simulations, AFAIK, vector finite elements from the H(curl) space are needed (e.g. Nedelec elements). I have no experience with such elements, and those are not implemented in SfePy. So unless you are sure that H1 elements are OK for your problem, SfePy might not be the right tool to use.
This is very clear, thanks. To be clear on my end, because this is a scattering problem, this actually reduces to solving the vector Helmholtz equation for a single field, either E or H. It's not the full Maxwell's equations, this is a standard approach for light scattering. The other field is easily determined from the first, as I mentioned above for H.
My understanding right now is that I can solve the vector Helmholtz equation within SfePy, and I believe I've implemented it correctly, but could be wrong.
Yes, the vector Helmholtz equation should be possible.
It is not currently possible, but it could be rather easily implemented, assuming your term evaluation function would be pure Python (or Cython). The reason of the current implementation is more or less historical - a number of term evaluation functions are implemented in C, so they expect just simple float64 arrays.
I'm happy to do pure Python or Cython. How would you suggest I proceed to get the complex numbers? Is there any way I can make this work with SfePy as it is now? I'd be willing to do PRs for SfePy in the near future if I am sure I can solve my problem within it -- which is why I'd like to work with it as it is now.
You would need to modify eval_complex() in sfepy/terms/terms.py so that it does not split the term arguments into the real and imaginary parts. For this, try using the term_mode argument. The call would then be something like:
val = problem.evaluate('your_term.i.Omega(E)', mode='eval', term_mode='complex')
and in eval_complex() you just check for term_mode == 'complex' and in that case just call
status = self.call_function(out, fargs)
I'll think of a cleaner way in the meantime (some other small changes in the term evaluation code might be needed - maybe adding 'ceval' to possible values of mode argument.
If this is meant just for post-processing, then the curl operator could be added. However, for solving the full coupled Maxwell equations, H(curl) elements would be required, as mentioned above. This is currently out of scope for me to add.
Yes, it is meant for just post-processing. Basically I solve for E, and get from H = const * curl(E). I'd also be interested in using the curl operator on parameter fields in the problem, but definitely not for test or trial fields. Again, is there a way to do this now with custom terms, or does it require changes to the SfePy source?
A new term similar to ev_grad would do it for the moment. It might be also a good idea to add curl as a possible mode to FieldVariable.evaluate().
I think we're on the same page -- I *believe* I can use SfePy for this because I'm solving the equations for light scattering (vector Helmholtz equation), rather than the full Maxwell equations. But please correct me if I'm wrong. :)
Yes, the Helmholtz equation should work :)
r.
On 08/28/2018 08:24 AM, Robert Cimrman wrote:
On 08/28/2018 12:43 AM, Irwin Zaid wrote:
I'm happy to do pure Python or Cython. How would you suggest I proceed to get the complex numbers? Is there any way I can make this work with SfePy as it is now? I'd be willing to do PRs for SfePy in the near future if I am sure I can solve my problem within it -- which is why I'd like to work with it as it is now.
You would need to modify eval_complex() in sfepy/terms/terms.py so that it does not split the term arguments into the real and imaginary parts. For this, try using the term_mode argument. The call would then be something like:
val = problem.evaluate('your_term.i.Omega(E)', mode='eval', term_mode='complex')
and in eval_complex() you just check for term_mode == 'complex' and in that case just call
status = self.call_function(out, fargs)
I'll think of a cleaner way in the meantime (some other small changes in the term evaluation code might be needed - maybe adding 'ceval' to possible values of mode argument.
Well, I totally forgot to mention the most simple solution: you can define your own .eval_complex() in your term subclass that would override the standard one.
r.
participants (2)
-
Irwin Zaid -
Robert Cimrman