From 8f9126c37266dac8a757b06e97a6ca4a0bdaa051 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Wed, 15 Jan 2025 16:11:51 -0600 Subject: [PATCH] feat: Multiphase poromechanics with contact (#3228) Co-authored-by: Matteo Cusini --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 + .../poromechanicsFractures/co2flash.txt | 1 + .../poromechanicsFractures/elevation.txt | 2 + .../initCO2CompFrac.txt | 2 + .../poromechanicsFractures/initTemp.txt | 2 + .../initWaterCompFrac.txt | 2 + ...ultiphasePoromechanics_FaultModel_base.xml | 207 +++++ ...ltiphasePoromechanics_FaultModel_smoke.xml | 130 +++ .../phaseVolFraction_gas.txt | 21 + .../phaseVolFraction_water.txt | 21 + inputFiles/poromechanicsFractures/pvtgas.txt | 2 + .../poromechanicsFractures/pvtliquid.txt | 2 + .../poromechanicsFractures/relPerm_gas.txt | 21 + .../poromechanicsFractures/relPerm_water.txt | 21 + .../fluidFlow/CompositionalMultiphaseBase.cpp | 5 +- .../fluidFlow/CompositionalMultiphaseBase.hpp | 71 ++ .../fluidFlow/CompositionalMultiphaseFVM.cpp | 84 +- .../fluidFlow/CompositionalMultiphaseFVM.hpp | 23 +- .../CompositionalMultiphaseHybridFVM.cpp | 2 +- .../CompositionalMultiphaseUtilities.hpp | 13 + .../fluidFlow/FlowSolverBase.cpp | 1 + .../fluidFlow/FlowSolverBase.hpp | 22 + .../fluidFlow/SinglePhaseBase.cpp | 3 - .../fluidFlow/SinglePhaseBase.hpp | 19 - .../fluidFlow/SinglePhaseHybridFVM.cpp | 19 - .../fluidFlow/SinglePhaseHybridFVM.hpp | 9 - .../kernels/compositional/C1PPUPhaseFlux.hpp | 4 +- .../compositional/FluxComputeKernel.hpp | 4 +- .../kernels/compositional/IHUPhaseFlux.hpp | 4 +- .../kernels/compositional/PPUPhaseFlux.hpp | 26 +- .../compositional/PhaseComponentFlux.hpp | 88 ++ .../kernels/compositional/PotGrad.hpp | 20 +- .../compositional/SolutionCheckKernel.hpp | 8 +- .../ThermalSolutionCheckKernel.hpp | 6 +- .../wells/CompositionalMultiphaseWell.cpp | 4 +- .../CompositionalMultiphaseWellKernels.hpp | 2 +- .../multiphysics/CMakeLists.txt | 3 + .../multiphysics/MultiphasePoromechanics.cpp | 4 +- ...iphasePoromechanicsConformingFractures.cpp | 833 ++++++++++++++++++ ...iphasePoromechanicsConformingFractures.hpp | 203 +++++ .../PoromechanicsInitialization.cpp | 3 + .../multiphysics/PoromechanicsSolver.hpp | 4 +- ...ePhasePoromechanicsConformingFractures.hpp | 3 - ...glePhasePoromechanicsEmbeddedFractures.cpp | 13 +- ...glePhasePoromechanicsEmbeddedFractures.hpp | 4 +- ...iphasePoromechanicsConformingFractures.hpp | 410 +++++++++ .../SinglePhasePoromechanics.hpp | 1 + 48 files changed, 2239 insertions(+), 119 deletions(-) create mode 100644 inputFiles/poromechanicsFractures/co2flash.txt create mode 100644 inputFiles/poromechanicsFractures/elevation.txt create mode 100644 inputFiles/poromechanicsFractures/initCO2CompFrac.txt create mode 100644 inputFiles/poromechanicsFractures/initTemp.txt create mode 100644 inputFiles/poromechanicsFractures/initWaterCompFrac.txt create mode 100644 inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_base.xml create mode 100644 inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_smoke.xml create mode 100644 inputFiles/poromechanicsFractures/phaseVolFraction_gas.txt create mode 100644 inputFiles/poromechanicsFractures/phaseVolFraction_water.txt create mode 100644 inputFiles/poromechanicsFractures/pvtgas.txt create mode 100644 inputFiles/poromechanicsFractures/pvtliquid.txt create mode 100644 inputFiles/poromechanicsFractures/relPerm_gas.txt create mode 100644 inputFiles/poromechanicsFractures/relPerm_water.txt create mode 100644 src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp create mode 100644 src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp create mode 100644 src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 746f223b649..2cb02daa192 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3495-9552-257c87e + baseline: integratedTests/baseline_integratedTests-pr3228-9676-61994fe allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index c5617ba8ed9..ea19f4ac5cd 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3228 (2024-01-15) +===================== +deltaVolume added in multiphase. + PR #3495 (2024-01-08) ===================== Add missing logic to support switching from fixed mass rate injection rate constraint to max injection pressure. diff --git a/inputFiles/poromechanicsFractures/co2flash.txt b/inputFiles/poromechanicsFractures/co2flash.txt new file mode 100644 index 00000000000..8de4352642b --- /dev/null +++ b/inputFiles/poromechanicsFractures/co2flash.txt @@ -0,0 +1 @@ +FlashModel CO2Solubility 1e6 10e7 5e4 367.15 369.15 1 0 diff --git a/inputFiles/poromechanicsFractures/elevation.txt b/inputFiles/poromechanicsFractures/elevation.txt new file mode 100644 index 00000000000..9f9ff9d4f67 --- /dev/null +++ b/inputFiles/poromechanicsFractures/elevation.txt @@ -0,0 +1,2 @@ +0 +1e4 \ No newline at end of file diff --git a/inputFiles/poromechanicsFractures/initCO2CompFrac.txt b/inputFiles/poromechanicsFractures/initCO2CompFrac.txt new file mode 100644 index 00000000000..f3c5e203d6f --- /dev/null +++ b/inputFiles/poromechanicsFractures/initCO2CompFrac.txt @@ -0,0 +1,2 @@ +0.001 +0.001 diff --git a/inputFiles/poromechanicsFractures/initTemp.txt b/inputFiles/poromechanicsFractures/initTemp.txt new file mode 100644 index 00000000000..da7b374fdb7 --- /dev/null +++ b/inputFiles/poromechanicsFractures/initTemp.txt @@ -0,0 +1,2 @@ +368.15 +368.15 diff --git a/inputFiles/poromechanicsFractures/initWaterCompFrac.txt b/inputFiles/poromechanicsFractures/initWaterCompFrac.txt new file mode 100644 index 00000000000..c88775171f6 --- /dev/null +++ b/inputFiles/poromechanicsFractures/initWaterCompFrac.txt @@ -0,0 +1,2 @@ +0.999 +0.999 diff --git a/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_base.xml b/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_base.xml new file mode 100644 index 00000000000..fc1798c6851 --- /dev/null +++ b/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_base.xml @@ -0,0 +1,207 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_smoke.xml b/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_smoke.xml new file mode 100644 index 00000000000..cb87c0a50e3 --- /dev/null +++ b/inputFiles/poromechanicsFractures/multiphasePoromechanics_FaultModel_smoke.xml @@ -0,0 +1,130 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/poromechanicsFractures/phaseVolFraction_gas.txt b/inputFiles/poromechanicsFractures/phaseVolFraction_gas.txt new file mode 100644 index 00000000000..e2d1da1f128 --- /dev/null +++ b/inputFiles/poromechanicsFractures/phaseVolFraction_gas.txt @@ -0,0 +1,21 @@ +0 +0.05 +0.1 +0.15 +0.2 +0.25 +0.3 +0.35 +0.4 +0.45 +0.5 +0.55 +0.6 +0.65 +0.7 +0.75 +0.8 +0.85 +0.9 +0.95 +1.0 diff --git a/inputFiles/poromechanicsFractures/phaseVolFraction_water.txt b/inputFiles/poromechanicsFractures/phaseVolFraction_water.txt new file mode 100644 index 00000000000..e2d1da1f128 --- /dev/null +++ b/inputFiles/poromechanicsFractures/phaseVolFraction_water.txt @@ -0,0 +1,21 @@ +0 +0.05 +0.1 +0.15 +0.2 +0.25 +0.3 +0.35 +0.4 +0.45 +0.5 +0.55 +0.6 +0.65 +0.7 +0.75 +0.8 +0.85 +0.9 +0.95 +1.0 diff --git a/inputFiles/poromechanicsFractures/pvtgas.txt b/inputFiles/poromechanicsFractures/pvtgas.txt new file mode 100644 index 00000000000..9da1273477c --- /dev/null +++ b/inputFiles/poromechanicsFractures/pvtgas.txt @@ -0,0 +1,2 @@ +DensityFun SpanWagnerCO2Density 1e6 10e7 5e4 367.15 369.15 1 +ViscosityFun FenghourCO2Viscosity 1e6 10e7 5e4 367.15 369.15 1 diff --git a/inputFiles/poromechanicsFractures/pvtliquid.txt b/inputFiles/poromechanicsFractures/pvtliquid.txt new file mode 100644 index 00000000000..b0bfb5fc76e --- /dev/null +++ b/inputFiles/poromechanicsFractures/pvtliquid.txt @@ -0,0 +1,2 @@ +DensityFun PhillipsBrineDensity 1e6 10e7 5e4 367.15 369.15 1 0 +ViscosityFun PhillipsBrineViscosity 0 diff --git a/inputFiles/poromechanicsFractures/relPerm_gas.txt b/inputFiles/poromechanicsFractures/relPerm_gas.txt new file mode 100644 index 00000000000..d6c4e7c035e --- /dev/null +++ b/inputFiles/poromechanicsFractures/relPerm_gas.txt @@ -0,0 +1,21 @@ +0 +0 +0.0118 +0.0333 +0.0612 +0.0943 +0.1318 +0.1732 +0.2183 +0.2667 +0.3182 +0.3727 +0.4300 +0.4899 +0.5524 +0.6173 +0.6847 +0.7542 +0.8261 +0.9000 +0.9760 diff --git a/inputFiles/poromechanicsFractures/relPerm_water.txt b/inputFiles/poromechanicsFractures/relPerm_water.txt new file mode 100644 index 00000000000..d6c4e7c035e --- /dev/null +++ b/inputFiles/poromechanicsFractures/relPerm_water.txt @@ -0,0 +1,21 @@ +0 +0 +0.0118 +0.0333 +0.0612 +0.0943 +0.1318 +0.1732 +0.2183 +0.2667 +0.3182 +0.3727 +0.4300 +0.4899 +0.5524 +0.6173 +0.6847 +0.7542 +0.8261 +0.9000 +0.9760 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 5bb8896898d..c13be912619 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -42,8 +42,6 @@ #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SourceFluxStatistics.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/GlobalComponentFractionKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/PhaseVolumeFractionKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseVolumeFractionKernel.hpp" @@ -418,7 +416,7 @@ void CompositionalMultiphaseBase::setConstitutiveNames( ElementSubRegionBase & s string & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); fluidName = getConstitutiveName< MultiFluidBase >( subRegion ); GEOS_THROW_IF( fluidName.empty(), - GEOS_FMT( "{}: Fluid model not found on subregion {}", + GEOS_FMT( "{}: multiphase fluid model not found on subregion {}", getDataContext(), subRegion.getDataContext() ), InputError ); @@ -1333,7 +1331,6 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom if( m_useSimpleAccumulation ) kernelFlags.set( KernelFlags::SimpleAccumulation ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, arrayView1d< string const > const & regionNames ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 318662d00fc..cb9c1fe33d3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -22,6 +22,10 @@ #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "constitutive/solid/CoupledSolidBase.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/AccumulationKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/ThermalAccumulationKernel.hpp" namespace geos { @@ -98,6 +102,21 @@ class CompositionalMultiphaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) override; + /** + * @brief assembles the accumulation terms for all cells of a spcefici subRegion. + * @tparam SUBREGION_TYPE the subRegion type + * @param dofManager degree-of-freedom manager associated with the linear system + * @param subRegion the subRegion + * @param localMatrix the system matrix + * @param localRhs the system right-hand side vector + * + */ + template< typename SUBREGION_TYPE > + void accumulationAssemblyLaunch( DofManager const & dofManager, + SUBREGION_TYPE const & subRegion, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + virtual void resetStateToBeginningOfStep( DomainPartition & domain ) override; @@ -541,6 +560,58 @@ void CompositionalMultiphaseBase::applyFieldValue( real64 const & time_n, } ); } +template< typename SUBREGION_TYPE > +void CompositionalMultiphaseBase::accumulationAssemblyLaunch( DofManager const & dofManager, + SUBREGION_TYPE const & subRegion, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + constitutive::MultiFluidBase const & fluid = + getConstitutiveModel< constitutive::MultiFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) ); + constitutive::CoupledSolidBase const & solid = + getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + + string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); + + using namespace isothermalCompositionalMultiphaseBaseKernels; + + BitFlags< KernelFlags > kernelFlags; + if( m_useTotalMassEquation ) + kernelFlags.set( KernelFlags::TotalMassEquation ); + if( m_useSimpleAccumulation ) + kernelFlags.set( KernelFlags::SimpleAccumulation ); + + if( m_isThermal ) + { + thermalCompositionalMultiphaseBaseKernels:: + AccumulationKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + kernelFlags, + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); + } + else + { + isothermalCompositionalMultiphaseBaseKernels:: + AccumulationKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + kernelFlags, + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); + } +} } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index d2dbcb2f3a0..aeb88ef68e3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -54,12 +54,14 @@ #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseMobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/AquiferBCKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp" namespace geos { using namespace dataRepository; using namespace constitutive; +using namespace compositionalMultiphaseUtilities; // for ScalingType CompositionalMultiphaseFVM::CompositionalMultiphaseFVM( const string & name, Group * const parent ) @@ -1264,9 +1266,89 @@ void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, } -real64 CompositionalMultiphaseFVM::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain ) +void CompositionalMultiphaseFVM::assembleHydrofracFluxTerms( 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 & elemDofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); + + using namespace isothermalCompositionalMultiphaseFVMKernels; + + BitFlags< KernelFlags > kernelFlags; + if( m_hasCapPressure ) + kernelFlags.set( KernelFlags::CapPressure ); + if( m_useTotalMassEquation ) + kernelFlags.set( KernelFlags::TotalMassEquation ); + if( m_gravityDensityScheme == GravityDensityScheme::PhasePresence ) + kernelFlags.set( KernelFlags::CheckPhasePresenceInGravity ); + if( m_isThermal ) + { + // should not end up here but just in case + GEOS_ERROR( "Thermal not yet supported in CompositionalMultiphaseFVM::assembleHydrofracFluxTerms" ); + } + if( fluxApprox.upwindingParams().upwindingScheme != UpwindingScheme::PPU ) + { + // a bit tricky to check in advance + GEOS_ERROR( "Only PPU upwinding is supported in CompositionalMultiphaseFVM::assembleHydrofracFluxTerms" ); + } + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & ) + { + fluxApprox.forStencils< CellElementStencilTPFA, FaceElementToCellStencil >( mesh, [&]( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + isothermalCompositionalMultiphaseFVMKernels:: + FluxComputeKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + kernelFlags, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + multiphasePoromechanicsConformingFracturesKernels:: + FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + kernelFlags, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView(), + dR_dAper ); + } ); + } ); +} + +real64 CompositionalMultiphaseFVM::setNextDt( const geos::real64 & currentDt, geos::DomainPartition & domain ) +{ if( m_targetFlowCFL < 0 ) return CompositionalMultiphaseBase::setNextDt( currentDt, domain ); else diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp index 896dc87b7f1..14d66667035 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp @@ -163,6 +163,14 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const override; + virtual void + assembleHydrofracFluxTerms( 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 ) override final; virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const override; @@ -205,15 +213,6 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase static constexpr char const * gravityDensitySchemeString() { return "gravityDensityScheme"; } }; - /** - * @brief Solution scaling type - */ - enum class ScalingType : integer - { - Global, ///< Scale the Newton update with a unique scaling factor - Local ///< Scale the Newton update locally (modifies the Newton direction) - }; - /** * @brief Storage for value and element location, used to determine global max + location */ @@ -255,7 +254,7 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase } m_dbcParams; /// Solution scaling type - ScalingType m_scalingType; + compositionalMultiphaseUtilities::ScalingType m_scalingType; /// scheme for density treatment in gravity GravityDensityScheme m_gravityDensityScheme; @@ -293,10 +292,6 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase }; -ENUM_STRINGS( CompositionalMultiphaseFVM::ScalingType, - "Global", - "Local" ); - } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index 021ce8fd96b..53415a507cd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -541,7 +541,7 @@ bool CompositionalMultiphaseHybridFVM::checkSystemSolution( DomainPartition & do SolutionCheckKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping, m_allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType::Global, + compositionalMultiphaseUtilities::ScalingType::Global, scalingFactor, pressure, compDens, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp index 3af76d4928d..cbac0f13f2d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp @@ -28,6 +28,19 @@ namespace geos namespace compositionalMultiphaseUtilities { +/** + * @brief Solution scaling type, used in CompositionalMultiphaseFVM + */ +enum class ScalingType : integer +{ + Global, ///< Scale the Newton update with a unique scaling factor + Local ///< Scale the Newton update locally (modifies the Newton direction) +}; + +ENUM_STRINGS( ScalingType, + "Global", + "Local" ); + /** * @brief In each block, shift the elements from 0 to numRowsToShift-1 one position ahead and replaces the first element * with the sum of all elements from 0 to numRowsToShift-1 of the input block one-dimensional array of values. diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index c081f20d4ce..e9a55b86c46 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -143,6 +143,7 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies ) [&]( localIndex const, ElementSubRegionBase & subRegion ) { + subRegion.registerField< fields::flow::deltaVolume >( getName() ); subRegion.registerField< fields::flow::gravityCoefficient >( getName() ). setApplyDefaultValue( 0.0 ); subRegion.registerField< fields::flow::netToGross >( getName() ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 4352e059405..b97716d0b00 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -146,6 +146,28 @@ class FlowSolverBase : public PhysicsSolverBase real64 const & timeAtBeginningOfStep, real64 const & dt ); + /** + * @brief assembles the flux terms for all cells for the hydrofracture case + * @param time_n previous time value + * @param dt time step + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localMatrix the system matrix + * @param localRhs the system right-hand side vector + * @param dR_dAper + */ + virtual void assembleHydrofracFluxTerms( 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 not yet implemented." ); + } + virtual void initializeFluidState( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } virtual void initializeThermalState( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 98934988d8d..94f3d46a640 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -93,15 +93,12 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies ) [&]( localIndex const, ElementSubRegionBase & subRegion ) { - subRegion.registerField< deltaVolume >( getName() ); - subRegion.registerField< mobility >( getName() ); subRegion.registerField< dMobility_dPressure >( getName() ); subRegion.registerField< fields::flow::mass >( getName() ); subRegion.registerField< fields::flow::mass_n >( getName() ); - if( m_isThermal ) { subRegion.registerField< dMobility_dTemperature >( getName() ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index db72beaa75a..a4049c28a87 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -205,25 +205,6 @@ class SinglePhaseBase : public FlowSolverBase arrayView1d< real64 > const & localRhs, string const & jumpDofKey ) = 0; - /** - * @brief assembles the flux terms for all cells for the hydrofracture case - * @param time_n previous time value - * @param dt time step - * @param domain the physical domain object - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localMatrix the system matrix - * @param localRhs the system right-hand side vector - * @param dR_dAper - */ - virtual void - assembleHydrofracFluxTerms( 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 ) = 0; - struct viewKeyStruct : FlowSolverBase::viewKeyStruct { static constexpr char const * elemDofFieldString() { return "singlePhaseVariables"; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp index e0b9732cec2..5884b6fc416 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp @@ -311,25 +311,6 @@ void SinglePhaseHybridFVM::assembleEDFMFluxTerms( real64 const GEOS_UNUSED_PARAM localRhs ); } -void SinglePhaseHybridFVM::assembleHydrofracFluxTerms( 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 ); - GEOS_UNUSED_VAR ( dt ); - GEOS_UNUSED_VAR ( domain ); - GEOS_UNUSED_VAR ( dofManager ); - GEOS_UNUSED_VAR ( localMatrix ); - GEOS_UNUSED_VAR ( localRhs ); - GEOS_UNUSED_VAR ( dR_dAper ); - - GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." ); -} - void SinglePhaseHybridFVM::applyBoundaryConditions( real64 const time_n, real64 const dt, DomainPartition & domain, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp index 7a804ccf7ee..b2910490508 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp @@ -142,15 +142,6 @@ class SinglePhaseHybridFVM : public SinglePhaseBase arrayView1d< real64 > const & localRhs, string const & jumpDofKey ) override final; - virtual void - assembleHydrofracFluxTerms( 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 ) override final; - virtual void applyAquiferBC( real64 const time, real64 const dt, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp index f92335dd726..3d79301623a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp @@ -107,6 +107,7 @@ struct C1PPUPhaseFlux real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp], real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] ) { + real64 dPotGrad_dTrans = 0.0; real64 dPresGrad_dP[numFluxSupportPoints]{}; real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; @@ -114,7 +115,8 @@ struct C1PPUPhaseFlux PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, checkPhasePresenceInGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, - phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP, + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + potGrad, dPotGrad_dTrans, dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); // gravity head diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp index 78271f54416..45882953307 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp @@ -252,6 +252,7 @@ class FluxComputeKernel : public FluxComputeKernelBase real64 compFlux[numComp]{}; real64 dCompFlux_dP[numFluxSupportPoints][numComp]{}; real64 dCompFlux_dC[numFluxSupportPoints][numComp][numComp]{}; + real64 dCompFlux_dTrans[numComp]{}; real64 const trans[numFluxSupportPoints] = { stack.transmissibility[connectionIndex][0], stack.transmissibility[connectionIndex][1] }; @@ -350,7 +351,8 @@ class FluxComputeKernel : public FluxComputeKernelBase dPhaseFlux_dC, compFlux, dCompFlux_dP, - dCompFlux_dC ); + dCompFlux_dC, + dCompFlux_dTrans ); } // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp index 9365a96c830..97f5bdd5a2e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp @@ -1760,7 +1760,7 @@ struct IHUPhaseFlux real64 dummy[numComp]; real64 dDummy_dP[numFluxSupportPoints][numComp]; real64 dDummy_dC[numFluxSupportPoints][numComp][numComp]; - + real64 dDummy_dTrans[numComp]; for( integer jp = 0; jp < numPhase; ++jp ) { @@ -1777,7 +1777,7 @@ struct IHUPhaseFlux phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, k_up_ppu, potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, - dummy, dDummy_dP, dDummy_dC ); + dummy, dDummy_dP, dDummy_dC, dDummy_dTrans ); totFlux += phaseFlux; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp index 468adabb4b1..14fc73f029d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp @@ -68,6 +68,10 @@ struct PPUPhaseFlux * @param phaseFlux phase flux * @param dPhaseFlux_dP derivative of phase flux wrt pressure * @param dPhaseFlux_dC derivative of phase flux wrt comp density + * @param compFlux component flux + * @param dCompFlux_dP derivative of phase flux wrt pressure + * @param dCompFlux_dC derivative of phase flux wrt comp density + * @param dCompFlux_dTrans derivative of phase flux wrt transmissibility */ template< integer numComp, integer numFluxSupportPoints > GEOS_HOST_DEVICE @@ -96,13 +100,15 @@ struct PPUPhaseFlux ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, localIndex & k_up, real64 & potGrad, - real64 ( &phaseFlux ), + real64 & phaseFlux, real64 ( & dPhaseFlux_dP )[numFluxSupportPoints], real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp], real64 ( & compFlux )[numComp], real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp], - real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] ) + real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp], + real64 ( & dCompFlux_dTrans )[numComp] ) { + real64 dPotGrad_dTrans = 0; real64 dPresGrad_dP[numFluxSupportPoints]{}; real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; @@ -111,8 +117,9 @@ struct PPUPhaseFlux seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, - phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, - dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); + phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, + potGrad, dPotGrad_dTrans, dPresGrad_dP, + dPresGrad_dC, dGravHead_dP, dGravHead_dC ); // *** upwinding *** @@ -125,6 +132,11 @@ struct PPUPhaseFlux real64 const mobility = phaseMob[er_up][esr_up][ei_up][ip]; + // compute phase flux using upwind mobility + phaseFlux = mobility * potGrad; + + real64 const dPhaseFlux_dTrans = mobility * dPotGrad_dTrans; + // pressure gradient depends on all points in the stencil for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { @@ -136,8 +148,6 @@ struct PPUPhaseFlux dPhaseFlux_dC[ke][jc] *= mobility; } } - // compute phase flux using upwind mobility. - phaseFlux = mobility * potGrad; real64 const dMob_dP = dPhaseMob[er_up][esr_up][ei_up][ip][Deriv::dP]; arraySlice1d< real64 const, compflow::USD_PHASE_DC - 2 > dPhaseMobSub = @@ -151,8 +161,8 @@ struct PPUPhaseFlux } //distribute on phaseComponentFlux here - PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei, phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, phaseFlux - , dPhaseFlux_dP, dPhaseFlux_dC, compFlux, dCompFlux_dP, dCompFlux_dC ); + PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei, phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans, compFlux, dCompFlux_dP, dCompFlux_dC, dCompFlux_dTrans ); } }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp index 3d7f1cba38b..f7810873e5e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp @@ -119,6 +119,94 @@ struct PhaseComponentFlux } } } + + /** + * @brief Compute the component flux for a given phase + * @tparam numComp number of components + * @tparam numFluxSupportPoints number of flux support points + * @param ip phase index + * @param k_up uptream index for this phase + * @param seri arraySlice of the stencil-implied element region index + * @param sesri arraySlice of the stencil-implied element subregion index + * @param sei arraySlice of the stencil-implied element index + * @param phaseCompFrac phase component fraction + * @param dPhaseCompFrac derivative of phase component fraction wrt pressure, temperature, component fraction + * @param dCompFrac_dCompDens derivative of component fraction wrt component density + * @param phaseFlux phase flux + * @param dPhaseFlux_dP derivative of phase flux wrt pressure + * @param dPhaseFlux_dC derivative of phase flux wrt comp density + * @param dPhaseFlux_dTrans derivative of phase flux wrt transmissibility + * @param compFlux component flux + * @param dCompFlux_dP derivative of phase flux wrt pressure + * @param dCompFlux_dC derivative of phase flux wrt comp density + * @param dCompFlux_dTrans derivative of phase flux wrt transmissibility + */ + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + compute( localIndex const ip, + localIndex const k_up, + localIndex const ( &seri )[numFluxSupportPoints], + localIndex const ( &sesri )[numFluxSupportPoints], + localIndex const ( &sei )[numFluxSupportPoints], + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, + ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, + ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, + real64 const & phaseFlux, + real64 const ( &dPhaseFlux_dP )[numFluxSupportPoints], + real64 const ( &dPhaseFlux_dC )[numFluxSupportPoints][numComp], + real64 const & dPhaseFlux_dTrans, + real64 ( & compFlux )[numComp], + real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp], + real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp], + real64 ( & dCompFlux_dTrans )[numComp] ) + { + localIndex const er_up = seri[k_up]; + localIndex const esr_up = sesri[k_up]; + localIndex const ei_up = sei[k_up]; + + real64 dProp_dC[numComp]{}; + + // slice some constitutive arrays to avoid too much indexing in component loop + arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE_COMP-3 > phaseCompFracSub = + phaseCompFrac[er_up][esr_up][ei_up][0][ip]; + arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC-3 > dPhaseCompFracSub = + dPhaseCompFrac[er_up][esr_up][ei_up][0][ip]; + + // compute component fluxes and derivatives using upstream cell composition + for( integer ic = 0; ic < numComp; ++ic ) + { + real64 const ycp = phaseCompFracSub[ic]; + compFlux[ic] += phaseFlux * ycp; + + // derivatives stemming from phase flux + for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) + { + dCompFlux_dP[ke][ic] += dPhaseFlux_dP[ke] * ycp; + for( integer jc = 0; jc < numComp; ++jc ) + { + dCompFlux_dC[ke][ic][jc] += dPhaseFlux_dC[ke][jc] * ycp; + } + } + + // additional derivatives stemming from upstream cell phase composition + dCompFlux_dP[k_up][ic] += phaseFlux * dPhaseCompFracSub[ic][Deriv::dP]; + + // convert derivatives of comp fraction w.r.t. comp fractions to derivatives w.r.t. comp densities + applyChainRule( numComp, + dCompFrac_dCompDens[er_up][esr_up][ei_up], + dPhaseCompFracSub[ic], + dProp_dC, + Deriv::dC ); + for( integer jc = 0; jc < numComp; ++jc ) + { + dCompFlux_dC[k_up][ic][jc] += phaseFlux * dProp_dC[jc]; + } + + dCompFlux_dTrans[ic] += dPhaseFlux_dTrans * ycp; + } + } + }; } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp index b783b061712..f9ca248b4dc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp @@ -62,6 +62,7 @@ struct PotGrad ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, real64 & potGrad, + real64 & dPotGrad_dTrans, real64 ( & dPresGrad_dP )[numFluxSupportPoints], real64 ( & dPresGrad_dC )[numFluxSupportPoints][numComp], real64 ( & dGravHead_dP )[numFluxSupportPoints], @@ -85,7 +86,9 @@ struct PotGrad real64 dDensMean_dC[numFluxSupportPoints][numComp]{}; real64 presGrad = 0.0; + real64 dPresGrad_dTrans = 0.0; real64 gravHead = 0.0; + real64 dGravHead_dTrans = 0.0; real64 dCapPressure_dC[numComp]{}; calculateMeanDensity( checkPhasePresenceInGravity, ip, @@ -126,22 +129,25 @@ struct PotGrad } } - presGrad += trans[i] * (pres[er][esr][ei] - capPressure); - dPresGrad_dP[i] += trans[i] * (1 - dCapPressure_dP) - + dTrans_dPres[i] * (pres[er][esr][ei] - capPressure); + real64 const dP = pres[er][esr][ei] - capPressure; + presGrad += trans[i] * dP; + dPresGrad_dTrans += dP; + dPresGrad_dP[i] += trans[i] * (1 - dCapPressure_dP) + dTrans_dPres[i] * dP; for( integer jc = 0; jc < numComp; ++jc ) { dPresGrad_dC[i][jc] += -trans[i] * dCapPressure_dC[jc]; } - real64 const gravD = trans[i] * gravCoef[er][esr][ei]; - real64 const dGravD_dP = dTrans_dPres[i] * gravCoef[er][esr][ei]; + real64 const gC = gravCoef[er][esr][ei]; + real64 const gravD = trans[i] * gC; + real64 const dGravD_dTrans = gC; + real64 const dGravD_dP = dTrans_dPres[i] * gC; // the density used in the potential difference is always a mass density // unlike the density used in the phase mobility, which is a mass density // if useMass == 1 and a molar density otherwise gravHead += densMean * gravD; - + dGravHead_dTrans += densMean * dGravD_dTrans; // need to add contributions from both cells the mean density depends on for( integer j = 0; j < numFluxSupportPoints; ++j ) { @@ -155,7 +161,7 @@ struct PotGrad // compute phase potential gradient potGrad = presGrad - gravHead; - + dPotGrad_dTrans = dPresGrad_dTrans - dGravHead_dTrans; } template< integer numComp, integer numFluxSupportPoints > diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp index 5cbf59b932d..89d681bb02c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp @@ -61,7 +61,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer */ SolutionCheckKernel( integer const allowCompDensChopping, integer const allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView2d< real64 const, compflow::USD_COMP > const compDens, @@ -217,7 +217,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer StackVariables & stack, FUNC && kernelOp = NoOpFunc{} ) const { - bool const localScaling = m_scalingType == CompositionalMultiphaseFVM::ScalingType::Local; + bool const localScaling = m_scalingType == compositionalMultiphaseUtilities::ScalingType::Local; real64 const newPres = m_pressure[ei] + (localScaling ? m_pressureScalingFactor[ei] : m_scalingFactor) * m_localSolution[stack.localRow]; if( newPres < 0 ) @@ -280,7 +280,7 @@ class SolutionCheckKernel : public SolutionScalingAndCheckingKernelBase< integer real64 const m_scalingFactor; /// scaling type (global or local) - CompositionalMultiphaseFVM::ScalingType const m_scalingType; + compositionalMultiphaseUtilities::ScalingType const m_scalingType; }; @@ -306,7 +306,7 @@ class SolutionCheckKernelFactory static SolutionCheckKernel::StackVariables createAndLaunch( integer const allowCompDensChopping, integer const allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView2d< real64 const, compflow::USD_COMP > const compDens, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp index 35c969bb891..b0447b33782 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalSolutionCheckKernel.hpp @@ -60,7 +60,7 @@ class SolutionCheckKernel : public isothermalCompositionalMultiphaseBaseKernels: */ SolutionCheckKernel( integer const allowCompDensChopping, integer const allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView1d< real64 const > const temperature, @@ -103,7 +103,7 @@ class SolutionCheckKernel : public isothermalCompositionalMultiphaseBaseKernels: { Base::computeSolutionCheck( ei, stack, [&] () { - bool const localScaling = m_scalingType == CompositionalMultiphaseFVM::ScalingType::Local; + bool const localScaling = m_scalingType == compositionalMultiphaseUtilities::ScalingType::Local; // compute the change in temperature real64 const newTemp = m_temperature[ei] + (localScaling ? m_temperatureScalingFactor[ei] : m_scalingFactor * m_localSolution[stack.localRow + m_temperatureOffset]); if( newTemp < minTemperature ) @@ -149,7 +149,7 @@ class SolutionCheckKernelFactory static SolutionCheckKernel::StackVariables createAndLaunch( integer const allowCompDensChopping, integer const allowNegativePressure, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView1d< real64 const > const temperature, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index b891467863d..e33e4ec0b17 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -1535,7 +1535,7 @@ CompositionalMultiphaseWell::checkSystemSolution( DomainPartition & domain, compositionalMultiphaseWellKernels:: SolutionCheckKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_allowCompDensChopping, - CompositionalMultiphaseFVM::ScalingType::Global, + compositionalMultiphaseUtilities::ScalingType::Global, scalingFactor, pressure, compDens, @@ -1573,7 +1573,7 @@ CompositionalMultiphaseWell::checkSystemSolution( DomainPartition & domain, { //integer const m_allowCompDensChopping(true); integer const m_allowNegativePressure( false ); - CompositionalMultiphaseFVM::ScalingType const m_scalingType( CompositionalMultiphaseFVM::ScalingType::Global ); + compositionalMultiphaseUtilities::ScalingType const m_scalingType( compositionalMultiphaseUtilities::ScalingType::Global ); arrayView1d< real64 const > const pressure = subRegion.getField< fields::well::pressure >(); arrayView1d< real64 const > const temperature = diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp index c2e9078ba47..5722613372d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/CompositionalMultiphaseWellKernels.hpp @@ -810,7 +810,7 @@ class SolutionCheckKernelFactory template< typename POLICY > static isothermalCompositionalMultiphaseBaseKernels::SolutionCheckKernel::StackVariables createAndLaunch( integer const allowCompDensChopping, - CompositionalMultiphaseFVM::ScalingType const scalingType, + compositionalMultiphaseUtilities::ScalingType const scalingType, real64 const scalingFactor, arrayView1d< real64 const > const pressure, arrayView2d< real64 const, compflow::USD_COMP > const compDens, diff --git a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt index 89a34cbc276..2332116db38 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt @@ -9,6 +9,7 @@ set( physicsSolvers_headers multiphysics/HydrofractureSolver.hpp multiphysics/HydrofractureSolverKernels.hpp multiphysics/MultiphasePoromechanics.hpp + multiphysics/MultiphasePoromechanicsConformingFractures.hpp multiphysics/PhaseFieldFractureSolver.hpp multiphysics/PoromechanicsInitialization.hpp multiphysics/PoromechanicsFields.hpp @@ -16,6 +17,7 @@ set( physicsSolvers_headers multiphysics/LogLevelsInfo.hpp multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp + multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp multiphysics/poromechanicsKernels/PoromechanicsBase.hpp multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp @@ -46,6 +48,7 @@ set( physicsSolvers_sources multiphysics/FlowProppantTransportSolver.cpp multiphysics/HydrofractureSolver.cpp multiphysics/MultiphasePoromechanics.cpp + multiphysics/MultiphasePoromechanicsConformingFractures.cpp multiphysics/PhaseFieldFractureSolver.cpp multiphysics/PoromechanicsInitialization.cpp multiphysics/SinglePhasePoromechanics.cpp diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp index 6b01f872f5c..0d11624cfc1 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp @@ -32,7 +32,7 @@ #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp" -//#include "physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp" +#include "physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp" //#include "physicsSolvers/contact/SolidMechanicsEmbeddedFractures.hpp" namespace geos @@ -361,7 +361,7 @@ void MultiphasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::updateBulkDensity } template class MultiphasePoromechanics<>; -//template class MultiphasePoromechanics< CompositionalMultiphaseBase, SolidMechanicsLagrangeContact >; +template class MultiphasePoromechanics< CompositionalMultiphaseBase, SolidMechanicsLagrangeContact >; //template class MultiphasePoromechanics< CompositionalMultiphaseBase, SolidMechanicsEmbeddedFractures >; template class MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> >; diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp new file mode 100644 index 00000000000..8ef737ad998 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp @@ -0,0 +1,833 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultiphasePoromechanicsConformingFractures.cpp + */ + +#include "MultiphasePoromechanicsConformingFractures.hpp" + +#include "constitutive/contact/HydraulicApertureRelationSelector.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" +#include "finiteVolume/FluxApproximationBase.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; +using namespace fields; + +template< typename FLOW_SOLVER > +MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::MultiphasePoromechanicsConformingFractures( const string & name, + Group * const parent ) + : Base( name, parent ) +{ + // TODO: MGR recipe + // LinearSolverParameters & params = this->m_linearSolverParameters.get(); + // params.mgr.strategy = LinearSolverParameters::MGR::StrategyType::multiphasePoromechanicsConformingFractures; + // params.mgr.separateComponents = false; + // params.mgr.displacementFieldName = solidMechanics::totalDisplacement::key(); + // params.dofsPerNode = 3; +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::postInputInitialization() +{ + Base::postInputInitialization(); + + if( this->flowSolver()->isThermal() ) + { + GEOS_ERROR( "Thermal flow is not yet supported for multiphase poromechanics conforming fractures" ); + } +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::setupCoupling( DomainPartition const & domain, + DofManager & dofManager ) const +{ + /// We need to add 2 coupling terms: + // 1. Poromechanical coupling in the bulk + Base::setupCoupling( domain, dofManager ); + + // 2. Traction - pressure coupling in the fracture + dofManager.addCoupling( m_flowDofKey, + fields::contact::traction::key(), + DofManager::Connector::Elem ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity ) +{ + GEOS_MARK_FUNCTION; + + GEOS_UNUSED_VAR( setSparsity ); + + /// 1. Add all coupling terms handled directly by the DofManager + dofManager.setDomain( domain ); + this->setupDofs( domain, dofManager ); + dofManager.reorderByRank(); + + /// 2. Add coupling terms not added by the DofManager. + localIndex const numLocalRows = dofManager.numLocalDofs(); + + SparsityPattern< globalIndex > patternOriginal; + dofManager.setSparsityPattern( patternOriginal ); + + // Get the original row lengths (diagonal blocks only) + array1d< localIndex > rowLengths( patternOriginal.numRows() ); + for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + { + rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); + } + + // Add the number of nonzeros induced by coupling + addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() ); + + // Create a new pattern with enough capacity for coupled matrix + SparsityPattern< globalIndex > pattern; + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), + patternOriginal.numColumns(), + rowLengths.data() ); + + // Copy the original nonzeros + for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) + { + globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) ); + } + + // Add the nonzeros from coupling + addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); + + localMatrix.setName( this->getName() + "/matrix" ); + localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); + + rhs.setName( this->getName() + "/rhs" ); + rhs.create( numLocalRows, MPI_COMM_GEOS ); + + solution.setName( this->getName() + "/solution" ); + solution.create( numLocalRows, MPI_COMM_GEOS ); + + setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); + + // if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + // { + // createPreconditioner( domain ); + // } +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::assembleSystem( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + this->solidMechanicsSolver()->synchronizeFractureState( domain ); + + assembleElementBasedContributions( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs ); + + // Assemble fluxes 3D/2D and get dFluidResidualDAperture + this->flowSolver()->assembleHydrofracFluxTerms( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs, + getDerivativeFluxResidual_dNormalJump() ); + + // This step must occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. + assembleCouplingTerms( time_n, + dt, + domain, + dofManager, + localMatrix, + localRhs ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::assembleElementBasedContributions( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_UNUSED_VAR( time_n, dt ); + + /// 3. assemble Force Residual w.r.t. pressure and Flow mass residual w.r.t. displacement + + Base::assembleElementBasedTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); + + // Flow accumulation for fractures + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + this->flowSolver()->accumulationAssemblyLaunch( dofManager, subRegion, localMatrix, localRhs ); + } ); + } ); + + this->solidMechanicsSolver()->assembleContact( domain, dofManager, localMatrix, localRhs ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::assembleCouplingTerms( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_UNUSED_VAR( time_n, dt ); + // These 2 steps need to occur after the fluxes are assembled because that's when DerivativeFluxResidual_dAperture is filled. + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & regionNames ) + { + /// 3. assemble Force Residual w.r.t. pressure and Flow mass residual w.r.t. displacement + assembleForceResidualDerivativeWrtPressure( mesh, regionNames, dofManager, localMatrix, localRhs ); + assembleFluidMassResidualDerivativeWrtDisplacement( mesh, regionNames, dofManager, localMatrix, localRhs ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +setUpDflux_dApertureMatrix( DomainPartition & domain, + DofManager const & GEOS_UNUSED_PARAM( dofManager ), + CRSMatrix< real64, globalIndex > & localMatrix ) +{ + integer const numComp = this->flowSolver()->numFluidComponents(); + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & regionNames ) + { + std::unique_ptr< CRSMatrix< real64, localIndex > > & derivativeFluxResidual_dAperture = this->getRefDerivativeFluxResidual_dAperture(); + + { + // calculate number of fracture elements + localIndex numRows = 0; + mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, FaceElementSubRegion const & subRegion ) + { + numRows += subRegion.size(); + } ); + // number of columns (derivatives) = number of fracture elements + localIndex numCol = numRows; + // number of rows (equations) = number of fracture elements * number of components + numRows *= numComp; + + derivativeFluxResidual_dAperture = std::make_unique< CRSMatrix< real64, localIndex > >( numRows, numCol ); + derivativeFluxResidual_dAperture->setName( this->getName() + "/derivativeFluxResidual_dAperture" ); + + derivativeFluxResidual_dAperture->reserveNonZeros( localMatrix.numNonZeros() ); + localIndex maxRowSize = -1; + for( localIndex row = 0; row < localMatrix.numRows(); ++row ) + { + localIndex const rowSize = localMatrix.numNonZeros( row ); + maxRowSize = maxRowSize > rowSize ? maxRowSize : rowSize; + } + // TODO This is way too much. The With the full system rowSize is not a good estimate for this. + for( localIndex row = 0; row < numRows * numComp; ++row ) + { + derivativeFluxResidual_dAperture->reserveNonZeros( row, maxRowSize ); + } + } + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn = 0; iconn < stencil.size(); ++iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + for( localIndex k0 = 0; k0 < numFluxElems; ++k0 ) + { + for( localIndex k1 = 0; k1 < numFluxElems; ++k1 ) + { + for( integer ic = 0; ic < numComp; ic++ ) + { + derivativeFluxResidual_dAperture->insertNonZero( sei[iconn][k0] * numComp + ic, sei[iconn][k1], 0.0 ); + } + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +addTransmissibilityCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + GEOS_MARK_FUNCTION; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, // meshBodyName, + MeshLevel const & mesh, + arrayView1d< string const > const & ) // regionNames + { + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const flowDofKey = dofManager.getKey( m_flowDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & stabilizationMethod = fvManager.getFluxApproximation( this->solidMechanicsSolver()->getStabilizationName() ); + + stabilizationMethod.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + for( localIndex iconn=0; iconn( sesri[iconn][0] ); + + ArrayOfArraysView< localIndex const > const elemsToNodes = elementSubRegion.nodeList().toViewConst(); + + arrayView1d< globalIndex const > const faceElementDofNumber = + elementSubRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + for( localIndex k0=0; k0= 0 && rowNumber < rowLengths.size() ) + { + for( localIndex k1=0; k1 +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +addTransmissibilityCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const +{ + GEOS_MARK_FUNCTION; + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const flowDofKey = dofManager.getKey( m_flowDofKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + + // Get the finite volume method used to compute the stabilization + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fvDiscretization = fvManager.getFluxApproximation( this->flowSolver()->getDiscretizationName() ); + + SurfaceElementRegion const & fractureRegion = + elemManager.getRegion< SurfaceElementRegion >( this->solidMechanicsSolver()->getUniqueFractureRegionName() ); + FaceElementSubRegion const & fractureSubRegion = + fractureRegion.getUniqueSubRegion< FaceElementSubRegion >(); + + GEOS_ERROR_IF( !fractureSubRegion.hasWrapper( flow::pressure::key() ), + this->getDataContext() << ": The fracture subregion must contain pressure field." ); + + arrayView2d< localIndex const > const elem2dToFaces = fractureSubRegion.faceList().toViewConst(); + + arrayView1d< globalIndex const > const & + flowDofNumber = fractureSubRegion.getReference< globalIndex_array >( flowDofKey ); + + globalIndex const rankOffset = dofManager.rankOffset(); + + fvDiscretization.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) + { + forAll< serialPolicy >( stencil.size(), [=] ( localIndex const iconn ) + { + localIndex const numFluxElems = stencil.stencilSize( iconn ); + + // A fracture connector has to be an edge shared by two faces + if( numFluxElems == 2 ) + { + typename SurfaceElementStencil::IndexContainerViewConstType const & sei = stencil.getElementIndices(); + + // First index: face element. Second index: node + for( localIndex kf = 0; kf < 2; ++kf ) + { + // Set row DOF index + // Note that the 1-kf index is intentional, as this is coupling the pressure of one face cell + // to the nodes of the adjacent cell + localIndex const rowIndex = flowDofNumber[sei[iconn][1-kf]] - rankOffset; + + if( rowIndex >= 0 && rowIndex < pattern.numRows() ) + { + + // Get fracture, face and region/subregion/element indices (for elements on both sides) + localIndex const fractureIndex = sei[iconn][kf]; + + // Get the number of nodes + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elem2dToFaces[fractureIndex][0] ); + + // Loop over the two sides of each fracture element + for( localIndex kf1 = 0; kf1 < 2; ++kf1 ) + { + localIndex const faceIndex = elem2dToFaces[fractureIndex][kf1]; + + // Save the list of DOF associated with nodes + for( localIndex a=0; a( i ); + pattern.insertNonZero( rowIndex, colIndex ); + } + } + } + } + } + } + } ); + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +assembleForceResidualDerivativeWrtPressure( MeshLevel const & mesh, + arrayView1d< string const > const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + EdgeManager const & edgeManager = mesh.getEdgeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst(); + arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst(); + arrayView2d< real64 const > faceCenters = faceManager.faceCenter(); + arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); + arrayView1d< real64 const > faceAreas = faceManager.faceArea(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( m_flowDofKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + // Get the coordinates for all nodes + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + arrayView1d< globalIndex const > const & + flowDofNumber = subRegion.getReference< globalIndex_array >( flowDofKey ); + arrayView1d< real64 const > const & pressure = subRegion.getReference< array1d< real64 > >( flow::pressure::key() ); + arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + + forAll< serialPolicy >( subRegion.size(), [=]( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + + real64 Nbar[3]; + Nbar[ 0 ] = faceNormal[elemsToFaces[kfe][0]][0] - faceNormal[elemsToFaces[kfe][1]][0]; + Nbar[ 1 ] = faceNormal[elemsToFaces[kfe][0]][1] - faceNormal[elemsToFaces[kfe][1]][1]; + Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2]; + LvArray::tensorOps::normalize< 3 >( Nbar ); + globalIndex rowDOF[3 * m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types + real64 nodeRHS[3 * m_maxFaceNodes]; + stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes ); + globalIndex colDOF[1]; + colDOF[0] = flowDofNumber[kfe]; // pressure is always first + + for( localIndex kf=0; kf<2; ++kf ) + { + localIndex const faceIndex = elemsToFaces[kfe][kf]; + + // Compute local area contribution for each node + stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf], + nodePosition, + faceToNodeMap, + faceToEdgeMap, + edgeToNodeMap, + faceCenters, + faceNormal, + faceAreas, + nodalArea ); + for( localIndex a=0; a( globalNodalForce, Nbar, nodalForceMag ); + + for( localIndex i=0; i<3; ++i ) + { + rowDOF[3*a+i] = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i ); + // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + nodeRHS[3*a+i] = +globalNodalForce[i] * pow( -1, kf ); + + // Opposite sign w.r.t. theory because of minus sign in stiffness matrix definition (K < 0) + dRdP( 3*a+i ) = -nodalArea[a] * Nbar[i] * pow( -1, kf ); + } + } + + for( localIndex idof = 0; idof < numNodesPerFace * 3; ++idof ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRow< parallelHostAtomic >( localRow, + colDOF, + &dRdP[idof], + 1 ); + RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], nodeRHS[idof] ); + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >:: +assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, + arrayView1d< string const > const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) ) +{ + GEOS_MARK_FUNCTION; + + integer const numComp = this->flowSolver()->numFluidComponents(); + + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + EdgeManager const & edgeManager = mesh.getEdgeManager(); + ElementRegionManager const & elemManager = mesh.getElemManager(); + + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + ArrayOfArraysView< localIndex const > const faceToEdgeMap = faceManager.edgeList().toViewConst(); + arrayView2d< localIndex const > const & edgeToNodeMap = edgeManager.nodeList().toViewConst(); + arrayView2d< real64 const > faceCenters = faceManager.faceCenter(); + arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); + arrayView1d< real64 const > faceAreas = faceManager.faceArea(); + + CRSMatrixView< real64 const, localIndex const > const & + dFluxResidual_dNormalJump = getDerivativeFluxResidual_dNormalJump().toViewConst(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & flowDofKey = dofManager.getKey( m_flowDofKey ); + + arrayView1d< globalIndex const > const & + dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + // Get the coordinates for all nodes + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion const & subRegion ) + { + arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >(); + + arrayView1d< globalIndex const > const & flowDofNumber = subRegion.getReference< array1d< globalIndex > >( flowDofKey ); + + arrayView2d< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + arrayView1d< real64 const > const & area = subRegion.getElementArea().toViewConst(); + + arrayView1d< integer const > const & fractureState = subRegion.getField< fields::contact::fractureState >(); + + forAll< serialPolicy >( subRegion.size(), [&]( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0], kf1 = elemsToFaces[kfe][1]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + globalIndex nodeDOF[2*3*m_maxFaceNodes]; + + real64 Nbar[3]; + Nbar[ 0 ] = faceNormal[kf0][0] - faceNormal[kf1][0]; + Nbar[ 1 ] = faceNormal[kf0][1] - faceNormal[kf1][1]; + Nbar[ 2 ] = faceNormal[kf0][2] - faceNormal[kf1][2]; + LvArray::tensorOps::normalize< 3 >( Nbar ); + + stackArray2d< real64, 2*3*m_maxFaceNodes * MultiFluidBase::MAX_NUM_COMPONENTS > dRdU( MultiFluidBase::MAX_NUM_COMPONENTS, 2*3*m_maxFaceNodes ); + + bool const isFractureOpen = ( fractureState[kfe] == fields::contact::FractureState::Open ); + + // Accumulation derivative + if( isFractureOpen ) + { + for( localIndex kf=0; kf<2; ++kf ) + { + // Compute local area contribution for each node + stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe][kf], + nodePosition, + faceToNodeMap, + faceToEdgeMap, + edgeToNodeMap, + faceCenters, + faceNormal, + faceAreas, + nodalArea ); + + // TODO: move to something like this plus a static method. + // localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elemsToFaces[kfe][kf] ); + // stackArray1d nodalArea( numNodesPerFace ); + + for( localIndex a=0; a( i ); + real64 const dAper_dU = -pow( -1, kf ) * Nbar[i]; + for( integer ic = 0; ic < numComp; ic++ ) + { + dRdU[ ic ][ kf*3*numNodesPerFace + 3*a + i ] = dVolume_dAperture * dAper_dU * compDens[kfe][ic]; // assuming poro=1 + } + } + } + } + + localIndex const localRow = LvArray::integerConversion< localIndex >( flowDofNumber[kfe] - rankOffset ); + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + integer const numRows = numComp; + for( integer i = 0; i < numRows; ++i ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow+i, + nodeDOF, + dRdU[i], + 2 * 3 * numNodesPerFace ); + } + } + } + + // flux derivative + bool skipAssembly = true; + for( integer ic = 0; ic < numComp; ic++ ) + { + localIndex const numColumns = dFluxResidual_dNormalJump.numNonZeros( kfe * numComp + ic ); + arraySlice1d< localIndex const > const & columns = dFluxResidual_dNormalJump.getColumns( kfe * numComp + ic ); + // TODO uncomment and fix build + arraySlice1d< real64 const > const & values = dFluxResidual_dNormalJump.getEntries( kfe * numComp + ic ); + + skipAssembly &= !isFractureOpen; + + for( localIndex kfe1 = 0; kfe1 < numColumns; ++kfe1 ) + { + real64 const dR_dAper = values[kfe1]; + localIndex const kfe2 = columns[kfe1]; + + bool const isOpen = ( fractureState[kfe2] == fields::contact::FractureState::Open ); + skipAssembly &= !isOpen; + + for( localIndex kf = 0; kf < 2; ++kf ) + { + //TODO: We should avoid allocating LvArrays inside kernel + stackArray1d< real64, FaceManager::maxFaceNodes() > nodalArea; + this->solidMechanicsSolver()->computeFaceNodalArea( elemsToFaces[kfe2][kf], + nodePosition, + faceToNodeMap, + faceToEdgeMap, + edgeToNodeMap, + faceCenters, + faceNormal, + faceAreas, + nodalArea ); + + for( localIndex a=0; a( i ); + real64 const dAper_dU = -pow( -1, kf ) * Nbar[i] * ( nodalArea[a] / area[kfe2] ); + dRdU[ ic ][ kf*3*numNodesPerFace + 3*a + i ] = dR_dAper * dAper_dU; + } + } + } + } + + if( !skipAssembly ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( flowDofNumber[kfe] - rankOffset ); + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + integer const numRows = numComp; + for( integer i = 0; i < numRows; ++i ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow+i, + nodeDOF, + dRdU[i], + 2 * 3 * numNodesPerFace ); + } + } + } + } + } ); + } ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::updateState( DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + + // call base poromechanics update + Base::updateState( domain ); + // need to call solid mechanics update separately to compute face displacement jump + this->solidMechanicsSolver()->updateState( domain ); + + // remove the contribution of the hydraulic aperture from the stencil weights + this->flowSolver()->prepareStencilWeights( domain ); + + updateHydraulicApertureAndFracturePermeability( domain ); + + // update the stencil weights using the updated hydraulic aperture + this->flowSolver()->updateStencilWeights( domain ); +} + +template< typename FLOW_SOLVER > +void MultiphasePoromechanicsConformingFractures< FLOW_SOLVER >::updateHydraulicApertureAndFracturePermeability( DomainPartition & domain ) +{ + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion & subRegion ) + { + arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + arrayView1d< real64 const > const area = subRegion.getElementArea(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView2d< real64 const > const fractureTraction = subRegion.getField< fields::contact::traction >(); + arrayView1d< real64 const > const pressure = subRegion.getField< fields::flow::pressure >(); + arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< fields::flow::aperture0 >(); + + arrayView1d< real64 > const aperture = subRegion.getElementAperture(); + arrayView1d< real64 > const hydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); + arrayView1d< real64 > const deltaVolume = subRegion.getField< flow::deltaVolume >(); + + string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ); + CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( porousSolidName ); + + string const & hydraulicApertureRelationName = subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() ); + HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName ); + + constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicAperture ) + { + using HydraulicApertureType = TYPEOFREF( castedHydraulicAperture ); + typename HydraulicApertureType::KernelWrapper hydraulicApertureWrapper = castedHydraulicAperture.createKernelWrapper(); + + constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) + { + typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates(); + + poromechanicsFracturesKernels::StateUpdateKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + porousMaterialWrapper, + hydraulicApertureWrapper, + dispJump, + pressure, + area, + volume, + deltaVolume, + aperture, + oldHydraulicAperture, + hydraulicAperture, + fractureTraction ); + + } ); + } ); + } ); + } ); +} + +template class MultiphasePoromechanicsConformingFractures<>; +//template class MultiphasePoromechanicsConformingFractures< MultiphaseReservoirAndWells<> >; + +namespace +{ +//typedef MultiphasePoromechanicsConformingFractures< MultiphaseReservoirAndWells<> > +// MultiphasePoromechanicsConformingFractures; +//REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphasePoromechanicsConformingFractures, string const &, Group * const ) +typedef MultiphasePoromechanicsConformingFractures<> MultiphasePoromechanicsConformingFractures; +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphasePoromechanicsConformingFractures, string const &, Group * const ) +} + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp new file mode 100644 index 00000000000..cd8d7c0073b --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp @@ -0,0 +1,203 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultiphasePoromechanicsConformingFractures.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ + +#include "physicsSolvers/multiphysics/MultiphasePoromechanics.hpp" +#include "physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp" + +namespace geos +{ + +template< typename FLOW_SOLVER = CompositionalMultiphaseBase > +class MultiphasePoromechanicsConformingFractures : public MultiphasePoromechanics< FLOW_SOLVER, SolidMechanicsLagrangeContact > +{ +public: + + using Base = MultiphasePoromechanics< FLOW_SOLVER, SolidMechanicsLagrangeContact >; + using Base::m_solvers; + using Base::m_dofManager; + using Base::m_localMatrix; + using Base::m_rhs; + using Base::m_solution; + + /// String used to form the solverName used to register solvers in CoupledSolver + static string coupledSolverAttributePrefix() { return "poromechanicsConformingFractures"; } + + /** + * @brief main constructor for MultiphasePoromechanicsConformingFractures objects + * @param name the name of this instantiation of MultiphasePoromechanicsConformingFractures in the repository + * @param parent the parent group of this instantiation of MultiphasePoromechanicsConformingFractures + */ + MultiphasePoromechanicsConformingFractures( const string & name, + dataRepository::Group * const parent ); + + /// Destructor for the class + ~MultiphasePoromechanicsConformingFractures() override {} + + /** + * @brief name of the node manager in the object catalog + * @return string that contains the catalog name to generate a new SinglePhasePoromechanicsConformingFractures object through the object + * catalog. + */ + static string catalogName() + { + if constexpr ( std::is_same_v< FLOW_SOLVER, CompositionalMultiphaseBase > ) + { + return "MultiphasePoromechanicsConformingFractures"; + } + else + { + return FLOW_SOLVER::catalogName() + "PoromechanicsConformingFractures"; + } + } + + /** + * @copydoc PhysicsSolverBase::getCatalogName() + */ + string getCatalogName() const override { return catalogName(); } + + /** + * @defgroup Solver Interface Functions + * + * These functions provide the primary interface that is required for derived classes + */ + /**@{*/ + + virtual void setupCoupling( DomainPartition const & domain, + DofManager & dofManager ) const override; + + virtual void setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity = true ) override; + + virtual void assembleSystem( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override final; + + virtual void updateState( DomainPartition & domain ) override final; + + /**@}*/ + +protected: + + virtual void postInputInitialization() override; + +private: + + struct viewKeyStruct : public Base::viewKeyStruct + {}; + + static const localIndex m_maxFaceNodes=11; // Maximum number of nodes on a contact face + + void assembleElementBasedContributions( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + virtual void assembleCouplingTerms( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override final; + + void assembleForceResidualDerivativeWrtPressure( MeshLevel const & mesh, + arrayView1d< string const > const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, + arrayView1d< string const > const & regionNames, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + /** + * @Brief add the nnz induced by the flux-aperture coupling + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param rowLenghts the nnz in each row + */ + void addTransmissibilityCouplingNNZ( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< localIndex > const & rowLengths ) const; + + /** + * @Brief add the sparsity pattern induced by the flux-aperture coupling + * @param domain the physical domain object + * @param dofManager degree-of-freedom manager associated with the linear system + * @param pattern the sparsity pattern + */ + void addTransmissibilityCouplingPattern( DomainPartition const & domain, + DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const; + + /** + * @brief Set up the Dflux_dApertureMatrix object + * + * @param domain + * @param dofManager + * @param localMatrix + */ + void setUpDflux_dApertureMatrix( DomainPartition & domain, + DofManager const & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix ); + + /** + * @brief + * + * @param domain + */ + void updateHydraulicApertureAndFracturePermeability( DomainPartition & domain ); + + + std::unique_ptr< CRSMatrix< real64, localIndex > > & getRefDerivativeFluxResidual_dAperture() + { + return m_derivativeFluxResidual_dAperture; + } + + CRSMatrixView< real64, localIndex const > getDerivativeFluxResidual_dNormalJump() + { + return m_derivativeFluxResidual_dAperture->toViewConstSizes(); + } + + CRSMatrixView< real64 const, localIndex const > getDerivativeFluxResidual_dNormalJump() const + { + return m_derivativeFluxResidual_dAperture->toViewConst(); + } + + std::unique_ptr< CRSMatrix< real64, localIndex > > m_derivativeFluxResidual_dAperture; + + string const m_flowDofKey = CompositionalMultiphaseBase::viewKeyStruct::elemDofFieldString(); + +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index 02fa02be59b..11ec7caa5d4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -24,6 +24,7 @@ #include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsStatistics.hpp" #include "physicsSolvers/multiphysics/MultiphasePoromechanics.hpp" +#include "physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp" #include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp" #include "physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp" @@ -136,6 +137,7 @@ execute( real64 const time_n, namespace { typedef PoromechanicsInitialization< MultiphasePoromechanics<> > MultiphasePoromechanicsInitialization; +typedef PoromechanicsInitialization< MultiphasePoromechanicsConformingFractures<> > MultiphasePoromechanicsConformingFracturesInitialization; typedef PoromechanicsInitialization< MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> > > MultiphaseReservoirPoromechanicsInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanics<> > SinglePhasePoromechanicsInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFractures<> > SinglePhasePoromechanicsConformingFracturesInitialization; @@ -144,6 +146,7 @@ typedef PoromechanicsInitialization< SinglePhasePoromechanicsEmbeddedFractures > typedef PoromechanicsInitialization< SinglePhasePoromechanics< SinglePhaseReservoirAndWells<> > > SinglePhaseReservoirPoromechanicsInitialization; typedef PoromechanicsInitialization< HydrofractureSolver< SinglePhasePoromechanics<> > > HydrofractureInitialization; REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsConformingFracturesInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, MultiphaseReservoirPoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsConformingFracturesInitialization, string const &, Group * const ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 62a772e0278..9db368c00b7 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -21,10 +21,10 @@ #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSSOLVER_HPP_ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSSOLVER_HPP_ -#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/multiphysics/CoupledSolver.hpp" -#include "physicsSolvers/multiphysics/PoromechanicsFields.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/multiphysics/PoromechanicsFields.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" #include "constitutive/solid/PorousSolid.hpp" #include "constitutive/contact/HydraulicApertureBase.hpp" diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp index 0e657441f11..1965dc5a457 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp @@ -21,10 +21,7 @@ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP_ #include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" -#include "physicsSolvers/multiphysics/CoupledSolver.hpp" #include "physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" -#include "dataRepository/Group.hpp" namespace geos { diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index 6c6e2be5bc4..f25329cfe14 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -95,17 +95,10 @@ void SinglePhasePoromechanicsEmbeddedFractures::initializePostInitialConditionsP updateState( this->getGroupByPath< DomainPartition >( "/Problem/domain" ) ); } -void SinglePhasePoromechanicsEmbeddedFractures::setupDofs( DomainPartition const & domain, - DofManager & dofManager ) const +void SinglePhasePoromechanicsEmbeddedFractures::setupCoupling( DomainPartition const & domain, + DofManager & dofManager ) const { - GEOS_MARK_FUNCTION; - solidMechanicsSolver()->setupDofs( domain, dofManager ); - flowSolver()->setupDofs( domain, dofManager ); - - // Add coupling between displacement and cell pressures - dofManager.addCoupling( fields::solidMechanics::totalDisplacement::key(), - SinglePhaseBase::viewKeyStruct::elemDofFieldString(), - DofManager::Connector::Elem ); + Base::setupCoupling( domain, dofManager ); map< std::pair< string, string >, array1d< string > > meshTargets; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp index 95c3e7b7f82..f90d4ffe0be 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp @@ -58,8 +58,8 @@ class SinglePhasePoromechanicsEmbeddedFractures : public SinglePhasePoromechanic bool const setSparsity = true ) override; virtual void - setupDofs( DomainPartition const & domain, - DofManager & dofManager ) const override; + setupCoupling( DomainPartition const & domain, + DofManager & dofManager ) const override; virtual void assembleSystem( real64 const time, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp new file mode 100644 index 00000000000..33b07aaaaa1 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanicsConformingFractures.hpp @@ -0,0 +1,410 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file MultiphasePoromechanicsConformingFractures.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP + +#include "physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp" +//#include "physicsSolvers/fluidFlow/FluxKernelsHelper.hpp" + +namespace geos +{ + +namespace multiphasePoromechanicsConformingFracturesKernels +{ + +template< integer NUM_EQN, integer NUM_DOF > +class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper > +{ +public: + + /** + * @brief The type for element-based data. Consists entirely of ArrayView's. + * + * Can be converted from ElementRegionManager::ElementViewConstAccessor + * by calling .toView() or .toViewConst() on an accessor instance + */ + template< typename VIEWTYPE > + using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + + using AbstractBase = isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernelBase; + using DofNumberAccessor = AbstractBase::DofNumberAccessor; + using CompFlowAccessors = AbstractBase::CompFlowAccessors; + using MultiFluidAccessors = AbstractBase::MultiFluidAccessors; + using CapPressureAccessors = AbstractBase::CapPressureAccessors; + using PermeabilityAccessors = AbstractBase::PermeabilityAccessors; + using FracturePermeabilityAccessors = StencilMaterialAccessors< constitutive::PermeabilityBase, + fields::permeability::dPerm_dDispJump >; + + using AbstractBase::m_dt; + using AbstractBase::m_dofNumber; + using AbstractBase::m_gravCoef; + using AbstractBase::m_pres; + + using Base = isothermalCompositionalMultiphaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper >; + using Base::numDof; + using Base::numEqn; + using Base::maxNumElems; + using Base::maxNumConns; + using Base::maxStencilSize; + using Base::m_stencilWrapper; + using Base::m_seri; + using Base::m_sesri; + using Base::m_sei; + using Base::numFluxSupportPoints; + using Base::numComp; + using Base::m_numPhases; + using Base::m_kernelFlags; + + using Base::m_permeability; + using Base::m_dPerm_dPres; + using Base::m_phaseMob; + using Base::m_dPhaseMob; + using Base::m_phaseVolFrac; + using Base::m_dPhaseVolFrac; + using Base::m_phaseCompFrac; + using Base::m_dPhaseCompFrac; + using Base::m_dCompFrac_dCompDens; + using Base::m_phaseMassDens; + using Base::m_dPhaseMassDens; + using Base::m_phaseCapPressure; + using Base::m_dPhaseCapPressure_dPhaseVolFrac; + + FluxComputeKernel( integer const numPhases, + globalIndex const rankOffset, + SurfaceElementStencilWrapper const & stencilWrapper, + DofNumberAccessor const & dofNumberAccessor, + CompFlowAccessors const & compFlowAccessors, + MultiFluidAccessors const & multiFluidAccessors, + CapPressureAccessors const & capPressureAccessors, + PermeabilityAccessors const & permeabilityAccessors, + FracturePermeabilityAccessors const & fracturePermeabilityAccessors, + real64 const dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + : Base( numPhases, + rankOffset, + stencilWrapper, + dofNumberAccessor, + compFlowAccessors, + multiFluidAccessors, + capPressureAccessors, + permeabilityAccessors, + dt, + localMatrix, + localRhs, + kernelFlags ), + m_dR_dAper( dR_dAper ), + m_dPerm_dDispJump( fracturePermeabilityAccessors.get( fields::permeability::dPerm_dDispJump {} ) ) + {} + + + /** + * @struct StackVariables + * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack + */ + struct StackVariables : public Base::StackVariables + { +public: + + /** + * @brief Constructor for the stack variables + * @param[in] size size of the stencil for this connection + * @param[in] numElems number of elements for this connection + */ + GEOS_HOST_DEVICE + StackVariables( localIndex const size, localIndex numElems ) + : Base::StackVariables( size, numElems ), + localColIndices( numElems ), + // TODO + dFlux_dAperture( numComp, numElems, size ) + {} + + stackArray1d< localIndex, maxNumElems > localColIndices; + + // TODO + stackArray3d< real64, numEqn * maxNumElems * maxStencilSize > dFlux_dAperture; + + /// Derivatives of transmissibility with respect to the dispJump + real64 dTrans_dDispJump[maxNumConns][2][3]{}; + }; + + /** + * @brief Performs the setup phase for the kernel. + * @param[in] iconn the connection index + * @param[in] stack the stack variables + */ + GEOS_HOST_DEVICE + void setup( localIndex const iconn, + StackVariables & stack ) const + { + // set degrees of freedom indices for this face + for( integer i = 0; i < stack.stencilSize; ++i ) + { + globalIndex const offset = m_dofNumber[m_seri( iconn, i )][m_sesri( iconn, i )][m_sei( iconn, i )]; + for( integer jdof = 0; jdof < numDof; ++jdof ) + { + stack.dofColIndices[i * numDof + jdof] = offset + jdof; + } + stack.localColIndices[ i ] = m_sei( iconn, i ); + } + } + + /** + * @brief Compute the local flux contributions to the residual and Jacobian + * @tparam FUNC the type of the function that can be used to customize the computation of the flux + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + * @param[in] NoOpFunc the function used to customize the computation of the flux + */ + template< typename FUNC = NoOpFunc > + GEOS_HOST_DEVICE + void computeFlux( localIndex const iconn, + StackVariables & stack, + FUNC && compFluxKernelOp = NoOpFunc{} ) const + { + m_stencilWrapper.computeWeights( iconn, + m_permeability, + m_dPerm_dPres, + m_dPerm_dDispJump, + stack.transmissibility, + stack.dTrans_dPres, + stack.dTrans_dDispJump ); + + localIndex k[numFluxSupportPoints]; + localIndex connectionIndex = 0; + for( k[0] = 0; k[0] < stack.numConnectedElems; ++k[0] ) + { + for( k[1] = k[0] + 1; k[1] < stack.numConnectedElems; ++k[1] ) + { + /// cell indices + localIndex const seri[numFluxSupportPoints] = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )}; + localIndex const sesri[numFluxSupportPoints] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )}; + localIndex const sei[numFluxSupportPoints] = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )}; + + // clear working arrays + real64 compFlux[numComp]{}; + real64 dCompFlux_dP[numFluxSupportPoints][numComp]{}; + real64 dCompFlux_dC[numFluxSupportPoints][numComp][numComp]{}; + real64 dCompFlux_dTrans[numComp]{}; + + real64 const trans[numFluxSupportPoints] = { stack.transmissibility[connectionIndex][0], + stack.transmissibility[connectionIndex][1] }; + + real64 const dTrans_dPres[numFluxSupportPoints] = { stack.dTrans_dPres[connectionIndex][0], + stack.dTrans_dPres[connectionIndex][1] }; + + //***** calculation of flux ***** + // loop over phases, compute and upwind phase flux and sum contributions to each component's flux + for( integer ip = 0; ip < m_numPhases; ++ip ) + { + // create local work arrays + real64 potGrad = 0.0; + real64 phaseFlux = 0.0; + real64 dPhaseFlux_dP[numFluxSupportPoints]{}; + real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{}; + + localIndex k_up = -1; + + + isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux::compute< numComp, numFluxSupportPoints > + ( m_numPhases, + ip, + m_kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CapPressure ), + m_kernelFlags.isSet( isothermalCompositionalMultiphaseFVMKernels::KernelFlags::CheckPhasePresenceInGravity ), + seri, sesri, sei, + trans, + dTrans_dPres, + m_pres, + m_gravCoef, + m_phaseMob, m_dPhaseMob, + m_phaseVolFrac, + m_dPhaseVolFrac, + m_phaseCompFrac, m_dPhaseCompFrac, + m_dCompFrac_dCompDens, + m_phaseMassDens, m_dPhaseMassDens, + m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac, + k_up, + potGrad, + phaseFlux, + dPhaseFlux_dP, + dPhaseFlux_dC, + compFlux, + dCompFlux_dP, + dCompFlux_dC, + dCompFlux_dTrans ); + + // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives + // possible use: assemble the derivatives wrt temperature, and the flux term of the energy equation for this phase + compFluxKernelOp( ip, k, seri, sesri, sei, connectionIndex, + k_up, seri[k_up], sesri[k_up], sei[k_up], potGrad, + phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + + } // loop over phases + + /// populate local flux vector and derivatives + for( integer ic = 0; ic < numComp; ++ic ) + { + integer const eqIndex0 = k[0] * numEqn + ic; + integer const eqIndex1 = k[1] * numEqn + ic; + + stack.localFlux[eqIndex0] += m_dt * compFlux[ic]; + stack.localFlux[eqIndex1] -= m_dt * compFlux[ic]; + + for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) + { + localIndex const localDofIndexPres = k[ke] * numDof; + stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dCompFlux_dP[ke][ic]; + stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dCompFlux_dP[ke][ic]; + + for( integer jc = 0; jc < numComp; ++jc ) + { + localIndex const localDofIndexComp = localDofIndexPres + jc + 1; + stack.localFluxJacobian[eqIndex0][localDofIndexComp] += m_dt * dCompFlux_dC[ke][ic][jc]; + stack.localFluxJacobian[eqIndex1][localDofIndexComp] -= m_dt * dCompFlux_dC[ke][ic][jc]; + } + } + } + + for( integer ic = 0; ic < numComp; ++ic ) + { + real64 dFlux_dAper[2] = {0.0, 0.0}; + dFlux_dAper[0] = m_dt * dCompFlux_dTrans[ic] * stack.dTrans_dDispJump[connectionIndex][0][0]; + dFlux_dAper[1] = -m_dt * dCompFlux_dTrans[ic] * stack.dTrans_dDispJump[connectionIndex][1][0]; + + stack.dFlux_dAperture[ic][k[0]][k[0]] += dFlux_dAper[0]; + stack.dFlux_dAperture[ic][k[0]][k[1]] += dFlux_dAper[1]; + stack.dFlux_dAperture[ic][k[1]][k[0]] -= dFlux_dAper[0]; + stack.dFlux_dAperture[ic][k[1]][k[1]] -= dFlux_dAper[1]; + } +// + connectionIndex++; + } + } + } + + /** + * @brief Performs the complete phase for the kernel. + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + */ + template< typename FUNC = NoOpFunc > + GEOS_HOST_DEVICE + void complete( localIndex const iconn, + StackVariables & stack, + FUNC && kernelOp = NoOpFunc{} ) const + { + // Call Base::complete to assemble the mass balance equations + // In the lambda, fill the dR_dAper matrix + Base::complete( iconn, stack, [&] ( integer const i, + localIndex const localRow ) + { + // TODO + localIndex const ei = LvArray::integerConversion< localIndex >( m_sei( iconn, i ) ); + for( integer ic = 0; ic < numComp; ++ic ) + { + localIndex const row = ei * numComp + ic; + + m_dR_dAper.addToRowBinarySearch< parallelDeviceAtomic >( row, + stack.localColIndices.data(), + stack.dFlux_dAperture[ic][i].dataIfContiguous(), + stack.stencilSize ); + } + + // call the lambda to assemble additional terms, such as thermal terms + kernelOp( i, localRow ); + } ); + } + +private: + + CRSMatrixView< real64, localIndex const > m_dR_dAper; + + ElementViewConst< arrayView4d< real64 const > > const m_dPerm_dDispJump; +}; + + + +/** + * @class FluxComputeKernelFactory + */ +class FluxComputeKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @tparam STENCILWRAPPER the type of the stencil wrapper + * @param[in] rankOffset the offset of my MPI rank + * @param[in] dofKey string to get the element degrees of freedom numbers + * @param[in] solverName name of the solver (to name accessors) + * @param[in] elemManager reference to the element region manager + * @param[in] stencilWrapper reference to the stencil wrapper + * @param[in] dt time step size + * @param[inout] localMatrix the local CRS matrix + * @param[inout] localRhs the local right-hand side vector + */ + template< typename POLICY > + static void + createAndLaunch( integer const numComps, + integer const numPhases, + globalIndex const rankOffset, + string const & dofKey, + BitFlags< isothermalCompositionalMultiphaseFVMKernels::KernelFlags > kernelFlags, + string const & solverName, + ElementRegionManager const & elemManager, + SurfaceElementStencilWrapper const & stencilWrapper, + real64 const dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + { + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&]( auto NC ) + { + integer constexpr NUM_COMP = NC(); + integer constexpr NUM_DOF = NC() + 1; + + ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > dofNumberAccessor = + elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); + dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + + using kernelType = FluxComputeKernel< NUM_COMP, NUM_DOF >; + typename kernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); + typename kernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); + typename kernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName ); + typename kernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); + typename kernelType::FracturePermeabilityAccessors fracPermAccessors( elemManager, solverName ); + + kernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor, + compFlowAccessors, multiFluidAccessors, capPressureAccessors, permeabilityAccessors, fracPermAccessors, + dt, localMatrix, localRhs, kernelFlags, dR_dAper ); + + kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); + } ); + } +}; + +} // namespace multiphasePoromechanicsConformingFracturesKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_MULTIPHASEPOROMECHANICSCONFORMINGFRACTURES_HPP diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp index 302c8113b3d..00350d1f905 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp @@ -22,6 +22,7 @@ #include "physicsSolvers/multiphysics/PoromechanicsFields.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsBase.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" namespace geos {