Create Generalized Laplace Operator

[bugrahantemur]

This PR introduces a new Laplace operator implementation that can be used by any module in ExaDG.

The requirements that this new operator fulfills were listed in #357. A document with the math behind this operator is also linked in that issue.

In short, the operator is templated for the scalar and vector-valued problems, as well as coefficients that couple or do not couple between the components of the gradient.

For a proof-of-concept, this PR also removes the diffusive operator in the convection-diffusion module and replaces it with the new generalized Laplace operator. Other than this, as few places as possible are touched in this PR to keep it manageable for reviews.

---

Here are some details about the implementation:

• Operators that use the implemented kernel (these are right now the operator itself, and the combined operator in the convection-diffusion module) are responsible for implementing their own cell-loop/face-loop functions to decide how to calculate coefficients, how to store them (like having cell-based face loops or not), and how and when to update them.

• do_boundary_integral_cell_based() is introduced in the base class operator OperatorBase fully analogous to the already present do_face_int_integral_cell_based() as these two functions need to access the cell-based face coefficients unlike their non-cell-based counterparts.

• fill_matrix_free_data() in IncNS::SpatialOperatorBase is made virtual and overriden in IncNS::OperatorCoupled, since the generalized Laplace operator requires to set additional flags if it is used as a preconditioner.

---

Note: The convection-diffusion module is currently the only place where the operator is used. However, the applications that test this module do not use variable coefficients as the diffusivity is assumed constant, i.e. the variable coefficient functionality is disabled for optimization. During development, these tests were also successfully run with variable coefficients, but this is not part of the PR to keep changes to a minimum.

Closes #357

[bugrahantemur]

I'd appreciate a review if you can @nfehn @peterrum :)

[nfehn]

Yes, I have it on my agenda, but I am currently struggling to get ExaDG running on my new machine, get tests passing, etc.

[nfehn]

@bugrahantemur Now I think I got a first overview of this PR. I very much appreciate the effort that you made so far.

In my opinion, the current design has the following problem

• we have duplicated code between the GeneralizedLaplace::Operator and the CombinedOperator (all the functions calculating the variable coefficient). The code duplication is somewhat substantial, so I believe we really have to avoid this.
• the implementation seems to be somewhat limited to variable coefficient fields provided via Functions. Extending this to other types of models will probably require changes at several places in the code, because of the code duplication item above.

