WIP - Constitutive Model Reorg for Poro-Solid-Damage

[rrsettgast]

WIP
This is a proposal to reorganize the constitutive models. Specifically the definition of the various permutations of solid/poro/damage model implementations.

Propose a porosity data class (artist also known as RockBase:

class PoroNull;

class PoroData
{
   array2d< real64 > m_newPorosity;
   array2d< real64 > m_oldPorosity;
   array2d< real64 > m_dPorosity_dPressure;
   array1d< real64 > m_referencePorosity;
   real64 m_defaultReferencePorosity;
   real64 m_compressibility;
   real64 m_grainBulkModulus;
   array2d< real64 > m_grainDensity;
   real64 m_defaultGrainDensity = 0;
};

Propose a family of stand-alone solid material and associated update classes...for example

class ElasticIsotropicUpdate : public SolidBaseUpdate
{
   void smallStrainUpdate(...)
}
class ElasticIsotropic : public SolidBase
{
  real64 m_defaultBulkModulus;
  real64 m_defaultShearModulus;
  array1d< real64 > m_bulkModulus;
  array1d< real64 > m_shearModulus;
};

Propose a family of damage models:

class DamageNull;
class DamageBase;
class DamageScalar : public DamageBase;
class DamageSpectral : public DamageBase;
class DamageVolDevSplit : public DamageBase;

Propose a set of aggregate models using composition for the stand alone models:

template< typename SOLID,
          typename DAMAGE = DamageNull, 
          typename PORO = PoroNull >
class MuchoUpdate
{
   void smallStrainUpdate(...)
...
SOLID::update m_solidUpdate;
};



template< typename SOLID,
          typename DAMAGE = DamageNull, 
          typename PORO = PoroNull >
class MuchoModel
{
SOLID m_solid;
DAMAGE m_damage;
PORO m_Poro;
}

Perhaps also the time to start thinking about better alternatives to the update classes.

[CusiniM]

I think both the multiple inheritance and the composition options should give us what we need.

At the end of the day the things to consider are the following:

• the damage model will always need a solid model for it to make sense.

• we can think of a "special" poroModel which is used as a standalone and employs a trivial relationship to update porosity. This is going to be useful to run flow only cases and compare results with existing codes. This model should basically only have a compressibility and a porosity.

• we will want to use solid models by themselves (I guess this is obvious).

• the way certain properties are computed will depend on the type of the solidmodel. E.g., not all solid models compute the biotCoefficient in the same way and even the type of the biotCoefficient can be different (tensor vs a scalar). Thus, we will either need to instruct the solid models about how to compute the biot coefficient or, we will have to define different functions to do so based on the type of solid model. I think the damage model would then simply modify the biotCoefficient based on the damage.

[rrsettgast]

There would be some "pollution" of solid models if we choose to embed (for example) the computation of the biotCoefficient in the solid model itself. On the other hand, if we wanted to keep it clean, then we would potentially have some intermediate helper classes to perform the computation...or embed it a specific instantiation of the MuchoModel.

[joshua-white]

The biotCoefficient can be computed generically with minimal info from the solid and poro models. We need a public method K_g = m_poro.getGrainCompressibility() from poro model, and a method K=m_solid.getBulkModulus() from the solid model. Together, they define a scalar biotCoefficient as b=1-K/K_g. If we have a method getElasticStiffness(), we can generically compute a Biot tensor as \vec{b}=\vec{1}-vec{1}:\mathbb{C}^3 / (3 K_g). See, e.g., equation 30 in https://doi.org/10.1016/j.cma.2020.113225.

[CusiniM]

Ok, this sounds good.
The only difficulty is the idea of updating all properties at the same time we assemble the system. It's going to be hard for the damageModel to modify the biotCoefficient as it needs to if we do this because the DamageSolver is sequentially coupled to the PoroealsticSolver and reuses the same assembly kernel.. However, in my understanding, the damage model only needs to modify the biotCoefficient and not the porosity directly. Thus, I propose that we store the biotCoefficient in the porosityModel and update it before assemblying the system. At that point, the damage poroealstic Solver can call it's dedicated update of the biotCoefficient.

This only works if:

• the damage does not need to directly modify the porosity.
• the biotCoefficient can be computed only from the grainCompressibility and the bulkModulus and neither of them depends on the strain (or the stress) which will only be computed and updated inside the Poroelastic kernel.

@andrembcosta, @joshua-white, @castelletto1, @rrsettgast comments?

[CusiniM]

There is one more thing that I noticed and we may want to use this refactoring to improve it. Most ConstitutiveUpdates of constitutive models have a method called update() which is marked const. However, the method then calls another internal method called compute which actually modifies the value of the member variables (e.g, m_property and m_dProperty_dSomething.

This is done for the relative permeability (e.g., https://github.com/GEOSX/GEOSX/blob/develop/src/coreComponents/constitutive/relativePermeability/BrooksCoreyBakerRelativePermeability.hpp#L68-84 ) and I reproduced the same pattern for the porosity and the permeability.

To me it looks like that this way we are defining as const a method which, at the end of the day, is not really const.

[rrsettgast]

They have to be marked const because the Update object is captured by value, which triggers the move of the member views. If something is captured by value, it is const in the context of the the lambda. So it is just like the ArrayView. All the accessors are const, but the actual ability to modify data is defined by the type (Array< double const> or Array< double > ).

[TotoGaz]

The issue with

class MuchoModel : public SOLID, public DAMAGE, public PORO

is that it implies a is a relationship.

The question therefore being "Do I want my MuchoModel to actually be a SOLID, a DAMAGE, a PORO?".
I would assume that MuchoModel is the precise and controlled combination of those 3 instead.
Getting public inheritance let any user use the SOLID, DAMAGE and PORO any old way.
Resulting in

• more bugs
• more expensive entry ticket to GEOSX
• more complicated refactoring (broader footprint)

Composition could let you decide with the finest precision what you want to expose.

Another possibility could be to use private inheritance.
It maybe end-up with multiple inheritance issues, but you'll control your public exposure.

[rrsettgast]

@TotoGaz I agree. Composition is clearly better here.

[rrsettgast]

This work was started in #1455