Post processing complex solutions


I’m trying to solve a simple curl curl problem with complex data NGSolve. But I’m getting some strange results. I’ve set up a simple test problem, to see what is going on. I supply a Dirichlet boundary condition which is also the solution of curl curl A=0 in the domain and A is divergence free. If I set the problem up such that the exact solution is A= e_1 \times x (e_1 unit vector in the x_1 direction, x coordinate vector) then all is fine. But if I set up a problem with complex solution e.g. A=(1+2j) e_1 \times x then I get strange results. I think the issues are in the post processing of the solution that I’m doing.

To post process I compute integrals of \int _sphereA * A d omega over a sphere of radius 1. I know the exact solution is

8\pi /15 if A= e_1 \times x


(1+2i)^2 8 pi/15= (-3 +4i )8pi/15 if A= (1+2i) e_1 \times x

I get the right result if I integrate the coefficient function against itself, but only get the correct solution for

Integrate(inout * InnerProduct(Theta,Theta),mesh)

If the solution is real. I know that I should be able to solve this exactly using the lowest order elements as the function is linear.

In the above inout is a flag which is 1 in the sphere and 0 outside, the domain is box with a sphere and the centre and I’m solving curl curl A =0 with regularisation in Omega with Dirichlet BCs on the outside.

Could you please advise me as to whether I am using this incorrectly or whether there may potentially be a bug with the integrate function. I attach a simple code to illustrate this.

Many thanks,



Hi Ben,

use the CG-Solver instead of GMRes, then it works

inverse= CGSolver(a.mat, c.mat, precision=1e-6, maxsteps=200)

CG can solve with complex-symmetric matrices.

Exact= (-5.0265482451562615+6.702064326875025j)
Calculated= (-5.026398383236931+6.701864510664208j)
HalfandHalf= (-5.0264730864946525+6.701964115326205j)


Thank you Joachim for your reply. I’ve applied the conjugate gradients and I’m getting good results for the test problem, as you described. I’ve also applied it to a eddy current problem with a known analytical solution. One of your colleagues (Christopher) previously suggested that it was possible to reduce the number of degrees of freedom by removing the gradients in the non-conducting region and this is also working well and leads to a considerable reduction in the number of degree of freedom. For the eddy current problem I’m having to use the CG solver as the GMRES led to similar post processing issues described before. This is OK, but the convergence of the residual is not monotonically decreasing. A further issue I’ve found is that if I try to include static condensation with the assembly of the matrix from the eddy current or magnetostatic bilinear form it does not give reliable results for orders above 2 either with or without gradients. Without static condensation I get good accuracy compared to the analytical result with and without gradients. I’ve attached a simple curl curl problem with a known analytical solution, which illustrates the issues I’ve found with the static condensation, even if the gradients are included.