On-the-fly material property evaluation

[lindsayad]

Reason

We have one place in MOOSE currently where we actually have wrong physics because of the limitations of our material system. It's in INSFV where we compute coefficients for Rhie-Chow interpolations for the current elem and neighbor of a FaceInfo object. For either elem or neighbor, we have to loop over all bounding FaceInfo objects in order to determine the elemental coefficient. In those loop computations we are supposed to use the viscosity evaluated at the face, but we do not have access to viscosity values at arbitrary points in space because of the material system design, in which we prequest references to material properties. So for now we are using the elemental viscosity all the time. As soon as start to use spatially varying viscosities, our current design will be wrong.

Design

Either a complete overhaul of the current material system or addition of a new parallel system would have to take place in order to allow evaluation of material properties at arbitrary points. It's hard to say what the design should look like exactly but the API should probably be something like:

ADReal getADMaterialProperty(const std::string & name, const & Point coord, const Real time);

Impact

Allow physics developers more flexibility (to be correct).

[dschwen]

This would be neat for line material sampler. I'm assuming stateful props are a no go here?

[lindsayad]

We could theoretically build a cache of material property evaluations and then shift then like we do for current material properties. This would require that you be performing evaluations at the same point in space at each time step.

[GiudGiud]

Adding to this issue, objects that inherit from (1)MaterialPropertyInterface could use the ability to compute material properties on a face. This can be done by interpolating from the neighbor value, or by calculating the material property using the variable values on the face.

This would allow the VolumetricFlowRate postprocessors in the navier stokes module to compute the proper inlet flow of energy (material property).

[lindsayad]

How about this kind of design? @rwcarlsen

#include <iostream>
#include <map>
#include <string>
#include <functional>

class Elem
{
};
class FaceInfo
{
};
class Variable
{
  double _x;

public:
  Variable(const double x) : _x(x) {}
  double operator()(const Elem & elem) const { return _x; }
  double operator()(const FaceInfo & fi) const { return _x; }
  double operator()(const unsigned int qp) const { return _x; }
};

class MaterialProperty
{
public:
  using ElemFn = std::function<double(const Elem &)>;
  using FiFn = std::function<double(const FaceInfo &)>;
  using QpFn = std::function<double(const unsigned int &)>;

  template <typename PolymorphicLambda>
  MaterialProperty & operator=(PolymorphicLambda my_lammy)
  {
    elem_functor = my_lammy;
    fi_functor = my_lammy;
    qp_functor = my_lammy;
    return *this;
  }

  double operator()(const Elem & elem) const { return elem_functor(elem); }

  double operator()(const FaceInfo & fi) const { return fi_functor(fi); }

  double operator()(const unsigned int qp) const { return qp_functor(qp); }

private:
  ElemFn elem_functor;
  FiFn fi_functor;
  QpFn qp_functor;
};

class Material
{
public:
  Material() : _u(2)
  {
    _mat_props["u2"] = [this](const auto & geom_entity) -> double {
      return _u(geom_entity) * _u(geom_entity);
    };
    _mat_props["u3"] = [this](const auto & geom_entity) -> double {
      return _u(geom_entity) * _u(geom_entity) * _u(geom_entity);
    };
  }

  const MaterialProperty & getMatProp(const std::string & name) const
  {
    return _mat_props.at(name);
  }

private:
  std::map<std::string, MaterialProperty> _mat_props;
  Variable _u;
};

class ResidualObject
{
public:
  void assignMatProp(const MaterialProperty & prop) { _mat_prop = &prop; }

protected:
  const MaterialProperty * _mat_prop;
};

class Kernel : public ResidualObject
{
public:
  void computeResidual(const Elem & elem) { std::cout << (*_mat_prop)(elem) << std::endl; }
};

class FVFluxKernel : public ResidualObject
{
public:
  void computeResidual(const FaceInfo & fi) { std::cout << (*_mat_prop)(fi) << std::endl; }
};

int
main()
{
  Material mat;
  Kernel kern;
  kern.assignMatProp(mat.getMatProp("u2"));
  FVFluxKernel fv;
  fv.assignMatProp(mat.getMatProp("u3"));

  kern.computeResidual(Elem());
  fv.computeResidual(FaceInfo());
}

