Request for input: Laplace-based linear and bilinear integrators

[IdoAkkerman]

I am interested in implementing LinearFormIntegrator and BilinearFormIntegrator involving the Laplacian.

Any heads up on potential conflicting or duplicate effort is highly appreciated.
Any tips or warning for pitfalls are also welcome.
Last but not least, if I am skipping some useful functionality for other use cases please let me know.

My personal interest for these terms are stabilized methods. So first would be the scalar convection-diffusion case which would require the following:

• $(\Delta w, f)$ Linear form for a given Coefficient f
• $(w,q\Delta \phi)$ Bilinear form for a given Coefficient q
• $({\bf a}\cdot \nabla w,\Delta \phi)$ Bilinear form for a given VectorCoefficient a
• $(\Delta w,q \phi)$ Bilinear form for a given Coefficient q
• $(\Delta w,{\bf a}\cdot \nabla \phi)$ Bilinear form for a given VectorCoefficient a
• $(\Delta w,q\phi)$ Bilinear form for a given Coefficient q

Next would be vector based case, such as navier-stokes. Which would require the same as above but all the scalar shape functions become vector shape function. The coefficients will in principle be of Matrix form, probably more efficient routines are available for scalar or vector (= diagonal of matrix). I will consider implementing all three.

To enable the pressure and PSPG terms two mixed terms are required.

• $(\nabla q, Q \Delta {\bf u})$ Bilinear form for a given MatrixCoefficient Q
• $(\Delta {\bf w}, Q \nabla p)$ Bilinear form for a given MatrixCoefficient Q

Hessian based integrators will not be considered in this PR.

[IdoAkkerman]

Implementation effort can be reduced by using the TransposeIntegrator I assume?

[IdoAkkerman]

Is there a need for a Trans.LaplaceOrder(fe) ? As place holder I am going with fe.GetOrder() - 2 for now.

[v-dobrev]

This sounds great to me!

Is there a need for a Trans.LaplaceOrder(fe) ? As place holder I am going with fe.GetOrder() - 2 for now.

I'm not sure we really need to target exact integration of Laplacian terms, even if that is possible in some cases. I think using the lowest quadrature order (e.g., ignoring the transformation order) that works robustly is best. It is hard to test all possible use cases for "works robustly", so go with what you find in your tests or what is recommended in the literature, if there are such recommendations.

How do you plan to compute the Laplacian on a general element? Do you plan to use exact computation at quadrature points using the Hessian and the Jacobian of the transformation? An approximate alternative that only requires the Jacobian is to project the gradient (e.g. using nodal interpolation) back into the same FE space and then take the divergence. I think @bslazarov and/or @vladotomov wanted to try something like that but I'm not 100% certain.

[IdoAkkerman]

The Laplacian/Hessian is already implemented - in principle.
The base class has the (some) interfaces defined and for NURBS elements it should work.
I do not know about other elements.
There is likely still room for extensions, improvements and simplifications.

The Laplacian is compute as the trace of the Hessian which sounds expensive but is the easiest way to have the correct effect of the mapping . There is a check of to see if the Transformation is affine,
in that case a faster, more direct, way of computing the Laplacian is available.

[IdoAkkerman]

@michi002 My apologies, I pushed it to an other local repo. I will push it tomorrow, unfortunately I can not ssh into my machine due to overly zealous security...

[michi002]

@IdoAkkerman thank you. I had a look on your commit.
I think you only need non-mixed spaces for your purpose. However I would also need
$$ \langle \Delta u, v \rangle $$
for v piecewise constant.
Is it ok for you, if I just add the mixed version in your branch? Or should I branch out from your branch for that?
I will adpat to your workflow, as you started it :)

[IdoAkkerman]

@michi002 Convection-diffusion only requires the scalar implementation. However, the SUPG/PSPG terms for Stokes/Navier-Stokes involve terms where the velocity and pressure spaces are combined, or is there a way to bypass that?

Regarding, the term you are interested in, that can be obtained by by using the 4th term, viz.
$(\Delta w,q \phi)$ and choosing a zero order L2 finite element for $\phi$ (v in your comment) .