These problems are kind of solved in the IncNS module, where we have variable viscosity models in the form of turbulence models and generalized Newtonian models (see the recent PR #381). The solution is to use the design principle "composition", i.e. a class that operates on / manipulates the variable coefficient data structure living in Kernel. It is important to not realize this functionality inside the Operator, because this limits flexibility and causes code duplication. In principle, we want to be able to use variable coefficient fields based on an arbitrary underlying function or model, which might depend on other dof-vectors or originate from any "multiphysics" context. This modularity can in my opinion be achieved in an elegant way by composition, see the class IncNS::SpatialOperatorBase with member variables TurbulenceModel and GeneralizedNewtonianModel. The viscous Kernel or Operator, however, should not know where the variable viscosity comes from (analytical function vs. a "complex" model). In principle, one could add another model like AnalyticalViscosityModel derived from ViscosityModelBase and implementing a variable viscosity field provided by the user via a dealii::Function. Note also that in order to be able to combine different models (assuming they are "additive"), the variable viscosity models in the IncNS module implement an add_viscosity() (not a calculate_viscosity()) logic.

If not(coefficient_is_variable), we need to make sure in my opinion that the VariableCoefficient data structures storing coefficients for all q-points is not used, but is instead by-passed using the constant coefficient that is a member variable of OperatorData. Otherwise, I have the fear that we make the standard case with constant coefficient too slow in terms of computational efficiency, increasing the data transfer from main memory significantly when reading q-point data unnecessarily in the constant coefficient case. This logic with "optimized code path" for the constant coefficient case should be handled by Kernel (owning the VariableCoefficient data structure) according to the single-responsibility principle. In this sense, the VariableCoefficient object is only a feature of the Kernel that is only used if coefficient_is_variable.

[bugrahantemur]

@nfehn Your proposed approach sounds quite good, and a separate object implementing the coefficient calculation can be done easily, but I think there is an important point which may have slipped your attention, which is the operators inside the multigrid preconditioners.

Namely, the preconditioners operate on a number of actual operators, like the IncNS::MomentumOperator or the ConvDiff::CombinedOperator which exist simultaneously on multiple levels. These operators each create and own their own kernel objects (this works fine as long as every needed information is present in the kernel). If the kernel has mesh-dependent physics, like calculating the penalty parameter in the viscous or the diffusive kernels, update_after_grid_motion is called on the operator, see:

exadg/include/exadg/incompressible_navier_stokes/preconditioners/multigrid_preconditioner_momentum.cpp

Lines 234 to 240 in 9aba126

	MultigridPreconditioner<dim, Number>::update_operators_after_grid_motion()
	{
	for(unsigned int level = this->coarse_level; level <= this->fine_level; ++level)
	{
	get_operator(level)->update_after_grid_motion();
	}
	}

which in turn might do something like:

exadg/include/exadg/incompressible_navier_stokes/spatial_discretization/operators/momentum_operator.cpp

Lines 132 to 137 in 9aba126

	MomentumOperator<dim, Number>::update_after_grid_motion()
	{
	if(operator_data.viscous_problem)
	viscous_kernel->calculate_penalty_parameter(this->get_matrix_free(),
	this->get_data().dof_index);
	}

However, this cannot include updating the coefficients if the operator or the kernel is not supposed to own the coefficient calculation object. Was this clear, else we can maybe discuss in person?

[bugrahantemur]

I pushed a new design trying to fit into our constraints. Here are some points:

• I implemented a new generic AnalyticalTensorCoefficientFunctionModel which is templated for the kernels. This object can be used to evaluate a dealii::Function and set the values at the quadrature points for any kernel that supports the operations set_coefficients_cell() etc.

• In the kernel operator I added an enum to select different coefficient models. Right now only ConstantFunction and VariableFunction are supported. This was the only way I found which allows the kernel data to stay portable between different instances of the kernel, which is very important due to the MultiGrid implementation (see my comment from above). Other models, e.g. computing the coefficients from a solution vector, may add the desired dof_index to the kernel data (analogous to the convective kernel design) and implement these models independently from the kernel.

This design allowed an elegant way of avoiding code duplication, keeping the FE operators clean (as @nfehn highlighted), and also keeping the kernels as small as possible in terms of implemented algorithms but also as detailed as necessary in the current framework.

[nfehn]

However, this cannot include updating the coefficients if the operator or the kernel is not supposed to own the coefficient calculation object. Was this clear, else we can maybe discuss in person?

As a general answer to this question: You are right that multigrid has its own design challenges. In ExaDG, the general/current strategy is to ask the fine-level operator for the data, and then interpolate the data to coarser multigrid levels. We do this currently e.g. to interpolate the velocity vector (IncNS) or displacement vector (Structure) to coarser levels. These vectors originate from the linearization inside the Newton solver. Generally, I would argue that one does not even have the information to recompute the variable coefficient on the coarser levels as it is done on the fine level. The variable coefficient might dependent on the solution vector, which is only known on the fine level and not on the coarser levels. If I take the example of a turbulence SGS model: Interpolating the fine-level turbulent viscosity to coarser levels is something different than interpolating the velocity to coarser levels and recomputing the turbulent viscosity on the coarser levels (using the viscosity model). However, I think the former is what one wants to have for multigrid?
In the Structure module, the operator owns the MaterialHandler, which could also be seen similar to the variable coefficient model. However, I think this operator is not used in a "combined" fashion with other operators. Since we want to write a generalized Laplace operator, I believe we can not follow the design chosen in the Structure module.

So in my opinion the procedure to enable multigrid should be:

1. Ask the fine-level operator for the data
2. Interpolate fine-level data to coarser levels e.g. in MyModule::MultigridPreconditioner::update() (this is an open problem in ExaDG, where implementations still have to be developed; questions: q-point vs dof-data layout for interpolation; which type of interpolation for structure-preservation)
3. Inject the coarse data into the coarse-level operators/kernels.
4. Ready for multigrid

The philosophy for our geometric multigrid algorithm implementation in ExaDG is:

The fine-level operator must be able to provide all the data to construct the coarse-level multigrid operators.

The general design of the PDE operator, however, should not be driven by the multigrid preconditioner. To solve the current problem with an analytical variable coefficient model regarding multigrid: I see no harm if the multigrid preconditioner for a certain module depends on the spatial operator of this module, which is the case for most multigrid preconditioners in ExaDG. In your case and to avoid the problem of interpolating fine-level data to coarse levels, one could easily pass the analytical variable coefficient model to this particular multigrid preconditioner during initialization. Just as the variable coefficient model is "module-dependent", the multigrid preconditioners may also be module-dependent. I see this as the burden/price of geometric multigrid.

[nfehn]

"Multiplication of what with coefficients of different ranks". The word "different" is misunderstandable. I would write "... from the left with a coefficient tensor of rank coefficient rank".

In my opinion, the variable names T, coeff_mult, x are too short. Also 2,0,3,1 is somewhat hard to understand.

[bugrahantemur]

I'll change the sentence, rename the variables and put a comment about the double contraction👍

[nfehn]

Let's also make the function clearer by using full names.

[nfehn]

See comment #414 (comment)

[nfehn]

"interior penalty" should be explained for this generalized Laplace operator. Does the operator support both continuous and discontinuous Galerkin discretizations?

[bugrahantemur]

I've done the implementation with only DG in my head. So I'm not sure if there are any intricacies apart from the interior penalty that I missed about continuous Galerkin. What do you think?

[nfehn]

I do not fully understand this. What about a 1-component solution field (gradient has dim components) with a diagonal 2-tensor as variable coefficient? The result of coeff_mult() would again be a vector with dim components. Would one call this a "coupling"?

[bugrahantemur]

See answer #414 (comment)

[nfehn]

Would the current implementation support a diagonal 2-tensor as variable coefficient? Please add some documentation if not already done.

[nfehn]

I would like to avoid this dependency of the analytical coefficient model on the kernel. The connection is actually the VariableCoefficients class and not the kernel in my opinion, just as you placed this functionality in the file variable_coefficients.h. I would not template VariableCoefficients here, but let the analytical model explicitly use the class VariableCoefficients.

[nfehn]

This code shows that kernel and the analytical model are strongly coupled from an implementation point of view (see arguments This as template and this; comment https://github.com/exadg/exadg/pull/414/files?diff=unified&w=0#r1206535671), while the general goal should be to minimize dependencies. I think the design would improve if kernel does not depend on the analytical model (and the other way round), however, both the kernel and the analytical model may depend on VariableCoefficients (as the low-level data structure).

[nfehn]

can one ask the kernel for these types (to avoid this duplication of code)?

[bugrahantemur]

I guess I could define them in the namespace as alias templates, so outside the classes instead of inside them. This is maybe even better. Not sure if one could ask the kernel, probably would work as well :)

[nfehn]

see comment #414 (comment)

[nfehn]

We should find a way to make this information unique, e.g. by referring to the kernel.

[nfehn]

Syntax: I find this syntax very complicated. Writing analytical_coefficient_function.evaluate(variable_coefficients) would be clearer / easier to understand (without need for documentation). Extracting that this line of code modifies variable_coefficients took me several minutes. Hence, this style of programming requires much more documentation in my opinion. I doubt this additional effort in documentation justifies the one line of code saved by this programming style.

Design: More importantly, this design does not solve my main criticsm from the last review, namely that GeneralizedLaplace::Kernel should in my opinion not depend on AnalyticalTensorCoefficientFunctionEvaluator.

Example: When adding a new type to CoefficientModel, the kernel/operator has to be recompiled. I consider this an unnecessary dependency.

If we do not want to provide access to Kernel::variable_coefficients via a public getter returning a reference to VariableCoefficients (in order to fill the variable coefficients from "outside"), a lambda function with argument VariableCoefficients & could be passed to the kernel during construction/initialization, which allows to compute the coefficient field.

What to do if one wants to realize an analytical coefficient function that is time dependent? I believe one has to touch again Kernel with the current design, but this is unnecessary. Would you adjust the constructor of AnalyticalTensorCoefficientFunctionEvaluator and make double time a member variable of that class? This is probably not a good design / unnecessary, because one rather wants to call the evaluate() function suggested above with a second argument time.

[nfehn]

For CoefficientModel::Constant, why should one provide a dealii::Function at all? Then, the dereferencing * will fail.

This highlights another disadvantage of the currently chosen "framework" design. It appears odd to add a separate member constant_coefficient to KernelData if this data structure already contains std::shared_ptr<dealii::Function<dim>> coefficient_function;. So to avoid such a redundancy, one forces "everyone" to use the dealii::Function interface.

However, the kernel does not require such information! Only the classes using the kernel or operator (which might originate from vastly different application scenarios) need this information.

[nfehn]

I am not a fan of such global definitions. To give an example: Why should one need to know outside the implementation of kernel/operator that the implementation is based on dealii::VectorizedArray?

[nfehn]

Can be considered error-prone when adding new options to CoefficientModel. Again, I see the origin of the problem in an unnecessary dependency.

[nfehn]

It is not clear to me why the documentation of these functions should explain "weighted by ...". Implementation of the weak form is in my opinion the responsibility of the Operator.

[bugrahantemur]

I wrote the documentation like this to make it clear where the name of the function comes from. Because that is our motivation to name this function get_value_flux(). The caller shall multiply the return value with the shape function values. The naming of the function is not obvious from its body, which for example also works with gradients here. I just wanted to prevent any confusion caused by the function name.

[nfehn]

where does the minus sign next to return come from?

[bugrahantemur]

That's because I didn't want to implement operator as the negative Laplace operator. But of course I can change it to be the negative Laplace as in other parts of ExaDG.

[nfehn]

I do not like that the present PR makes such subtle sign changes.

[nfehn]

I think explanations/motivation for these numbers would be helpful here.

[bugrahantemur]

As of now, there is a new version of the implementation for the Generalized Laplace operator.

In the new implementation,

• the kernel and the coefficient calculation is completely decoupled,
• the type aliases are put in the kernel, the operator inherits these from the kernel,
• the kernel holds a variant (an actual std::variant) of either a constant coefficient or a VariableCoefficients object,
• there are Get***Coefficient visitors for the variant for efficient access to the coefficients from the variant in a namespace called VarCoeffUtils,
• in the same namespace, an AnalyticalTensorCoefficientFunctionEvaluator is implemented, which can fill a VariableCoefficients object based on a dealii::Function,
• the ConvDiff::Operator implements initialize_diffusivity() and calculate_diffusivity() to manipulate the coefficients in the diffusive kernel,
• these two functions implement the logic for different diffusivity calculation models (a module enum parameter),
  • Constant, AnalyticalTimeFunction, AnalyticalTimeAndSpaceFunction, only the last one requires VariableCoefficients
• the Multigrid preconditioner of ConvDiff gets a constant diffusivity and runs with this, even in the case of variable diffusivity,
• the pressure convection diffusion operator in IncNS can now actually work with variable viscosities.
  • as a side effect, it was necessary to make the viscosity a tensor of rank 0 in the viscous kernel to realize this.

[bugrahantemur]

I'm wondering if std::visit to get coefficients on std::variant<VariableCoefficients<Coeff>, Coeff> in the innermost loop is going to ruin performance as it won't be inlined...

[nfehn]

I think one should ask the pde_operator for the diffusivity instead of extending the interface.

[nfehn]

Why returning a non-const reference and then overwrite the value? I think I would then prefer a set-function as we usually do in ExaDG.

[nfehn]

The problem is that this code is not placed correctly in update_after_grid_motion(). Also as a user, I can not expect that for a time-dependent problem on a non-moving mesh, the calculation of a time-dependent diffusivity will be done in a function that explicitly relates to grid-motion (ALE) problems.

[nfehn]

what is get doing here? Could you explain and/or add a comment?

[nfehn]

Let me raise a general question: I see the novelty of the present PR in the introduction of a generalized Laplace operator. However, it does not really seem to fit optimally that this is new implementation is then used in the present PR to replace the diffusive operator in the ConvDiff module. So far, the ConvDiff module does not support variable diffusivities, while we have variable diffusivities/viscosities in the IncNS module, which is therefore much more suitable as an application example for a generalized, variable coefficient Laplace operator. As a consequence, the present PR has to introduce new functionality for the ConvDiff module (functions like calculate_viscosisty() in the spatial operator, changes to the multigrid preconditioner), even though there seems to be no motivation for such functionality currently from an application perspective. As a consequence, I believe the present PR mixes too many topics, that make code review unnecessarily complex and somewhat hinder the progress in the sense of merging new functionality step by step. So I want to raise the question whether it would be better to first introduce the generalized Laplace operator before using it in a concrete module in ExaDG.

[nfehn]

Related to the comment #357 (comment), I believe the present PR needs to decide whether to support for now variable coefficients for the generalized Laplace operator or not. I mainly see the following two options:

• If not, using the new generalized Laplace operator for the diffusive operator in the ConvDiff already in the present PR is reasonable since this operator does so far not support variable coefficients. Then, no redundancies are introduced by the present PR. The present PR replaces ConvDiff::DiffusiveOperator by a more general implementation. The present PR does, however, not extend the functionality of the ConvDiff module.
• If yes, the implementation of the generalized Laplace operator introduced here probably needs to care more about the DG methodology for the generalized Laplace operator, and this should probably be introduced stand-alone, i.e. the present PR should really focus on such aspects and not on changes in the ConvDiff module. This would imply that we have for now some redundancy with IncNS::ViscousOperator, but this is acceptable since we have the plan to replace this functionality later by the new variant. If no one wants to do this in the future, one would simply remove the new generalized Laplace operator again from ExaDG, which is easy since it is not used so far. Of course, this is not intended, but I just want to describe the general procedure to avoid code duplications in ExaDG over a longer period. In any case, I think this is better than introducing another variable coefficient operator in ConvDiff actually not needed right now.

[nfehn]

Could you explain why this is needed? Why didn't we need it before, or has it been a bug before?

[nfehn]

I find it un-intuitive to see from this code whether just a reference is copied or the whole data. I mentioned somewhere else that getter functions returning non-const references are also more difficult to understand in my opinion.

Shouldn't such code look like this:

if(this->viscous_kernel->viscosity_is_variable())
  pressure_conv_diff_operator->set_diffusivity(this->viscous_kernel->get_variable_viscosity());
else
  pressure_conv_diff_operator->set_diffusivity(this->viscous_kernel->get_const_viscosity());

[nfehn]

This should not be realized inside variable_coefficients.h.

[nfehn]

To be honest, I do not see the advantages of the variant/visitor design patter here. In my opinion, one needs to write more lines of code, and the code is more difficult to understand than a very simple if-else logic. Since we are here in the generalized Laplace operator, there is no concern right now that this if-else logic needs to be duplicated at multiple places. And even in such a case, one could also introduce free functions in the future, but only once we need this complexity. So where is the advantage of the present design? Defining two functions operator() with the same name is essentially a complicated way of implementing an if-else logic.

[nfehn]

In terms of code design, I am not convinced that it is a good idea that the most complex case (with face, q variables) dictates the "interface" also for the simplest case of a constant coefficient (not even mentioning that it is quite hard to identify an "interface" from this struct definition). While I do not see many other variants to be realized in the future, I also do not understand the motivation for the present design if we had more variants in the future. If all these variants have different interfaces, the interface of this struct will grow according to the number of all variants, just as the simple if-elseif logic.

[bugrahantemur]

I think this PR is stuck right now between trying to have a use for the generalized Laplace operator, and the features it needs to have due to its original motivation, which was to be used in the context of the Darcy equation. I am strongly against merging such a big contribution without using it anywhere. On the other hand, trying to use it with existing modules, as pointed out, then causes problems, for example by introducing features that are not application motivated. The interfaces with the modules become a problem, as well.

In light of this, I think it is wiser that I further develop this operator together with the Darcy module in my own repository, use and test it there. Also concerning external constraints like development time, I think this is the best way to go. After the operator has been used there successfully with the Darcy module, we can later decide what parts of it we want to have in ExaDG and how we want it specifically. We can the import parts here step by step. This would be easier for everybody, I guess.

So, I would actually suggest closing this PR for now. What do you think @nfehn ?

[nfehn]

Overall, I think we made progress in the course of this PR and developed ideas that I consider useful, so it would be somewhat sad if we do not finalize things, fully or in parts. I expressed my ideas regarding the current status and content-related aspects here: #414 (comment) and #414 (comment).

Realistically, the developments in ExaDG will go on and it might happen that codes then start to more and more diverge. I see ExaDG also as a platform for collaboration and exchange. It might be that this will become more difficult when working in separate repositories. My very personal assessment is to continue in ExaDG. I think it will pay off.

I want to emphasize that, in the end, it is totally your own decision how to proceed, because only you know your constraints. I can offer help/advice regarding DG discretizations / variable coefficients (and, potentially, also software design), but due to constraints on my side I can not invest much time if things are not going to be brought into ExaDG in the near future.

[bugrahantemur]

Yes, we have come a long way due to our continuous efforts and that's great. What I'm talking about is just going to be transient (I'm planning something like 2-3 months maximum). In the meantime, I will keep up with the developments in the main repo, but I consider this a necessary step for the speed of development I need right now (conference mid-September 😅) As soon as I have things working and got some results, I want to put it all into the main repo, just as was our plan all along. Of course, it will all build on the state (of this PR) we have achieved until now.