EDIT: SOLVED; I defined my trial and test functions incorrectly. Please delete this thread if possible.
Summary: I am solving for two variables in coupled PDEs simultaneously, but one of them is zero everywhere and I do not understand why. Can you see why?
Hello,
I am very new to Netgen/NGSolve, and have been playing with documentation and tutorials for the last week.
Now, I am trying to simultaneously solve a system of coupled time-harmonic Maxwell’s equations in a cube, for both the electric field E and magnetic field H. This is a very simple implementation, I feel - but my lack of experience with NGSolve is at fault. Here is the system of PDEs:
j w (eps E) - curl(H) = - J_E
j w (mu H) + curl(E) = 0
where j is the imaginary unit, w is frequency, eps/mu are material coefficients and J_E is the source term. j/w were introduced to replace time derivatives, and the variational form used in my code (attached) was the one following integration by parts (from which the resulting boundary terms of the curl terms was set to 0).
In order to test my NGSolve implementation, I manufactured a solution for the electric field E and magnetic field H and used it to compare the numerical solution.
For sanity, I started by solving for the source-free case J_E = 0. However, this resulted in there being a zero solution for H. This does not match with how the exact solution looks like (both when values are output using print(), and on Netgen) when interpolated over the mesh.
The same issue was seen when I analytically solved for J_E using the differential form above and included it as a source term in my variational formulation. Still, I end up with a zero solution for the numerical result of H. I have tried debugging this issue for a few days now but cannot see what’s wrong. So, could someone help me out by pointing out my mistake?
I would also appreciate any input in how I am implementing this problem because this is the first one I am trying on my own.
I have attached a minimal implementation. Note that the issue remains when I use CGSolver with various preconditioners. So I am suspecting that my issue is with how I setup my problem.
Thanks in advance!
EDIT: I just saw the code and noticed an error in the last line that prints output, it should be:
print(u.components[0].vec)
print(u.components[1].vec) # This outputs (0, 0) everywhere.
instead of
print(uu.components[0].vec)
print(uu.components[1].vec) # This outputs (0, 0) everywhere.