[michi002]

@IdoAkkerman
Yes, this is the bilinear form I need. But we also need to also implement the member-routine for mixed spaces
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat);
or not?
I have only seen the routine for the same ansatz and test-function
/// Support for use in BilinearForm. Can be used only when appropriate. virtual void AssembleElementMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &elmat)

Concerning the coefficient:
I guess you are right, If the coefficients are differentiable one can just use the chain rule.

[IdoAkkerman]

Implementing AssembleElementMatrix2 should be minimal effort, so worth to include in the PR I guess.
However, the use case for AssembleElementMatrix2 is not clear to me. By itself not so much of a concern, but this would make it a bit harder to debug and test (unit_test).

Do you have a clear application in mind?

[michi002]

@IdoAkkerman
Yes, I want to apply it in distributed optimal control. I want to solve the equation
$-\Delta u = z$
for $z$, where $u$ is given. For implementing control constraints we want a simple connection between $z$ and $u$.
If we choose $u$ to be a $C^1$ spline and $z$ to be piecewise constant we can easily do that by
$\underline{z}_h = M_h^{-1} K_h \underline{u}_h$,
where $M_h$ is the (diagonal) mass for piecewise constants and $K_h$ is defined via
$\langle -\Delta u_h, p_h \rangle$
and here $p_h$ is also piecewise constant.

I do not know if it is an important case for the general userbase, but I think a mixed version (AssembleElementMatrix2) might come in handy for some users.

[IdoAkkerman]

So you are converting
$-\Delta u = z$
into a discrete weak form
$\langle-\Delta u^h, p_h \rangle =\langle z, p_h \rangle$
where $u^h$ are $C^1$ spline and $p_h$ are piecewise constants.
I was contemplating this as test case :)
But I do not see the logic of doing this... but perhaps it is something application specific.

Why not test with $C_1$ splines similar as used for $u^h$ -- and use integration by parts.
$\langle\nabla u^h, \nabla p_h \rangle - \langle\nabla u^h, n p_h \rangle_\Gamma=\langle z, p_h \rangle$
This will yield the standard optimal Galerkin formulation for the Poisson problem.

[michi002]

@IdoAkkerman
I am considering this for distributed optimal control, where $z$ (the control) is the quantity of interest.
We want to approximate the control $z_h$ via piecewise constants. The Laplace problem for the state, I solve as you suggest with $C^1$ test and ansatzfunctions.
So suppose the state $u_h$ is given (via a minimization problem), I want to reconstruct the right hand side, that gives this state.
Then I have to find $z_h \in Z_h$ such that
$\langle z_h, p_h \rangle = \langle -\Delta u_h, p_h\rangle \quad \forall p_h \in Y_h.$
If I now take z_h to be piecewise constant, I either have to take $u_h$ to be a $C^1$ function, to apply the strong Laplacian and $Y_h=Z_h$, or $Y_h$ to be piecewise linear continuous functions, to apply integration by parts, yielding a rectangular left hand side. In the second case I would have to solve some kind of saddle point problem, which I would like to avoid.

If you are really interested why we do that, you can look at our preprint: https://arxiv.org/abs/2404.10350
But I tried to summarize it. In the preprint we also did not consider control constraints yet, where this form is also advantageous.

I hope that answers your question.

[michi002]

@IdoAkkerman I have now implemented a mixed version of your LaplaceIntegrator called MixedLaplaceIntegrator.
I have pushed it to your branch (I hope that is ok).
I think in your implementation of it, there is a small error.
In line 1559 I think you want to calculate the laplacian of the trial functions, instead of its shape:

el.CalcPhysShape(Trans, shape);
el.CalcPhysShape(Trans, laplace);

I would like to test my routine, by projection the laplacian of a $C^1$ spline to piecewise constants defined on the same nurbs mesh. Then I would refine the mesh and look at the convergence order.
But I do not know how you do it at MFEM? I just saw big files for unit tests for bilinear form integrators in dimensions one to three, which were a little confusing for me.

[IdoAkkerman]

