[wip] refact: Fit ALM as LM in the Poromechanics hierarchy

inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ThermoPoroElastic_ALM_conforming_base.xml

@@ -0,0 +1,200 @@
+<?xml version="1.0" ?>
+
+<Problem>
+  <Included>
+    <File
+      name="../../../thermoPoromechanicsFractures/ThermoPoroElastic_base.xml"/>
+  </Included>
+
+  <Solvers
+    gravityVector="{ 0.0, 0.0, 0.0 }">
+    <SinglePhasePoromechanicsConformingFracturesALM
+      name="poroSolver"
+      targetRegions="{ RockMatrix, Fracture }"
+      initialDt="1e-3"
+      flowSolverName="flowSolver"
+      solidSolverName="fractureMechSolver"
+      logLevel="1"
+      isThermal="1">
+      <NonlinearSolverParameters
+        newtonTol="1.0e-5"
+        newtonMaxIter="20"
+        maxTimeStepCuts="4"/>
+      <LinearSolverParameters
+        directParallel="0"
+        solverType="direct"/>
+    </SinglePhasePoromechanicsConformingFracturesALM>
+
+    <SolidMechanicsAugmentedLagrangianContact
+      name="fractureMechSolver"
+      logLevel="1"
+      simultaneous="1"
+      symmetric="1"
+      iterPenaltyN="10.0"
+      iterPenaltyT="10.0"
+      tolJumpN="1.e-05"
+      tolJumpT="1.e-04"
+      tolTauLimit="1.e-02"
+      timeIntegrationOption="QuasiStatic"
+      discretization="FE1"
+      targetRegions="{ RockMatrix, Fracture }">
+      <NonlinearSolverParameters
+        newtonTol="1.0e-5"
+        logLevel="1"
+        maxNumConfigurationAttempts="50"
+        newtonMaxIter="20"
+        lineSearchAction="None"
+        lineSearchMaxCuts="4"/>
+      <LinearSolverParameters
+        directParallel="0"
+        solverType="direct"/>
+    </SolidMechanicsAugmentedLagrangianContact>
+
+    <SinglePhaseFVM
+      name="flowSolver"
+      discretization="singlePhaseTPFA"
+      targetRegions="{ RockMatrix, Fracture }"
+      isThermal="1">
+      <NonlinearSolverParameters
+        normType="L2"/>
+    </SinglePhaseFVM>
+
+    <SurfaceGenerator
+      logLevel="0"
+      name="SurfaceGenerator"
+      targetRegions="{ RockMatrix }"
+      fractureRegion="Fracture"
+      initialRockToughness="1e6"
+      mpiCommOrder="1"/>
+  </Solvers>
+
+  <NumericalMethods>
+    <FiniteElements>
+      <FiniteElementSpace
+        name="FE1"
+        order="1"/>
+    </FiniteElements>
+
+    <FiniteVolume>
+      <TwoPointFluxApproximation
+        name="singlePhaseTPFA"/>
+    </FiniteVolume>
+  </NumericalMethods>
+
+  <ElementRegions>
+    <CellElementRegion
+      name="RockMatrix"
+      cellBlocks="{ * }"
+      materialList="{ porousRock, water, thermalCond }"/>
+
+    <SurfaceElementRegion
+      name="Fracture"
+      materialList="{ water, fractureFilling, fractureContact, thermalCond, hApertureModel }"
+      defaultAperture="1.e-4"/>
+  </ElementRegions>
+
+  <FieldSpecifications>
+    <FieldSpecification
+      name="frac"
+      initialCondition="1"
+      setNames="{ fracture }"
+      objectPath="faceManager"
+      fieldName="ruptureState"
+      scale="1"/>
+
+    <FieldSpecification
+      name="separableFace"
+      initialCondition="1"
+      setNames="{ fracture }"
+      objectPath="faceManager"
+      fieldName="isFaceSeparable"
+      scale="1"/>
+
+    <FieldSpecification
+      name="fracPorosity"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/Fracture"
+      fieldName="fracturePorosity_porosity"
+      scale="1.00"/>
+
+    <FieldSpecification
+      name="initialPressure"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/RockMatrix/cb1"
+      fieldName="pressure"
+      scale="0.0"/>
+
+    <FieldSpecification
+      name="initialPressureFracture"
+      setNames="{ all }"
+      objectPath="ElementRegions/Fracture"
+      fieldName="pressure"
+      scale="1.0e5"/>
+
+    <FieldSpecification
+      name="initialTemperature"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/RockMatrix/cb1"
+      fieldName="temperature"
+      scale="0.0"/>
+
+    <FieldSpecification
+      name="initialTemperatureFracture"
+      setNames="{ all }"
+      objectPath="ElementRegions/Fracture"
+      fieldName="temperature"
+      scale="300"/>
+
+    <FieldSpecification
+      name="xnegconstraint"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="0"
+      scale="0.0"
+      setNames="{ xneg, xpos }"/>
+
+    <FieldSpecification
+      name="yposconstraint"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="1"
+      scale="0.0"
+      setNames="{ ypos, yneg }"/>
+
+    <FieldSpecification
+      name="zconstraint"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="2"
+      scale="0.0"
+      setNames="{ zneg, zpos }"/>
+
+    <FieldSpecification
+      name="facePressure"
+      objectPath="faceManager"
+      fieldName="pressure"
+      scale="0.0"
+      setNames="{ xneg, xpos }"/>
+
+    <FieldSpecification
+      name="faceTemperature"
+      objectPath="faceManager"
+      fieldName="temperature"
+      scale="0.0"
+      setNames="{ xneg, xpos }"/>
+  </FieldSpecifications>
+
+  <Outputs>
+    <VTK
+      name="vtkOutput"
+      plotFileRoot="ThermoPoroElastic_ALM_conforming"
+      plotLevel="2"
+      format="binary"/>
+
+    <Restart
+      name="restartOutput"/>
+  </Outputs>
+</Problem>

