I ran into an apparent accuracy problem with NGSolve, I wonder if there is some mistake in the way I am using NGSolve or else if there is a workaround for this.
To debug this, I made a simple test case with a known analytic solution: capacitance calculation of a parallel plate capacitor. The attached .py file creates two parallel plates 1 distance unit apart with an area 1x1 unit^2 each. The solution of the Laplace equation is a simple linear potential distribution, the field (-grad u) is a constant 1 between the plates and 0 outside.
The potential distribution calculated by NGSolve looks ok, the field defined as gfE.Set( -grad(gfu)) is less clean (a bit noisy) but otherwise doesn’t look wrong. However, when I calculate surface integrals of normal field at the two electrodes I get ~0.5 instead of 1. This is how I calculate the integrals in the attached .py:
charge = (Integrate(sqrt((gfE*normals)**2), mesh, definedon=mesh.Boundaries(electrodes[i])))
Any idea what is going on and how to improve the accuracy of this calculation?
Thanks and regards,
It seems that the problem is with the way integration is carried out. Integrate(sqrt((gfE*normals)**2), mesh, definedon=mesh.Boundaries(electrodes[i])) calculates the normal field integral along the interface between the electrode (Dirichlet boundary condition) and dielectric. Of course, inside the electrode the electric field is zero, so the interface integral seems to use the average of electric field 0 inside the electrode and non-zero field E on the dielectric side. Overall the effect is to lower the calculated field integral by 2x.
Is there a way to limit integration of normal field to just the dielectric side of the boundary and skip the electrode side? If not, I guess just multiplying the integral by 2x should produce the same result to reverse the apparent averaging of field outside E and inside 0 of the electrodes.
it makes perfect sense what you get. The integration does the right thing.
You choose the H(div) - space (i.e. face elements) for the electric field, which admits a unique normal-component En on every face. The interpolation gfE.Set (grad(u)) does an arithmetic averaging of the normal component from both sides, the dielectric, and the electrode (where grad(u) is zero). From that you get half the value.
(A) Define the space for the E-field only in the dielectric, see attachment
(B) use only the dielectric region for the interpolation to the E-field. Then no averaging (with 0) happens.
I also gave a material name to the dielectric in you example file.
Thanks very much for your help, your solution A makes sense and works well.