@michi002 Thanks for the input, I was to lazy to change the name. I think the functionality is okay, but you are right for a final commit this needs to be rectified.
I do not see the need for a a mixed class, doesn't the LaplaceIntegrator::AssembleElementMatrix2 member function not take care of that?

[IdoAkkerman]

@michi002 Sorry you are correct. I did make a typo! Thanks!!

Perhaps I am also mistaken about the mixed stuff, but you might need to convince me...

[michi002]

@michi002 Thanks for the input, I was to lazy to change the name. I think the functionality is okay, but you are right for a final commit this needs to be rectified. I do not see the need for a a mixed class, doesn't the LaplaceIntegrator::AssembleElementMatrix2 member function not take care of that?

@IdoAkkerman
Yes. But the AssembleElementMatrix2 is not implemented for this class. Also I thought the virtual function for that is also only available in the mixed base class, but it seems it is actually there. This also confuses me. @tzanio can you please explain why there are two base classes then?

[IdoAkkerman]

I would not be surprised if bilininteg and lininteg could benefit from some pruning and rationalization.
However, this might be risky in terms of legacy stuff. Perhaps something for mfem 5.0.

My guess is that some of the mixed implementations have become defacto deprecated.

[michi002]

@IdoAkkerman
But you did not implement the mixed method. Should I combine my routine with your class or was this on purpose?

[IdoAkkerman]

The AssembleElementMatrix2 routines starting on lines 1596, 1650 and 1731 of the cpp-file, should provide the integrator for the mixed cases.

[michi002]

@IdoAkkerman
Can it be, that it is not yet pushed to the server? I do not see any AssembleElementMatrix2 routines in bilininteg.cpp for the new classes on the dev-stab-mini branch.
Sry for the spam, but there seems no possiblity to send private messages here on github.

[Rickbude]

I am very interested in this integrator (I would like to use it as a BoundaryIntegrator, for implementing third-order asymptotic boundary conditions for Poisson's equation, see also #4281). So I tried out this LaplaceLaplaceIntegrator, but I did not get the results I expected.

I have not yet completely ruled out implementation issues or derivation issues on my side, but one thing I noticed is that this integrator seems to produce very large matrix entries compared to, for example, the DiffusionIntegrator. This in turn destroys the system matrix condition number (verified with MATLABs condest). This seems to be mainly due to the values returned by CalcPhysLaplacian scaling with 1/(element size)^2, as opposed to the values returned by CalcPhysDShape that grow with 1/(element size). Is this really what you would expect? Below some measurements:

Element size (m)	Order	Normlinf(Laplacian)	MaxMaxNorm(Dhshape)
10 $\times$ 10	3	0.46	0.23
10 $\times$ 10	4	0.91	0.45
1 $\times$ 1	3	46	2.3
1 $\times$ 1	4	91	4.5
0.1 $\times$ 0.1	3	4651	23
0.1 $\times$ 0.1	4	9142	45
0.01 $\times$ 0.01	3	465166	233
0.01 $\times$ 0.01	4	1.21721e+06	449

Furthermore, I modified example 0 to use the LaplaceLaplaceIntegrator instead of the DiffusionIntegrator, but found that any refinement (such that you have more than 1 cell in the mesh) makes the problem practically insolvable. The system matrix will have horrible condition numbers, and not be full rank anymore. Looking online, it seems like H1 elements are not suitable for this problem, because C1-continuity is required. But then I don't understand why the miniapp in miniapps/stabilized/biharm.cpp doés work, what is the key difference in that code?

I added my code as attachment: ex0_biharmonic.zip

[IdoAkkerman]

@michi002 Apologies, my own workflows are confusing me (I had only pushed to my own local repo).
The mixed routines should be in now.

[IdoAkkerman]

@Rickbude Good to know there is more interest.

1. I have not considered whether the integrators can be used on the boundary, or not... I am not directly seeing an issue, but just as a general word of caution.

2. The element-size scaling you report seems correct to me. The shape itself does not scale with the element-size. These are direct consequences of the 0, 1-st and 2-nd derivative (think about what an fd formula would look like).

3. As a consequence of point 2, a direct C1 discretisation of a biharmonic problem has an unfavorable condition number when compared with the equivalent mixed H1 discretization. I am not a preconditioning
   specialist, but I have heard multiple people make this comment.

4. The LaplaceLaplaceIntegrator is not a drop in replacement of DiffusionIntegrator.
   One solves the biharmonic problem, the other solves the Laplace problem.
   That being said... you should be able to reuse the code and get something. However, I think there are 2 things to keep in mind:
   a)The biharmonic problem requires, on paper, other BCs than the Laplace problem.
   In practice this can still be okay if you enforce a Dirichlet BC everywhere (Laplace would be okay to have it on only a portion of the boundary - perhaps biharm can also live with it??). For the biharmonic problem a second BC all be implied in this case (normal of third derivative set te zero, I beleive).
   b) Because of the higher order derivatives in the biharmonic problem, H1 elements are indeed not sufficient.
   You need H2 (or something slightly less restricted) which a single patch C1 NURBS discretization would provide.

