Calculating residual in finite element space with dirichlet boundary condition

Main Forum Old-Forum - do not post here
1

Hi!

I am working on a program that can calculate the induced eddy current given a periodic changing magnetic field in the frequency domain.
For the domain where eddy durrents can occur I use the following finite element space:

fes = HCurl(mesh, order = 3, nograds = True, complex = True, definedon = “iron”, dirichlet = “gamma_iron”)

I need the dirichlet boundary condition to make sure that the eddy current cannot escape the iron domain.
When using the inverse or the CG solver the system of equations converges and the result is physically correct.
To make sure the result is correct I want to calculate the residual of the problem: r = a*sol - f
r := residual, a := bilinear form, sol := solution, f := linear form

But the residual is not as low as expected although the solver says otherwise.
This only occurs in this finite element space when there is a dirichlet boundary condition.
I tried deleting all non-free dofs from the matrix and vectors using fes.FreeDofs(True) but that did not change the result.
Is there a way to calculate the residual just like the CG solver does?

Here is a short summary of all settings I set:

fes = HCurl(mesh, order = 3, nograds = True, complex = True, definedon = “iron”, dirichlet = “gamma_iron”)
a = BilinearForm(fes, symmetric = True)
a += SymbolicBFI(…)
a.Assemble()
a_inv = a.mat.Inverse(freedofs = fes.FreeDofs())
f = LinearForm(fes)
f += SymbolicLFI(…, definedon = mesh.Materials(“iron”))
f.Assemble()
sol = GridFunction(fes)
sol.vec.data = a.invf.vec
res = GridFunction(fes)
res.vec.data = a.matsol.vec - f.vec

2

Hi Nils,

you gave the answer yourself: The residual of the linear equation should be 0 only for the non-Dirichlet rows.

You can define a projector
proj = Projector(fes.FreeDofs())
which is a diagonal matrix with diagonal entries 1/0 corresponding to the BitArray.

r2.data = proj * r

Conjugate gradients computes the error via
sqrt ( InnerProductd (r, pre *r) )
where pre is a preconditioner. The range of the preconditioner is exactly the sub-space corresponding to non-Dirichlet dofs.

Best,
Joachim

3

Hi Joachim.

Thank You for the quick help. I did not know you could create a projector like that.
I solved my problem by calculating the residual like this:

fes = HCurl(mesh, order = 3, nograds = True, complex = True, definedon = “iron”, dirichlet = “gamma_iron”)
a = BilinearForm(fes, symmetric = True)
a += SymbolicBFI(…)
a.Assemble()
a_inv = a.mat.Inverse(freedofs = fes.FreeDofs())
f = LinearForm(fes)
f += SymbolicLFI(…, definedon = mesh.Materials(“iron”))
f.Assemble()
sol = GridFunction(fes)
sol.vec.data = a.invf.vec
res = sol.vec.CreateVector()
res.data = a.matsol.vec - f.vec
res2 = sol.vec.CreateVector()
proj = Projector(fes.FreeDofs(), True)
res2.data = proj*res
res_value = numpy.absolute(sqrt(InnerProduct(res2,res2)))

Best regards,
Nils