inputFiles/poromechanicsFractures/Contact/AugmentedLagrangianMultipliers/ThermoPoroElastic_ALM_conforming_smoke.xml

@@ -0,0 +1,54 @@
+<?xml version="1.0" ?>
+
+<Problem>
+  <Included>
+    <File
+      name="./ThermoPoroElastic_ALM_conforming_base.xml"/>
+  </Included>
+
+  <Mesh>
+    <InternalMesh
+      name="mesh1"
+      elementTypes="{ C3D8 }"
+      xCoords="{ 0, 20 }"
+      yCoords="{ 0, 20 }"
+      zCoords="{ 0, 1 }"
+      nx="{ 10 }"
+      ny="{ 11 }"
+      nz="{ 2 }"
+      cellBlockNames="{ cb1 }"/>
+  </Mesh>
+
+  <Geometry>
+    <Box
+      name="fracture"
+      xMin="{ 9.5, -0.01, -0.01 }"
+      xMax="{ 10.5, 20.01, 1.01 }"/>
+  </Geometry>
+
+  <Events
+    maxTime="1.0e6">
+    <SoloEvent
+      name="preFracture"
+      target="/Solvers/SurfaceGenerator"/>
+
+    <SoloEvent
+      name="preFracPlot"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="outputs"
+      cycleFrequency="1"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="solverApplication0"
+      forceDt="1.0e5"
+      target="/Solvers/poroSolver"/>
+
+    <PeriodicEvent
+      name="restarts"
+      timeFrequency="5.0e5"
+      target="/Outputs/restartOutput"/>
+  </Events>
+</Problem>

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp

@@ -176,6 +176,18 @@ class FlowSolverBase : public PhysicsSolverBase
     GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." );
   }
 
+  virtual void assembleHydrofracFluxTermsALM( real64 const time_n,
+                                              real64 const dt,
+                                              DomainPartition const & domain,
+                                              DofManager const & dofManager,
+                                              CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                              arrayView1d< real64 > const & localRhs,
+                                              CRSMatrixView< real64, localIndex const > const & dR_dAper )
+  {
+    GEOS_UNUSED_VAR ( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper );
+    GEOS_ERROR( "Poroelastic fluxes with conforming fractures ALM not yet implemented." );
+  }
+
   void initializeState( DomainPartition & domain );
 
   virtual void initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); }

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp

@@ -268,6 +268,12 @@ class SinglePhaseBase : public FlowSolverBase
   virtual void
   updateState ( DomainPartition & domain ) override;
 
+  /**
+   * @brief Getter for the number of fluid components (species)
+   * @return the number of components
+   */
+  integer numFluidComponents() const { return 1; }
+
   /**
    * @brief Function to update all constitutive state and dependent variables
    * @param subRegion subregion that contains the fields

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp

@@ -45,6 +45,8 @@
 #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp"
 #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp"
 #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp"
+#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp"
+#include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp"
 
 /**
  * @namespace the geos namespace that encapsulates the majority of the code
@@ -669,6 +671,92 @@ void SinglePhaseFVM< BASE >::assembleHydrofracFluxTerms( real64 const GEOS_UNUSE
 
 }
 
+template< typename BASE >
+void SinglePhaseFVM< BASE >::assembleHydrofracFluxTermsALM( real64 const GEOS_UNUSED_PARAM ( time_n ),
+                                                            real64 const dt,
+                                                            DomainPartition const & domain,
+                                                            DofManager const & dofManager,
+                                                            CRSMatrixView< real64, globalIndex const > const & localMatrix,
+                                                            arrayView1d< real64 > const & localRhs,
+                                                            CRSMatrixView< real64, localIndex const > const & dR_dAper )
+{
+  GEOS_MARK_FUNCTION;
+
+  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
+  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
+  FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName );
+
+  string const & dofKey = dofManager.getKey( SinglePhaseBase::viewKeyStruct::elemDofFieldString() );
+
+
+  this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
+                                                                      MeshLevel const & mesh,
+                                                                      string_array const & )
+  {
+    fluxApprox.forStencils< CellElementStencilTPFA, FaceElementToCellStencil >( mesh, [&]( auto & stencil )
+    {
+      typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper();
+
+      if( m_isThermal )
+      {
+        thermalSinglePhaseFVMKernels::
+          FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(),
+                                                                               dofKey,
+                                                                               this->getName(),
+                                                                               mesh.getElemManager(),
+                                                                               stencilWrapper,
+                                                                               dt,
+                                                                               localMatrix.toViewConstSizes(),
+                                                                               localRhs.toView() );
+      }
+      else
+      {
+        singlePhaseFVMKernels::
+          FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(),
+                                                                               dofKey,
+                                                                               this->getName(),
+                                                                               mesh.getElemManager(),
+                                                                               stencilWrapper,
+                                                                               dt,
+                                                                               localMatrix.toViewConstSizes(),
+                                                                               localRhs.toView() );
+      }
+    } );
+
+    fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( auto & stencil )
+    {
+      typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper();
+
+      if( m_isThermal )
+      {
+        thermalSinglePhasePoromechanicsConformingFracturesALMKernels::
+          ConnectorBasedAssemblyKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(),
+                                                                                          dofKey,
+                                                                                          this->getName(),
+                                                                                          mesh.getElemManager(),
+                                                                                          stencilWrapper,
+                                                                                          dt,
+                                                                                          localMatrix.toViewConstSizes(),
+                                                                                          localRhs.toView(),
+                                                                                          dR_dAper );
+      }
+      else
+      {
+        singlePhasePoromechanicsConformingFracturesALMKernels::
+          ConnectorBasedAssemblyKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(),
+                                                                                          dofKey,
+                                                                                          this->getName(),
+                                                                                          mesh.getElemManager(),
+                                                                                          stencilWrapper,
+                                                                                          dt,
+                                                                                          localMatrix.toViewConstSizes(),
+                                                                                          localRhs.toView(),
+                                                                                          dR_dAper );
+      }
+    } );
+  } );
+}
+
 template< typename BASE >
 void
 SinglePhaseFVM< BASE >::applyBoundaryConditions( real64 const time_n,