Outputs:

(moose-libmesh) lindad@localhost:~/programming/cpp/functor-materials(master)$ ./a.out
4
8

[GiudGiud]

I like it.
Having face values for materials can reduce our use of neighborMaterialProperties for interpolation purposes too.
Would there be a way to handle gradients of material properties? Thinking of porosity.

Would you just skip computing material properties during reinit then?
Or compute and store at reinit and fetch when computing residuals?

This would be a new material property interface right? Like we would change the interface the kernels inherit

[dschwen]

Gradients would still require shape functions, which materials don't have. You could compute a gradient property by the material. Would that help?

The design is slick, but we'd replacing accessing a reference to a fixed bit of memory with a virtual dispatch and a function call.

[rwcarlsen]

@lindsayad - not bad. Definitely would be an improvement. @dschwen - So what you're saying is we'd be replacing a poor design with a better one :P

[lindsayad]

The design is slick, but we'd replacing accessing a reference to a fixed bit of memory with a virtual dispatch and a function call.

I don't envision my design replacing the current material system. I think that for CFEM, what we have now is pretty optimal. However, what we have now is really bad for finite volume calculations where each degree of freedom can involve a pretty large stencil of elements' degrees of freedom. I would like to name this new material/property something like FunctorMaterial/FunctorMaterialProperty and it could work for finite element calculations, and even be essential if a finite element calculation involves some non-local information. But in practice these new objects could be probably be just as well named FVMaterial/FVMaterialProperty as that would indicate its primary use case.

[rwcarlsen]

My concern is the kitchen sink problem. If we just keep adding features to scratch one or two itches everytime a need comes up, and don't figure out how to integrate/evolve moose cohesively - we will just end up with SCALE for multi-physics - which is about as far from complementary as code-metaphors can get.

[lindsayad]

But people often complain about MOOSE being slow. The more general we are, usually the slower we will be. I understand your concern, but when a lot of milestones are about optimization and I imagine we will continue to have milestones around optimization, then I think it makes sense, although it may be distasteful, to have parallel systems. If our users were to say that this kind of FunctorMaterial design/concept was good enough/fast enough for them, then I would be all for removing our current index-into-precomputed-reference-to-array design.

[lindsayad]

@idaholab/moose-ccb should FunctorMaterial replace/become Material or should the two systems live side-by-side? The former would require a lot of up front work (for which perhaps we do not have appropriate funding right now) followed by less maintenance. The latter would require about a week of up front work followed by more maintenance. I like the latter. I think it's fair for me to say that @rwcarlsen likes the former. What are other opinions?

[rwcarlsen]

I think we should make sure we are addressing as many use-cases as possible. A while back I vaguely remember Ben was saying they had some funky material prop location needs/desires for xfem stuff.

[lindsayad]

It seems like only @rwcarlsen and I care about this 😄

[loganharbour]

I know that proper material caching was also on the roadmap at some point in time... how would that play in to our choice here?

[WilkAndy]

We have "nodal Materials" in PorousFlow that evaluate Material Properties at the nodes of the element. This is to enable numerical stabilization. We store the information inside the current Material datastructures, which is painful because we have to be careful when the number of quadpoints != the number of nodes. In most cases it's also inefficient because for a single node we perform the same calculation number_of_connected_elements times when we're calculating those nodal Material properties (only in one case - porosity - is the nodal value different depending on which element you're in).

We use DerivativeMaterial (not AD)

We also have stateful.

We have and "nodal Material" props evaluated at faces for BCs

We also have "nodal Material" evaluated at DiracKernel points (which has never properly worked, i believe, #10471).

[lindsayad]

I know that proper material caching was also on the roadmap at some point in time... how would that play in to our choice here?

Caching could cut-down CPU time but would increase memory consumption. I think would have to profile to understand which resource is affected most critically. If the CPU decreases dramatically with a corresponding dramatic increase in memory, then perhaps you could create a configure option or even a run-time option that enables/disables caching depending on what you care about.

[rwcarlsen]

We could even have a fixed/automatically sized cache that drops/retains as necessary to keep an appropriate memory footprint.

[lindsayad]

I believe this could be helpful for computing quantities of interest like energy dissipation from friction (based off a conversation with @recuero )