Refactoring interface with contact relations

[CusiniM]

Goal

I am trying to think about how to restructure the interface with the contact relations and the permeability models to allow for the following:

• homogenize the various contact mechanics solvers interface with the contact models.
• use any contact model with any fracture solver. Currently, e.g., the Coulomb contact model can only be used in the LagrangianContactSolver.
• allow for a more flexible interface to permeability update in faults and fractures. The current design basically requires the definition of the hydraulicAperture which kind of assumes a lubrication theory type of model. This works well for open fractures but makes it hard to use models like those that are being implemented in [Feature] implement fault permeability models #1564 by @jhuang2601.

Proposal

I am thinking of doing the following things:

• define a field dispJump in all contact solvers (2 out of 3 already explicitly have one). This will always contain the displacement jump vector in a basis formed by the surface normal and two tangential vectors.
• define a field traction in all contact mechanics solvers written in the same basis as the dsipJump.
• define a member variable, arrayView2d<real64> m_traction, in ContactBase (a view to the traction field registered on the subRegion by the solvers). I think it would be nice for ContactBase to actually own the data but it may be difficult because the LagrangianSolver uses the traction as a primary variable and I am not sure that the DofManager will work properly if we do that (e.g., I am thinking of methods like addField() ). @klevzoff , @Andrea, let me know if I am wrong about this point.
• define the field hydraulicAperture in ContactBase (with proper getters if needed).
• remove the field hydraulicAperture registered by the flow solver and only have it registered by poromechanics solvers with fractures.
• remove array1d< real64 > m_elementAperture from the SurfaceElementSubRegion since it has more to do with the physics than with the geometry and it carries the same info as the dispJump. This way the subRegion only has geometrical information.
• Get rid of the following calls in FlowSolverBase: https://github.com/GEOSX/GEOSX/blob/449701f1aaf7d0a2c235e2060760ab0c4fca3340/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp#L53 and https://github.com/GEOSX/GEOSX/blob/449701f1aaf7d0a2c235e2060760ab0c4fca3340/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp#L272. The permeability in FlowSolverBase can be assumed to only be a function of flow variables so we can get rid of this dependency of the aperture.
• PoromechanicsSolvers with fractures (i.e., HydrofractureSolver, SinglePhasePoromechanicsSolverEmbeddedFractures) have and updateState method in which they launch a kernel to update the fracture permeability and the hydraulic aperture and compute derivatives w.r.t. to the dispJump. The HydrofracSolver can then chain rule to get derivatives w.r.t. the nodal displacements.
• Define a new interface in CoupledSolid and in the permeability models:

virtual void update(  localIndex const k,  localIndex const q,
                                 arrayView2d< real64 const > const & dispJump,
                                 arrayView2d< real64 const > const & traction, 
                                 arrayView1d<real64 const> const & hydraulicAperture ) const
{
// each permeability model does whatever it needs.
}

@andrea-franceschini @rrsettgast @joshua-white @castelletto1 @klevzoff @francoishamon @TotoGaz I don't think the solution I m proposing is perfect yet. Thoughts? Suggestions for improvement?

[klevzoff]

• define a member variable, arrayView2d<real64> m_traction, in ContactBase (a view to the traction field registered on the subRegion by the solvers). I think it would be nice for ContactBase to actually own the data but it may be difficult because the LagrangianSolver uses the traction as a primary variable and I am not sure that the DofManager will work properly if we do that (e.g., I am think of methods like addField() ). @klevzoff , @Andrea, let me know if I am wrong about this point.

addField() would be fine, I'm assuming you're referring to things like addVectorToField(), which indeed might be impacted in two ways:

• You'll have to pass a destination field name as e.g. "fractureContact_traction" where "fractureContact" is the name of the constitutive model, to allow DofManager to find the field by it's subregion wrapper name
• The field cannot be higher than 2D since DofManager doesn't currently know how to map DoF component indices to multidimensional slice indices, although I think this can be fixed. Assuming the traction field doesn't have quadrature point dimension, this is not a problem.

[TotoGaz]

It would surely be interesting to have some addField or addVectorField... if we have the Field/ScalarField/VectorField concepts too! 🎉

Any options?

[rrsettgast]

• allow for a more flexible interface to permeability update in faults and fractures. The current design basically requires the definition of the hydraulicAperture which kind of assumes a lubrication theory type of model. This works well for open fractures but makes it hard to use models like those that are being implemented in [Feature] implement fault permeability models #1564 by @jhuang2601.

The hydraulic aperture should be a function of what? Normal stress? Some dilation slip state variable?

Proposal

I am thinking of doing the following things:

• define a field dispJump in all contact solvers (2 out of 3 already explicitly have one). This will always contain the displacement jump vector in a basis formed by the surface normal and two tangential vectors.

I would argue that the dispJump is derived from the primary displacement variables. Also, if you wanted to implement a Nitche's method, then you would want a full integration on contact surface.

• define a member variable, arrayView2d<real64> m_traction, in ContactBase (a view to the traction field registered on the subRegion by the solvers). I think it would be nice for ContactBase to actually own the data but it may be difficult because the LagrangianSolver uses the traction as a primary variable and I am not sure that the DofManager will work properly if we do that (e.g., I am thinking of methods like addField() ). @klevzoff , @Andrea, let me know if I am wrong about this point.

As long as it is proper state variable. What is the state variable for contact constitutive relations? I would think it would be some sort of plastic slip?

• define the field hydraulicAperture in ContactBase (with proper getters if needed).

