PDE IntegrandLib

[TTitscher]

We want to provide some default integrands to the user. They should probably all have a similar design regarding their API or their history data.

My first idea: Provide a combination of static functions

static Vector IntegrandInLib::Gradient(CellData, CellIpData, HistoryData)
{    // calculate gradient using the history data }
static HistoryData IntegrandInLib::History(CellData, CellIpData, OldHistoryData)
{   // calculate history data }

and template helpers:

//! @brief only provides a gradient method that automatically updates the history data upon evaluation
template <typename TStaticIntegrand>
struct ExplicitIntegrand
{
   std::vector<TSTaticItnegrand::HistoryDataType> mData;
   Vector Gradient(...)
  {
    auto history = TStaticIntegrand::History(...);
    mData[...] = history;
    return TStaticIntegrand::Gradient(..., history);
  }
};

The benefit: The physics/math is defined in the static functions. You can use them to build your own integrands. For the default cases (explicit and implicit), there will be templates as the one above to do that for you. But you are also free to combine multiple of these static functions into a mixed integrand.

Do you like that approach? For the momentum balance equation, we can even define constitutive laws that provide the stress-strain relationship so we can use multiple laws with the same integrand. How will this fit in here?

[joergfunger]

In my opinion, the "old" implementation in NuTo for explicit time integration with history variables is not fully correct (Don't think that anyone has used it so far). What we currently (in the old NuTo) do in the explicit Runge-Kutta steps is an update while evaluating the function. I think that this is not correct. We should calculate the final solution as a linear combination of the intermediate steps (each step integrates its history data from the previous converged time step and calculates the residual), and only then we can perform an update. This means that we should not do update and stress calculation simultaneously.

Furthermore, I also thought about decoupling history data evaluation from gradient calculation. However, I think for many problems (e.g. plasticity) this is not optimal, since the return mapping algorithm directly provides the stresses required to calculate the gradient. Furthermore (and that is even more important), the algorithmic stiffness is not only a function of the current plastic strain, but it requires also the plastic multipliers (gamma), which is usually not part of the history data. Thus, we would have to store these multipliers as well.

As for the handling of the history data, I would prefer the solution from @vhirtham in the current pull request class CreepLaw : public Laws::MechanicsInterface<1>, public HistoryDataContiguousMemory<CreepHistoryData> where the user can just derive its integrand from a templated history data class that handles the history data.

[vhirtham]

Don't know much (nothing) about the requirements for history data of explicit schemes, but I think the history data should get its own class. The additional functionality is only needed by some laws and only they should provide this interface IMHO. If we use inheritance to plug the history data into the integrand and all our default integrands use the default history data, we also automatically get the same interfaces. If inheritance is too much of a pain for us, we can also just use the history data class as a member, but then we have to provide pass through functions (for example for initialization) with a standard interface. Personally I don't think that inheritance is so super bad, especially if it is just of dept 1.

Additionally I think we should agree on the interfaces of an integrand first. Then we can discuss details of the implementation. What do we need? Gradient, Jacobian_dt 0-X (previously know as Hessian :stuck_out_tongue: ) and the history data stuff. What else? How do we handle the fact that not every law needs a Jacobian1 or higher? Do they have to implement a variant that throws an exception or can I simply just skip the implementation if the law does not need it? Do I need an info routine how many time derivatives my law can handle or for other stuff? What about members - Are they always set in the constructor and can't be accessed anymore or do we need getters (setters)?

Addressing your suggestions: Well this split is an interesting approach, but I think that is more an implementation detail to reduce code for a specific group of integrands that have a lot in common. I would say, the basic requirements for our integrand lib should be a common and intuitive interface. The implementation details can be discussed in a pull request. The "user" don't care about them. He only cares about a comfortable interface.

[TTitscher]

Additionally I think we should agree on the interfaces of an integrand first.

The result of about one year of discussion was to not have a common interface. :) Because if we omit this common interface, all your other questions are not relevant anymore. Big advantage.

I might have opened this issue too early. New suggestion: How about we delay this discussion until we have some integrands to talk about. Like 4 or 5. And then we try to couple them. We may find things to refactor. And we may find patterns that makes using those integrands feely simple/hard.

[vhirtham]

It was not about a general interface for an integrand, just for the integrand lib. As you stated:

They should probably all have a similar design regarding their API or their history data.