[Rickbude]

@IdoAkkerman thanks for your reply! The element scaling the way it does due to it being equivalent to the 2nd derivative is so logical, I don't know why this did not occur to me... But thanks for confirming this is correct behavior nonetheless. I can now scratch of a suspect for why my boundary condition might not be working..

I realize that the LaplaceLaplaceIntegrator is not a drop-in replacement for the DiffusionIntegrator. I was indeed trying to solve the biharmonic equation as a simple test case for this integrator, with homogeneous dirichlet BCs as you suggest. But with H1 elements, which are thus insufficient. I am not familiar with NURBS discretization, so I will have to read up on that. Thanks for these pointers!

[IdoAkkerman]

@Rickbude Perhaps look in the miniapps/nurbs directory for how to initialize a NURBS FESpace.
I will add a biharmonic problem to that directory in the dev-stab-mini branch. I implemented it to check the integrators.

[michi002]

@IdoAkkerman

Now I have the latest commit. However I think there was maybe some mismatch between somewhere.
CMake now complains

CMake Error at config/cmake/modules/MfemCmakeUtilities.cmake:57 (add_library):
  Cannot find source file:

    /home/mreichelt/Desktop/git_repos/mfem/linalg/tadvector.hpp

  Tried extensions .c .C .c++ .cc .cpp .cxx .cu .mpp .m .M .mm .ixx .cppm
  .ccm .cxxm .c++m .h .hh .h++ .hm .hpp .hxx .in .txx .f .F .for .f77 .f90
  .f95 .f03 .hip .ispc

and this tadvector.hpp is also not in master branch.

Nevertheless I looked at your routine for the laplace integrator. First there is again the same bug with the shape instead of the laplacian. Further I think there might be a problem, as we calculate the laplacian of the trial function. However the assembly routine takes the transformation of the test-function. I am not sure if this works, if I take for example the L2 space of piecewise constants over the mesh. At least in some test I had problems. This is why I implemented it via the adjoint operator in the first place. Maybe the transformation is the same also in e.g. the L2 case?

[IdoAkkerman]

@michi002 My apologies the branch is under development... I am working on several things, and I did not consider other users. I removed the tadvector.hpp file from the CMakelist.txt file. That should fix the issue...

I am not seeing the bug, its looks fine to me.
Lines 1556, 1557 1590 1591 1635 1636 1678 1679 1720 1756 & 1757 all look fine to me.
Please could you point me to the explicit line, if you think there is still an issue.

I do not understand the quote:

However the assembly routine takes the transformation of the test-function.

The transformations are taken via the Trans entity, which is in principle its own thing.
Depending on how the routine is used this could be related to the test or trail space (in case of an isoparametric approach).

[Rickbude]

d085026 removed a few cpp files from minapps/stabilized (laplace.cpp, biharm.cpp), but they are still referenced in miniapps/stabilized/CMakeLists.txt, causing configuration errors. Should these miniapps be removed from the CMakeLists, or should the cpp files be added back? Also I noticed that tadvector.hpp (and taddensemat.hpp) is still referenced in linalg/autodiff.hpp, causing compile errors (should these files be moved over from miniapps/autodiff?)