I don't think we would want the concept of hydraulicAperture in the contact model....not that I have a better suggestion.

• remove array1d< real64 > m_elementAperture from the SurfaceElementSubRegion since it has more to do with the physics than with the geometry and it carries the same info as the dispJump. This way the subRegion only has geometrical information.

The aperture does define the volume of the surface element in some respects. But the jump would be more descriptive

[CusiniM]

• allow for a more flexible interface to permeability update in faults and fractures. The current design basically requires the definition of the hydraulicAperture which kind of assumes a lubrication theory type of model. This works well for open fractures but makes it hard to use models like those that are being implemented in [Feature] implement fault permeability models #1564 by @jhuang2601.

The hydraulic aperture should be a function of what? Normal stress? Some dilation slip state variable?

So, I was more referring to the fact that currently we only have an interface to compute the permeability as a function of the hydraulicAperture using the cubic law and that some of the models we are now implementing compute the permeability based on something else (e.g., the slip magnitude). As for the hydraulicAperture itself I think we can make it a generic function of the dispJump?

Proposal

I am thinking of doing the following things:

• define a field dispJump in all contact solvers (2 out of 3 already explicitly have one). This will always contain the displacement jump vector in a basis formed by the surface normal and two tangential vectors.

I would argue that the dispJump is derived from the primary displacement variables. Also, if you wanted to implement a Nitche's method, then you would want a full integration on contact surface.

I think that depends on the formulation but irrespective of how the dispJump is computed (whether it is the main unknown or it is computed from the nodal displacements) I feel that it is very convenient to express hydromechanical properties (mainly thinking of fracture/fault permeability and storativity) as functions of it. At least for the cases I have in mind. Then, we can recover the other dependencies in the physics solver when assembling derivatives.

• define a member variable, arrayView2d<real64> m_traction, in ContactBase (a view to the traction field registered on the subRegion by the solvers). I think it would be nice for ContactBase to actually own the data but it may be difficult because the LagrangianSolver uses the traction as a primary variable and I am not sure that the DofManager will work properly if we do that (e.g., I am thinking of methods like addField() ). @klevzoff , @Andrea, let me know if I am wrong about this point.

As long as it is proper state variable. What is the state variable for contact constitutive relations? I would think it would be some sort of plastic slip?

• define the field hydraulicAperture in ContactBase (with proper getters if needed).

I don't think we would want the concept of hydraulicAperture in the contact model....not that I have a better suggestion.

Yeah, couldn't really think of a better location... Unless it's a field on the SurfaceElementSubRegion but it will have to be registered by each solver dealing with fractures because we don't have a FracturePoromechanicsSolverBase (maybe we should at some point).

• remove array1d< real64 > m_elementAperture from the SurfaceElementSubRegion since it has more to do with the physics than with the geometry and it carries the same info as the dispJump. This way the subRegion only has geometrical information.

The aperture does define the volume of the surface element in some respects. But the jump would be more descriptive

It's just that in my mind these are geometrically surfaces and so they should not even have a volume. The volume is more of a modelling choice because it depends on how we choose to compute the hydraulic aperture since that's the volume that really matters. This point is honestly not so important and there are arguments in favour of keeping it as it is because it's convenient for a number of other things.

[jhuang2601]

• allow for a more flexible interface to permeability update in faults and fractures. The current design basically requires the definition of the hydraulicAperture which kind of assumes a lubrication theory type of model. This works well for open fractures but makes it hard to use models like those that are being implemented in [Feature] implement fault permeability models #1564 by @jhuang2601.

The hydraulic aperture should be a function of what? Normal stress? Some dilation slip state variable?

For the planned implementation, the hydraulic aperture would be a function of

• relative shear displacement (slip) in the shear displacement dependent permeability model
• normal stress and slip in the Willis-Richards model (Barton-Bandis type correlation)
• some state variables/correlations (e.g., hydraulic aperture vs. mechanical opening) with the input table functions (more generalized model)

I agree with @CusiniM that the existing interface to compute the permeability in the PermeabilityBase.hpp should be updated to allow more flexibility.

Proposal

I am thinking of doing the following things:

• define a field dispJump in all contact solvers (2 out of 3 already explicitly have one). This will always contain the displacement jump vector in a basis formed by the surface normal and two tangential vectors.

I would argue that the dispJump is derived from the primary displacement variables. Also, if you wanted to implement a Nitche's method, then you would want a full integration on contact surface.

• define a member variable, arrayView2d<real64> m_traction, in ContactBase (a view to the traction field registered on the subRegion by the solvers). I think it would be nice for ContactBase to actually own the data but it may be difficult because the LagrangianSolver uses the traction as a primary variable and I am not sure that the DofManager will work properly if we do that (e.g., I am thinking of methods like addField() ). @klevzoff , @Andrea, let me know if I am wrong about this point.

As long as it is proper state variable. What is the state variable for contact constitutive relations? I would think it would be some sort of plastic slip?

As of now, only the elastic slip (Coulomb failure criterion) is considered and implemented for the Lagrange contact model. You're proposing to develop plastic slip model for contact region?

• define the field hydraulicAperture in ContactBase (with proper getters if needed).

I don't think we would want the concept of hydraulicAperture in the contact model....not that I have a better suggestion.

• remove array1d< real64 > m_elementAperture from the SurfaceElementSubRegion since it has more to do with the physics than with the geometry and it carries the same info as the dispJump. This way the subRegion only has geometrical information.

The aperture does define the volume of the surface element in some respects. But the jump would be more descriptive