From c6ef7ca3fb4fa3e069ad54263ac2408a5d04e82e Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 18 Jun 2026 16:04:38 -0500 Subject: [PATCH] enble thermal effect for ALM solver --- .../fluidFlow/FlowSolverBase.hpp | 12 + .../fluidFlow/SinglePhaseFVM.cpp | 88 ++++ .../fluidFlow/SinglePhaseFVM.hpp | 9 + .../multiphysics/CMakeLists.txt | 2 + ...asePoromechanicsConformingFracturesALM.cpp | 156 ++++++- ...asePoromechanicsConformingFracturesALM.hpp | 9 + ...asePoromechanicsConformingFracturesALM.hpp | 330 +++++++++++++++ ...asePoromechanicsConformingFracturesALM.hpp | 391 ++++++++++++++++++ 8 files changed, 984 insertions(+), 13 deletions(-) create mode 100644 src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp create mode 100644 src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 75afb4fcd65..1328bacc12a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -169,6 +169,18 @@ class FlowSolverBase : public PhysicsSolverBase GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." ); } + virtual void assembleHydrofracFluxTermsALM( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + { + GEOS_UNUSED_VAR ( time_n, dt, domain, dofManager, localMatrix, localRhs, dR_dAper ); + GEOS_ERROR( "Poroelastic fluxes with conforming fractures ALM not yet implemented." ); + } + void initializeState( DomainPartition & domain ); virtual void initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index b507ea1b22c..f8d57fcf269 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -45,6 +45,8 @@ #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp" /** * @namespace the geos namespace that encapsulates the majority of the code @@ -651,6 +653,92 @@ void SinglePhaseFVM< BASE >::assembleHydrofracFluxTerms( real64 const GEOS_UNUSE } +template< typename BASE > +void SinglePhaseFVM< BASE >::assembleHydrofracFluxTermsALM( real64 const GEOS_UNUSED_PARAM ( time_n ), + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) +{ + GEOS_MARK_FUNCTION; + + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName ); + + string const & dofKey = dofManager.getKey( SinglePhaseBase::viewKeyStruct::elemDofFieldString() ); + + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel const & mesh, + string_array const & ) + { + fluxApprox.forStencils< CellElementStencilTPFA, FaceElementToCellStencil >( mesh, [&]( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + if( m_isThermal ) + { + thermalSinglePhaseFVMKernels:: + FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), + dofKey, + this->getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } + else + { + singlePhaseFVMKernels:: + FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), + dofKey, + this->getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } + } ); + + fluxApprox.forStencils< SurfaceElementStencil >( mesh, [&]( auto & stencil ) + { + typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); + + if( m_isThermal ) + { + thermalSinglePhasePoromechanicsConformingFracturesALMKernels:: + ConnectorBasedAssemblyKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), + dofKey, + this->getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView(), + dR_dAper ); + } + else + { + singlePhasePoromechanicsConformingFracturesALMKernels:: + ConnectorBasedAssemblyKernelFactory::createAndLaunch< parallelDevicePolicy<> >( dofManager.rankOffset(), + dofKey, + this->getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView(), + dR_dAper ); + } + } ); + } ); +} + template< typename BASE > void SinglePhaseFVM< BASE >::applyBoundaryConditions( real64 const time_n, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp index a85545c0c43..aeaeae61bce 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp @@ -180,6 +180,15 @@ class SinglePhaseFVM : public BASE arrayView1d< real64 > const & localRhs, CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; + virtual void + assembleHydrofracFluxTermsALM( real64 const time_n, + real64 const dt, + DomainPartition const & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; + /**@}*/ virtual void diff --git a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt index 36e8a23c9c0..684ab99233a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt @@ -45,6 +45,7 @@ set( multiPhysicsSolvers_headers poromechanicsKernels/SinglePhasePoromechanics.hpp poromechanicsKernels/SinglePhasePoromechanics_impl.hpp poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp + poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp poromechanicsKernels/SinglePhasePoromechanicsDamage.hpp poromechanicsKernels/SinglePhasePoromechanicsDamage_impl.hpp poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp @@ -58,6 +59,7 @@ set( multiPhysicsSolvers_headers poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFractures.hpp + poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp SinglePhasePoromechanics.hpp SinglePhasePoromechanicsEmbeddedFractures.hpp diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index ff7d9662a60..ebc5881ee2f 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -19,7 +19,9 @@ #include "SinglePhasePoromechanicsConformingFracturesALM.hpp" +#include "constitutive/contact/HydraulicApertureRelationSelector.hpp" #include "finiteVolume/FluxApproximationBase.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" namespace geos { @@ -34,15 +36,28 @@ SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::SinglePhasePorome : Base( name, parent ) {} +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::postInputInitialization() +{ + Base::postInputInitialization(); + + GEOS_WARNING_IF( this->getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::FullyImplicit, + "FullyImplicit coupling not implemented for this solver. A sequential coupling approach will be used." ); + + this->getNonlinearSolverParameters().m_couplingType = NonlinearSolverParameters::CouplingType::Sequential; +} + template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupCoupling( DomainPartition const & domain, DofManager & dofManager ) const { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain, dofManager ); + Base::setupCoupling( domain, dofManager ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + dofManager.addCoupling( this->getFlowDofKey(), + fields::contact::traction::key(), + DofManager::Connector::Elem ); } @@ -60,7 +75,10 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( GEOS_UNUSED_VAR( domain, dofManager, localMatrix, rhs, solution, setSparsity ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } } @@ -77,7 +95,10 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleSyst GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } } template< typename FLOW_SOLVER > @@ -92,7 +113,10 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleElem GEOS_UNUSED_VAR( time_n, dt, domain, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } } @@ -107,15 +131,23 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleCoup GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( domain, dofManager, localMatrix, localRhs, time_n, dt ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } } template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::updateState( DomainPartition & domain ) { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( domain ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + + Base::updateState( domain ); + this->solidMechanicsSolver()->updateState( domain ); + + this->flowSolver()->prepareStencilWeights( domain ); + updateHydraulicApertureAndFracturePermeability( domain ); + this->flowSolver()->updateStencilWeights( domain ); } template< typename FLOW_SOLVER > @@ -125,7 +157,11 @@ setUpDflux_dApertureMatrix( DomainPartition & domain, CRSMatrix< real64, globalIndex > & localMatrix ) { GEOS_UNUSED_VAR( domain, localMatrix ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } } template< typename FLOW_SOLVER > @@ -137,7 +173,11 @@ addTransmissibilityCouplingNNZ( DomainPartition const & domain, GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( domain, dofManager, rowLengths ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } } @@ -150,7 +190,11 @@ addTransmissibilityCouplingPattern( DomainPartition const & domain, GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( domain, dofManager, pattern ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } } template< typename FLOW_SOLVER > @@ -165,7 +209,11 @@ assembleForceResidualDerivativeWrtPressure( string const & meshName, GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( meshName, mesh, regionNames, dofManager, localMatrix, localRhs ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } } template< typename FLOW_SOLVER > @@ -179,7 +227,89 @@ assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( mesh, regionNames, dofManager, localMatrix ); - GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + + if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) + { + GEOS_ERROR( "SinglePhasePoromechanicsConformingFracturesALM does not support FullyImplicit coupling type." ); + } + +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::mapSolutionBetweenSolvers( DomainPartition & domain, + integer const solverType ) +{ + GEOS_MARK_FUNCTION; + + if( solverType == static_cast< integer >( Base::SolverType::SolidMechanics ) + && !this->m_performStressInitialization ) + { + this->flowSolver()->prepareStencilWeights( domain ); + updateHydraulicApertureAndFracturePermeability( domain ); + this->flowSolver()->updateStencilWeights( domain ); + } + + Base::mapSolutionBetweenSolvers( domain, solverType ); +} + +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::updateHydraulicApertureAndFracturePermeability( DomainPartition & domain ) +{ + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, + [&]( localIndex const, + FaceElementSubRegion & subRegion ) + { + arrayView2d< real64 const > const dispJump = subRegion.getField< fields::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< fields::flow::hydraulicAperture >(); + arrayView1d< real64 > const deltaVolume = subRegion.getField< fields::flow::deltaVolume >(); + arrayView1d< integer > const & fractureState = subRegion.getField< fields::contact::fractureState >(); + + string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ); + CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( porousSolidName ); + + string const & hydraulicApertureRelationName = subRegion.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(); + + 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, + fractureState ); + } ); + } ); + } ); + } ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index 15f24d838d9..4267619df0a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -183,6 +183,11 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec DofManager const & dofManager, CRSMatrix< real64, globalIndex > & localMatrix ); + virtual void mapSolutionBetweenSolvers( DomainPartition & domain, + integer const solverType ) override; + + void updateHydraulicApertureAndFracturePermeability( DomainPartition & domain ); + std::unique_ptr< CRSMatrix< real64, localIndex > > & getRefDerivativeFluxResidual_dAperture() { return m_derivativeFluxResidual_dAperture; @@ -202,6 +207,10 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec string const m_pressureKey = SinglePhaseBase::viewKeyStruct::elemDofFieldString(); +protected: + + virtual void postInputInitialization() override final; + }; } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp new file mode 100644 index 00000000000..d4ea015b384 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -0,0 +1,330 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SinglePhasePoromechanicsConformingFracturesALM.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP + +#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp" +#include "codingUtilities/Utilities.hpp" + +namespace geos +{ + +namespace singlePhasePoromechanicsConformingFracturesALMKernels +{ + +template< integer NUM_EQN, integer NUM_DOF > +class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::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 = singlePhaseFVMKernels::FluxComputeKernelBase; + using DofNumberAccessor = AbstractBase::DofNumberAccessor; + using SinglePhaseFlowAccessors = AbstractBase::SinglePhaseFlowAccessors; + using SinglePhaseFluidAccessors = AbstractBase::SinglePhaseFluidAccessors; + using PermeabilityAccessors = AbstractBase::PermeabilityAccessors; + using FracturePermeabilityAccessors = StencilMaterialAccessors< constitutive::PermeabilityBase, + fields::permeability::dPerm_dDispJump >; + + using AbstractBase::m_dt; + using AbstractBase::m_rankOffset; + using AbstractBase::m_dofNumber; + using AbstractBase::m_permeability; + using AbstractBase::m_dPerm_dPres; + using AbstractBase::m_gravCoef; + using AbstractBase::m_pres; + using AbstractBase::m_mob; + using AbstractBase::m_dMob; + using AbstractBase::m_dMob_dPres; + using AbstractBase::m_dens; + using AbstractBase::m_dDens; + + + using Base = singlePhaseFVMKernels::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; + + ConnectorBasedAssemblyKernel( globalIndex const rankOffset, + SurfaceElementStencilWrapper const & stencilWrapper, + DofNumberAccessor const & flowDofNumberAccessor, + SinglePhaseFlowAccessors const & singlePhaseFlowAccessors, + SinglePhaseFluidAccessors const & singlePhaseFluidAccessors, + PermeabilityAccessors const & permeabilityAccessors, + FracturePermeabilityAccessors const & fracturePermeabilityAccessors, + real64 const & dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + : Base( rankOffset, + stencilWrapper, + flowDofNumberAccessor, + singlePhaseFlowAccessors, + singlePhaseFluidAccessors, + permeabilityAccessors, + dt, + localMatrix, + localRhs ), + 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 ), + dFlux_dAperture( numElems, size ) + {} + + stackArray1d< localIndex, maxNumElems > localColIndices; + + stackArray2d< real64, 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 && kernelOp = NoOpFunc{} ) const + + { + + m_stencilWrapper.computeWeights( iconn, + m_permeability, + m_dPerm_dPres, + m_dPerm_dDispJump, + stack.transmissibility, + stack.dTrans_dPres, + stack.dTrans_dDispJump ); + + + localIndex k[2]; + localIndex connectionIndex = 0; + for( k[0]=0; k[0] + 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 ) + { + + localIndex const row = LvArray::integerConversion< localIndex >( m_sei( iconn, i ) ); + + m_dR_dAper.addToRowBinarySearch< parallelDeviceAtomic >( row, + stack.localColIndices.data(), + stack.dFlux_dAperture[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 ConnectorBasedAssemblyKernelFactory + */ +class ConnectorBasedAssemblyKernelFactory +{ +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( globalIndex const rankOffset, + string const & dofKey, + 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 ) + { + integer constexpr NUM_DOF = 1; // pressure + integer constexpr NUM_EQN = 1; + + ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > flowDofNumberAccessor = + elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); + flowDofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + + using kernelType = ConnectorBasedAssemblyKernel< NUM_EQN, NUM_DOF >; + typename kernelType::SinglePhaseFlowAccessors flowAccessors( elemManager, solverName ); + typename kernelType::SinglePhaseFluidAccessors fluidAccessors( elemManager, solverName ); + typename kernelType::PermeabilityAccessors permAccessors( elemManager, solverName ); + typename kernelType::FracturePermeabilityAccessors fracPermAccessors( elemManager, solverName ); + + kernelType kernel( rankOffset, stencilWrapper, flowDofNumberAccessor, + flowAccessors, fluidAccessors, permAccessors, fracPermAccessors, + dt, localMatrix, localRhs, dR_dAper ); + + kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); + } +}; + +} // namespace SinglePhasePoromechanicsConformingFracturesALMKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp new file mode 100644 index 00000000000..66f90840ec1 --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsConformingFracturesALM.hpp @@ -0,0 +1,391 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SinglePhasePoromechanicsConformingFracturesALM.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP + +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFracturesALM.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernelBase.hpp" + +namespace geos +{ + +namespace thermalSinglePhasePoromechanicsConformingFracturesALMKernels +{ + +template< integer NUM_EQN, integer NUM_DOF > +class ConnectorBasedAssemblyKernel : public singlePhasePoromechanicsConformingFracturesALMKernels::ConnectorBasedAssemblyKernel< NUM_EQN, NUM_DOF > +{ +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 SinglePhaseFVMAbstractBase = singlePhaseFVMKernels::FluxComputeKernelBase; + using DofNumberAccessor = SinglePhaseFVMAbstractBase::DofNumberAccessor; + using SinglePhaseFlowAccessors = SinglePhaseFVMAbstractBase::SinglePhaseFlowAccessors; + using SinglePhaseFluidAccessors = SinglePhaseFVMAbstractBase::SinglePhaseFluidAccessors; + using PermeabilityAccessors = SinglePhaseFVMAbstractBase::PermeabilityAccessors; + using FracturePermeabilityAccessors = StencilMaterialAccessors< constitutive::PermeabilityBase, + fields::permeability::dPerm_dDispJump >; + using SinglePhaseFVMAbstractBase::m_dt; + using SinglePhaseFVMAbstractBase::m_rankOffset; + using SinglePhaseFVMAbstractBase::m_dofNumber; + using SinglePhaseFVMAbstractBase::m_gravCoef; + using SinglePhaseFVMAbstractBase::m_mob; + using SinglePhaseFVMAbstractBase::m_dens; + using SinglePhaseFVMAbstractBase::m_dDens; + using SinglePhaseFVMAbstractBase::m_dMob; + + using SinglePhaseFVMBase = singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper >; + using SinglePhaseFVMBase::numDof; + using SinglePhaseFVMBase::numEqn; + using SinglePhaseFVMBase::maxNumElems; + using SinglePhaseFVMBase::maxNumConns; + using SinglePhaseFVMBase::maxStencilSize; + using SinglePhaseFVMBase::m_stencilWrapper; + using SinglePhaseFVMBase::m_seri; + using SinglePhaseFVMBase::m_sesri; + using SinglePhaseFVMBase::m_sei; + using SinglePhaseFVMBase::m_ghostRank; + + using Base = singlePhasePoromechanicsConformingFracturesALMKernels::ConnectorBasedAssemblyKernel< NUM_EQN, NUM_DOF >; + + using ThermalSinglePhaseFlowAccessors = + StencilAccessors< fields::flow::temperature >; + + using ThermalSinglePhaseFluidAccessors = + StencilMaterialAccessors< constitutive::SingleFluidBase, + fields::singlefluid::enthalpy, + fields::singlefluid::dEnthalpy >; + + using ThermalConductivityAccessors = + StencilMaterialAccessors< constitutive::SinglePhaseThermalConductivityBase, + fields::thermalconductivity::effectiveConductivity >; + + + + ConnectorBasedAssemblyKernel( globalIndex const rankOffset, + SurfaceElementStencilWrapper const & stencilWrapper, + DofNumberAccessor const & flowDofNumberAccessor, + SinglePhaseFlowAccessors const & singlePhaseFlowAccessors, + ThermalSinglePhaseFlowAccessors const & thermalSinglePhaseFlowAccessors, + SinglePhaseFluidAccessors const & singlePhaseFluidAccessors, + ThermalSinglePhaseFluidAccessors const & thermalSinglePhaseFluidAccessors, + PermeabilityAccessors const & permeabilityAccessors, + FracturePermeabilityAccessors const & edfmPermeabilityAccessors, + ThermalConductivityAccessors const & thermalConductivityAccessors, + real64 const & dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs, + CRSMatrixView< real64, localIndex const > const & dR_dAper ) + : Base( rankOffset, + stencilWrapper, + flowDofNumberAccessor, + singlePhaseFlowAccessors, + singlePhaseFluidAccessors, + permeabilityAccessors, + edfmPermeabilityAccessors, + dt, + localMatrix, + localRhs, + dR_dAper ), + m_temp( thermalSinglePhaseFlowAccessors.get( fields::flow::temperature {} ) ), + m_enthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::enthalpy {} ) ), + m_dEnthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) ), + m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) ) + {} + + + /** + * @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 ), + energyFlux( 0.0 ), + dEnergyFlux_dTrans( 0.0 ), + dEnergyFlux_dP( size ), + dEnergyFlux_dT( size ), + dEnergyFlux_dDispJump( size, 3 ) + {} + using SinglePhaseFVMBase::StackVariables::stencilSize; + using SinglePhaseFVMBase::StackVariables::numFluxElems; + using SinglePhaseFVMBase::StackVariables::transmissibility; + using SinglePhaseFVMBase::StackVariables::dTrans_dPres; + using SinglePhaseFVMBase::StackVariables::dofColIndices; + using SinglePhaseFVMBase::StackVariables::localFlux; + using SinglePhaseFVMBase::StackVariables::localFluxJacobian; + + // Thermal transmissibility (for now, no derivatives) + + real64 thermalTransmissibility[maxNumConns][2]{}; + + // Energy fluxes and derivatives + + /// Energy fluxes + real64 energyFlux; + /// Derivative of the Energy fluxes wrt transmissibility + real64 dEnergyFlux_dTrans; + /// Derivatives of energy fluxes wrt pressure + stackArray1d< real64, maxStencilSize > dEnergyFlux_dP; + /// Derivatives of energy fluxes wrt temperature + stackArray1d< real64, maxStencilSize > dEnergyFlux_dT; + /// Derivatives of energy fluxes wrt dispJump + stackArray2d< real64, maxStencilSize *3 > dEnergyFlux_dDispJump{}; + + }; + + /** + * @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 + */ + GEOS_HOST_DEVICE + void computeFlux( localIndex const iconn, + StackVariables & stack ) const + + { + // *********************************************** + // First, we call the base computeFlux to compute: + // 1) compFlux and its derivatives (including derivatives wrt temperature), + // 2) enthalpy part of energyFlux and its derivatives (including derivatives wrt temperature) + // + // Computing dFlux_dT and the enthalpy flux requires quantities already computed in the base computeFlux, + // such as potGrad, fluxVal, and the indices of the upwind cell + // We use the lambda below (called **inside** the phase loop of the base computeFlux) to access these variables + Base::computeFlux( iconn, stack, [&] ( localIndex const (&k)[2], + localIndex const (&seri)[2], + localIndex const (&sesri)[2], + localIndex const (&sei)[2], + localIndex const, + real64 const & alpha, + real64 const & mobility, + real64 const & potGrad, + real64 const & massFlux, + real64 const & dMassFlux_dTrans, + real64 const (&dMassFlux_dP)[2] ) + { + real64 trans[2] = {stack.transmissibility[0][0], stack.transmissibility[0][1]}; + real64 dMassFlux_dT[2]{}; + + singlePhaseFluxKernelsHelper::computeEnthalpyFlux( seri, sesri, sei, + trans, + m_enthalpy, + m_dEnthalpy, + m_gravCoef, + m_dDens, + m_dMob, + alpha, + mobility, + potGrad, + massFlux, + dMassFlux_dTrans, + dMassFlux_dP, + dMassFlux_dT, + stack.energyFlux, + stack.dEnergyFlux_dTrans, + stack.dEnergyFlux_dP, + stack.dEnergyFlux_dT ); + + // add dMassFlux_dT to localFluxJacobian + for( integer ke = 0; ke < 2; ++ke ) + { + localIndex const localDofIndexTemp = k[ke] * numDof + numDof - 1; + stack.localFluxJacobian[k[0]*numEqn][localDofIndexTemp] += m_dt * dMassFlux_dT[ke]; + stack.localFluxJacobian[k[1]*numEqn][localDofIndexTemp] -= m_dt * dMassFlux_dT[ke]; + } + } ); + + // ***************************************************** + // Computation of the conduction term in the energy flux + // Note that the enthalpy term in the energy was computed above + // Note that this term is computed using an explicit treatment of conductivity for now + + // Step 1: compute the thermal transmissibilities at this face + m_stencilWrapper.computeWeights( iconn, + m_thermalConductivity, + m_thermalConductivity, // we have to pass something here, so we just use thermal conductivity + stack.thermalTransmissibility, + stack.dTrans_dPres ); // again, we have to pass something here, but this is unused for now + + localIndex k[2]; + localIndex connectionIndex = 0; + + for( k[0] = 0; k[0] < stack.numFluxElems; ++k[0] ) + { + for( k[1] = k[0] + 1; k[1] < stack.numFluxElems; ++k[1] ) + { + real64 const thermalTrans[2] = { stack.thermalTransmissibility[connectionIndex][0], stack.thermalTransmissibility[connectionIndex][1] }; + + localIndex const seri[2] = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )}; + localIndex const sesri[2] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )}; + localIndex const sei[2] = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )}; + + // Step 2: compute temperature difference at the interface + singlePhaseFluxKernelsHelper::computeConductiveFlux( seri, sesri, sei, m_temp, thermalTrans, stack.energyFlux, stack.dEnergyFlux_dT ); + + // add energyFlux and its derivatives to localFlux and localFluxJacobian + stack.localFlux[k[0]*numEqn + numEqn - 1] += m_dt * stack.energyFlux; + stack.localFlux[k[1]*numEqn + numEqn - 1] -= m_dt * stack.energyFlux; + + for( integer ke = 0; ke < 2; ++ke ) + { + integer const localDofIndexPres = k[ke] * numDof; + stack.localFluxJacobian[k[0]*numEqn + numEqn - 1][localDofIndexPres] = m_dt * stack.dEnergyFlux_dP[ke]; + stack.localFluxJacobian[k[1]*numEqn + numEqn - 1][localDofIndexPres] = -m_dt * stack.dEnergyFlux_dP[ke]; + integer const localDofIndexTemp = localDofIndexPres + 1; + stack.localFluxJacobian[k[0]*numEqn + numEqn - 1][localDofIndexTemp] = m_dt * stack.dEnergyFlux_dT[ke]; + stack.localFluxJacobian[k[1]*numEqn + numEqn - 1][localDofIndexTemp] = -m_dt * stack.dEnergyFlux_dT[ke]; + } + + connectionIndex++; + } + } + } + + /** + * @brief Performs the complete phase for the kernel. + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + */ + GEOS_HOST_DEVICE + void complete( localIndex const iconn, + StackVariables & stack ) const + { + // Call Base::complete to assemble the mass balance equations + // In the lambda, add contribution to residual and jacobian into the energy balance equation + Base::complete( iconn, stack, [&] ( integer const i, + localIndex const localRow ) + { + // The no. of fluxes is equal to the no. of equations in m_localRhs and m_localMatrix + RAJA::atomicAdd( parallelDeviceAtomic{}, &SinglePhaseFVMAbstractBase::m_localRhs[localRow + numEqn-1], stack.localFlux[i * numEqn + numEqn-1] ); + + SinglePhaseFVMAbstractBase::m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow + numEqn-1, + stack.dofColIndices.data(), + stack.localFluxJacobian[i * numEqn + numEqn-1].dataIfContiguous(), + stack.stencilSize * numDof ); + + } ); + } + + +private: + + /// Views on temperature + ElementViewConst< arrayView1d< real64 const > > const m_temp; + + /// Views on enthalpies + ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_enthalpy; + + /// Views on enthalpies + ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const m_dEnthalpy; + + /// View on thermal conductivity + ElementViewConst< arrayView3d< real64 const > > m_thermalConductivity; + +}; + + + +/** + * @class ConnectorBasedAssemblyKernelFactory + */ +class ConnectorBasedAssemblyKernelFactory +{ +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( globalIndex const rankOffset, + string const & flowDofKey, + 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 ) + { + integer constexpr NUM_DOF = 2; // pressure + temperature + integer constexpr NUM_EQN = 2; // mass balance + energy balance + + + ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > flowDofNumberAccessor = + elemManager.constructArrayViewAccessor< globalIndex, 1 >( flowDofKey ); + flowDofNumberAccessor.setName( solverName + "/accessors/" + flowDofKey ); + + using kernelType = ConnectorBasedAssemblyKernel< NUM_EQN, NUM_DOF >; + typename kernelType::SinglePhaseFlowAccessors flowAccessors( elemManager, solverName ); + typename kernelType::ThermalSinglePhaseFlowAccessors thermalFlowAccessors( elemManager, solverName ); + + typename kernelType::SinglePhaseFluidAccessors fluidAccessors( elemManager, solverName ); + typename kernelType::ThermalSinglePhaseFluidAccessors thermalFluidAccessors( elemManager, solverName ); + + typename kernelType::PermeabilityAccessors permAccessors( elemManager, solverName ); + typename kernelType::FracturePermeabilityAccessors edfmPermAccessors( elemManager, solverName ); + typename kernelType::ThermalConductivityAccessors thermalConductivityAccessors( elemManager, solverName ); + + kernelType kernel( rankOffset, stencilWrapper, + flowDofNumberAccessor, + flowAccessors, thermalFlowAccessors, fluidAccessors, thermalFluidAccessors, + permAccessors, edfmPermAccessors, thermalConductivityAccessors, + dt, localMatrix, localRhs, dR_dAper ); + + kernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); + } +}; + + + +} // namespace SinglePhaseProppantFluxKernels + +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_THERMALSINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP