Discussion: Design of face loops (acoustics module)

[nfehn]

Follow up to PR #587, related to discussions https://github.com/exadg/exadg/pull/587/files#r1389376858 and https://github.com/exadg/exadg/pull/587/files#r1389390122. I also added a .cpp file, because the code was otherwise not built.

The second commit shows how to get rid of the template parameters. By using default arguments for pressure_p, velocity_p in face_kernel(..., pressure_p = pressure_m, velocity_p = velocity_m), we could even get rid of the unused function parameters passed to the face-kernel in case of boundary face loops. The template parameter weight_neighbor still exists, but could now actually be removed when using the design suggested here.

Note that the new design needs essentially the same number of lines of code. The new version should however be more flexible according to discussion https://github.com/exadg/exadg/pull/587/files#r1389376858. It also does not need template parameters. Hence, I would suggest to change the design.

@jh66637 FYI

[nfehn]

here, one can now easily call e.g. a free function calculate_interior_value() taking into account homogeneous/inhomogeneous/full operators, see discussion https://github.com/exadg/exadg/pull/587/files#r1389376858

[jh66637]

Thank you for proposing another design. Can we postpone this discussion weather the current design is changed in this way until I filed a PR, how to address non-matching discretizations? This will have some implications.

[nfehn]

Could you explain what is different for the non-matching case in terms of the required interfaces? From previous discussions I thought that the get_value() interface of BoundaryFaceIntegrator is the relevant one to also treat the non-matching case, where another type of (Boundary)FaceIntegrator will be needed/used. I tried to propose a design that is in line with this requirement/interface, but is more flexible w.r.t. the evaluation of interior values (homogeneous, inhomogeneous, full operator). I would like to continue the discussion in this direction.

[jh66637]

Dependent on the non-matching methodology, pressure_m and pressure_p will be of type FEPointEvaluation.

I am currently working on a tutorial how to make use of Nitsche-type mortaring in deal.II. I really think we should postpone this discussion until the respective PR is filed and discussed in deal.II.

[nfehn]

If pressure_m is of type FEPointEvaluation, does this imply that even more template parameters will be necessary for the face_loop() as compared to the present code version with 3 template parameters?

[jh66637]

The design I would choose has additional template parameters. At this point, I want to emphasize that using template parameters also has advantages and that compilation times can still be kept low by shifting the relevant parts to source files.

In any way, we can discuss options on how to avoid the template parameters, but I think it would be good to postpone this discussion.

[nfehn]

In PR #587, we discussed to introduce template parameters because they will be needed later for the non-matching case (even though they are not needed right now for the functionality currently available, i.e. matching case only). If I understand you correctly, we would need even more template parameters to realize the non-matching case in the future.

Then, however, it is not clear to me why we introduced template parameters at all in PR #587 (and I would like to refer to those discussions where we exchanged pros and cons or different perspectives on the use template parameters). I think the dis-agreement currently is whether we base discussion on the present code, or on future code. In my opinion, the standard procedure would be the former variant and to keep the code as simple as possible for now (according to the current requirements), and add complexity once it is needed or find solutions to problems once they occur.

The implication of the latter variant is in my opinion that, realistically, we already spent more time with PR #587 (and references to upcoming PRs) than following - for now - the standard implementation procedure in ExaDG for matching problems and introducing the non-matching case (with potential design changes) once we address this concretely.

[jh66637]

I think the current design is superior to the design in the rest of ExaDG since the code that calculates the face integrals is only written once. This was made possible with the template parameter. Personally, I don't like the changes in this PR. The selling point here is that all template parameters are removed. For this the interface is changed. Personally, I prefer the interface in the current code, and I don't see a lot of advantages removing the template parameters here. Removing the boolean even introduces another if statement in the innermost loop which might again have performance implications (we should at least keep in mind to quantify the performance difference), as discussed similarly before for other template parameters (where we also agreed to take a closer look at the performance differences). I know we disagree on the general use of template parameters in this part of the code, and I can also live with the proposed interface in the meantime.

I was continuously asking to discuss this later, since we will have the same discussion again when introducing Nitsche-type mortaring. I don't think the presented design can be extended to mortaring (the current design is easy to extend), but I might be mistaken. I know this is referencing upcoming code, but I wanted to point out that I see implications later to save us some time.

If you want to, you can go with proposed changes. I will most certainly have to undo these changes, and we can discuss the design of the function at this point again.

[nfehn]

@jh66637 What I somewhat miss in this discussion is that you did not comment on the arguments provided here: https://github.com/exadg/exadg/pull/587/files#r1389376858, https://github.com/exadg/exadg/pull/587/files#r1389390122, #597 (comment).

Last week, I agreed to your proposed template version, but now I think it seems to have shortcomings (see links above), which is why I came up with the present PR and with the goal to have a discussion about these shortcomings.

As a sidenote/comment and independently of how we proceed with this PR: It is only partly true that "the code that calculates the face integrals is only written once". Looking at

exadg/include/exadg/acoustic_conservation_equations/spatial_discretization/operators/operator.h

Lines 338 to 378 in 51a6a7e

	if(data.formulation == Formulation::Weak)
	{
	scalar const flux_momentum_weak = rhocc * flux_momentum * n;
	vector const flux_mass_weak = rho_inv * flux_mass * n;
	pressure_m.submit_value(flux_momentum_weak, q);
	velocity_m.submit_value(flux_mass_weak, q);
	if constexpr(weight_neighbor)
	{
	// minus signs since n⁺ = - n⁻
	pressure_p.submit_value(-flux_momentum_weak, q);
	velocity_p.submit_value(-flux_mass_weak, q);
	}
	}
	else if(data.formulation == Formulation::Strong)
	{
	pressure_m.submit_value(rhocc * (flux_momentum - um) * n, q);
	velocity_m.submit_value(rho_inv * (flux_mass - pm) * n, q);
	if constexpr(weight_neighbor)
	{
	// minus signs since n⁺ = - n⁻
	pressure_p.submit_value(-rhocc * (flux_momentum - up) * n, q);
	velocity_p.submit_value(-rho_inv * (flux_mass - pp) * n, q);
	}
	}
	else if(data.formulation == Formulation::SkewSymmetric)
	{
	scalar const flux_momentum_weak = rhocc * flux_momentum * n;
	pressure_m.submit_value(flux_momentum_weak, q);
	velocity_m.submit_value(rho_inv * (flux_mass - pm) * n, q);
	if constexpr(weight_neighbor)
	{
	// minus signs since n⁺ = - n⁻
	pressure_p.submit_value(-flux_momentum_weak, q);
	velocity_p.submit_value(-rho_inv * (flux_mass - pp) * n, q);
	}
	}

one can see that the logic

   pressure_m.submit_value(..., q); 
   velocity_m.submit_value(..., q); 
  
   if constexpr(weight_neighbor) 
   { 
     pressure_p.submit_value(..., q); 
     velocity_p.submit_value(..., q); 
   } 

is implemented three times. Some implementations in ExaDG address this by computing tuples of fluxes, see

exadg/include/exadg/convection_diffusion/spatial_discretization/operators/convective_operator.h

Line 413 in 51a6a7e

calculate_flux_interior_and_neighbor(unsigned int const q,

. So far, I did not address this aspect for the acoustics module, but the version proposed in this PR goes into this direction. It goes into the direction of writing a function that receives values (and gradients) and returns fluxes. The advantage of such an implementation of face integrals is that it is independent of the data types that call get_value() and those that call submit_value(). One could argue such an implementation needs to implement the face integral only once as well. Of course, the get_value() part is written more than once for such an implementation, but this is exactly the part that is different for different "FEEvaluation types". Implementing classes like BoundaryFaceIntegratorP/U (as a consequence of the template design) can be seen as another way of implementing the get_value() part more than once.

[nfehn]

Let's set this PR (or, better, discussion) "on hold".

[jh66637]

@nfehn I hope this answers all the questions:

To the question of how to deal with interior values. As already suggested, we could have an additional class (similar to the class that provides the boundary values).

I don't think the duplication of submit_value() is the key problem. If we would have only one formulation there is no duplication. As you already suggested we could eliminate this by a function that returns a std::tuple.

While looking at the ConvDiff module I think I get the problem and understand your design suggestion. The problem we face is that OperatorBase is purely based on virtual functions.

I am hereby proposing an alternative solution that tries to fit in with the current OperatorBase. We can still keep the template parameters in the actual kernel and extend OperatorBase. Let me know if this would be a possible solution we can compromise on:

struct Test
{
  int get_value(){return 1.0;}
};

struct TestExterior
{
  int get_value(){return 1.0;}
};

struct OperatorBase
{
  // use different overloads to get the correct function call
  // keep all the virtual interface functions to be able not to change a lot in OperatorBase
  virtual void
  do_face_loop_internal(Test class1, Test class2);

  virtual void
  do_face_loop_boundary(Test class1, TestExterior class2);

};

// intermediate class only used to call the correct version of the templated function
// It would be possible to get rid of this class. However, this would require that OperatorBase 
// has to be completely in the header or that each operator type has to be registered manually 
// in the source file. I think there are ways to keep this registration short but I am not sure atm. 
// If you want to I can look it up.
template<typename Operator>
struct OperatorCRTP : OperatorBase
{
  // only needed if do_face_integral can not be implemented as static function
  auto
  this_op() 
  {
    return *static_cast<Operator*>(this);
  }

  // call the function in the derived class with the correct template parameters
  void
  do_face_loop_internal(Test class1, Test class2) final
  {
    this_op().template do_face_integral<true>(class1,class2);
  }
  void
  do_face_loop_boundary(Test class1, TestExterior class2) final
  {
    this_op().template do_face_integral<false>(class1, class2);
  }
};


struct ConvDiffOp : public OperatorCRTP<ConvDiffOp>
{
  // implement the logic what is happening on a face / boundary face / internal face / external face once
  template<bool weight_neighbor, typename Class1, typename Class2>
  void
  do_face_integral(Class1 class1, Class2 class2)
  {
    //...
  }
};

A short sidenote: I have been working on the non-matching implementation during the last few days and I think it is safe to say, that a function that automatically does the correct thing needs

template<InterorIntegrator, ExteriorIntegrator>
do_face_int(InteriorIntegrator int_m, ExteriorIntegrator int_p)
{
}

Where InteriorIntegrator can be of type FEFaceEval, FEPointEval and ExteriorIntegrator can be of type FEFaceEval BoundaryFaceEval and (possibly two different types) of FERemoteEval. So I think we can not avoid templates if we want to have one function that implements the logic that is happening at any face.

[nfehn]

In commit 41d3a33, I removed the ExteriorBoundaryIntegrator classes (which are currently not needed) in order to have a more direct comparison between (i) the template-version and (ii) the standard version of ExaDG (currently).

To make it simple, I suggest to discuss certain things independently from OperatorBase, because the acoustics operator does currently not derive from OperatorBase. Let me summarize the main properties of the two different versions from my perspective:

Properties of current standard version:

• need to write another face/boundary-face loop for the non-matching case
• no code duplication of physics since a kernel function is called that receives values and computes fluxes
• code is more explicit and easy to understand (subjective, of course)

Properties of template version:

• need to implement several classes (e.g. providing get_value()) by following the interface defined by the templated function
• no code duplication of physics since a templated function is called that operates on "Integrator" objects that follow a common interface like get_value()
• definition of interfaces via templates is more implicit and more difficult to understand in my opinion (again, subjective, of course)

In my opinion, it is far from obvious that the advantages of the templated version really outperform its disadvantages (as compared to the standard version). The potential increase in compile times for the templated version is not my main concern or the only concern actually (in particular not for the explicit acoustic operator addressed currently). (However, I would indeed like to avoid a pure header implementation of the base class OperatorBase or instantiating for every single operator for the reason of compile times.) My main concern is that the template version complicates the code a lot, and in the end does actually not reduce the amount of code that has to be written. Also, I did not fully understand so far how it is actually possible to use the same face/boundary-face loop implementation for the matching and the non-matching case when using the templated version. What I can offer is that I make again a one-to-one comparison between the templated version and the standard version once we have the non-matching case for acoustics. Then, the decision is made for the more simple solution.

In the end and when looking just at the code in ExaDG, one can probably argue for one or the other version (as the two of us currently do ;)) as long as we are not concrete/explicit about the non-matching case. However, we should also incorporate other aspects. In my opinion, the present discussion raised different questions:

1. Why is the Integrator/FEFaceEvaluation/FEEvaluation for the get()-part (solution function) and the submit()-part (test function) the same object in a "standard" dealii::matrixfree implementation?
2. In application code writing physics: Do we want the "kernel" to include the loop over quadrature points or just the operation on a single q-point?
3. Based on / related to point 2.: should the kernel function operate on (vectorized) dealii::Tensor data objects or on Integrator/FEFaceEvaluation/FEEvaluation objects.

[jh66637]

Sorry for the late reply. I finished the tutorial on how to deal with non-matching interfaces dealii/dealii#16299. With this, we have a better base for the discussion in the non-matching point.

I think we should discuss this with OperatorBase and non-matching in mind since we should aim for a unified procedure in designing the face loop in ÈxaDG. I am quite happy with the aproach in the tutorial since the Integrator objects can have a lot of types and the face_kernel has to be written only once. Note, that the functions face_kernel and inhomogenous_face_kernel from the tutorial could be easily unified.