Divergence free vector field on ref-cube.mesh

[markusrenoldner]

Dear community,

clearly the vector field $\boldsymbol{v}\in H(curl)$ is divergence free

$$ \boldsymbol{v} := \begin{pmatrix} y\\ -x\\0 \end{pmatrix}$$

It doesnt seem to be divergence free in MFEM. To check, i start with implementing the vectorfield

void vec_func (const mfem::Vector &x, mfem::Vector &returnvalue) {    
    returnvalue(0) = x(1);
    returnvalue(1) = -x(0);
    returnvalue(2) = 0;
}

i project this onto a GridFunction

mfem::GridFunction vec(&ND); 
mfem::VectorFunctionCoefficient vec_coeff(mesh.Dimension(), vec_func);
vec.ProjectCoefficient(vec_coeff);

to find the (weak) divergence i use the MixedVectorGradientIntegrator. It generates a matrix $M$ with components
$$M_{ij}:=(\nabla \phi_i, \boldsymbol{\psi}_j), \quad \phi_i \in L^2, \quad \boldsymbol{\psi}_j \in H(curl)$$

...which I then transpose and multiply with the vectorfield:

mfem::MixedBilinearForm blf_div(&CG, &ND);
blf_div.AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator()); 
blf_div.Assemble();
blf_div.Finalize();
mfem::SparseMatrix mat_div(blf_div.SpMat());
mfem::SparseMatrix *mat_div_T;
mat_div_T = Transpose(mat_div);
mat_div.Finalize();
mat_div_T->Finalize();

mfem::Vector div_vec (CG.GetNDofs());
mat_div_T->Mult(vec,div_vec);
std::cout <<div_vec.Norml2() << "\n";

The output is smaller than $1$, but not even close to machine precision on the ref-cube.mesh.
Somehow on the periodic-cube.mesh the divergence was zero up to machine precision. Can somebody explain this?

[pazner]

Hello @markusrenoldner,

When you perform integration by parts on the MixedVectorGradientIntegrator to get the divergence, you also pick up a boundary term. The periodic cube has no boundaries, so this term doesn't appear.

[markusrenoldner]

Thank you for your quick response.
I have to apologize - of course you are right and I did not give you the complete expression for $\boldsymbol{v}$.

The vector field is actually:

$$ \boldsymbol{v} = \begin{cases} \operatorname{cos}^2(K(x^2+y^2+z^2)) \cdot\begin{pmatrix} y \\-x\\0 \end{pmatrix} & \text{if }x^2+y^2+z^2 < R^2\\ \quad \text{ }0 & \text{else} \end{cases} $$

Now $\boldsymbol{v}$ is still divergence-free, it is actually twice continuously differentiable, and it is zero on the boundary, if one choses the constant $K$ and the radius $R$ of the support of $\boldsymbol{v}$ correctly.

The divergence is still not 0, as I desribed in the first post.
I computed the divergence on $\sim10^7$ elements, and the divergence ended up to be around $10^{-7}$. Is this to be expected?

[pazner]

The function $\boldsymbol{v}$ cannot be represented exactly in the finite element space, and so you are really performing computations on its interpolant, which might not be divergence-free (differentiation and interpolation do not commute in general).

[markusrenoldner]

Thanks a lot for elaborating!

[mlstowell]

Hello, @markusrenoldner,

I agree with everything @pazner wrote but I was curious to visualize the results and see where the errors appear. First, I should point out that the data stored in div_vec do not correspond directly to divergence values. This vector actually stores a LinearForm rather than a GridFunction. To visualize the errors we need a GridFunction representation which is computed from div_vec as div_gf = M^{-1} div_vec where M is a CG mass matrix in this case.

After computing the necessary GridFunction we can see that the errors follow the transition where x^2+y^2+z^2 = R^2 (I chose R=0.4):

Uniformly refining the mesh reduces the errors by a factor of 2 as we might expect.

Using higher order basis functions also improves the errors. I was also able to reproduce similar errors on the periodic-cube.mesh. Note however that the domains are different in reference-cube.mesh and periodic-cube.mesh. They domains have different sizes and are centered differently with respect to the origin. Computing errors of virtually zero on the periodic mesh may have been a coincidence. These things do sometimes happen when meshes and the functions interpolated on them share certain symmetries.

If accuracy in such a weak divergence is critical you may need to align the mesh with any interesting features of your expected solutions. In this example a mesh with edges aligned, or nearly aligned, with the sphere x^2+y^2+z^2 = R^2 would probably produce significantly lower errors.

Best wishes,
Mark

[markusrenoldner]

Ah thats a very interesting answer thanks a lot for doing the visualizations!

Just to be clear, you chose $R=0.4$ and therefore
$$K=\frac{\pi}{2R^2} = \frac{\pi}{0.32}=9.81748...$$
as $R= 0.5\sqrt{2\pi/K}$, right?

[mlstowell]

Actually, I misread the argument in the cosine. I thought it was $\cos(K R)$ rather than $\cos(K R^2)$ so I used $K = \frac{\pi}{2R}$.

Rerunning with $\cos(K R^2)$ and $K = \frac{\pi}{2R^2}$ produces a similar pattern but with a larger magnitude error. This is not surprising since the transition caused by the cosine become sharper with $R^2$.

[markusrenoldner]

Ah alright, thanks for clarifying!

[mlstowell]

Incidentally, producing H(curl) fields with zero divergence is a common issue in computational electromagnetics (often goes by the name "divergence cleaning"). MFEM has a class designed to accomplish this. It's located amongst the miniapp codes, see DivergenceFreeProjector.

Unfortunately, I just noticed that there is a small issue with this code. I'll try to submit a fix, in the next day or two, to improve its effectiveness.

[mlstowell]

I just initiated PR #3533 to improve the performance of the IrrotationalProjector and also the DivergenceFreeProjector which is closely related.