diff --git a/inputFiles/thermalMultiphaseFlow/co2_thermal_2d.xml b/inputFiles/thermalMultiphaseFlow/co2_thermal_2d.xml new file mode 100644 index 00000000000..d23b8900942 --- /dev/null +++ b/inputFiles/thermalMultiphaseFlow/co2_thermal_2d.xml @@ -0,0 +1,223 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/thermalMultiphaseFlow/co2flash.txt b/inputFiles/thermalMultiphaseFlow/co2flash.txt new file mode 100644 index 00000000000..4c5a766d852 --- /dev/null +++ b/inputFiles/thermalMultiphaseFlow/co2flash.txt @@ -0,0 +1 @@ +FlashModel CO2Solubility 1e6 7.5e7 5e5 299.15 369.15 10 0 diff --git a/inputFiles/thermalMultiphaseFlow/pvtgas.txt b/inputFiles/thermalMultiphaseFlow/pvtgas.txt new file mode 100644 index 00000000000..ccd4e817f67 --- /dev/null +++ b/inputFiles/thermalMultiphaseFlow/pvtgas.txt @@ -0,0 +1,3 @@ +DensityFun SpanWagnerCO2Density 1e6 7.5e7 5e5 299.15 369.15 10 +ViscosityFun FenghourCO2Viscosity 1e6 7.5e7 5e5 299.15 369.15 10 +EnthalpyFun CO2Enthalpy 1e6 7.5e7 5e5 299.15 369.15 10 diff --git a/inputFiles/thermalMultiphaseFlow/pvtliquid.txt b/inputFiles/thermalMultiphaseFlow/pvtliquid.txt new file mode 100644 index 00000000000..4df0c25f042 --- /dev/null +++ b/inputFiles/thermalMultiphaseFlow/pvtliquid.txt @@ -0,0 +1,3 @@ +DensityFun PhillipsBrineDensity 1e6 7.5e7 5e5 299.15 369.15 10 0 +ViscosityFun PhillipsBrineViscosity 0 +EnthalpyFun BrineEnthalpy 1e6 7.5e7 5e5 299.15 369.15 10 0 diff --git a/integratedTests b/integratedTests index a1d96eb348f..16dd843b5b3 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit a1d96eb348f1169ec06481a0891c06b04ceabd26 +Subproject commit 16dd843b5b3d53a437b02b492b4a032b6dbdbb06 diff --git a/src/coreComponents/constitutive/fluid/CO2BrineFluid.cpp b/src/coreComponents/constitutive/fluid/CO2BrineFluid.cpp index 1567cf15be2..48b20435d01 100644 --- a/src/coreComponents/constitutive/fluid/CO2BrineFluid.cpp +++ b/src/coreComponents/constitutive/fluid/CO2BrineFluid.cpp @@ -185,7 +185,7 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::createPVTModels() { string_array const strs = stringutilities::tokenize( str, " " ); - if( strs[0] == toString( SubModelInputNames::DENSITY ) ) + if( strs[0] == "DensityFun" ) { if( strs[1] == PHASE1::Density::catalogName() ) { @@ -196,7 +196,7 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::createPVTModels() phase2InputParams[PHASE2::InputParamOrder::DENSITY] = strs; } } - else if( strs[0] == toString( SubModelInputNames::VISCOSITY ) ) + else if( strs[0] == "ViscosityFun" ) { if( strs[1] == PHASE1::Viscosity::catalogName() ) { @@ -207,7 +207,7 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::createPVTModels() phase2InputParams[PHASE2::InputParamOrder::VISCOSITY] = strs; } } - else if( strs[0] == toString( SubModelInputNames::ENTHALPY ) ) + else if( strs[0] == "EnthalpyFun" ) { if( strs[1] == PHASE1::Enthalpy::catalogName() ) { diff --git a/src/coreComponents/constitutive/fluid/CO2BrineFluid.hpp b/src/coreComponents/constitutive/fluid/CO2BrineFluid.hpp index 875b3f2085b..c2752fdd8e6 100644 --- a/src/coreComponents/constitutive/fluid/CO2BrineFluid.hpp +++ b/src/coreComponents/constitutive/fluid/CO2BrineFluid.hpp @@ -569,25 +569,6 @@ CO2BrineFluid< PHASE1, PHASE2, FLASH >::KernelWrapper:: m_totalDensity( k, q ) ); } -/// Declare strings associated with enumeration values -/// Needed for now, because we don't use the catalogNames for input (yet) -ENUM_STRINGS( CO2BrinePhillipsFluid::SubModelInputNames, - "DensityFun", - "ViscosityFun", - "EnthalpyFun" ); -ENUM_STRINGS( CO2BrinePhillipsThermalFluid::SubModelInputNames, - "DensityFun", - "ViscosityFun", - "EnthalpyFun" ); -ENUM_STRINGS( CO2BrineEzrokhiFluid::SubModelInputNames, - "DensityFun", - "ViscosityFun", - "EnthalpyFun" ); -ENUM_STRINGS( CO2BrineEzrokhiThermalFluid::SubModelInputNames, - "DensityFun", - "ViscosityFun", - "EnthalpyFun" ); - } // namespace constitutive diff --git a/src/coreComponents/constitutive/fluid/MultiFluidExtrinsicData.hpp b/src/coreComponents/constitutive/fluid/MultiFluidExtrinsicData.hpp index 3ab287e1ec1..a589e4c2661 100644 --- a/src/coreComponents/constitutive/fluid/MultiFluidExtrinsicData.hpp +++ b/src/coreComponents/constitutive/fluid/MultiFluidExtrinsicData.hpp @@ -156,7 +156,7 @@ EXTRINSIC_MESH_DATA_TRAIT( dPhaseInternalEnergy, 0, NOPLOT, NO_WRITE, - "Derivative of phase internal energy with respect to pressure, temperature, and global component fraction" ); + "Derivative of phase internal energy with respect to pressure, temperature, and global component fractions" ); EXTRINSIC_MESH_DATA_TRAIT( phaseCompFraction, "phaseCompFraction", diff --git a/src/coreComponents/constitutive/fluid/multiFluidSelector.hpp b/src/coreComponents/constitutive/fluid/multiFluidSelector.hpp index 96321059c35..6977c38f81f 100644 --- a/src/coreComponents/constitutive/fluid/multiFluidSelector.hpp +++ b/src/coreComponents/constitutive/fluid/multiFluidSelector.hpp @@ -44,7 +44,10 @@ void constitutiveUpdatePassThru( MultiFluidBase const & fluid, CompositionalMultiphaseFluid, #endif CO2BrinePhillipsFluid, - CO2BrineEzrokhiFluid >::execute( fluid, std::forward< LAMBDA >( lambda ) ); + CO2BrineEzrokhiFluid, + CO2BrinePhillipsThermalFluid /*, // if I uncomment the two models at the same time, the compiler segfaults on + Lassen! + CO2BrineEzrokhiThermalFluid*/>::execute( fluid, std::forward< LAMBDA >( lambda ) ); } template< typename LAMBDA > @@ -57,7 +60,10 @@ void constitutiveUpdatePassThru( MultiFluidBase & fluid, CompositionalMultiphaseFluid, #endif CO2BrinePhillipsFluid, - CO2BrineEzrokhiFluid >::execute( fluid, std::forward< LAMBDA >( lambda ) ); + CO2BrineEzrokhiFluid, + CO2BrinePhillipsThermalFluid /*, // if I uncomment the two models at the same time, the compiler segfaults on + Lassen! + CO2BrineEzrokhiThermalFluid*/>::execute( fluid, std::forward< LAMBDA >( lambda ) ); } } // namespace constitutive diff --git a/src/coreComponents/constitutive/solid/SolidInternalEnergy.cpp b/src/coreComponents/constitutive/solid/SolidInternalEnergy.cpp index 618574f28d3..add71a381c9 100644 --- a/src/coreComponents/constitutive/solid/SolidInternalEnergy.cpp +++ b/src/coreComponents/constitutive/solid/SolidInternalEnergy.cpp @@ -30,34 +30,34 @@ SolidInternalEnergy::SolidInternalEnergy( string const & name, Group * const par ConstitutiveBase( name, parent ), m_internalEnergy(), m_dInternalEnergy_dTemperature(), - m_specificHeatCapacity(), + m_volumetricHeatCapacity(), m_referenceTemperature(), m_referenceInternalEnergy() { registerWrapper( viewKeyStruct::internalEnergyString(), &m_internalEnergy ). setPlotLevel( PlotLevel::LEVEL_0 ). setApplyDefaultValue( 0.0 ). - setDescription( "Internal energy of the solid" ); + setDescription( "Internal energy of the solid per unit volume [J/m^3]" ); registerWrapper( viewKeyStruct::oldInternalEnergyString(), &m_oldInternalEnergy ). setApplyDefaultValue( 0.0 ). - setDescription( "Internal energy of the solid at the previous time-step" ); + setDescription( "Internal energy of the solid per unit volume at the previous time-step [J/m^3]" ); registerWrapper( viewKeyStruct::dInternalEnergy_dTemperatureString(), &m_dInternalEnergy_dTemperature ). setApplyDefaultValue( 0.0 ). - setDescription( "Derivative of the solid internal energy w.r.t. temperature" ); + setDescription( "Derivative of the solid internal energy w.r.t. temperature [J/(m^3.K)]" ); - registerWrapper( viewKeyStruct::specificHeatCapacityString(), &m_specificHeatCapacity ). - setApplyDefaultValue( 0.0 ). - setDescription( "Solid specific heat capacity" ); + registerWrapper( viewKeyStruct::volumetricHeatCapacityString(), &m_volumetricHeatCapacity ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Solid volumetric heat capacity [J/(kg.K)]" ); registerWrapper( viewKeyStruct::referenceTemperatureString(), &m_referenceTemperature ). - setApplyDefaultValue( 0.0 ). - setDescription( "Reference temperature" ); + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Reference temperature [K]" ); registerWrapper( viewKeyStruct::referenceInternalEnergyString(), &m_referenceInternalEnergy ). - setApplyDefaultValue( 0.0 ). - setDescription( "Internal energy at the reference temperature" ); + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Internal energy at the reference temperature [J/kg]" ); } void SolidInternalEnergy::allocateConstitutiveData( Group & parent, @@ -81,6 +81,8 @@ void SolidInternalEnergy::saveConvergedState() const } ); } -} +REGISTER_CATALOG_ENTRY( ConstitutiveBase, SolidInternalEnergy, string const &, Group * const ) -} +} /* namespace constitutive */ + +} /* namespace geosx */ diff --git a/src/coreComponents/constitutive/solid/SolidInternalEnergy.hpp b/src/coreComponents/constitutive/solid/SolidInternalEnergy.hpp index aff61795562..ada374344bd 100644 --- a/src/coreComponents/constitutive/solid/SolidInternalEnergy.hpp +++ b/src/coreComponents/constitutive/solid/SolidInternalEnergy.hpp @@ -35,12 +35,12 @@ class SolidInternalEnergyUpdates SolidInternalEnergyUpdates( arrayView2d< real64 > const & internalEnergy, arrayView2d< real64 > const & dInternalEnergy_dTemperature, - real64 const & specificHeatCapacity, + real64 const & volumetricHeatCapacity, real64 const & referenceTemperature, real64 const & referenceInternalEnergy ): m_internalEnergy( internalEnergy ), m_dInternalEnergy_dTemperature( dInternalEnergy_dTemperature ), - m_specificHeatCapacity( specificHeatCapacity ), + m_volumetricHeatCapacity( volumetricHeatCapacity ), m_referenceTemperature( referenceTemperature ), m_referenceInternalEnergy( referenceInternalEnergy ) {} @@ -59,8 +59,8 @@ class SolidInternalEnergyUpdates real64 & internalEnergy, real64 & dInternalEnergy_dTemperature ) const { - internalEnergy = m_referenceInternalEnergy + m_specificHeatCapacity * ( temperature - m_referenceTemperature ); - dInternalEnergy_dTemperature = m_specificHeatCapacity; + internalEnergy = m_referenceInternalEnergy + m_volumetricHeatCapacity * ( temperature - m_referenceTemperature ); + dInternalEnergy_dTemperature = m_volumetricHeatCapacity; } private: @@ -71,8 +71,8 @@ class SolidInternalEnergyUpdates /// Derivative of the solid internal energy w.r.t. the temperature arrayView2d< real64 > m_dInternalEnergy_dTemperature; - /// Solid specific heat capacity - real64 m_specificHeatCapacity; + /// Solid volumetric heat capacity + real64 m_volumetricHeatCapacity; /// Reference temperature real64 m_referenceTemperature; @@ -99,7 +99,7 @@ class SolidInternalEnergy : public ConstitutiveBase static constexpr char const * internalEnergyString() { return "internalEnergy"; } static constexpr char const * oldInternalEnergyString() { return "oldInternalEnergy"; } static constexpr char const * dInternalEnergy_dTemperatureString() { return "dInternalEnergy_dTemperature"; } - static constexpr char const * specificHeatCapacityString() { return "specificHeatCapacity"; } + static constexpr char const * volumetricHeatCapacityString() { return "volumetricHeatCapacity"; } static constexpr char const * referenceTemperatureString() { return "referenceTemperature"; } static constexpr char const * referenceInternalEnergyString() { return "referenceInternalEnergy"; } } viewKeys; @@ -114,7 +114,7 @@ class SolidInternalEnergy : public ConstitutiveBase { return KernelWrapper( m_internalEnergy, m_dInternalEnergy_dTemperature, - m_specificHeatCapacity, + m_volumetricHeatCapacity, m_referenceTemperature, m_referenceInternalEnergy ); } @@ -152,8 +152,8 @@ class SolidInternalEnergy : public ConstitutiveBase /// Derivative of the solid internal energy w.r.t. the temperature array2d< real64 > m_dInternalEnergy_dTemperature; - /// Solid specific heat capacity - real64 m_specificHeatCapacity; + /// Solid volumetric heat capacity + real64 m_volumetricHeatCapacity; /// Reference temperature real64 m_referenceTemperature; diff --git a/src/coreComponents/constitutive/thermalConductivity/ConstantThermalConductivity.cpp b/src/coreComponents/constitutive/thermalConductivity/ConstantThermalConductivity.cpp index 2170d1dff33..11fd59865a0 100644 --- a/src/coreComponents/constitutive/thermalConductivity/ConstantThermalConductivity.cpp +++ b/src/coreComponents/constitutive/thermalConductivity/ConstantThermalConductivity.cpp @@ -32,7 +32,7 @@ ConstantThermalConductivity::ConstantThermalConductivity( string const & name, G registerWrapper( viewKeyStruct::thermalConductivityComponentsString(), &m_thermalConductivityComponents ). setInputFlag( InputFlags::REQUIRED ). setRestartFlags( RestartFlags::NO_WRITE ). - setDescription( "xx, yy, and zz components of a diagonal thermal conductivity tensor [W/(m.K)]" ); + setDescription( "xx, yy, and zz components of a diagonal thermal conductivity tensor [J/(s.m.K)]" ); } std::unique_ptr< ConstitutiveBase > diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt index de4e7e4bf1b..89712b95be1 100644 --- a/src/coreComponents/physicsSolvers/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/CMakeLists.txt @@ -16,9 +16,11 @@ set( physicsSolvers_headers contact/SolidMechanicsEmbeddedFractures.hpp fluidFlow/CompositionalMultiphaseBase.hpp fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp - fluidFlow/CompositionalMultiphaseBaseKernels.hpp + fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp + fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp fluidFlow/CompositionalMultiphaseFVM.hpp - fluidFlow/CompositionalMultiphaseFVMKernels.hpp + fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp + fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp fluidFlow/CompositionalMultiphaseHybridFVM.hpp fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp fluidFlow/CompositionalMultiphaseUtilities.hpp @@ -95,8 +97,8 @@ set( physicsSolvers_sources contact/LagrangianContactSolver.cpp contact/SolidMechanicsEmbeddedFractures.cpp fluidFlow/CompositionalMultiphaseBase.cpp - fluidFlow/CompositionalMultiphaseFVM.cpp - fluidFlow/CompositionalMultiphaseFVMKernels.cpp + fluidFlow/CompositionalMultiphaseFVM.cpp + fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp fluidFlow/CompositionalMultiphaseHybridFVM.cpp fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp fluidFlow/FlowSolverBase.cpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 6291b214332..f3ceca70708 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -28,6 +28,8 @@ #include "constitutive/fluid/multiFluidSelector.hpp" #include "constitutive/relativePermeability/RelativePermeabilityExtrinsicData.hpp" #include "constitutive/relativePermeability/relativePermeabilitySelector.hpp" +#include "constitutive/solid/SolidBase.hpp" +#include "constitutive/solid/SolidInternalEnergy.hpp" #include "constitutive/thermalConductivity/thermalConductivitySelector.hpp" #include "constitutive/permeability/PermeabilityExtrinsicData.hpp" #include "fieldSpecification/AquiferBoundaryCondition.hpp" @@ -36,9 +38,10 @@ #include "fieldSpecification/SourceFluxBoundaryCondition.hpp" #include "mesh/DomainPartition.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" -#include "physicsSolvers/fluidFlow/FlowSolverBaseExtrinsicData.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseExtrinsicData.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" +#include "physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp" #if defined( __INTEL_COMPILER ) #pragma GCC optimize "O0" @@ -49,7 +52,6 @@ namespace geosx using namespace dataRepository; using namespace constitutive; -using namespace compositionalMultiphaseBaseKernels; CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, Group * const parent ) @@ -74,6 +76,11 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Use mass formulation instead of molar" ); + this->registerWrapper( viewKeyStruct::isThermalString(), &m_isThermal ). + setApplyDefaultValue( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Flag indicating whether the problem is thermal or not." ); + this->registerWrapper( viewKeyStruct::computeCFLNumbersString(), &m_computeCFLNumbers ). setApplyDefaultValue( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -143,7 +150,8 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) m_numPhases = referenceFluid.numFluidPhases(); m_numComponents = referenceFluid.numFluidComponents(); } - m_numDofPerCell = m_numComponents + 1; + // n_c components + one pressure ( + one temperature if needed ) + m_numDofPerCell = m_isThermal ? m_numComponents + 2 : m_numComponents + 1; // 2. Register and resize all fields as necessary forMeshTargets( meshBodies, [&]( string const &, @@ -179,9 +187,13 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerExtrinsicData< initialPressure >( getName() ); subRegion.registerExtrinsicData< deltaPressure >( getName() ); - subRegion.registerExtrinsicData< bcPressure >( getName() ); + subRegion.registerExtrinsicData< bcPressure >( getName() ); // needed for the application of boundary conditions + // these fields are always registered for the evaluation of the fluid properties subRegion.registerExtrinsicData< temperature >( getName() ); + subRegion.registerExtrinsicData< deltaTemperature >( getName() ); + + subRegion.registerExtrinsicData< bcTemperature >( getName() ); // needed for the application of boundary conditions // The resizing of the arrays needs to happen here, before the call to initializePreSubGroups, // to make sure that the dimensions are properly set before the timeHistoryOutput starts its initialization. @@ -214,6 +226,15 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerExtrinsicData< dPhaseMobility_dGlobalCompDensity >( getName() ). reference().resizeDimension< 1, 2 >( m_numPhases, m_numComponents ); + if( m_isThermal ) + { + subRegion.registerExtrinsicData< dPhaseVolumeFraction_dTemperature >( getName() ). + reference().resizeDimension< 1 >( m_numPhases ); + + subRegion.registerExtrinsicData< dPhaseMobility_dTemperature >( getName() ). + reference().resizeDimension< 1 >( m_numPhases ); + } + if( m_computeCFLNumbers ) { subRegion.registerExtrinsicData< phaseOutflux >( getName() ). @@ -277,6 +298,18 @@ void CompositionalMultiphaseBase::setConstitutiveNames( ElementSubRegionBase & s if( m_isThermal ) { + string & solidInternalEnergyName = subRegion.registerWrapper< string >( viewKeyStruct::solidInternalEnergyNamesString() ). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setSizedFromParent( 0 ). + setDescription( "Name of the solid internal energy constitutive model to use" ). + reference(); + + solidInternalEnergyName = getConstitutiveName< SolidInternalEnergy >( subRegion ); + GEOSX_THROW_IF( solidInternalEnergyName.empty(), + GEOSX_FMT( "Solid internal energy model not found on subregion {}", subRegion.getName() ), + InputError ); + string & thermalConductivityName = subRegion.registerWrapper< string >( viewKeyStruct::thermalConductivityNamesString() ). setPlotLevel( PlotLevel::NOPLOT ). setRestartFlags( RestartFlags::NO_WRITE ). @@ -285,12 +318,10 @@ void CompositionalMultiphaseBase::setConstitutiveNames( ElementSubRegionBase & s reference(); thermalConductivityName = getConstitutiveName< ThermalConductivityBase >( subRegion ); - GEOSX_THROW_IF( relPermName.empty(), - GEOSX_FMT( "Thermal Conductivity model not found on subregion {}", subRegion.getName() ), + GEOSX_THROW_IF( thermalConductivityName.empty(), + GEOSX_FMT( "Thermal conductivity model not found on subregion {}", subRegion.getName() ), InputError ); } - - } @@ -345,8 +376,6 @@ void CompositionalMultiphaseBase::initializeAquiferBC( ConstitutiveManager const integer const waterPhaseIndex = fluid0.getWaterPhaseIndex(); bc.setWaterPhaseIndex( waterPhaseIndex ); - - arrayView1d< real64 const > const & aquiferWaterPhaseCompFrac = bc.getWaterPhaseComponentFraction(); arrayView1d< string const > const & aquiferWaterPhaseCompNames = bc.getWaterPhaseComponentNames(); @@ -454,7 +483,8 @@ void CompositionalMultiphaseBase::updateComponentFraction( ObjectManagerBase & d { GEOSX_MARK_FUNCTION; - ComponentFractionKernelFactory:: + isothermalCompositionalMultiphaseBaseKernels:: + ComponentFractionKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, dataGroup ); @@ -467,12 +497,24 @@ void CompositionalMultiphaseBase::updatePhaseVolumeFraction( ObjectManagerBase & string const & fluidName = dataGroup.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( dataGroup, fluidName ); - PhaseVolumeFractionKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dataGroup, - fluid ); - + if( m_isThermal ) + { + thermalCompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid ); + } + else + { + isothermalCompositionalMultiphaseBaseKernels:: + PhaseVolumeFractionKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid ); + } } void CompositionalMultiphaseBase::updateFluidModel( ObjectManagerBase & dataGroup ) const @@ -482,6 +524,7 @@ void CompositionalMultiphaseBase::updateFluidModel( ObjectManagerBase & dataGrou arrayView1d< real64 const > const pres = dataGroup.getExtrinsicData< extrinsicMeshData::flow::pressure >(); arrayView1d< real64 const > const dPres = dataGroup.getExtrinsicData< extrinsicMeshData::flow::deltaPressure >(); arrayView1d< real64 const > const temp = dataGroup.getExtrinsicData< extrinsicMeshData::flow::temperature >(); + arrayView1d< real64 const > const dTemp = dataGroup.getExtrinsicData< extrinsicMeshData::flow::deltaTemperature >(); arrayView2d< real64 const, compflow::USD_COMP > const compFrac = dataGroup.getExtrinsicData< extrinsicMeshData::flow::globalCompFraction >(); @@ -494,12 +537,15 @@ void CompositionalMultiphaseBase::updateFluidModel( ObjectManagerBase & dataGrou using ExecPolicy = typename FluidType::exec_policy; typename FluidType::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); - FluidUpdateKernel::launch< ExecPolicy >( dataGroup.size(), - fluidWrapper, - pres, - dPres, - temp, - compFrac ); + thermalCompositionalMultiphaseBaseKernels:: + FluidUpdateKernel:: + launch< ExecPolicy >( dataGroup.size(), + fluidWrapper, + pres, + dPres, + temp, + dTemp, + compFrac ); } ); } @@ -517,9 +563,11 @@ void CompositionalMultiphaseBase::updateRelPermModel( ObjectManagerBase & dataGr { typename TYPEOFREF( castedRelPerm ) ::KernelWrapper relPermWrapper = castedRelPerm.createKernelWrapper(); - RelativePermeabilityUpdateKernel::launch< parallelDevicePolicy<> >( dataGroup.size(), - relPermWrapper, - phaseVolFrac ); + isothermalCompositionalMultiphaseBaseKernels:: + RelativePermeabilityUpdateKernel:: + launch< parallelDevicePolicy<> >( dataGroup.size(), + relPermWrapper, + phaseVolFrac ); } ); } @@ -539,13 +587,35 @@ void CompositionalMultiphaseBase::updateCapPressureModel( ObjectManagerBase & da { typename TYPEOFREF( castedCapPres ) ::KernelWrapper capPresWrapper = castedCapPres.createKernelWrapper(); - CapillaryPressureUpdateKernel::launch< parallelDevicePolicy<> >( dataGroup.size(), - capPresWrapper, - phaseVolFrac ); + isothermalCompositionalMultiphaseBaseKernels:: + CapillaryPressureUpdateKernel:: + launch< parallelDevicePolicy<> >( dataGroup.size(), + capPresWrapper, + phaseVolFrac ); } ); } } +void CompositionalMultiphaseBase::updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const +{ + arrayView1d< real64 const > const temp = dataGroup.getExtrinsicData< extrinsicMeshData::flow::temperature >(); + arrayView1d< real64 const > const dTemp = dataGroup.getExtrinsicData< extrinsicMeshData::flow::deltaTemperature >(); + + string const & solidInternalEnergyName = dataGroup.getReference< string >( viewKeyStruct::solidInternalEnergyNamesString() ); + SolidInternalEnergy & solidInternalEnergy = getConstitutiveModel< SolidInternalEnergy >( dataGroup, solidInternalEnergyName ); + + SolidInternalEnergy::KernelWrapper solidInternalEnergyWrapper = solidInternalEnergy.createKernelUpdates(); + + // TODO: this should go somewhere, handle the case of flow in fracture, etc + + thermalCompositionalMultiphaseBaseKernels:: + SolidInternalEnergyUpdateKernel:: + launch< parallelDevicePolicy<> >( dataGroup.size(), + solidInternalEnergyWrapper, + temp, + dTemp ); +} + void CompositionalMultiphaseBase::updateFluidState( ObjectManagerBase & subRegion ) const { GEOSX_MARK_FUNCTION; @@ -677,12 +747,11 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, // - This step depends phaseRelPerm updatePhaseMobility( subRegion ); - // 4.6 Finally, we initialize the thermal conductivity + // 4.6 Finally, we initialize the rock thermal quantities: conductivity and solid internal energy // // Note: // - This must be called after updatePorosityAndPermeability and updatePhaseVolumeFraction // - This step depends on porosity and phaseVolFraction - // - Energy balance is not supported yet, so the following flag is always false for now if( m_isThermal ) { // initialized porosity @@ -693,6 +762,13 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, getConstitutiveModel< ThermalConductivityBase >( subRegion, thermalConductivityName ); conductivityMaterial.initializeRockFluidState( porosity, phaseVolFrac ); // note that there is nothing to update here because thermal conductivity is explicit for now + + updateSolidInternalEnergyModel( subRegion ); + string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString() ); + SolidInternalEnergy const & solidInternalEnergyMaterial = + getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); + solidInternalEnergyMaterial.saveConvergedState(); + } } ); @@ -871,32 +947,34 @@ void CompositionalMultiphaseBase::computeHydrostaticEquilibrium() typename FluidType::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); // note: inside this kernel, serialPolicy is used, and elevation/pressure values don't go to the GPU - HydrostaticPressureKernel::ReturnType const returnValue = - HydrostaticPressureKernel::launch( numPointsInTable, - numComps, - numPhases, - ipInit, - maxNumEquilIterations, - equilTolerance, - gravVector, - minElevation, - elevationIncrement, - datumElevation, - datumPressure, - fluidWrapper, - compFracTableWrappers.toViewConst(), - tempTableWrapper, - elevationValues.toNestedView(), - pressureValues.toView() ); - - GEOSX_THROW_IF( returnValue == HydrostaticPressureKernel::ReturnType::FAILED_TO_CONVERGE, + isothermalCompositionalMultiphaseBaseKernels:: + HydrostaticPressureKernel::ReturnType const returnValue = + isothermalCompositionalMultiphaseBaseKernels:: + HydrostaticPressureKernel::launch( numPointsInTable, + numComps, + numPhases, + ipInit, + maxNumEquilIterations, + equilTolerance, + gravVector, + minElevation, + elevationIncrement, + datumElevation, + datumPressure, + fluidWrapper, + compFracTableWrappers.toViewConst(), + tempTableWrapper, + elevationValues.toNestedView(), + pressureValues.toView() ); + + GEOSX_THROW_IF( returnValue == isothermalCompositionalMultiphaseBaseKernels::HydrostaticPressureKernel::ReturnType::FAILED_TO_CONVERGE, CompositionalMultiphaseBase::catalogName() << " " << getName() << ": hydrostatic pressure initialization failed to converge in region " << region.getName() << "! \n" << "Try to loosen the equilibration tolerance, or increase the number of equilibration iterations. \n" << "If nothing works, something may be wrong in the fluid model, see ", std::runtime_error ); - GEOSX_LOG_RANK_0_IF( returnValue == HydrostaticPressureKernel::ReturnType::DETECTED_MULTIPHASE_FLOW, + GEOSX_LOG_RANK_0_IF( returnValue == isothermalCompositionalMultiphaseBaseKernels::HydrostaticPressureKernel::ReturnType::DETECTED_MULTIPHASE_FLOW, CompositionalMultiphaseBase::catalogName() << " " << getName() << ": currently, GEOSX assumes that there is only one mobile phase when computing the hydrostatic pressure. \n" << "We detected multiple phases using the provided datum pressure, temperature, and component fractions. \n" @@ -1037,7 +1115,9 @@ void CompositionalMultiphaseBase::backupFields( MeshLevel & mesh, forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { if( elemGhostRank[ei] >= 0 ) + { return; + } for( integer ip = 0; ip < numPhase; ++ip ) { @@ -1045,6 +1125,7 @@ void CompositionalMultiphaseBase::backupFields( MeshLevel & mesh, phaseMobOld[ei][ip] = phaseMob[ei][ip]; } } ); + } ); } @@ -1085,8 +1166,6 @@ void CompositionalMultiphaseBase::assembleSystem( real64 const GEOSX_UNUSED_PARA dofManager, localMatrix, localRhs ); - - } void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( DomainPartition & domain, @@ -1111,16 +1190,34 @@ void CompositionalMultiphaseBase::assembleAccumulationAndVolumeBalanceTerms( Dom MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); CoupledSolidBase const & solid = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - ElementBasedAssemblyKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - dofKey, - subRegion, - fluid, - solid, - localMatrix, - localRhs ); + if( m_isThermal ) + { + thermalCompositionalMultiphaseBaseKernels:: + ElementBasedAssemblyKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); + } + else + { + isothermalCompositionalMultiphaseBaseKernels:: + ElementBasedAssemblyKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + dofKey, + subRegion, + fluid, + solid, + localMatrix, + localRhs ); + } } ); } ); } @@ -1252,13 +1349,16 @@ namespace { bool validateDirichletBC( DomainPartition & domain, + integer const isThermal, integer const numComp, real64 const time ) { constexpr integer MAX_NC = MultiFluidBase::MAX_NUM_COMPONENTS; FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); - map< string, map< string, map< string, ComponentMask< MAX_NC > > > > bcStatusMap; // map to check consistent application of BC + // maps to check consistent application of BC + map< string, map< string, map< string, ComponentMask< MAX_NC > > > > bcStatusMap; + map< string, map< string, set< string > > > bcTempStatusMap; bool bcConsistent = true; // 1. Check pressure Dirichlet BCs @@ -1272,7 +1372,7 @@ bool validateDirichletBC( DomainPartition & domain, Group & subRegion, string const & ) { - // 1.0. Check whether pressure has already been applied to this set + // Check whether pressure has already been applied to this set string const & subRegionName = subRegion.getName(); string const & regionName = subRegion.getParent().getParent().getName(); @@ -1285,7 +1385,44 @@ bool validateDirichletBC( DomainPartition & domain, subRegionSetMap[setName].setNumComp( numComp ); } ); - // 2. Check composition BC (global component fraction) + // 2. Check temperature Dirichlet BCs + if( isThermal ) + { + fsManager.apply( time, + domain.getMeshBody( 0 ).getMeshLevel( 0 ), + "ElementRegions", + extrinsicMeshData::flow::temperature::key(), + [&]( FieldSpecificationBase const &, + string const & setName, + SortedArrayView< localIndex const > const &, + Group & subRegion, + string const & ) + { + string const & subRegionName = subRegion.getName(); + string const & regionName = subRegion.getParent().getParent().getName(); + + // 2.1 Check whether temperature has already been applied to this set + auto & tempSubRegionSetMap = bcTempStatusMap[regionName][subRegionName]; + if( tempSubRegionSetMap.count( setName ) > 0 ) + { + bcConsistent = false; + GEOSX_WARNING( GEOSX_FMT( "Conflicting pressure boundary conditions on set {}/{}/{}", regionName, subRegionName, setName ) ); + } + tempSubRegionSetMap.insert( setName ); + + // 2.2 Check that there is pressure bc applied to this set + auto & presSubRegionSetMap = bcStatusMap[regionName][subRegionName]; + if( presSubRegionSetMap.count( setName ) == 0 ) + { + bcConsistent = false; + GEOSX_WARNING( GEOSX_FMT( "Pressure boundary condition not prescribed on set {}/{}/{}", regionName, subRegionName, setName ) ); + } + + // no need to set the number of components here, it was done while checking pressure + } ); + } + + // 3. Check composition BC (global component fraction) fsManager.apply( time, domain.getMeshBody( 0 ).getMeshLevel( 0 ), "ElementRegions", @@ -1296,7 +1433,7 @@ bool validateDirichletBC( DomainPartition & domain, Group & subRegion, string const & ) { - // 2.0. Check pressure and record composition bc application + // 3.1 Check pressure, temperature, and record composition bc application string const & subRegionName = subRegion.getName(); string const & regionName = subRegion.getParent().getParent().getName(); integer const comp = fs.getComponent(); @@ -1307,6 +1444,15 @@ bool validateDirichletBC( DomainPartition & domain, bcConsistent = false; GEOSX_WARNING( GEOSX_FMT( "Pressure boundary condition not prescribed on set {}/{}/{}", regionName, subRegionName, setName ) ); } + if( isThermal ) + { + auto & tempSubRegionSetMap = bcTempStatusMap[regionName][subRegionName]; + if( tempSubRegionSetMap.count( setName ) == 0 ) + { + bcConsistent = false; + GEOSX_WARNING( GEOSX_FMT( "Temperature boundary condition not prescribed on set {}/{}/{}", regionName, subRegionName, setName ) ); + } + } if( comp < 0 || comp >= numComp ) { bcConsistent = false; @@ -1323,7 +1469,7 @@ bool validateDirichletBC( DomainPartition & domain, compMask.set( comp ); } ); - // 2.3 Check consistency between composition BC applied to sets + // 3.2 Check consistency between composition BC applied to sets for( auto const & regionEntry : bcStatusMap ) { for( auto const & subRegionEntry : regionEntry.second ) @@ -1362,7 +1508,7 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, // Only validate BC at the beginning of Newton loop if( m_nonlinearSolverParameters.m_numNewtonIterations == 0 ) { - bool const bcConsistent = validateDirichletBC( domain, m_numComponents, time + dt ); + bool const bcConsistent = validateDirichletBC( domain, m_isThermal, m_numComponents, time + dt ); GEOSX_ERROR_IF( !bcConsistent, GEOSX_FMT( "CompositionalMultiphaseBase {}: inconsistent boundary conditions", getName() ) ); } @@ -1393,7 +1539,27 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, extrinsicMeshData::flow::bcPressure::key() ); } ); - // 2. Apply composition BC (global component fraction) and store them for constitutive call + // 2. Apply temperature Dirichlet BCs, store in a separate field + if( m_isThermal ) + { + fsManager.apply( time + dt, + domain.getMeshBody( 0 ).getMeshLevel( 0 ), + "ElementRegions", + extrinsicMeshData::flow::temperature::key(), + [&]( FieldSpecificationBase const & fs, + string const &, + SortedArrayView< localIndex const > const & targetSet, + Group & subRegion, + string const & ) + { + fs.applyFieldValue< FieldSpecificationEqual, parallelDevicePolicy<> >( targetSet, + time + dt, + subRegion, + extrinsicMeshData::flow::bcTemperature::key() ); + } ); + } + + // 3. Apply composition BC (global component fraction) and store them for constitutive call fsManager.apply( time + dt, domain.getMeshBody( 0 ).getMeshLevel( 0 ), "ElementRegions", @@ -1413,7 +1579,7 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, globalIndex const rankOffset = dofManager.rankOffset(); string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); - // 3. Call constitutive update, back-calculate target global component densities and apply to the system + // 4. Call constitutive update, back-calculate target global component densities and apply to the system fsManager.apply( time + dt, domain.getMeshBody( 0 ).getMeshLevel( 0 ), "ElementRegions", @@ -1427,10 +1593,13 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); + // in the isothermal case, we use the reservoir temperature to enforce the boundary condition + string const temperatureKey = m_isThermal ? extrinsicMeshData::flow::bcTemperature::key() : extrinsicMeshData::flow::temperature::key(); + arrayView1d< real64 const > const bcPres = subRegion.getReference< array1d< real64 > >( extrinsicMeshData::flow::bcPressure::key() ); - arrayView1d< real64 const > const temp = - subRegion.getReference< array1d< real64 > >( extrinsicMeshData::flow::temperature::key() ); + arrayView1d< real64 const > const bcTemp = + subRegion.getReference< array1d< real64 > >( temperatureKey ); arrayView2d< real64 const, compflow::USD_COMP > const compFrac = subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( extrinsicMeshData::flow::globalCompFraction::key() ); @@ -1440,11 +1609,13 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, using ExecPolicy = typename FluidType::exec_policy; typename FluidType::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); - FluidUpdateKernel::launch< ExecPolicy >( targetSet, - fluidWrapper, - bcPres, - temp, - compFrac ); + thermalCompositionalMultiphaseBaseKernels:: + FluidUpdateKernel:: + launch< ExecPolicy >( targetSet, + fluidWrapper, + bcPres, + bcTemp, + compFrac ); } ); arrayView1d< integer const > const ghostRank = @@ -1455,6 +1626,10 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, subRegion.getReference< array1d< real64 > >( extrinsicMeshData::flow::pressure::key() ); arrayView1d< real64 const > const dPres = subRegion.getReference< array1d< real64 > >( extrinsicMeshData::flow::deltaPressure::key() ); + arrayView1d< real64 const > const temp = + subRegion.getReference< array1d< real64 > >( extrinsicMeshData::flow::temperature::key() ); + arrayView1d< real64 const > const dTemp = + subRegion.getReference< array1d< real64 > >( extrinsicMeshData::flow::deltaTemperature::key() ); arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( extrinsicMeshData::flow::globalCompDensity::key() ); arrayView2d< real64 const, compflow::USD_COMP > const dCompDens = @@ -1462,6 +1637,7 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity(); integer const numComp = m_numComponents; + integer const isThermal = m_isThermal; forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const a ) { localIndex const ei = targetSet[a]; @@ -1474,7 +1650,7 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, localIndex const localRow = dofIndex - rankOffset; real64 rhsValue; - // 3.1. Apply pressure value to the matrix/rhs + // 4.1. Apply pressure value to the matrix/rhs FieldSpecificationEqual::SpecifyFieldValue( dofIndex, rankOffset, localMatrix, @@ -1483,7 +1659,19 @@ void CompositionalMultiphaseBase::applyDirichletBC( real64 const time, pres[ei] + dPres[ei] ); localRhs[localRow] = rhsValue; - // 3.2. For each component, apply target global density value + // 4.2. Apply temperature value to the matrix/rhs + if( isThermal ) + { + FieldSpecificationEqual::SpecifyFieldValue( dofIndex + numComp + 1, + rankOffset, + localMatrix, + rhsValue, + bcTemp[ei], + temp[ei] + dTemp[ei] ); + localRhs[localRow + numComp + 1] = rhsValue; + } + + // 4.3. For each component, apply target global density value for( integer ic = 0; ic < numComp; ++ic ) { FieldSpecificationEqual::SpecifyFieldValue( dofIndex + ic + 1, @@ -1516,6 +1704,8 @@ void CompositionalMultiphaseBase::chopNegativeDensities( DomainPartition & domai { GEOSX_MARK_FUNCTION; + using namespace isothermalCompositionalMultiphaseBaseKernels; + integer const numComp = m_numComponents; forMeshTargets( domain.getMeshBodies(), [&]( string const &, @@ -1572,10 +1762,23 @@ void CompositionalMultiphaseBase::resetStateToBeginningOfStep( DomainPartition & dPres.zero(); dCompDens.zero(); - // update porosity and permeability + if( m_isThermal ) + { + arrayView1d< real64 > const & dTemp = + subRegion.template getExtrinsicData< extrinsicMeshData::flow::deltaTemperature >(); + dTemp.zero(); + } + + // update porosity, permeability updatePorosityAndPermeability( subRegion ); // update all fluid properties updateFluidState( subRegion ); + // for thermal simulations, update solid internal energy + if( m_isThermal ) + { + updateSolidInternalEnergyModel( subRegion ); + } + } ); } ); } @@ -1614,12 +1817,25 @@ void CompositionalMultiphaseBase::implicitStepComplete( real64 const & time, forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { pres[ei] += dPres[ei]; - for( localIndex ic = 0; ic < numComp; ++ic ) + for( integer ic = 0; ic < numComp; ++ic ) { compDens[ei][ic] += dCompDens[ei][ic]; } } ); + if( m_isThermal ) + { + arrayView1d< real64 const > const dTemp = + subRegion.getExtrinsicData< extrinsicMeshData::flow::deltaTemperature >(); + arrayView1d< real64 > const temp = + subRegion.getExtrinsicData< extrinsicMeshData::flow::temperature >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) + { + temp[ei] += dTemp[ei]; + } ); + } + // Step 3: save the converged fluid state string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase const & fluidMaterial = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); @@ -1680,10 +1896,15 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain ) SurfaceElementSubRegion >( regionNames, [&]( localIndex const, auto & subRegion ) { - // update porosity and permeability + // update porosity, permeability, and solid internal energy updatePorosityAndPermeability( subRegion ); // update all fluid properties updateFluidState( subRegion ); + // for thermal, update solid internal energy + if( m_isThermal ) + { + updateSolidInternalEnergyModel( subRegion ); + } } ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index b35550e9d56..bcdd57de806 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -138,19 +138,25 @@ class CompositionalMultiphaseBase : public FlowSolverBase /** * @brief Update all relevant relperm models using current values of phase volume fraction - * @param castedRelPerm the group storing the required fields + * @param dataGroup the group storing the required fields */ - void updateRelPermModel( ObjectManagerBase & castedRelPerm ) const; + void updateRelPermModel( ObjectManagerBase & dataGroup ) const; /** * @brief Update all relevant capillary pressure models using current values of phase volume fraction - * @param castedCapPres the group storing the required fields + * @param dataGroup the group storing the required fields */ - void updateCapPressureModel( ObjectManagerBase & castedCapPres ) const; + void updateCapPressureModel( ObjectManagerBase & dataGroup ) const; + + /** + * @brief Update all relevant solid internal energy models using current values of temperature + * @param dataGroup the group storing the required fields + */ + void updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const; /** * @brief Recompute phase mobility from constitutive and primary variables - * @param domain the domain containing the mesh and fields + * @param dataGroup the group storing the required field */ virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const = 0; @@ -162,13 +168,13 @@ class CompositionalMultiphaseBase : public FlowSolverBase * @brief Getter for the number of fluid components (species) * @return the number of components */ - localIndex numFluidComponents() const { return m_numComponents; } + integer numFluidComponents() const { return m_numComponents; } /** * @brief Getter for the number of fluid phases * @return the number of phases */ - localIndex numFluidPhases() const { return m_numPhases; } + integer numFluidPhases() const { return m_numPhases; } /** * @brief Getter for the name of the reference fluid model name @@ -219,6 +225,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase static constexpr char const * useMassFlagString() { return "useMass"; } + static constexpr char const * isThermalString() { return "isThermal"; } + static constexpr char const * computeCFLNumbersString() { return "computeCFLNumbers"; } static constexpr char const * relPermNamesString() { return "relPermNames"; } @@ -227,6 +235,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase static constexpr char const * thermalConductivityNamesString() { return "thermalConductivityNames"; } + static constexpr char const * solidInternalEnergyNamesString() { return "solidInternalEnergyNames"; } + static constexpr char const * maxCompFracChangeString() { return "maxCompFractionChange"; } static constexpr char const * allowLocalCompDensChoppingString() { return "allowLocalCompDensityChopping"; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp index 363b6a818c3..a4d54c5f472 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp @@ -207,6 +207,15 @@ EXTRINSIC_MESH_DATA_TRAIT( bcPressure, WRITE_AND_READ, "Boundary condition pressure" ); +EXTRINSIC_MESH_DATA_TRAIT( bcTemperature, + "bcTemperature", + array1d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Boundary temperature" ); + + } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 3176ef0b34d..a9cdc2da7ca 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -34,17 +34,17 @@ #include "mesh/DomainPartition.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseExtrinsicData.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" +#include "physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp" +#include "physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp" namespace geosx { using namespace dataRepository; using namespace constitutive; -using namespace compositionalMultiphaseFVMKernels; -using namespace compositionalMultiphaseBaseKernels; CompositionalMultiphaseFVM::CompositionalMultiphaseFVM( const string & name, Group * const parent ) @@ -108,18 +108,38 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, { typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); - FaceBasedAssemblyKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dofManager.rankOffset(), - elemDofKey, - m_hasCapPressure, - getName(), - mesh.getElemManager(), - stencilWrapper, - dt, - localMatrix.toViewConstSizes(), - localRhs.toView() ); + if( m_isThermal ) + { + thermalCompositionalMultiphaseFVMKernels:: + FaceBasedAssemblyKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + m_hasCapPressure, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } + else + { + isothermalCompositionalMultiphaseFVMKernels:: + FaceBasedAssemblyKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dofManager.rankOffset(), + elemDofKey, + m_hasCapPressure, + getName(), + mesh.getElemManager(), + stencilWrapper, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); + } } ); } ); } @@ -163,10 +183,14 @@ void CompositionalMultiphaseFVM::computeCFLNumbers( real64 const & dt, FiniteVolumeManager & fvManager = numericalMethodManager.getFiniteVolumeManager(); FluxApproximationBase & fluxApprox = fvManager.getFluxApproximation( m_discretizationName ); - CFLFluxKernel::CompFlowAccessors compFlowAccessors( mesh.getElemManager(), getName() ); - CFLFluxKernel::MultiFluidAccessors multiFluidAccessors( mesh.getElemManager(), getName() ); - CFLFluxKernel::PermeabilityAccessors permeabilityAccessors( mesh.getElemManager(), getName() ); - CFLFluxKernel::RelPermAccessors relPermAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + CFLFluxKernel::CompFlowAccessors compFlowAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + CFLFluxKernel::MultiFluidAccessors multiFluidAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + CFLFluxKernel::PermeabilityAccessors permeabilityAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + CFLFluxKernel::RelPermAccessors relPermAccessors( mesh.getElemManager(), getName() ); // TODO: find a way to compile with this modifiable accessors in CompFlowAccessors, and remove them from here ElementRegionManager::ElementViewAccessor< arrayView2d< real64, compflow::USD_PHASE > > const phaseOutfluxAccessor = @@ -177,28 +201,30 @@ void CompositionalMultiphaseFVM::computeCFLNumbers( real64 const & dt, mesh.getElemManager().constructViewAccessor< array2d< real64, compflow::LAYOUT_COMP >, arrayView2d< real64, compflow::USD_COMP > >( extrinsicMeshData::flow::componentOutflux::key() ); + fluxApprox.forAllStencils( mesh, [&] ( auto & stencil ) { typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); // While this kernel is waiting for a factory class, pass all the accessors here - KernelLaunchSelector1< CFLFluxKernel >( m_numComponents, - m_numPhases, - dt, - stencilWrapper, - compFlowAccessors.get( extrinsicMeshData::flow::pressure{} ), - compFlowAccessors.get( extrinsicMeshData::flow::gravityCoefficient{} ), - compFlowAccessors.get( extrinsicMeshData::flow::phaseVolumeFraction{} ), - permeabilityAccessors.get( extrinsicMeshData::permeability::permeability{} ), - permeabilityAccessors.get( extrinsicMeshData::permeability::dPerm_dPressure{} ), - relPermAccessors.get( extrinsicMeshData::relperm::phaseRelPerm{} ), - multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseViscosity{} ), - multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseDensity{} ), - multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseMassDensity{} ), - multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseCompFraction{} ), - phaseOutfluxAccessor.toNestedView(), - compOutfluxAccessor.toNestedView() ); + isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector1 + < isothermalCompositionalMultiphaseFVMKernels::CFLFluxKernel >( m_numComponents, + m_numPhases, + dt, + stencilWrapper, + compFlowAccessors.get( extrinsicMeshData::flow::pressure{} ), + compFlowAccessors.get( extrinsicMeshData::flow::gravityCoefficient{} ), + compFlowAccessors.get( extrinsicMeshData::flow::phaseVolumeFraction{} ), + permeabilityAccessors.get( extrinsicMeshData::permeability::permeability{} ), + permeabilityAccessors.get( extrinsicMeshData::permeability::dPerm_dPressure{} ), + relPermAccessors.get( extrinsicMeshData::relperm::phaseRelPerm{} ), + multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseViscosity{} ), + multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseDensity{} ), + multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseMassDensity{} ), + multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseCompFraction{} ), + phaseOutfluxAccessor.toNestedView(), + compOutfluxAccessor.toNestedView() ); } ); } ); @@ -248,22 +274,23 @@ void CompositionalMultiphaseFVM::computeCFLNumbers( real64 const & dt, arrayView2d< real64 const > const & porosity = solidModel.getPorosity(); - KernelLaunchSelector2< CFLKernel >( m_numComponents, m_numPhases, - subRegion.size(), - volume, - porosity, - compDens, - compFrac, - phaseVolFrac, - phaseRelPerm, - dPhaseRelPerm_dPhaseVolFrac, - phaseVisc, - phaseOutflux, - compOutflux, - phaseCFLNumber, - compCFLNumber, - subRegionMaxPhaseCFLNumber, - subRegionMaxCompCFLNumber ); + isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector2 + < isothermalCompositionalMultiphaseFVMKernels::CFLKernel >( m_numComponents, m_numPhases, + subRegion.size(), + volume, + porosity, + compDens, + compFrac, + phaseVolFrac, + phaseRelPerm, + dPhaseRelPerm_dPhaseVolFrac, + phaseVisc, + phaseOutflux, + compOutflux, + phaseCFLNumber, + compCFLNumber, + subRegionMaxPhaseCFLNumber, + subRegionMaxCompCFLNumber ); localMaxPhaseCFLNumber = LvArray::math::max( localMaxPhaseCFLNumber, subRegionMaxPhaseCFLNumber ); localMaxCompCFLNumber = LvArray::math::max( localMaxCompCFLNumber, subRegionMaxCompCFLNumber ); @@ -284,7 +311,8 @@ real64 CompositionalMultiphaseFVM::calculateResidualNorm( DomainPartition const { GEOSX_MARK_FUNCTION; - real64 localResidualNorm = 0.0; + real64 localFlowResidualNorm = 0.0; + real64 localEnergyResidualNorm = 0.0; globalIndex const rankOffset = dofManager.rankOffset(); string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); @@ -297,6 +325,7 @@ real64 CompositionalMultiphaseFVM::calculateResidualNorm( DomainPartition const [&]( localIndex const, ElementSubRegionBase const & subRegion ) { + arrayView1d< globalIndex const > dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank(); arrayView1d< real64 const > const volume = subRegion.getElementVolume(); @@ -306,32 +335,80 @@ real64 CompositionalMultiphaseFVM::calculateResidualNorm( DomainPartition const arrayView2d< real64 const, multifluid::USD_FLUID > const totalDensOld = fluid.totalDensityOld(); string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() ); - CoupledSolidBase const & solid = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - arrayView1d< real64 const > const referencePorosity = solid.getReferencePorosity(); - - real64 subRegionResidualNorm = 0.0; - ResidualNormKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localRhs, - rankOffset, - numFluidComponents(), - dofNumber, - elemGhostRank, - referencePorosity, - volume, - totalDensOld, - subRegionResidualNorm ); - localResidualNorm += subRegionResidualNorm; + CoupledSolidBase const & solidModel = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); + arrayView1d< real64 const > const referencePorosity = solidModel.getReferencePorosity(); + + real64 subRegionFlowResidualNorm = 0.0; + real64 subRegionEnergyResidualNorm = 0.0; + + if( m_isThermal ) + { + arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFracOld = + subRegion.getExtrinsicData< extrinsicMeshData::flow::phaseVolumeFractionOld >(); + arrayView3d< real64 const, multifluid::USD_PHASE > const phaseDensOld = fluid.phaseDensityOld(); + arrayView3d< real64 const, multifluid::USD_PHASE > const phaseInternalEnergyOld = fluid.phaseInternalEnergyOld(); + + string const & solidInternalEnergyName = subRegion.getReference< string >( viewKeyStruct::solidInternalEnergyNamesString() ); + SolidInternalEnergy const & solidInternalEnergy = getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); + arrayView2d< real64 const > const solidInternalEnergyOld = solidInternalEnergy.getOldInternalEnergy(); + + thermalCompositionalMultiphaseBaseKernels:: + ResidualNormKernel:: + launch< parallelDevicePolicy<> >( localRhs, + rankOffset, + numFluidPhases(), + numFluidComponents(), + dofNumber, + elemGhostRank, + referencePorosity, + volume, + solidInternalEnergyOld, + phaseVolFracOld, + totalDensOld, + phaseDensOld, + phaseInternalEnergyOld, + subRegionFlowResidualNorm, + subRegionEnergyResidualNorm ); + } + else + { + isothermalCompositionalMultiphaseBaseKernels:: + ResidualNormKernel:: + launch< parallelDevicePolicy<> >( localRhs, + rankOffset, + numFluidComponents(), + dofNumber, + elemGhostRank, + referencePorosity, + volume, + totalDensOld, + subRegionFlowResidualNorm ); + } + localFlowResidualNorm += subRegionFlowResidualNorm; + localEnergyResidualNorm += subRegionEnergyResidualNorm; } ); } ); - // compute global residual norm - real64 const residual = std::sqrt( MpiWrapper::sum( localResidualNorm ) ); - - if( getLogLevel() >= 1 && logger::internal::rank == 0 ) + // compute global residual norms + real64 residual = 0.0; + if( m_isThermal ) { - std::cout << GEOSX_FMT( " ( Rfluid ) = ( {:4.2e} ) ;", residual ); + real64 const flowResidual = std::sqrt( MpiWrapper::sum( localFlowResidualNorm ) ); + real64 const energyResidual = std::sqrt( MpiWrapper::sum( localEnergyResidualNorm ) ); + residual = std::sqrt( flowResidual*flowResidual + energyResidual*energyResidual ); + if( getLogLevel() >= 1 && logger::internal::rank == 0 ) + { + std::cout << GEOSX_FMT( " ( Rfluid ) = ( {:4.2e} ) ; ( Renergy ) = ( {:4.2e} ) ; ", flowResidual, energyResidual ); + } + } + else + { + residual = std::sqrt( MpiWrapper::sum( localFlowResidualNorm ) ); + if( getLogLevel() >= 1 && logger::internal::rank == 0 ) + { + std::cout << GEOSX_FMT( " ( Rfluid ) = ( {:4.2e} ) ; ", residual ); + } } - return residual; } @@ -348,7 +425,7 @@ real64 CompositionalMultiphaseFVM::scalingForSystemSolution( DomainPartition con return 1.0; } - real64 constexpr eps = compositionalMultiphaseBaseKernels::minDensForDivision; + real64 constexpr eps = isothermalCompositionalMultiphaseBaseKernels::minDensForDivision; real64 const maxCompFracChange = m_maxCompFracChange; localIndex const NC = m_numComponents; @@ -448,18 +525,18 @@ bool CompositionalMultiphaseFVM::checkSystemSolution( DomainPartition const & do subRegion.getExtrinsicData< extrinsicMeshData::flow::deltaGlobalCompDensity >(); localIndex const subRegionSolutionCheck = - SolutionCheckKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localSolution, - dofManager.rankOffset(), - numFluidComponents(), - dofNumber, - elemGhostRank, - pres, - dPres, - compDens, - dCompDens, - m_allowCompDensChopping, - scalingFactor ); + isothermalCompositionalMultiphaseBaseKernels:: + SolutionCheckKernel::launch< parallelDevicePolicy<> >( localSolution, + dofManager.rankOffset(), + numFluidComponents(), + dofNumber, + elemGhostRank, + pres, + dPres, + compDens, + dCompDens, + m_allowCompDensChopping, + scalingFactor ); localCheck = std::min( localCheck, subRegionSolutionCheck ); } ); @@ -476,6 +553,7 @@ void CompositionalMultiphaseFVM::applySystemSolution( DofManager const & dofMana GEOSX_MARK_FUNCTION; DofManager::CompMask pressureMask( m_numDofPerCell, 0, 1 ); + DofManager::CompMask componentMask( m_numDofPerCell, 1, m_numComponents+1 ); dofManager.addVectorToField( localSolution, viewKeyStruct::elemDofFieldString(), @@ -487,7 +565,17 @@ void CompositionalMultiphaseFVM::applySystemSolution( DofManager const & dofMana viewKeyStruct::elemDofFieldString(), extrinsicMeshData::flow::deltaGlobalCompDensity::key(), scalingFactor, - ~pressureMask ); + componentMask ); + + if( m_isThermal ) + { + DofManager::CompMask temperatureMask( m_numDofPerCell, m_numComponents+1, m_numComponents+2 ); + dofManager.addVectorToField( localSolution, + viewKeyStruct::elemDofFieldString(), + extrinsicMeshData::flow::deltaTemperature::key(), + scalingFactor, + temperatureMask ); + } // if component density chopping is allowed, some component densities may be negative after the update // these negative component densities are set to zero in this function @@ -495,6 +583,7 @@ void CompositionalMultiphaseFVM::applySystemSolution( DofManager const & dofMana { chopNegativeDensities( domain ); } + forMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, arrayView1d< string const > const & ) @@ -502,6 +591,12 @@ void CompositionalMultiphaseFVM::applySystemSolution( DofManager const & dofMana std::map< string, string_array > fieldNames; fieldNames["elems"].emplace_back( extrinsicMeshData::flow::deltaPressure::key() ); fieldNames["elems"].emplace_back( extrinsicMeshData::flow::deltaGlobalCompDensity::key() ); + + if( m_isThermal ) + { + fieldNames["elems"].emplace_back( extrinsicMeshData::flow::deltaTemperature::key() ); + } + CommunicationTools::getInstance().synchronizeFields( fieldNames, mesh, domain.getNeighbors(), true ); } ); } @@ -517,12 +612,26 @@ void CompositionalMultiphaseFVM::updatePhaseMobility( ObjectManagerBase & dataGr string const & relpermName = dataGroup.getReference< string >( viewKeyStruct::relPermNamesString() ); RelativePermeabilityBase const & relperm = getConstitutiveModel< RelativePermeabilityBase >( dataGroup, relpermName ); - PhaseMobilityKernelFactory:: - createAndLaunch< parallelDevicePolicy<> >( m_numComponents, - m_numPhases, - dataGroup, - fluid, - relperm ); + if( m_isThermal ) + { + thermalCompositionalMultiphaseFVMKernels:: + PhaseMobilityKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid, + relperm ); + } + else + { + isothermalCompositionalMultiphaseFVMKernels:: + PhaseMobilityKernelFactory:: + createAndLaunch< parallelDevicePolicy<> >( m_numComponents, + m_numPhases, + dataGroup, + fluid, + relperm ); + } } void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, @@ -535,6 +644,7 @@ void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, GEOSX_MARK_FUNCTION; FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); + forMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, arrayView1d< string const > const & ) @@ -548,8 +658,10 @@ void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, mesh.getElemManager().constructArrayViewAccessor< globalIndex, 1 >( elemDofKey ); elemDofNumber.setName( getName() + "/accessors/" + elemDofKey ); - AquiferBCKernel::CompFlowAccessors compFlowAccessors( mesh.getElemManager(), getName() ); - AquiferBCKernel::MultiFluidAccessors multiFluidAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + AquiferBCKernel::CompFlowAccessors compFlowAccessors( mesh.getElemManager(), getName() ); + isothermalCompositionalMultiphaseFVMKernels:: + AquiferBCKernel::MultiFluidAccessors multiFluidAccessors( mesh.getElemManager(), getName() ); fsManager.apply< AquiferBoundaryCondition >( time + dt, mesh, @@ -586,32 +698,33 @@ void CompositionalMultiphaseFVM::applyAquiferBC( real64 const time, arrayView1d< real64 const > const & aquiferWaterPhaseCompFrac = bc.getWaterPhaseComponentFraction(); // While this kernel is waiting for a factory class, pass all the accessors here - KernelLaunchSelector1< AquiferBCKernel >( m_numComponents, - m_numPhases, - waterPhaseIndex, - allowAllPhasesIntoAquifer, - stencil, - dofManager.rankOffset(), - elemDofNumber.toNestedViewConst(), - aquiferBCWrapper, - aquiferWaterPhaseDens, - aquiferWaterPhaseCompFrac, - compFlowAccessors.get( extrinsicMeshData::ghostRank{} ), - compFlowAccessors.get( extrinsicMeshData::flow::pressure{} ), - compFlowAccessors.get( extrinsicMeshData::flow::deltaPressure{} ), - compFlowAccessors.get( extrinsicMeshData::flow::gravityCoefficient{} ), - compFlowAccessors.get( extrinsicMeshData::flow::phaseVolumeFraction{} ), - compFlowAccessors.get( extrinsicMeshData::flow::dPhaseVolumeFraction_dPressure{} ), - compFlowAccessors.get( extrinsicMeshData::flow::dPhaseVolumeFraction_dGlobalCompDensity{} ), - compFlowAccessors.get( extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity{} ), - multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseDensity{} ), - multiFluidAccessors.get( extrinsicMeshData::multifluid::dPhaseDensity{} ), - multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseCompFraction{} ), - multiFluidAccessors.get( extrinsicMeshData::multifluid::dPhaseCompFraction{} ), - time, - dt, - localMatrix.toViewConstSizes(), - localRhs.toView() ); + isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector1 + < isothermalCompositionalMultiphaseFVMKernels::AquiferBCKernel >( m_numComponents, + m_numPhases, + waterPhaseIndex, + allowAllPhasesIntoAquifer, + stencil, + dofManager.rankOffset(), + elemDofNumber.toNestedViewConst(), + aquiferBCWrapper, + aquiferWaterPhaseDens, + aquiferWaterPhaseCompFrac, + compFlowAccessors.get( extrinsicMeshData::ghostRank{} ), + compFlowAccessors.get( extrinsicMeshData::flow::pressure{} ), + compFlowAccessors.get( extrinsicMeshData::flow::deltaPressure{} ), + compFlowAccessors.get( extrinsicMeshData::flow::gravityCoefficient{} ), + compFlowAccessors.get( extrinsicMeshData::flow::phaseVolumeFraction{} ), + compFlowAccessors.get( extrinsicMeshData::flow::dPhaseVolumeFraction_dPressure{} ), + compFlowAccessors.get( extrinsicMeshData::flow::dPhaseVolumeFraction_dGlobalCompDensity{} ), + compFlowAccessors.get( extrinsicMeshData::flow::dGlobalCompFraction_dGlobalCompDensity{} ), + multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseDensity{} ), + multiFluidAccessors.get( extrinsicMeshData::multifluid::dPhaseDensity{} ), + multiFluidAccessors.get( extrinsicMeshData::multifluid::phaseCompFraction{} ), + multiFluidAccessors.get( extrinsicMeshData::multifluid::dPhaseCompFraction{} ), + time, + dt, + localMatrix.toViewConstSizes(), + localRhs.toView() ); } ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index ed1fed23a12..215ab9690ca 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -28,9 +28,9 @@ #include "finiteVolume/MimeticInnerProductDispatch.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseExtrinsicData.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp" /** @@ -41,7 +41,7 @@ namespace geosx using namespace dataRepository; using namespace constitutive; -using namespace compositionalMultiphaseBaseKernels; +using namespace isothermalCompositionalMultiphaseBaseKernels; using namespace compositionalMultiphaseHybridFVMKernels; using namespace mimeticInnerProduct; @@ -513,7 +513,7 @@ real64 CompositionalMultiphaseHybridFVM::scalingForSystemSolution( DomainPartiti return 1.0; } - real64 constexpr eps = compositionalMultiphaseBaseKernels::minDensForDivision; + real64 constexpr eps = isothermalCompositionalMultiphaseBaseKernels::minDensForDivision; real64 const maxCompFracChange = m_maxCompFracChange; real64 const maxRelativePresChange = m_maxRelativePresChange; @@ -675,18 +675,18 @@ bool CompositionalMultiphaseHybridFVM::checkSystemSolution( DomainPartition cons subRegion.getExtrinsicData< extrinsicMeshData::flow::deltaGlobalCompDensity >(); localIndex const subRegionSolutionCheck = - compositionalMultiphaseBaseKernels::SolutionCheckKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localSolution, - dofManager.rankOffset(), - numFluidComponents(), - elemDofNumber, - elemGhostRank, - elemPres, - dElemPres, - compDens, - dCompDens, - m_allowCompDensChopping, - scalingFactor ); + isothermalCompositionalMultiphaseBaseKernels:: + SolutionCheckKernel::launch< parallelDevicePolicy<> >( localSolution, + dofManager.rankOffset(), + numFluidComponents(), + elemDofNumber, + elemGhostRank, + elemPres, + dElemPres, + compDens, + dCompDens, + m_allowCompDensChopping, + scalingFactor ); localCheck = std::min( localCheck, subRegionSolutionCheck ); } ); @@ -703,14 +703,14 @@ bool CompositionalMultiphaseHybridFVM::checkSystemSolution( DomainPartition cons faceManager.getExtrinsicData< extrinsicMeshData::flow::deltaFacePressure >(); localIndex const faceSolutionCheck = - compositionalMultiphaseHybridFVMKernels::SolutionCheckKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localSolution, - dofManager.rankOffset(), - faceDofNumber, - faceGhostRank, - facePres, - dFacePres, - scalingFactor ); + compositionalMultiphaseHybridFVMKernels:: + SolutionCheckKernel::launch< parallelDevicePolicy<> >( localSolution, + dofManager.rankOffset(), + faceDofNumber, + faceGhostRank, + facePres, + dFacePres, + scalingFactor ); localCheck = std::min( localCheck, faceSolutionCheck ); } ); @@ -807,17 +807,16 @@ real64 CompositionalMultiphaseHybridFVM::calculateResidualNorm( DomainPartition arrayView1d< real64 const > const & referencePorosity = solid.getReferencePorosity(); real64 subRegionResidualNorm = 0.0; - compositionalMultiphaseBaseKernels:: - ResidualNormKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localRhs, - rankOffset, - numFluidComponents(), - elemDofNumber, - elemGhostRank, - referencePorosity, - volume, - totalDensOld, - subRegionResidualNorm ); + isothermalCompositionalMultiphaseBaseKernels:: + ResidualNormKernel::launch< parallelDevicePolicy<> >( localRhs, + rankOffset, + numFluidComponents(), + elemDofNumber, + elemGhostRank, + referencePorosity, + volume, + totalDensOld, + subRegionResidualNorm ); localResidualNorm += subRegionResidualNorm; } ); @@ -832,19 +831,18 @@ real64 CompositionalMultiphaseHybridFVM::calculateResidualNorm( DomainPartition // 2. Compute the residual for the face-based constraints real64 faceResidualNorm = 0.0; compositionalMultiphaseHybridFVMKernels:: - ResidualNormKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localRhs, - rankOffset, - numFluidPhases(), - faceDofNumber.toNestedViewConst(), - faceGhostRank.toNestedViewConst(), - m_regionFilter.toViewConst(), - elemRegionList.toNestedViewConst(), - elemSubRegionList.toNestedViewConst(), - elemList.toNestedViewConst(), - compFlowAccessors.get( extrinsicMeshData::elementVolume{} ), - compFlowAccessors.get( extrinsicMeshData::flow::phaseMobilityOld{} ), - faceResidualNorm ); + ResidualNormKernel::launch< parallelDevicePolicy<> >( localRhs, + rankOffset, + numFluidPhases(), + faceDofNumber.toNestedViewConst(), + faceGhostRank.toNestedViewConst(), + m_regionFilter.toViewConst(), + elemRegionList.toNestedViewConst(), + elemSubRegionList.toNestedViewConst(), + elemList.toNestedViewConst(), + compFlowAccessors.get( extrinsicMeshData::elementVolume{} ), + compFlowAccessors.get( extrinsicMeshData::flow::phaseMobilityOld{} ), + faceResidualNorm ); localResidualNorm += faceResidualNorm; } ); @@ -867,6 +865,7 @@ void CompositionalMultiphaseHybridFVM::applySystemSolution( DofManager const & d DomainPartition & domain ) { GEOSX_MARK_FUNCTION; + // 1. apply the elem-based update DofManager::CompMask pressureMask( m_numDofPerCell, 0, 1 ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp index ed2d10731ef..22830ba8096 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels.hpp @@ -27,9 +27,10 @@ #include "mesh/ElementRegionManager.hpp" #include "mesh/ObjectManagerBase.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/StencilAccessors.hpp" + namespace geosx { @@ -670,11 +671,11 @@ struct FluxKernel * @brief Define the interface for the property kernel in charge of computing the phase mobilities */ template< integer NUM_COMP, integer NUM_PHASE > -class PhaseMobilityKernel : public compositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > +class PhaseMobilityKernel : public isothermalCompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { public: - using Base = compositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; + using Base = isothermalCompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; using Base::numComp; /// Compile time value for the number of phases @@ -709,10 +710,10 @@ class PhaseMobilityKernel : public compositionalMultiphaseBaseKernels::PropertyK * @param[in] ei the element index * @param[in] phaseMobilityKernelOp the function used to customize the kernel */ - template< typename FUNC = compositionalMultiphaseBaseKernels::NoOpFunc > + template< typename FUNC = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc > GEOSX_HOST_DEVICE void compute( localIndex const ei, - FUNC && phaseMobilityKernelOp = compositionalMultiphaseBaseKernels::NoOpFunc{} ) const + FUNC && phaseMobilityKernelOp = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc{} ) const { using Deriv = multifluid::DerivativeOffset; @@ -841,7 +842,7 @@ class PhaseMobilityKernelFactory { if( numPhase == 2 ) { - compositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { integer constexpr NUM_COMP = NC(); PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); @@ -850,7 +851,7 @@ class PhaseMobilityKernelFactory } else if( numPhase == 3 ) { - compositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { integer constexpr NUM_COMP = NC(); PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); @@ -867,7 +868,7 @@ struct ResidualNormKernel template< typename VIEWTYPE > using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; - template< typename POLICY, typename REDUCE_POLICY > + template< typename POLICY > static void launch( arrayView1d< real64 const > const & localResidual, globalIndex const rankOffset, @@ -882,7 +883,7 @@ struct ResidualNormKernel ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMobOld, real64 & localResidualNorm ) { - RAJA::ReduceSum< REDUCE_POLICY, real64 > sumScaled( 0.0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > sumScaled( 0.0 ); forAll< POLICY >( facePresDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iface ) { @@ -932,7 +933,7 @@ struct ResidualNormKernel struct SolutionCheckKernel { - template< typename POLICY, typename REDUCE_POLICY > + template< typename POLICY > static localIndex launch( arrayView1d< real64 const > const & localSolution, globalIndex const rankOffset, @@ -942,7 +943,7 @@ struct SolutionCheckKernel arrayView1d< real64 const > const & dFacePres, real64 const scalingFactor ) { - RAJA::ReduceMin< REDUCE_POLICY, integer > check( 1 ); + RAJA::ReduceMin< ReducePolicy< POLICY >, integer > check( 1 ); forAll< POLICY >( dofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iface ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp index 8b28e24f829..d4c53d24842 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp @@ -28,31 +28,52 @@ namespace compositionalMultiphaseUtilities { /** - * @brief Shifts all elements one position ahead and replaces the first element - * with the sum of all elements in each block of the input block one-dimensional - * array of values. + * @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. * * @tparam VEC type of one-dimensional array of values - * @param N block size - * @param NB number of blocks + * @param numRowsToShift number of rows to shift in each block + * @param numRowsInBlock block size + * @param numBlocks number of blocks * @param v block one-dimensional array of values + * + * @detail numRowsToShift is different from numRowsInBlock if there is an equation that we do *not* want to + * modify, like the energy flux for thermal simulations. In that specific case, numRowsToShift + * is set to numRowsInBlock-1 to make sure that we don't modify the last row of each block (the energy flux) + * + * Let's consider the following blocks arising when computing the localFlux in the thermal compositional solver (for two-component flow) + * + * + * (mass_comp_1)_1 \ + * (mass_comp_2)_1 |=> block 1 + * (energy)_1 / + * + * (mass_comp_1)_2 \ + * (mass_comp_2)_2 |=> block 2 + * (energy)_2 / + * + * In this case: + * - numRowsToShift = 2 + * - numRowsInBlock = 3 + * - numBlocks = 2 */ template< typename VEC > GEOSX_HOST_DEVICE -void shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( integer const N, - integer const NB, +void shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( integer const numRowsToShift, + integer const numRowsInBlock, + integer const numBlocks, VEC && v ) { - for( integer i = 0; i < NB; ++i ) + for( integer i = 0; i < numBlocks; ++i ) { - integer const ind = i * N + N - 1; + integer const ind = i * numRowsInBlock + numRowsToShift - 1; real64 tmp = v[ind]; - for( int j = ind - 1; j >= i * N; --j ) + for( int j = ind - 1; j >= i * numRowsInBlock; --j ) { v[j+1] = v[j]; tmp += v[j]; } - v[i*N] = tmp; + v[i*numRowsInBlock] = tmp; } } @@ -61,56 +82,81 @@ void shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( integer const N, * with the sum of all elements in the input one-dimensional array of values. * * @tparam VEC type of one-dimensional array of values - * @param N vector size + * @param numRows vector size * @param v one-dimensional array of values */ template< typename VEC > GEOSX_HOST_DEVICE -void shiftElementsAheadByOneAndReplaceFirstElementWithSum( integer const N, +void shiftElementsAheadByOneAndReplaceFirstElementWithSum( integer const numRows, VEC && v ) { - shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( N, 1, v ); + shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numRows, numRows, 1, v ); } /** - * @brief Shifts all rows one position ahead and replaces the first row - * with the sum of all rows in each block of the input block row + * @brief In each block, shift the elements from 0 to numRowsToShift-1, shifts all rows one position ahead + * and replaces the first row with the sum of all rows from 0 to numRowsToShift-1 of the input block row * two-dimensional array of values. * * @tparam MATRIX type of two-dimensional array of values * @tparam VEC type of one-dimensional array of values - * @param M block number of rows - * @param N block number of columns - * @param NB number of row blocks + * @param numRowsToShift number of rows to shift in each block + * @param numRowsInBlock block number of rows + * @param numColsInBlock block number of columns + * @param numBlocks number of row blocks * @param mat block row two-dimensional array of values - * @param work one-dimensional working array of values of size N + * @param work one-dimensional working array of values of size numColsInBlock + * + * @detail numRowsToShift is different from numRowsInBlock if there is an equation that we do *not* want to + * modify, like the energy flux for thermal simulations. In that specific case, numRowsToShift + * is set to numRowsInBlock-1 to make sure that we don't modify the last row of each block (the energy flux) + * + * Let's consider the following blocks arising when computing the localFluxJacobian in the thermal compositional solver (for two-component + * flow) + * + * p \rho_1 \rho_2 T + * (mass_comp_1)_1 . . . . \ + * (mass_comp_2)_1 . . . . |=> block 1 + * (energy)_1 . . . . / + * + * (mass_comp_1)_2 . . . . \ + * (mass_comp_2)_2 . . . . |=> block 2 + * (energy)_2 . . . . / + * + * In this case: + * - numRowsToShift = 2 + * - numRowsInBlock = 3 + * - numColsInBlock = 4 + * - numBlocks = 2 + * */ template< typename MATRIX, typename VEC > GEOSX_HOST_DEVICE -void shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( integer const M, - integer const N, - integer const NB, +void shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( integer const numRowsToShift, + integer const numRowsInBlock, + integer const numColsInBlock, + integer const numBlocks, MATRIX && mat, VEC && work ) { - for( integer k = 0; k < NB; ++k ) + for( integer k = 0; k < numBlocks; ++k ) { - integer const ind = k * M + M - 1; - for( integer j = 0; j < N; ++j ) + integer const ind = k * numRowsInBlock + numRowsToShift - 1; + for( integer j = 0; j < numColsInBlock; ++j ) { work[j] = mat[ind][j]; } - for( integer i = ind - 1; i >= k * M; --i ) + for( integer i = ind - 1; i >= k * numRowsInBlock; --i ) { - for( integer j = 0; j < N; ++j ) + for( integer j = 0; j < numColsInBlock; ++j ) { mat[i+1][j] = mat[i][j]; work[j] += mat[i][j]; } } - for( integer j = 0; j < N; ++j ) + for( integer j = 0; j < numColsInBlock; ++j ) { - mat[k*M][j] = work[j]; + mat[k*numRowsInBlock][j] = work[j]; } } } @@ -122,19 +168,19 @@ void shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( integer const M, * * @tparam MATRIX type of two-dimensional array of values * @tparam VEC type of one-dimensional array of values - * @param M number of rows - * @param N number of columns + * @param numRowsInBlock number of rows + * @param numColsInBlock number of columns * @param mat block row two-dimensional array of values - * @param work one-dimensional working array of values of size N + * @param work one-dimensional working array of values of size numColsInBlock */ template< typename MATRIX, typename VEC > GEOSX_HOST_DEVICE -void shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( integer const M, - integer const N, +void shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( integer const numRowsInBlock, + integer const numColsInBlock, MATRIX && mat, VEC && work ) { - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( M, N, 1, mat, work ); + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numRowsInBlock, numRowsInBlock, numColsInBlock, 1, mat, work ); } } // namespace compositionalMultiphaseUtilities diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp similarity index 93% rename from src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp rename to src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp index 990d20f21fa..d1e6004547f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp @@ -13,11 +13,11 @@ */ /** - * @file CompositionalMultiphaseBaseKernels.hpp + * @file IsothermalCompositionalMultiphaseBaseKernels.hpp */ -#ifndef GEOSX_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEBASEKERNELS_HPP -#define GEOSX_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEBASEKERNELS_HPP +#ifndef GEOSX_PHYSICSSOLVERS_FLUIDFLOW_ISOTHERMALCOMPOSITIONALMULTIPHASEBASEKERNELS_HPP +#define GEOSX_PHYSICSSOLVERS_FLUIDFLOW_ISOTHERMALCOMPOSITIONALMULTIPHASEBASEKERNELS_HPP #include "common/DataLayouts.hpp" #include "common/DataTypes.hpp" @@ -33,7 +33,7 @@ namespace geosx { -namespace compositionalMultiphaseBaseKernels +namespace isothermalCompositionalMultiphaseBaseKernels { using namespace constitutive; @@ -355,7 +355,7 @@ class PhaseVolumeFractionKernel : public PropertyKernelBase< NUM_COMP > // call the lambda in the phase loop to allow the reuse of the phaseVolFrac and totalDensity // possible use: assemble the derivatives wrt temperature - phaseVolFractionKernelOp( ip, phaseVolFrac[ip], totalDensity ); + phaseVolFractionKernelOp( ip, phaseVolFrac[ip], phaseDensInv, totalDensity ); // now finalize the computation by multiplying by total density for( integer jc = 0; jc < numComp; ++jc ) @@ -438,82 +438,6 @@ class PhaseVolumeFractionKernelFactory } }; -/******************************** FluidUpdateKernel ********************************/ - -struct FluidUpdateKernel -{ - template< typename POLICY, typename FLUID_WRAPPER > - static void - launch( localIndex const size, - FLUID_WRAPPER const & fluidWrapper, - arrayView1d< real64 const > const & pres, - arrayView1d< real64 const > const & temp, - arrayView2d< real64 const, compflow::USD_COMP > const & compFrac ) - { - forAll< POLICY >( size, [=] GEOSX_HOST_DEVICE ( localIndex const k ) - { - for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q ) - { - fluidWrapper.update( k, q, pres[k], temp[k], compFrac[k] ); - } - } ); - } - - template< typename POLICY, typename FLUID_WRAPPER > - static void - launch( localIndex const size, - FLUID_WRAPPER const & fluidWrapper, - arrayView1d< real64 const > const & pres, - arrayView1d< real64 const > const & dPres, - arrayView1d< real64 const > const & temp, - arrayView2d< real64 const, compflow::USD_COMP > const & compFrac ) - { - forAll< POLICY >( size, [=] GEOSX_HOST_DEVICE ( localIndex const k ) - { - for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q ) - { - fluidWrapper.update( k, q, pres[k] + dPres[k], temp[k], compFrac[k] ); - } - } ); - } - - template< typename POLICY, typename FLUID_WRAPPER > - static void - launch( SortedArrayView< localIndex const > const & targetSet, - FLUID_WRAPPER const & fluidWrapper, - arrayView1d< real64 const > const & pres, - arrayView1d< real64 const > const & dPres, - arrayView1d< real64 const > const & temp, - arrayView2d< real64 const, compflow::USD_COMP > const & compFrac ) - { - forAll< POLICY >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const a ) - { - localIndex const k = targetSet[a]; - for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q ) - { - fluidWrapper.update( k, q, pres[k] + dPres[k], temp[k], compFrac[k] ); - } - } ); - } - - template< typename POLICY, typename FLUID_WRAPPER > - static void - launch( SortedArrayView< localIndex const > const & targetSet, - FLUID_WRAPPER const & fluidWrapper, - arrayView1d< real64 const > const & pres, - arrayView1d< real64 const > const & temp, - arrayView2d< real64 const, compflow::USD_COMP > const & compFrac ) - { - forAll< POLICY >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const a ) - { - localIndex const k = targetSet[a]; - for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q ) - { - fluidWrapper.update( k, q, pres[k], temp[k], compFrac[k] ); - } - } ); - } -}; /******************************** RelativePermeabilityUpdateKernel ********************************/ @@ -796,6 +720,7 @@ class ElementBasedAssemblyKernel + phaseCompFrac[ip][ic] * dPhaseAmount_dC[jc]; stack.localJacobian[ic][jc + 1] += dPhaseCompAmount_dC; } + } // call the lambda in the phase loop to allow the reuse of the phase amounts and their derivatives @@ -806,7 +731,7 @@ class ElementBasedAssemblyKernel } /** - * @brief Compute the local accumulation contributions to the residual and Jacobian + * @brief Compute the local volume balance contributions to the residual and Jacobian * @tparam FUNC the type of the function that can be used to customize the kernel * @param[in] ei the element index * @param[inout] stack the stack variables @@ -836,17 +761,17 @@ class ElementBasedAssemblyKernel } } + // call the lambda in the phase loop to allow the reuse of the phase amounts and their derivatives + // possible use: assemble the derivatives wrt temperature, and use oneMinusPhaseVolFracSum if poreVolumeNew depends on temperature + phaseVolFractionSumKernelOp( oneMinusPhaseVolFracSum ); + // scale saturation-based volume balance by pore volume (for better scaling w.r.t. other equations) stack.localResidual[numComp] = stack.poreVolumeNew * oneMinusPhaseVolFracSum; - for( integer idof = 0; idof < numComp+1; ++idof ) + for( integer idof = 0; idof < numDof; ++idof ) { stack.localJacobian[numComp][idof] *= stack.poreVolumeNew; } stack.localJacobian[numComp][0] += stack.dPoreVolume_dPres * oneMinusPhaseVolFracSum; - - // call the lambda in the phase loop to allow the reuse of the phase amounts and their derivatives - // possible use: assemble the derivatives wrt temperature, and use oneMinusPhaseVolFracSum if poreVolumeNew depends on temperature - phaseVolFractionSumKernelOp( oneMinusPhaseVolFracSum ); } /** @@ -866,8 +791,9 @@ class ElementBasedAssemblyKernel shiftElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.localResidual ); // add contribution to residual and jacobian into: - // - the component mass balance equations - // - the volume balance equations + // - the component mass balance equations (i = 0 to i = numComp-1) + // - the volume balance equations (i = numComp) + // note that numDof includes derivatives wrt temperature if this class is derived in ThermalKernels for( integer i = 0; i < numComp+1; ++i ) { m_localRhs[stack.localRow + i] += stack.localResidual[i]; @@ -1005,7 +931,7 @@ class ElementBasedAssemblyKernelFactory struct ResidualNormKernel { - template< typename POLICY, typename REDUCE_POLICY > + template< typename POLICY > static void launch( arrayView1d< real64 const > const & localResidual, globalIndex const rankOffset, integer const numComponents, @@ -1016,7 +942,7 @@ struct ResidualNormKernel arrayView2d< real64 const, multifluid::USD_FLUID > const & totalDensOld, real64 & localResidualNorm ) { - RAJA::ReduceSum< REDUCE_POLICY, real64 > localSum( 0.0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > localSum( 0.0 ); forAll< POLICY >( dofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { @@ -1042,7 +968,7 @@ struct ResidualNormKernel struct SolutionCheckKernel { - template< typename POLICY, typename REDUCE_POLICY > + template< typename POLICY > static localIndex launch( arrayView1d< real64 const > const & localSolution, globalIndex const rankOffset, @@ -1058,7 +984,7 @@ struct SolutionCheckKernel { real64 constexpr eps = minDensForDivision; - RAJA::ReduceMin< REDUCE_POLICY, integer > check( 1 ); + RAJA::ReduceMin< ReducePolicy< POLICY >, integer > check( 1 ); forAll< POLICY >( dofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { @@ -1447,9 +1373,9 @@ void KernelLaunchSelector2( integer const numComp, integer const numPhase, ARGS } } -} // namespace compositionalMultiphaseBaseKernels +} // namespace isothermalCompositionalMultiphaseBaseKernels } // namespace geosx -#endif //GEOSX_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEBASEKERNELS_HPP +#endif //GEOSX_PHYSICSSOLVERS_FLUIDFLOW_ISOTHERMALCOMPOSITIONALMULTIPHASEBASEKERNELS_HPP diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp similarity index 98% rename from src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp rename to src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp index 51399c52c8b..43b0e3ca8ad 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp @@ -13,10 +13,10 @@ */ /** - * @file CompositionalMultiphaseFVMKernels.cpp + * @file IsothermalCompositionalMultiphaseFVMKernels.cpp */ -#include "CompositionalMultiphaseFVMKernels.hpp" +#include "IsothermalCompositionalMultiphaseFVMKernels.hpp" #include "CompositionalMultiphaseUtilities.hpp" #include "finiteVolume/CellElementStencilTPFA.hpp" @@ -28,7 +28,7 @@ namespace geosx { -namespace compositionalMultiphaseFVMKernels +namespace isothermalCompositionalMultiphaseFVMKernels { /******************************** FaceBasedAssemblyKernel ********************************/ @@ -152,7 +152,7 @@ CFLFluxKernel:: real64 const mobility = phaseRelPerm[er_up][esr_up][ei_up][0][ip] / phaseVisc[er_up][esr_up][ei_up][0][ip]; // increment the phase (volumetric) outflux of the upstream cell - real64 const absPhaseFlux = fabs( dt * mobility * potGrad ); + real64 const absPhaseFlux = LvArray::math::abs( dt * mobility * potGrad ); RAJA::atomicAdd( parallelDeviceAtomic{}, &phaseOutflux[er_up][esr_up][ei_up][ip], absPhaseFlux ); // increment the component (mass/molar) outflux of the upstream cell @@ -320,7 +320,7 @@ CFLKernel:: real64 const coef0 = denom * mob[ip1] / mob[ip0] * dMob_dVolFrac[ip0]; real64 const coef1 = -denom * mob[ip0] / mob[ip1] * dMob_dVolFrac[ip1]; - phaseCFLNumber = fabs( coef0*phaseOutflux[ip0] + coef1*phaseOutflux[ip1] ); + phaseCFLNumber = LvArray::math::abs( coef0*phaseOutflux[ip0] + coef1*phaseOutflux[ip1] ); } // three-phase flow regime else if( numMobilePhases == 3 ) @@ -350,7 +350,7 @@ CFLKernel:: } phaseCFLNumber = f[0][0] + f[1][1]; phaseCFLNumber += sqrt( phaseCFLNumber*phaseCFLNumber - 4 * ( f[0][0]*f[1][1] - f[1][0]*f[0][1] ) ); - phaseCFLNumber = 0.5 * fabs( phaseCFLNumber ) / poreVol; + phaseCFLNumber = 0.5 * LvArray::math::abs( phaseCFLNumber ) / poreVol; } } @@ -702,7 +702,6 @@ INST_AquiferBCKernel( 5 ); #undef INST_AquiferBCKernel - -} // namespace compositionalMultiphaseFVMKernels +} // namespace isothermalCompositionalMultiphaseFVMKernels } // namespace geosx diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp similarity index 93% rename from src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp rename to src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp index 956cdc87631..9c4c33566a8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp @@ -13,11 +13,11 @@ */ /** - * @file CompositionalMultiphaseFVMKernels.hpp + * @file IsothermalCompositionalMultiphaseFVMKernels.hpp */ -#ifndef GEOSX_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEFVMKERNELS_HPP -#define GEOSX_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEFVMKERNELS_HPP +#ifndef GEOSX_PHYSICSSOLVERS_FLUIDFLOW_ISOTHERMALCOMPOSITIONALMULTIPHASEFVMKERNELS_HPP +#define GEOSX_PHYSICSSOLVERS_FLUIDFLOW_ISOTHERMALCOMPOSITIONALMULTIPHASEFVMKERNELS_HPP #include "common/DataLayouts.hpp" #include "common/DataTypes.hpp" @@ -35,14 +35,14 @@ #include "mesh/utilities/MeshMapUtilities.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseExtrinsicData.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/StencilAccessors.hpp" namespace geosx { -namespace compositionalMultiphaseFVMKernels +namespace isothermalCompositionalMultiphaseFVMKernels { using namespace constitutive; @@ -56,11 +56,11 @@ using namespace constitutive; * @brief Define the interface for the property kernel in charge of computing the phase mobilities */ template< integer NUM_COMP, integer NUM_PHASE > -class PhaseMobilityKernel : public compositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > +class PhaseMobilityKernel : public isothermalCompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { public: - using Base = compositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; + using Base = isothermalCompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; using Base::numComp; /// Compile time value for the number of phases @@ -97,10 +97,10 @@ class PhaseMobilityKernel : public compositionalMultiphaseBaseKernels::PropertyK * @param[in] ei the element index * @param[in] phaseMobilityKernelOp the function used to customize the kernel */ - template< typename FUNC = compositionalMultiphaseBaseKernels::NoOpFunc > + template< typename FUNC = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc > GEOSX_HOST_DEVICE void compute( localIndex const ei, - FUNC && phaseMobilityKernelOp = compositionalMultiphaseBaseKernels::NoOpFunc{} ) const + FUNC && phaseMobilityKernelOp = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc{} ) const { using Deriv = multifluid::DerivativeOffset; @@ -240,7 +240,7 @@ class PhaseMobilityKernelFactory { if( numPhase == 2 ) { - compositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { integer constexpr NUM_COMP = NC(); PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); @@ -249,7 +249,7 @@ class PhaseMobilityKernelFactory } else if( numPhase == 3 ) { - compositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { integer constexpr NUM_COMP = NC(); PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); @@ -259,6 +259,7 @@ class PhaseMobilityKernelFactory } }; + /******************************** FaceBasedAssemblyKernel ********************************/ /** @@ -313,11 +314,11 @@ class FaceBasedAssemblyKernelBase * @param[in] numPhases the number of fluid phases * @param[in] rankOffset the offset of my MPI rank * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not - * @param[in] dofNumberAccessor - * @param[in] compFlowAccessors - * @param[in] multiFluidAccessors - * @param[in] capPressureAccessors - * @param[in] permeabilityAccessors + * @param[in] dofNumberAccessor accessor for the dof numbers + * @param[in] compFlowAccessors accessor for wrappers registered by the solver + * @param[in] multiFluidAccessors accessor for wrappers registered by the multifluid model + * @param[in] capPressureAccessors accessor for wrappers registered by the cap pressure model + * @param[in] permeabilityAccessors accessor for wrappers registered by the permeability model * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector @@ -573,20 +574,16 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase /** * @brief Compute the local flux contributions to the residual and Jacobian - * @tparam FUNC1 the type of the function that can be used to customize the computation of the phase fluxes - * @tparam FUNC2 the type of the function that can be used to customize the assembly into the local Jacobian + * @tparam FUNC the type of the function that can be used to customize the computation of the phase fluxes * @param[in] iconn the connection index * @param[inout] stack the stack variables - * @param[in] phaseFluxKernelOp the function used to customize the computation of the phase fluxes - * @param[in] localFluxJacobianKernelOp the function used to customize the computation of the assembly into the local Jacobian + * @param[in] compFluxKernelOp the function used to customize the computation of the component fluxes */ - template< typename FUNC1 = compositionalMultiphaseBaseKernels::NoOpFunc, - typename FUNC2 = compositionalMultiphaseBaseKernels::NoOpFunc > + template< typename FUNC = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc > GEOSX_HOST_DEVICE void computeFlux( localIndex const iconn, StackVariables & stack, - FUNC1 && phaseFluxKernelOp = compositionalMultiphaseBaseKernels::NoOpFunc{}, - FUNC2 && localFluxJacobianKernelOp = compositionalMultiphaseBaseKernels::NoOpFunc{} ) const + FUNC && compFluxKernelOp = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc{} ) const { using Deriv = multifluid::DerivativeOffset; @@ -811,8 +808,7 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase // 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 - // note: more variables may need to be passed, but it is hard to tell which ones will be needed for now - phaseFluxKernelOp( ip, er_up, esr_up, ei_up, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); + compFluxKernelOp( ip, k_up, er_up, esr_up, ei_up, potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); } @@ -821,30 +817,23 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase // populate local flux vector and derivatives for( integer ic = 0; ic < numComp; ++ic ) { - stack.localFlux[ic] = m_dt * stack.compFlux[ic]; - stack.localFlux[numComp + ic] = -m_dt * stack.compFlux[ic]; + stack.localFlux[ic] = m_dt * stack.compFlux[ic]; + stack.localFlux[numEqn + ic] = -m_dt * stack.compFlux[ic]; for( integer ke = 0; ke < stack.stencilSize; ++ke ) { localIndex const localDofIndexPres = ke * numDof; - stack.localFluxJacobian[ic][localDofIndexPres] = m_dt * stack.dCompFlux_dP[ke][ic]; - stack.localFluxJacobian[numComp + ic][localDofIndexPres] = -m_dt * stack.dCompFlux_dP[ke][ic]; + stack.localFluxJacobian[ic][localDofIndexPres] = m_dt * stack.dCompFlux_dP[ke][ic]; + stack.localFluxJacobian[numEqn + ic][localDofIndexPres] = -m_dt * stack.dCompFlux_dP[ke][ic]; for( integer jc = 0; jc < numComp; ++jc ) { localIndex const localDofIndexComp = localDofIndexPres + jc + 1; - stack.localFluxJacobian[ic][localDofIndexComp] = m_dt * stack.dCompFlux_dC[ke][ic][jc]; - stack.localFluxJacobian[numComp + ic][localDofIndexComp] = -m_dt * stack.dCompFlux_dC[ke][ic][jc]; + stack.localFluxJacobian[ic][localDofIndexComp] = m_dt * stack.dCompFlux_dC[ke][ic][jc]; + stack.localFluxJacobian[numEqn + ic][localDofIndexComp] = -m_dt * stack.dCompFlux_dC[ke][ic][jc]; } } } - - // call the lambda to allow assembly into the localFluxJacobian - // possible uses: - assemble the derivatives of compFlux wrt temperature into localFluxJacobian; - // - assemble energyFlux into localFlux and localFluxJacobian - // TODO: this second kernelOp becomes unnecessary if we increment the residual/jacobian directly in the phase loop - localFluxJacobianKernelOp( stack.localFlux, stack.localFluxJacobian ); - } /** @@ -852,20 +841,24 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase * @param[in] iconn the connection index * @param[inout] stack the stack variables */ + template< typename FUNC = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc > GEOSX_HOST_DEVICE void complete( localIndex const iconn, - StackVariables & stack ) const + StackVariables & stack, + FUNC && assemblyKernelOp = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc{} ) const { using namespace compositionalMultiphaseUtilities; // Apply equation/variable change transformation(s) stackArray1d< real64, maxStencilSize * numDof > work( stack.stencilSize * numDof ); - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numDof*stack.stencilSize, stack.numFluxElems, + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComp, numEqn, numDof*stack.stencilSize, stack.numFluxElems, stack.localFluxJacobian, work ); - shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, stack.numFluxElems, + shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComp, numEqn, stack.numFluxElems, stack.localFlux ); - // Add to residual/jacobian + // add contribution to residual and jacobian into: + // - the component mass balance equations (i = 0 to i = numComp-1) + // note that numDof includes derivatives wrt temperature if this class is derived in ThermalKernels for( integer i = 0; i < stack.numFluxElems; ++i ) { if( m_ghostRank[m_seri( iconn, i )][m_sesri( iconn, i )][m_sei( iconn, i )] < 0 ) @@ -877,13 +870,16 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase for( integer ic = 0; ic < numComp; ++ic ) { - RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[localRow + ic], stack.localFlux[i * numComp + ic] ); + RAJA::atomicAdd( parallelDeviceAtomic{}, &m_localRhs[localRow + ic], stack.localFlux[i * numEqn + ic] ); m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic > ( localRow + ic, stack.dofColIndices.data(), - stack.localFluxJacobian[i * numComp + ic].dataIfContiguous(), + stack.localFluxJacobian[i * numEqn + ic].dataIfContiguous(), stack.stencilSize * numDof ); } + + // call the lambda to assemble additional terms, such as thermal terms + assemblyKernelOp( i, localRow ); } } } @@ -913,7 +909,7 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase } ); } -private: +protected: // Stencil information @@ -963,7 +959,7 @@ class FaceBasedAssemblyKernelFactory CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { - compositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&] ( auto NC ) + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComps, [&] ( auto NC ) { integer constexpr NUM_COMP = NC(); integer constexpr NUM_DOF = NC()+1; @@ -1206,9 +1202,9 @@ struct AquiferBCKernel }; -} // namespace compositionalMultiphaseFVMKernels +} // namespace isothermalCompositionalMultiphaseFVMKernels } // namespace geosx -#endif //GEOSX_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASEFVMKERNELS_HPP +#endif //GEOSX_PHYSICSSOLVERS_FLUIDFLOW_ISOTHERMALCOMPOSITIONALMULTIPHASEFVMKERNELS_HPP diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp index 4f940401491..fe0316ad1ef 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp @@ -274,8 +274,8 @@ struct FluidUpdateKernel struct ResidualNormKernel { - template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > - static void launch( LOCAL_VECTOR const localResidual, + template< typename POLICY > + static void launch( arrayView1d< real64 const > const & localResidual, globalIndex const rankOffset, arrayView1d< globalIndex const > const & presDofNumber, arrayView1d< integer const > const & ghostRank, @@ -284,9 +284,9 @@ struct ResidualNormKernel arrayView2d< real64 const > const & poroOld, real64 * localResidualNorm ) { - RAJA::ReduceSum< REDUCE_POLICY, real64 > localSum( 0.0 ); - RAJA::ReduceSum< REDUCE_POLICY, real64 > normSum( 0.0 ); - RAJA::ReduceSum< REDUCE_POLICY, localIndex > count( 0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > localSum( 0.0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > normSum( 0.0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, localIndex > count( 0 ); forAll< POLICY >( presDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const a ) { @@ -310,8 +310,8 @@ struct ResidualNormKernel struct SolutionCheckKernel { - template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > - static localIndex launch( LOCAL_VECTOR const localSolution, + template< typename POLICY > + static localIndex launch( arrayView1d< real64 const > const & localSolution, globalIndex const rankOffset, arrayView1d< globalIndex const > const & presDofNumber, arrayView1d< integer const > const & ghostRank, @@ -319,7 +319,7 @@ struct SolutionCheckKernel arrayView1d< real64 const > const & dPres, real64 const scalingFactor ) { - RAJA::ReduceMin< REDUCE_POLICY, localIndex > minVal( 1 ); + RAJA::ReduceMin< ReducePolicy< POLICY >, localIndex > minVal( 1 ); forAll< POLICY >( presDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index d088e1122f5..bc079617f06 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -137,14 +137,14 @@ real64 SinglePhaseFVM< BASE >::calculateResidualNorm( DomainPartition const & do arrayView2d< real64 const > const & porosityOld = solidModel.getOldPorosity(); - ResidualNormKernel::launch< parallelDevicePolicy<>, parallelDeviceReduce >( localRhs, - rankOffset, - dofNumber, - elemGhostRank, - volume, - densOld, - porosityOld, - localResidualNorm ); + ResidualNormKernel::launch< parallelDevicePolicy<> >( localRhs, + rankOffset, + dofNumber, + elemGhostRank, + volume, + densOld, + porosityOld, + localResidualNorm ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp index e27d1a8a3e2..32b33381d26 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp @@ -478,15 +478,15 @@ real64 SinglePhaseHybridFVM::calculateResidualNorm( DomainPartition const & doma { arrayView2d< real64 const > const & porosityOld = castedSolidModel.getOldPorosity(); - singlePhaseBaseKernels::ResidualNormKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localRhs, - rankOffset, - elemDofNumber, - elemGhostRank, - volume, - densOld, - porosityOld, - localResidualNorm ); + singlePhaseBaseKernels:: + ResidualNormKernel::launch< parallelDevicePolicy<> >( localRhs, + rankOffset, + elemDofNumber, + elemGhostRank, + volume, + densOld, + porosityOld, + localResidualNorm ); } ); string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); @@ -506,17 +506,17 @@ real64 SinglePhaseHybridFVM::calculateResidualNorm( DomainPartition const & doma defaultViscosity /= subRegionCounter; // 2. Compute the residual for the face-based constraints - singlePhaseHybridFVMKernels::ResidualNormKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localRhs, - rankOffset, - faceDofNumber.toNestedViewConst(), - faceGhostRank.toNestedViewConst(), - elemRegionList.toNestedViewConst(), - elemSubRegionList.toNestedViewConst(), - elemList.toNestedViewConst(), - flowAccessors.get( extrinsicMeshData::elementVolume{} ), - defaultViscosity, - &localResidualNorm[3] ); + singlePhaseHybridFVMKernels:: + ResidualNormKernel::launch< parallelDevicePolicy<> >( localRhs, + rankOffset, + faceDofNumber.toNestedViewConst(), + faceGhostRank.toNestedViewConst(), + elemRegionList.toNestedViewConst(), + elemSubRegionList.toNestedViewConst(), + elemList.toNestedViewConst(), + flowAccessors.get( extrinsicMeshData::elementVolume{} ), + defaultViscosity, + &localResidualNorm[3] ); } ); @@ -572,14 +572,14 @@ bool SinglePhaseHybridFVM::checkSystemSolution( DomainPartition const & domain, subRegion.getExtrinsicData< extrinsicMeshData::flow::deltaPressure >(); localIndex const subRegionSolutionCheck = - singlePhaseBaseKernels::SolutionCheckKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localSolution, - rankOffset, - elemDofNumber, - elemGhostRank, - pres, - dPres, - scalingFactor ); + singlePhaseBaseKernels:: + SolutionCheckKernel::launch< parallelDevicePolicy<> >( localSolution, + rankOffset, + elemDofNumber, + elemGhostRank, + pres, + dPres, + scalingFactor ); if( subRegionSolutionCheck == 0 ) { @@ -599,14 +599,14 @@ bool SinglePhaseHybridFVM::checkSystemSolution( DomainPartition const & domain, localIndex const faceSolutionCheck = - singlePhaseBaseKernels::SolutionCheckKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localSolution, - rankOffset, - faceDofNumber, - faceGhostRank, - facePres, - dFacePres, - scalingFactor ); + singlePhaseBaseKernels:: + SolutionCheckKernel::launch< parallelDevicePolicy<> >( localSolution, + rankOffset, + faceDofNumber, + faceGhostRank, + facePres, + dFacePres, + scalingFactor ); if( faceSolutionCheck == 0 ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp index 56e7108b826..02b7159afe4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp @@ -638,9 +638,9 @@ struct ResidualNormKernel template< typename VIEWTYPE > using ElementViewConst = typename ElementRegionManager::ElementViewConst< VIEWTYPE >; - template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > + template< typename POLICY > static void - launch( LOCAL_VECTOR const localResidual, + launch( arrayView1d< real64 const > const & localResidual, globalIndex const rankOffset, arrayView1d< globalIndex const > const & facePresDofNumber, arrayView1d< integer const > const & faceGhostRank, @@ -652,7 +652,7 @@ struct ResidualNormKernel real64 * localResidualNorm ) { - RAJA::ReduceSum< REDUCE_POLICY, real64 > sumScaled( 0.0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > sumScaled( 0.0 ); forAll< POLICY >( facePresDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iface ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp new file mode 100644 index 00000000000..83c2021756e --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp @@ -0,0 +1,668 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ThermalCompositionalMultiphaseBaseKernels.hpp + */ + +#ifndef GEOSX_PHYSICSSOLVERS_FLUIDFLOW_THERMALCOMPOSITIONALMULTIPHASEBASEKERNELS_HPP +#define GEOSX_PHYSICSSOLVERS_FLUIDFLOW_THERMALCOMPOSITIONALMULTIPHASEBASEKERNELS_HPP + +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" + +namespace geosx +{ + +namespace thermalCompositionalMultiphaseBaseKernels +{ + +using namespace constitutive; + + +/******************************** PhaseVolumeFractionKernel ********************************/ + +/** + * @class PhaseVolumeFractionKernel + * @tparam NUM_COMP number of fluid components + * @tparam NUM_PHASE number of fluid phases + * @brief Define the interface for the property kernel in charge of computing the phase volume fractions + */ +template< integer NUM_COMP, integer NUM_PHASE > +class PhaseVolumeFractionKernel : public isothermalCompositionalMultiphaseBaseKernels::PhaseVolumeFractionKernel< NUM_COMP, NUM_PHASE > +{ +public: + + using Base = isothermalCompositionalMultiphaseBaseKernels::PhaseVolumeFractionKernel< NUM_COMP, NUM_PHASE >; + using Base::m_dPhaseDens; + using Base::m_dPhaseFrac; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + PhaseVolumeFractionKernel( ObjectManagerBase & subRegion, + MultiFluidBase const & fluid ) + : Base( subRegion, fluid ), + m_dPhaseVolFrac_dTemp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dTemperature >() ) + {} + + /** + * @brief Compute the phase volume fractions in an element + * @param[in] ei the element index + */ + GEOSX_HOST_DEVICE + void compute( localIndex const ei ) const + { + using Deriv = multifluid::DerivativeOffset; + + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseDens = m_dPhaseDens[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseFrac = m_dPhaseFrac[ei][0]; + + arraySlice1d< real64, compflow::USD_PHASE - 1 > const dPhaseVolFrac_dTemp = m_dPhaseVolFrac_dTemp[ei]; + LvArray::forValuesInSlice( dPhaseVolFrac_dTemp, []( real64 & val ){ val = 0.0; } ); + + // Call the base compute the compute the phase volume fraction + Base::compute( ei, [&] ( localIndex const ip, + real64 const & phaseVolFrac, + real64 const & phaseDensInv, + real64 const & totalDensity ) + { + // when this lambda is called, we are in the phase loop + // for each phase ip, compute the derivative of phase volume fraction wrt temperature + dPhaseVolFrac_dTemp[ip] = (dPhaseFrac[ip][Deriv::dT] - phaseVolFrac * dPhaseDens[ip][Deriv::dT]) * phaseDensInv; + dPhaseVolFrac_dTemp[ip] *= totalDensity; + } ); + } + +protected: + + // outputs + + /// Views on thermal derivatives of phase volume fractions + arrayView2d< real64, compflow::USD_PHASE > m_dPhaseVolFrac_dTemp; + +}; + +/** + * @class PhaseVolumeFractionKernelFactory + */ +class PhaseVolumeFractionKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numComp the number of fluid components + * @param[in] numPhase the number of fluid phases + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + */ + template< typename POLICY > + static void + createAndLaunch( integer const numComp, + integer const numPhase, + ObjectManagerBase & subRegion, + MultiFluidBase const & fluid ) + { + if( numPhase == 2 ) + { + isothermalCompositionalMultiphaseBaseKernels:: + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseVolumeFractionKernel< NUM_COMP, 2 > kernel( subRegion, fluid ); + PhaseVolumeFractionKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + isothermalCompositionalMultiphaseBaseKernels:: + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseVolumeFractionKernel< NUM_COMP, 3 > kernel( subRegion, fluid ); + PhaseVolumeFractionKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + } +}; + + +/******************************** ElementBasedAssemblyKernel ********************************/ + +/** + * @class ElementBasedAssemblyKernel + * @tparam NUM_COMP number of fluid components + * @tparam NUM_DOF number of degrees of freedom + * @brief Define the interface for the assembly kernel in charge of thermal accumulation and volume balance + */ +template< localIndex NUM_COMP, localIndex NUM_DOF > +class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseKernels::ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > +{ +public: + + using Base = isothermalCompositionalMultiphaseBaseKernels::ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF >; + using Base::numComp; + using Base::numDof; + using Base::numEqn; + using Base::m_numPhases; + using Base::m_rankOffset; + using Base::m_dofNumber; + using Base::m_elemGhostRank; + using Base::m_volume; + using Base::m_porosityOld; + using Base::m_porosityNew; + using Base::m_dPoro_dPres; + using Base::m_dCompFrac_dCompDens; + using Base::m_phaseVolFracOld; + using Base::m_phaseVolFrac; + using Base::m_dPhaseVolFrac_dPres; + using Base::m_dPhaseVolFrac_dCompDens; + using Base::m_phaseDensOld; + using Base::m_phaseDens; + using Base::m_dPhaseDens; + using Base::m_phaseCompFracOld; + using Base::m_phaseCompFrac; + using Base::m_dPhaseCompFrac; + using Base::m_localMatrix; + using Base::m_localRhs; + + /** + * @brief Constructor + * @param[in] numPhases the number of fluid phases + * @param[in] rankOffset the offset of my MPI rank + * @param[in] dofKey the string key to retrieve the degress of freedom numbers + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] solid the solid model + * @param[inout] localMatrix the local CRS matrix + * @param[inout] localRhs the local right-hand side vector + */ + ElementBasedAssemblyKernel( localIndex const numPhases, + globalIndex const rankOffset, + string const dofKey, + ElementSubRegionBase const & subRegion, + MultiFluidBase const & fluid, + CoupledSolidBase const & solid, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + : Base( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ), + m_dPhaseVolFrac_dTemp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dTemperature >() ), + m_phaseInternalEnergyOld( fluid.phaseInternalEnergyOld() ), + m_phaseInternalEnergy( fluid.phaseInternalEnergy() ), + m_dPhaseInternalEnergy( fluid.dPhaseInternalEnergy() ), + m_rockInternalEnergyOld( solid.getOldInternalEnergy() ), + m_rockInternalEnergy( solid.getInternalEnergy() ), + m_dRockInternalEnergy_dTemp( solid.getDinternalEnergy_dTemperature() ) + {} + + struct StackVariables : public Base::StackVariables + { +public: + + GEOSX_HOST_DEVICE + StackVariables() + : Base::StackVariables() + {} + + using Base::StackVariables::poreVolumeNew; + using Base::StackVariables::poreVolumeOld; + using Base::StackVariables::dPoreVolume_dPres; + using Base::StackVariables::localRow; + using Base::StackVariables::dofIndices; + using Base::StackVariables::localResidual; + using Base::StackVariables::localJacobian; + + // Solid energy + + /// Solid energy at time n+1 + real64 solidEnergyNew = 0.0; + + /// Solid energy at the previous converged time step + real64 solidEnergyOld = 0.0; + + /// Derivative of solid internal energy with respect to pressure + real64 dSolidEnergy_dPres = 0.0; + + /// Derivative of solid internal energy with respect to temperature + real64 dSolidEnergy_dTemp = 0.0; + + }; + + /** + * @brief Performs the setup phase for the kernel. + * @param[in] ei the element index + * @param[in] stack the stack variables + */ + GEOSX_HOST_DEVICE + void setup( localIndex const ei, + StackVariables & stack ) const + { + Base::setup( ei, stack ); + + // initialize the solid volume + real64 const solidVolumeNew = m_volume[ei] * ( 1.0 - m_porosityNew[ei][0] ); + real64 const solidVolumeOld = m_volume[ei] * ( 1.0 - m_porosityOld[ei][0] ); + real64 const dSolidVolume_dPres = -m_volume[ei] * m_dPoro_dPres[ei][0]; + + // initialize the solid internal energy + stack.solidEnergyNew = solidVolumeNew * m_rockInternalEnergy[ei][0]; + stack.solidEnergyOld = solidVolumeOld * m_rockInternalEnergyOld[ei][0]; + stack.dSolidEnergy_dPres = dSolidVolume_dPres * m_rockInternalEnergy[ei][0]; + stack.dSolidEnergy_dTemp = solidVolumeNew * m_dRockInternalEnergy_dTemp[ei][0]; + } + + /** + * @brief Compute the local accumulation contributions to the residual and Jacobian + * @tparam FUNC the type of the function that can be used to customize the kernel + * @param[in] ei the element index + * @param[inout] stack the stack variables + */ + GEOSX_HOST_DEVICE + void computeAccumulation( localIndex const ei, + StackVariables & stack ) const + { + using Deriv = multifluid::DerivativeOffset; + + Base::computeAccumulation( ei, stack, [&] ( integer const ip, + real64 const & phaseAmountNew, + real64 const & phaseAmountOld, + real64 const & dPhaseAmount_dP, + real64 const (&dPhaseAmount_dC)[numComp] ) + { + // We are in the loop over phases, ip provides the current phase index. + // We have to do two things: + // 1- Assemble the derivatives of the component mass balance equations with respect to temperature + // 2- Assemble the phase-dependent part of the accumulation term of the energy equation + + real64 dPhaseInternalEnergy_dC[numComp]{}; + + // construct the slices + arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei]; + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac = m_phaseVolFrac[ei]; + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > dPhaseVolFrac_dTemp = m_dPhaseVolFrac_dTemp[ei]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseDens = m_phaseDens[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseDens = m_dPhaseDens[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_COMP - 2 > phaseCompFrac = m_phaseCompFrac[ei][0]; + arraySlice3d< real64 const, multifluid::USD_PHASE_COMP_DC - 2 > dPhaseCompFrac = m_dPhaseCompFrac[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseInternalEnergyOld = m_phaseInternalEnergyOld[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseInternalEnergy = m_phaseInternalEnergy[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseInternalEnergy = m_dPhaseInternalEnergy[ei][0]; + + // Step 1: assemble the derivatives of the component mass balance equations with respect to temperature + + real64 const dPhaseAmount_dT = stack.poreVolumeNew + * (dPhaseVolFrac_dTemp[ip] * phaseDens[ip] + phaseVolFrac[ip] * dPhaseDens[ip][Deriv::dT] ); + for( integer ic = 0; ic < numComp; ++ic ) + { + stack.localJacobian[ic][numDof-1] += dPhaseAmount_dT * phaseCompFrac[ip][ic] + + phaseAmountNew * dPhaseCompFrac[ip][ic][Deriv::dT]; + } + + // Step 2: assemble the phase-dependent part of the accumulation term of the energy equation + + real64 const phaseEnergyNew = phaseAmountNew * phaseInternalEnergy[ip]; + real64 const phaseEnergyOld = phaseAmountOld * phaseInternalEnergyOld[ip]; + real64 const dPhaseEnergy_dP = dPhaseAmount_dP * phaseInternalEnergy[ip] + + phaseAmountNew * dPhaseInternalEnergy[ip][Deriv::dP]; + real64 const dPhaseEnergy_dT = dPhaseAmount_dT * phaseInternalEnergy[ip] + + phaseAmountNew * dPhaseInternalEnergy[ip][Deriv::dT]; + + // local accumulation + stack.localResidual[numEqn-1] += phaseEnergyNew - phaseEnergyOld; + + // derivatives w.r.t. pressure and temperature + stack.localJacobian[numEqn-1][0] += dPhaseEnergy_dP; + stack.localJacobian[numEqn-1][numDof-1] += dPhaseEnergy_dT; + + // derivatives w.r.t. component densities + applyChainRule( numComp, dCompFrac_dCompDens, dPhaseInternalEnergy[ip], dPhaseInternalEnergy_dC, Deriv::dC ); + for( integer jc = 0; jc < numComp; ++jc ) + { + stack.localJacobian[numEqn-1][jc + 1] += phaseInternalEnergy[ip] * dPhaseAmount_dC[jc] + + dPhaseInternalEnergy_dC[jc] * phaseAmountNew; + } + } ); + + // Step 3: assemble the solid part of the accumulation term + + // local accumulation and derivatives w.r.t. pressure and temperature + stack.localResidual[numEqn-1] += stack.solidEnergyNew - stack.solidEnergyOld; + stack.localJacobian[numEqn-1][0] += stack.dSolidEnergy_dPres; + stack.localJacobian[numEqn-1][numDof-1] += stack.dSolidEnergy_dTemp; + + } + + /** + * @brief Compute the local volume balance contributions to the residual and Jacobian + * @tparam FUNC the type of the function that can be used to customize the kernel + * @param[in] ei the element index + * @param[inout] stack the stack variables + */ + GEOSX_HOST_DEVICE + void computeVolumeBalance( localIndex const ei, + StackVariables & stack ) const + { + Base::computeVolumeBalance( ei, stack, [&] ( real64 const & oneMinusPhaseVolFraction ) + { + GEOSX_UNUSED_VAR( oneMinusPhaseVolFraction ); + + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > dPhaseVolFrac_dTemp = m_dPhaseVolFrac_dTemp[ei]; + + for( integer ip = 0; ip < m_numPhases; ++ip ) + { + stack.localJacobian[numEqn-2][numDof-1] -= dPhaseVolFrac_dTemp[ip]; + } + } ); + } + + GEOSX_HOST_DEVICE + void complete( localIndex const ei, + StackVariables & stack ) const + { + // Step 1: assemble the component mass balance equations and volume balance equations + Base::complete( ei, stack ); + + // Step 2: assemble the energy equation + m_localRhs[stack.localRow + numEqn-1] += stack.localResidual[numEqn-1]; + m_localMatrix.template addToRow< serialAtomic >( stack.localRow + numEqn-1, + stack.dofIndices, + stack.localJacobian[numEqn-1], + numDof ); + } + +protected: + + /// Views on derivatives wrt to temperature for phase volume fraction + arrayView2d< real64 const, compflow::USD_PHASE > m_dPhaseVolFrac_dTemp; + + /// Views on phase internal energy + arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseInternalEnergyOld; + arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseInternalEnergy; + arrayView4d< real64 const, multifluid::USD_PHASE_DC > m_dPhaseInternalEnergy; + + /// Views on rock internal energy + arrayView2d< real64 const > m_rockInternalEnergyOld; + arrayView2d< real64 const > m_rockInternalEnergy; + arrayView2d< real64 const > m_dRockInternalEnergy_dTemp; + +}; + +/** + * @class ElementBasedAssemblyKernelFactory + */ +class ElementBasedAssemblyKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numComps the number of fluid components + * @param[in] numPhases the number of fluid phases + * @param[in] rankOffset the offset of my MPI rank + * @param[in] dofKey the string key to retrieve the degress of freedom numbers + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] solid the solid model + * @param[inout] localMatrix the local CRS matrix + * @param[inout] localRhs the local right-hand side vector + */ + template< typename POLICY > + static void + createAndLaunch( localIndex const numComps, + localIndex const numPhases, + globalIndex const rankOffset, + string const dofKey, + ElementSubRegionBase const & subRegion, + MultiFluidBase const & fluid, + CoupledSolidBase const & solid, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + isothermalCompositionalMultiphaseBaseKernels:: + internal::kernelLaunchSelectorCompSwitch( numComps, [&] ( auto NC ) + { + localIndex constexpr NUM_COMP = NC(); + localIndex constexpr NUM_DOF = NC()+2; + ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > + kernel( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ); + ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF >::template + launch< POLICY, ElementBasedAssemblyKernel< NUM_COMP, NUM_DOF > >( subRegion.size(), kernel ); + } ); + } + +}; + +/******************************** FluidUpdateKernel ********************************/ + +struct FluidUpdateKernel +{ + template< typename POLICY, typename FLUID_WRAPPER > + static void + launch( localIndex const size, + FLUID_WRAPPER const & fluidWrapper, + arrayView1d< real64 const > const & pres, + arrayView1d< real64 const > const & temp, + arrayView2d< real64 const, compflow::USD_COMP > const & compFrac ) + { + forAll< POLICY >( size, [=] GEOSX_HOST_DEVICE ( localIndex const k ) + { + for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q ) + { + fluidWrapper.update( k, q, pres[k], temp[k], compFrac[k] ); + } + } ); + } + + template< typename POLICY, typename FLUID_WRAPPER > + static void + launch( localIndex const size, + FLUID_WRAPPER const & fluidWrapper, + arrayView1d< real64 const > const & pres, + arrayView1d< real64 const > const & dPres, + arrayView1d< real64 const > const & temp, + arrayView1d< real64 const > const & dTemp, + arrayView2d< real64 const, compflow::USD_COMP > const & compFrac ) + { + forAll< POLICY >( size, [=] GEOSX_HOST_DEVICE ( localIndex const k ) + { + for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q ) + { + fluidWrapper.update( k, q, pres[k] + dPres[k], temp[k] + dTemp[k], compFrac[k] ); + } + } ); + } + + template< typename POLICY, typename FLUID_WRAPPER > + static void + launch( SortedArrayView< localIndex const > const & targetSet, + FLUID_WRAPPER const & fluidWrapper, + arrayView1d< real64 const > const & pres, + arrayView1d< real64 const > const & dPres, + arrayView1d< real64 const > const & temp, + arrayView1d< real64 const > const & dTemp, + arrayView2d< real64 const, compflow::USD_COMP > const & compFrac ) + { + forAll< POLICY >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const a ) + { + localIndex const k = targetSet[a]; + for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q ) + { + fluidWrapper.update( k, q, pres[k] + dPres[k], temp[k] + dTemp[k], compFrac[k] ); + } + } ); + } + + template< typename POLICY, typename FLUID_WRAPPER > + static void + launch( SortedArrayView< localIndex const > const & targetSet, + FLUID_WRAPPER const & fluidWrapper, + arrayView1d< real64 const > const & pres, + arrayView1d< real64 const > const & temp, + arrayView2d< real64 const, compflow::USD_COMP > const & compFrac ) + { + forAll< POLICY >( targetSet.size(), [=] GEOSX_HOST_DEVICE ( localIndex const a ) + { + localIndex const k = targetSet[a]; + for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q ) + { + fluidWrapper.update( k, q, pres[k], temp[k], compFrac[k] ); + } + } ); + } +}; + +/******************************** SolidInternalEnergyUpdateKernel ********************************/ + +struct SolidInternalEnergyUpdateKernel +{ + + template< typename POLICY, typename SOLID_INTERNAL_ENERGY_WRAPPER > + static void + launch( localIndex const size, + SOLID_INTERNAL_ENERGY_WRAPPER const & solidInternalEnergyWrapper, + arrayView1d< real64 const > const & temp, + arrayView1d< real64 const > const & dTemp ) + { + forAll< POLICY >( size, [=] GEOSX_HOST_DEVICE ( localIndex const k ) + { + solidInternalEnergyWrapper.update( k, temp[k] + dTemp[k] ); + } ); + } +}; + +/******************************** ResidualNormKernel ********************************/ + +struct ResidualNormKernel +{ + + template< typename POLICY > + static void launch( arrayView1d< real64 const > const & localResidual, + globalIndex const rankOffset, + integer const numComponents, + integer const numPhases, + arrayView1d< globalIndex const > const & dofNumber, + arrayView1d< integer const > const & ghostRank, + arrayView1d< real64 const > const & refPoro, + arrayView1d< real64 const > const & volume, + arrayView2d< real64 const > const & solidInternalEnergyOld, + arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFracOld, + arrayView2d< real64 const, multifluid::USD_FLUID > const & totalDensOld, + arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDensOld, + arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseInternalEnergyOld, + real64 & localFlowResidualNorm, + real64 & localEnergyResidualNorm ) + { + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > localFlowSum( 0.0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > localEnergySum( 0.0 ); + + forAll< POLICY >( dofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) + { + if( ghostRank[ei] < 0 ) + { + localIndex const localRow = dofNumber[ei] - rankOffset; + + // first, compute the normalizers for the component mass balance and energy balance equations + // TODO: apply a separate treatment to the volume balance equation + real64 const poreVolume = refPoro[ei] * volume[ei]; + real64 const flowNormalizer = totalDensOld[ei][0] * poreVolume; + real64 energyNormalizer = solidInternalEnergyOld[ei][0] * ( 1.0 - refPoro[ei] ) * volume[ei]; + for( integer ip = 0; ip < numPhases; ++ip ) + { + energyNormalizer += phaseInternalEnergyOld[ei][0][ip] * phaseDensOld[ei][0][ip] * phaseVolFracOld[ei][ip] * poreVolume; + } + + // then, compute the normalized residual + for( integer idof = 0; idof < numComponents + 1; ++idof ) + { + real64 const val = localResidual[localRow + idof] / flowNormalizer; + localFlowSum += val * val; + } + real64 const val = localResidual[localRow + numComponents + 1] / energyNormalizer; + localEnergySum += val * val; + } + } ); + localFlowResidualNorm += localFlowSum.get(); + localEnergyResidualNorm += localEnergySum.get(); + } + +}; + + +/******************************** Kernel launch machinery ********************************/ + +// TODO: remove, move, avoid duplication + +namespace internal +{ + +template< typename T, typename LAMBDA > +void KernelLaunchSelectorCompSwitch( T value, LAMBDA && lambda ) +{ + static_assert( std::is_integral< T >::value, "KernelLaunchSelectorCompSwitch: type should be integral" ); + + switch( value ) + { + case 1: + { lambda( std::integral_constant< T, 1 >() ); return; } + case 2: + { lambda( std::integral_constant< T, 2 >() ); return; } + case 3: + { lambda( std::integral_constant< T, 3 >() ); return; } + case 4: + { lambda( std::integral_constant< T, 4 >() ); return; } + case 5: + { lambda( std::integral_constant< T, 5 >() ); return; } + default: + { GEOSX_ERROR( "Unsupported number of components: " << value ); } + } +} + +} // namespace helpers + +template< typename KERNELWRAPPER, typename ... ARGS > +void KernelLaunchSelector1( localIndex const numComp, ARGS && ... args ) +{ + internal::KernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + KERNELWRAPPER::template launch< NC() >( std::forward< ARGS >( args )... ); + } ); +} + +template< typename KERNELWRAPPER, typename ... ARGS > +void KernelLaunchSelector2( localIndex const numComp, localIndex const numPhase, ARGS && ... args ) +{ + internal::KernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + switch( numPhase ) + { + case 2: + { KERNELWRAPPER::template launch< NC(), 2 >( std::forward< ARGS >( args )... ); return; } + case 3: + { KERNELWRAPPER::template launch< NC(), 3 >( std::forward< ARGS >( args )... ); return; } + default: + { GEOSX_ERROR( "Unsupported number of phases: " << numPhase ); } + } + } ); +} + +} // namespace thermalCompositionalMultiphaseBaseKernels + +} // namespace geosx + + +#endif //GEOSX_PHYSICSSOLVERS_FLUIDFLOW_THERMALCOMPOSITIONALMULTIPHASEBASEKERNELS_HPP diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp new file mode 100644 index 00000000000..0006f5919fa --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseFVMKernels.hpp @@ -0,0 +1,673 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ThermalCompositionalMultiphaseFVMKernels.hpp + */ + +#ifndef GEOSX_PHYSICSSOLVERS_FLUIDFLOW_THERMALCOMPOSITIONALMULTIPHASEFVMKERNELS_HPP +#define GEOSX_PHYSICSSOLVERS_FLUIDFLOW_THERMALCOMPOSITIONALMULTIPHASEFVMKERNELS_HPP + +#include "constitutive/thermalConductivity/ThermalConductivityBase.hpp" +#include "constitutive/thermalConductivity/ThermalConductivityExtrinsicData.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp" + +namespace geosx +{ + +namespace thermalCompositionalMultiphaseFVMKernels +{ + +using namespace constitutive; + +/******************************** PhaseMobilityKernel ********************************/ + +/** + * @class PhaseMobilityKernel + * @tparam NUM_COMP number of fluid components + * @tparam NUM_PHASE number of fluid phases + * @brief Define the interface for the property kernel in charge of computing the phase mobilities + */ +template< integer NUM_COMP, integer NUM_PHASE > +class PhaseMobilityKernel : public isothermalCompositionalMultiphaseFVMKernels::PhaseMobilityKernel< NUM_COMP, NUM_PHASE > +{ +public: + + using Base = isothermalCompositionalMultiphaseFVMKernels::PhaseMobilityKernel< NUM_COMP, NUM_PHASE >; + using Base::numPhase; + using Base::m_phaseDens; + using Base::m_dPhaseDens; + using Base::m_phaseVisc; + using Base::m_dPhaseVisc; + using Base::m_dPhaseRelPerm_dPhaseVolFrac; + + /** + * @brief Constructor + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] relperm the relperm model + */ + PhaseMobilityKernel( ObjectManagerBase & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + : Base( subRegion, fluid, relperm ), + m_dPhaseVolFrac_dTemp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseVolumeFraction_dTemperature >() ), + m_dPhaseMob_dTemp( subRegion.getExtrinsicData< extrinsicMeshData::flow::dPhaseMobility_dTemperature >() ) + {} + + /** + * @brief Compute the phase mobilities in an element + * @param[in] ei the element index + */ + GEOSX_HOST_DEVICE + void compute( localIndex const ei ) const + { + using Deriv = multifluid::DerivativeOffset; + + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const phaseDens = m_phaseDens[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseDens = m_dPhaseDens[ei][0]; + arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > const phaseVisc = m_phaseVisc[ei][0]; + arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > const dPhaseVisc = m_dPhaseVisc[ei][0]; + arraySlice2d< real64 const, relperm::USD_RELPERM_DS - 2 > const dPhaseRelPerm_dPhaseVolFrac = m_dPhaseRelPerm_dPhaseVolFrac[ei][0]; + arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const dPhaseVolFrac_dTemp = m_dPhaseVolFrac_dTemp[ei]; + + arraySlice1d< real64, compflow::USD_PHASE - 1 > const dPhaseMob_dTemp = m_dPhaseMob_dTemp[ei]; + LvArray::forValuesInSlice( dPhaseMob_dTemp, []( real64 & val ){ val = 0.0; } ); + + Base::compute( ei, [&] ( localIndex const ip, + real64 const & phaseMob, + real64 const & GEOSX_UNUSED_PARAM( dPhaseMob_dPres ), + arraySlice1d< real64 const, compflow::USD_PHASE_DC - 2 > const & GEOSX_UNUSED_PARAM( dPhaseMob_dComp ) ) + { + // Step 1: compute the derivative of relPerm[ip] wrt temperature + real64 dRelPerm_dT = 0.0; + for( integer jp = 0; jp < numPhase; ++jp ) + { + dRelPerm_dT += dPhaseRelPerm_dPhaseVolFrac[ip][jp] * dPhaseVolFrac_dTemp[jp]; + } + + // Step 2: compute the derivative of phaseMob[ip] wrt temperature + dPhaseMob_dTemp[ip] = dRelPerm_dT * phaseDens[ip] / phaseVisc[ip] + + phaseMob * (dPhaseDens[ip][Deriv::dT] / phaseDens[ip] - dPhaseVisc[ip][Deriv::dT] / phaseVisc[ip] ); + } ); + } + +protected: + + // inputs + + /// Views on thermal derivatives of phase volume fractions + arrayView2d< real64 const, compflow::USD_PHASE > m_dPhaseVolFrac_dTemp; + + // outputs + + /// Views on thermal derivatives of phase mobilities + arrayView2d< real64, compflow::USD_PHASE > m_dPhaseMob_dTemp; + +}; + +/** + * @class PhaseMobilityKernelFactory + */ +class PhaseMobilityKernelFactory +{ +public: + + /** + * @brief Create a new kernel and launch + * @tparam POLICY the policy used in the RAJA kernel + * @param[in] numComp the number of fluid components + * @param[in] numPhase the number of fluid phases + * @param[in] subRegion the element subregion + * @param[in] fluid the fluid model + * @param[in] relperm the relperm model + */ + template< typename POLICY > + static void + createAndLaunch( integer const numComp, + integer const numPhase, + ObjectManagerBase & subRegion, + MultiFluidBase const & fluid, + RelativePermeabilityBase const & relperm ) + { + if( numPhase == 2 ) + { + isothermalCompositionalMultiphaseBaseKernels:: + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 2 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 2 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + else if( numPhase == 3 ) + { + isothermalCompositionalMultiphaseBaseKernels:: + internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + PhaseMobilityKernel< NUM_COMP, 3 > kernel( subRegion, fluid, relperm ); + PhaseMobilityKernel< NUM_COMP, 3 >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } + } +}; + + +/******************************** FaceBasedAssemblyKernel ********************************/ + +/** + * @class FaceBasedAssemblyKernel + * @tparam NUM_COMP number of fluid components + * @tparam NUM_DOF number of degrees of freedom + * @tparam STENCILWRAPPER the type of the stencil wrapper + * @brief Define the interface for the assembly kernel in charge of flux terms + */ +template< integer NUM_COMP, integer NUM_DOF, typename STENCILWRAPPER > +class FaceBasedAssemblyKernel : public isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER > +{ +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::FaceBasedAssemblyKernelBase; + using DofNumberAccessor = AbstractBase::DofNumberAccessor; + using CompFlowAccessors = AbstractBase::CompFlowAccessors; + using MultiFluidAccessors = AbstractBase::MultiFluidAccessors; + using CapPressureAccessors = AbstractBase::CapPressureAccessors; + using PermeabilityAccessors = AbstractBase::PermeabilityAccessors; + + using AbstractBase::m_dt; + using AbstractBase::m_numPhases; + using AbstractBase::m_hasCapPressure; + using AbstractBase::m_rankOffset; + using AbstractBase::m_dofNumber; + using AbstractBase::m_gravCoef; + using AbstractBase::m_phaseMob; + using AbstractBase::m_dPhaseMassDens; + using AbstractBase::m_phaseCompFrac; + using AbstractBase::m_dPhaseCompFrac; + using AbstractBase::m_dCompFrac_dCompDens; + using AbstractBase::m_dPhaseCapPressure_dPhaseVolFrac; + + using Base = isothermalCompositionalMultiphaseFVMKernels::FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; + using Base::numComp; + 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 ThermalCompFlowAccessors = + StencilAccessors< extrinsicMeshData::flow::temperature, + extrinsicMeshData::flow::deltaTemperature, + extrinsicMeshData::flow::dPhaseMobility_dTemperature, + extrinsicMeshData::flow::dPhaseVolumeFraction_dTemperature >; + + using ThermalMultiFluidAccessors = + StencilMaterialAccessors< MultiFluidBase, + extrinsicMeshData::multifluid::phaseEnthalpy, + extrinsicMeshData::multifluid::dPhaseEnthalpy >; + + using ThermalConductivityAccessors = + StencilMaterialAccessors< ThermalConductivityBase, + extrinsicMeshData::thermalconductivity::effectiveConductivity >; + // for now, we treat thermal conductivity explicitly + + /** + * @brief Constructor for the kernel interface + * @param[in] numPhases the number of fluid phases + * @param[in] rankOffset the offset of my MPI rank + * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not + * @param[in] stencilWrapper reference to the stencil wrapper + * @param[in] dofNumberAccessor accessor for the dofs numbers + * @param[in] compFlowAccessor accessor for wrappers registered by the solver + * @param[in] thermalCompFlowAccessors accessor for *thermal* wrappers registered by the solver + * @param[in] multiFluidAccessor accessor for wrappers registered by the multifluid model + * @param[in] thermalMultiFluidAccessors accessor for *thermal* wrappers registered by the multifluid model + * @param[in] capPressureAccessors accessor for wrappers registered by the cap pressure model + * @param[in] permeabilityAccessors accessor for wrappers registered by the permeability model + * @param[in] thermalConductivityAccessors accessor for wrappers registered by the thermal conductivity model + * @param[in] dt time step size + * @param[inout] localMatrix the local CRS matrix + * @param[inout] localRhs the local right-hand side vector + */ + FaceBasedAssemblyKernel( integer const numPhases, + globalIndex const rankOffset, + integer const hasCapPressure, + STENCILWRAPPER const & stencilWrapper, + DofNumberAccessor const & dofNumberAccessor, + CompFlowAccessors const & compFlowAccessors, + ThermalCompFlowAccessors const & thermalCompFlowAccessors, + MultiFluidAccessors const & multiFluidAccessors, + ThermalMultiFluidAccessors const & thermalMultiFluidAccessors, + CapPressureAccessors const & capPressureAccessors, + PermeabilityAccessors const & permeabilityAccessors, + ThermalConductivityAccessors const & thermalConductivityAccessors, + real64 const & dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + : Base( numPhases, + rankOffset, + hasCapPressure, + stencilWrapper, + dofNumberAccessor, + compFlowAccessors, + multiFluidAccessors, + capPressureAccessors, + permeabilityAccessors, + dt, + localMatrix, + localRhs ), + m_temp( thermalCompFlowAccessors.get( extrinsicMeshData::flow::temperature {} ) ), + m_dTemp( thermalCompFlowAccessors.get( extrinsicMeshData::flow::deltaTemperature {} ) ), + m_dPhaseMob_dTemp( thermalCompFlowAccessors.get( extrinsicMeshData::flow::dPhaseMobility_dTemperature {} ) ), + m_dPhaseVolFrac_dTemp( thermalCompFlowAccessors.get( extrinsicMeshData::flow::dPhaseVolumeFraction_dTemperature {} ) ), + m_phaseEnthalpy( thermalMultiFluidAccessors.get( extrinsicMeshData::multifluid::phaseEnthalpy {} ) ), + m_dPhaseEnthalpy( thermalMultiFluidAccessors.get( extrinsicMeshData::multifluid::dPhaseEnthalpy {} ) ), + m_thermalConductivity( thermalConductivityAccessors.get( extrinsicMeshData::thermalconductivity::effectiveConductivity {} ) ) + {} + + struct StackVariables : public Base::StackVariables + { +public: + + GEOSX_HOST_DEVICE + StackVariables( localIndex const size, localIndex numElems ) + : Base::StackVariables( size, numElems ), + dCompFlux_dT( size, numComp ), + energyFlux( 0.0 ), + dEnergyFlux_dP( size ), + dEnergyFlux_dT( size ), + dEnergyFlux_dC( size, numComp ) + {} + + using Base::StackVariables::stencilSize; + using Base::StackVariables::numFluxElems; + using Base::StackVariables::transmissibility; + using Base::StackVariables::dTrans_dPres; + using Base::StackVariables::dofColIndices; + using Base::StackVariables::localFlux; + using Base::StackVariables::localFluxJacobian; + + // Component fluxes and derivatives + + /// Derivatives of component fluxes wrt temperature + stackArray2d< real64, maxStencilSize * numComp > dCompFlux_dT; + + // Thermal transmissibility (for now, no derivatives) + + real64 thermalTransmissibility[maxNumConns][2]{}; + + // Energy fluxes and derivatives + + /// Energy fluxes + real64 energyFlux; + /// 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 component densities + stackArray2d< real64, maxStencilSize * numComp > dEnergyFlux_dC; + + }; + + /** + * @brief Compute the local flux contributions to the residual and Jacobian + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + */ + GEOSX_HOST_DEVICE + void computeFlux( localIndex const iconn, + StackVariables & stack ) const + { + using Deriv = multifluid::DerivativeOffset; + + // *********************************************** + // 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 dCompFlux_dT and the enthalpy flux requires quantities already computed in the base computeFlux, + // such as potGrad, phaseFlux, 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, [&] ( integer const ip, + localIndex const k_up, + localIndex const er_up, + localIndex const esr_up, + localIndex const ei_up, + real64 const & potGrad, + real64 const & phaseFlux, + real64 const (&dPhaseFlux_dP)[maxStencilSize], + real64 const (&dPhaseFlux_dC)[maxStencilSize][numComp] ) + { + GEOSX_UNUSED_VAR( dPhaseFlux_dP, dPhaseFlux_dC ); + + // We are in the loop over phases, ip provides the current phase index. + + // Step 1: compute the derivatives of the mean density at the interface wrt temperature + + stackArray1d< real64, maxNumElems > dDensMean_dT( stack.numFluxElems ); + + for( integer i = 0; i < stack.numFluxElems; ++i ) + { + localIndex const er = m_seri( iconn, i ); + localIndex const esr = m_sesri( iconn, i ); + localIndex const ei = m_sei( iconn, i ); + + real64 const dDens_dT = m_dPhaseMassDens[er][esr][ei][0][ip][Deriv::dT]; + dDensMean_dT[i] = 0.5 * dDens_dT; + } + + // Step 2: compute the derivatives of the phase potential difference wrt temperature + //***** calculation of flux ***** + + stackArray1d< real64, maxStencilSize > dPresGrad_dT( stack.stencilSize ); + stackArray1d< real64, maxStencilSize > dGravHead_dT( stack.numFluxElems ); + + // compute potential difference MPFA-style + for( integer i = 0; i < stack.stencilSize; ++i ) + { + localIndex const er = m_seri( iconn, i ); + localIndex const esr = m_sesri( iconn, i ); + localIndex const ei = m_sei( iconn, i ); + + // Step 2.1: compute derivative of capillary pressure wrt temperature + real64 dCapPressure_dT = 0.0; + if( m_hasCapPressure ) + { + for( integer jp = 0; jp < m_numPhases; ++jp ) + { + real64 const dCapPressure_dS = m_dPhaseCapPressure_dPhaseVolFrac[er][esr][ei][0][ip][jp]; + dCapPressure_dT += dCapPressure_dS * m_dPhaseVolFrac_dTemp[er][esr][ei][jp]; + } + } + + // Step 2.2: compute derivative of phase pressure difference wrt temperature + dPresGrad_dT[i] -= stack.transmissibility[0][i] * dCapPressure_dT; + real64 const gravD = stack.transmissibility[0][i] * m_gravCoef[er][esr][ei]; + + // Step 2.3: compute derivative of gravity potential difference wrt temperature + for( integer j = 0; j < stack.numFluxElems; ++j ) + { + dGravHead_dT[j] += dDensMean_dT[j] * gravD; + } + } + + // Step 3: compute the derivatives of the (upwinded) compFlux wrt temperature + // *** upwinding *** + + // note: the upwinding is done in the base class, which is in charge of + // computing the following quantities: potGrad, phaseFlux, k_up, er_up, esr_up, ei_up + + real64 dPhaseFlux_dT[maxStencilSize]{}; + + // Step 3.1: compute the derivative of phase flux wrt temperature + for( integer ke = 0; ke < stack.stencilSize; ++ke ) + { + dPhaseFlux_dT[ke] += dPresGrad_dT[ke]; + } + for( integer ke = 0; ke < stack.numFluxElems; ++ke ) + { + dPhaseFlux_dT[ke] -= dGravHead_dT[ke]; + } + for( integer ke = 0; ke < stack.stencilSize; ++ke ) + { + dPhaseFlux_dT[ke] *= m_phaseMob[er_up][esr_up][ei_up][ip]; + } + dPhaseFlux_dT[k_up] += m_dPhaseMob_dTemp[er_up][esr_up][ei_up][ip] * potGrad; + + // Step 3.2: compute the derivative of component flux wrt temperature + + // slice some constitutive arrays to avoid too much indexing in component loop + arraySlice1d< real64 const, multifluid::USD_PHASE_COMP - 3 > phaseCompFracSub = + m_phaseCompFrac[er_up][esr_up][ei_up][0][ip]; + arraySlice2d< real64 const, multifluid::USD_PHASE_COMP_DC - 3 > dPhaseCompFracSub = + m_dPhaseCompFrac[er_up][esr_up][ei_up][0][ip]; + + for( integer ic = 0; ic < numComp; ++ic ) + { + real64 const ycp = phaseCompFracSub[ic]; + for( integer ke = 0; ke < stack.stencilSize; ++ke ) + { + stack.dCompFlux_dT[ke][ic] += dPhaseFlux_dT[ke] * ycp; + } + stack.dCompFlux_dT[k_up][ic] += phaseFlux * dPhaseCompFracSub[ic][Deriv::dT]; + } + + // Step 4: compute the enthalpy flux + + real64 const enthalpy = m_phaseEnthalpy[er_up][esr_up][ei_up][0][ip]; + stack.energyFlux += phaseFlux * enthalpy; + + for( integer ke = 0; ke < stack.stencilSize; ++ke ) + { + stack.dEnergyFlux_dP[ke] += dPhaseFlux_dP[ke] * enthalpy; + stack.dEnergyFlux_dT[ke] += dPhaseFlux_dT[ke] * enthalpy; + + for( integer jc = 0; jc < numComp; ++jc ) + { + stack.dEnergyFlux_dC[ke][jc] += dPhaseFlux_dC[ke][jc] * enthalpy; + } + } + + stack.dEnergyFlux_dP[k_up] += phaseFlux * m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip][Deriv::dP]; + stack.dEnergyFlux_dT[k_up] += phaseFlux * m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip][Deriv::dT]; + + real64 dProp_dC[numComp]{}; + applyChainRule( numComp, + m_dCompFrac_dCompDens[er_up][esr_up][ei_up], + m_dPhaseEnthalpy[er_up][esr_up][ei_up][0][ip], + dProp_dC, + Deriv::dC ); + for( integer jc = 0; jc < numComp; ++jc ) + { + stack.dEnergyFlux_dC[k_up][jc] += phaseFlux * dProp_dC[jc]; + } + } ); + + // ***************************************************** + // Computation of the conduction term in the energy flux + // Note that the phase 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 + // Below, the thermal conductivity used to compute (explicitly) the thermal conducivity + // To avoid modifying the signature of the "computeWeights" function for now, we pass m_thermalConductivity twice + // TODO: modify computeWeights to accomodate explicit coefficients + 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 + + // Step 2: compute temperature difference at the interface + for( integer i = 0; i < stack.stencilSize; ++i ) + { + localIndex const er = m_seri( iconn, i ); + localIndex const esr = m_sesri( iconn, i ); + localIndex const ei = m_sei( iconn, i ); + + stack.energyFlux += stack.thermalTransmissibility[0][i] * ( m_temp[er][esr][ei] + m_dTemp[er][esr][ei] ); + stack.dEnergyFlux_dT[i] += stack.thermalTransmissibility[0][i]; + } + + // ********************************************************************************** + // At this point, we have computed the energyFlux and the compFlux for all components + // We have to do two things here: + // 1) Add dCompFlux_dTemp to the localFluxJacobian of the component mass balance equations + // 2) Add energyFlux and its derivatives to the localFlux(Jacobian) of the energy balance equation + + // Step 1: add dCompFlux_dTemp to localFluxJacobian + for( integer ic = 0; ic < numComp; ++ic ) + { + for( integer ke = 0; ke < stack.stencilSize; ++ke ) + { + integer const localDofIndexTemp = ke * numDof + numDof - 1; + stack.localFluxJacobian[ic][localDofIndexTemp] = m_dt * stack.dCompFlux_dT[ke][ic]; + stack.localFluxJacobian[numEqn + ic][localDofIndexTemp] = -m_dt * stack.dCompFlux_dT[ke][ic]; + } + } + + // Step 2: add energyFlux and its derivatives to localFlux and localFluxJacobian + integer const localRowIndexEnergy = numEqn-1; + stack.localFlux[localRowIndexEnergy] = m_dt * stack.energyFlux; + stack.localFlux[numEqn + localRowIndexEnergy] = -m_dt * stack.energyFlux; + + for( integer ke = 0; ke < stack.stencilSize; ++ke ) + { + integer const localDofIndexPres = ke * numDof; + stack.localFluxJacobian[localRowIndexEnergy][localDofIndexPres] = m_dt * stack.dEnergyFlux_dP[ke]; + stack.localFluxJacobian[numEqn + localRowIndexEnergy][localDofIndexPres] = -m_dt * stack.dEnergyFlux_dP[ke]; + integer const localDofIndexTemp = localDofIndexPres + numDof - 1; + stack.localFluxJacobian[localRowIndexEnergy][localDofIndexTemp] = m_dt * stack.dEnergyFlux_dT[ke]; + stack.localFluxJacobian[numEqn + localRowIndexEnergy][localDofIndexTemp] = -m_dt * stack.dEnergyFlux_dT[ke]; + + for( integer jc = 0; jc < numComp; ++jc ) + { + integer const localDofIndexComp = localDofIndexPres + jc + 1; + stack.localFluxJacobian[localRowIndexEnergy][localDofIndexComp] = m_dt * stack.dEnergyFlux_dC[ke][jc]; + stack.localFluxJacobian[numEqn + localRowIndexEnergy][localDofIndexComp] = -m_dt * stack.dEnergyFlux_dC[ke][jc]; + } + } + } + + /** + * @brief Performs the complete phase for the kernel. + * @param[in] iconn the connection index + * @param[inout] stack the stack variables + */ + GEOSX_HOST_DEVICE + void complete( localIndex const iconn, + StackVariables & stack ) const + { + // Call Case::complete to assemble the component mass balance equations (i = 0 to i = numDof-2) + // In the lambda, add contribution to residual and jacobian into the energy balance equation + Base::complete( iconn, stack, [&] ( integer const i, + localIndex const localRow ) + { + // beware, there is volume balance eqn in m_localRhs and m_localMatrix! + RAJA::atomicAdd( parallelDeviceAtomic{}, &AbstractBase::m_localRhs[localRow + numEqn], stack.localFlux[i * numEqn + numEqn-1] ); + AbstractBase::m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic > + ( localRow + numEqn, + stack.dofColIndices.data(), + stack.localFluxJacobian[i * numEqn + numEqn-1].dataIfContiguous(), + stack.stencilSize * numDof ); + + } ); + } + +protected: + + /// Views on temperature + ElementViewConst< arrayView1d< real64 const > > const m_temp; + ElementViewConst< arrayView1d< real64 const > > const m_dTemp; + + /// Views on derivatives of phase mobilities, volume fractions, mass densities and phase comp fractions + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_dPhaseMob_dTemp; + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const m_dPhaseVolFrac_dTemp; + + /// Views on phase enthalpies + ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const m_phaseEnthalpy; + ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const m_dPhaseEnthalpy; + + /// View on thermal conductivity + ElementViewConst< arrayView3d< real64 const > > m_thermalConductivity; + // for now, we treat thermal conductivity explicitly + +}; + +/** + * @class FaceBasedAssemblyKernelFactory + */ +class FaceBasedAssemblyKernelFactory +{ +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] numComps the number of fluid components + * @param[in] numPhases the number of fluid phases + * @param[in] rankOffset the offset of my MPI rank + * @param[in] dofKey string to get the element degrees of freedom numbers + * @param[in] hasCapPressure flag specifying whether capillary pressure is used or not + * @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, typename STENCILWRAPPER > + static void + createAndLaunch( integer const numComps, + integer const numPhases, + globalIndex const rankOffset, + string const & dofKey, + integer const hasCapPressure, + string const & solverName, + ElementRegionManager const & elemManager, + STENCILWRAPPER const & stencilWrapper, + real64 const & dt, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + isothermalCompositionalMultiphaseBaseKernels:: + internal::kernelLaunchSelectorCompSwitch( numComps, [&] ( auto NC ) + { + integer constexpr NUM_COMP = NC(); + integer constexpr NUM_DOF = NC()+2; + + ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > dofNumberAccessor = + elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey ); + dofNumberAccessor.setName( solverName + "/accessors/" + dofKey ); + + using KernelType = FaceBasedAssemblyKernel< NUM_COMP, NUM_DOF, STENCILWRAPPER >; + typename KernelType::CompFlowAccessors compFlowAccessors( elemManager, solverName ); + typename KernelType::ThermalCompFlowAccessors thermalCompFlowAccessors( elemManager, solverName ); + typename KernelType::MultiFluidAccessors multiFluidAccessors( elemManager, solverName ); + typename KernelType::ThermalMultiFluidAccessors thermalMultiFluidAccessors( elemManager, solverName ); + typename KernelType::CapPressureAccessors capPressureAccessors( elemManager, solverName ); + typename KernelType::PermeabilityAccessors permeabilityAccessors( elemManager, solverName ); + typename KernelType::ThermalConductivityAccessors thermalConductivityAccessors( elemManager, solverName ); + + KernelType kernel( numPhases, rankOffset, hasCapPressure, stencilWrapper, dofNumberAccessor, + compFlowAccessors, thermalCompFlowAccessors, multiFluidAccessors, thermalMultiFluidAccessors, + capPressureAccessors, permeabilityAccessors, thermalConductivityAccessors, + dt, localMatrix, localRhs ); + KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); + } ); + } +}; + + +} // namespace thermalCompositionalMultiphaseFVMKernels + +} // namespace geosx + + +#endif //GEOSX_PHYSICSSOLVERS_FLUIDFLOW_THERMALCOMPOSITIONALMULTIPHASEFVMKERNELS_HPP diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 0304ae0823a..d6a4421f16c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -30,10 +30,11 @@ #include "constitutive/relativePermeability/RelativePermeabilityExtrinsicData.hpp" #include "dataRepository/Group.hpp" #include "mesh/DomainPartition.hpp" -#include "mesh/WellElementSubRegion.hpp" #include "mesh/PerforationData.hpp" +#include "mesh/WellElementSubRegion.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" +#include "physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellExtrinsicData.hpp" #include "physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp" #include "physicsSolvers/fluidFlow/wells/SinglePhaseWellKernels.hpp" @@ -154,6 +155,7 @@ void CompositionalMultiphaseWell::registerDataOnMesh( Group & meshBodies ) subRegion.registerExtrinsicData< extrinsicMeshData::well::deltaPressure >( getName() ); subRegion.registerExtrinsicData< extrinsicMeshData::well::temperature >( getName() ); + subRegion.registerExtrinsicData< extrinsicMeshData::well::deltaTemperature >( getName() ); // The resizing of the arrays needs to happen here, before the call to initializePreSubGroups, // to make sure that the dimensions are properly set before the timeHistoryOutput starts its initialization. @@ -475,7 +477,7 @@ void CompositionalMultiphaseWell::updateComponentFraction( WellElementSubRegion { GEOSX_MARK_FUNCTION; - compositionalMultiphaseBaseKernels:: + isothermalCompositionalMultiphaseBaseKernels:: ComponentFractionKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, subRegion ); @@ -745,6 +747,7 @@ void CompositionalMultiphaseWell::updateFluidModel( WellElementSubRegion & subRe arrayView1d< real64 const > const & pres = subRegion.getExtrinsicData< extrinsicMeshData::well::pressure >(); arrayView1d< real64 const > const & dPres = subRegion.getExtrinsicData< extrinsicMeshData::well::deltaPressure >(); arrayView1d< real64 const > const & temp = subRegion.getExtrinsicData< extrinsicMeshData::well::temperature >(); + arrayView1d< real64 const > const dTemp = subRegion.getExtrinsicData< extrinsicMeshData::well::deltaTemperature >(); arrayView2d< real64 const, compflow::USD_COMP > const & compFrac = subRegion.getExtrinsicData< extrinsicMeshData::well::globalCompFraction >(); @@ -757,12 +760,15 @@ void CompositionalMultiphaseWell::updateFluidModel( WellElementSubRegion & subRe using ExecPolicy = typename FluidType::exec_policy; typename FluidType::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); - compositionalMultiphaseBaseKernels::FluidUpdateKernel::launch< ExecPolicy >( subRegion.size(), - fluidWrapper, - pres, - dPres, - temp, - compFrac ); + thermalCompositionalMultiphaseBaseKernels:: + FluidUpdateKernel:: + launch< ExecPolicy >( subRegion.size(), + fluidWrapper, + pres, + dPres, + temp, + dTemp, + compFrac ); } ); } @@ -772,7 +778,7 @@ void CompositionalMultiphaseWell::updatePhaseVolumeFraction( WellElementSubRegio string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase & fluid = subRegion.getConstitutiveModel< MultiFluidBase >( fluidName ); - compositionalMultiphaseBaseKernels:: + isothermalCompositionalMultiphaseBaseKernels:: PhaseVolumeFractionKernelFactory:: createAndLaunch< parallelDevicePolicy<> >( m_numComponents, m_numPhases, @@ -793,7 +799,6 @@ void CompositionalMultiphaseWell::updateTotalMassDensity( WellElementSubRegion & m_numPhases, subRegion, fluid ); - } void CompositionalMultiphaseWell::updateSubRegionState( MeshLevel const & meshLevel, @@ -875,25 +880,26 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain ) // 1) Loop over all perforations to compute an average mixture density and component fraction // 2) Initialize the reference pressure // 3) Estimate the pressures in the well elements using the average density - PresTempCompFracInitializationKernel::launch( perforationData.size(), - subRegion.size(), - numComp, - numPhase, - perforationData.getNumPerforationsGlobal(), - wellControls, - 0.0, // initialization done at t = 0 - resCompFlowAccessors.get( extrinsicMeshData::flow::pressure{} ), - resCompFlowAccessors.get( extrinsicMeshData::flow::temperature{} ), - resCompFlowAccessors.get( extrinsicMeshData::flow::globalCompDensity{} ), - resCompFlowAccessors.get( extrinsicMeshData::flow::phaseVolumeFraction{} ), - resMultiFluidAccessors.get( extrinsicMeshData::multifluid::phaseMassDensity{} ), - resElementRegion, - resElementSubRegion, - resElementIndex, - wellElemGravCoef, - wellElemPressure, - wellElemTemp, - wellElemCompFrac ); + PresTempCompFracInitializationKernel:: + launch( perforationData.size(), + subRegion.size(), + numComp, + numPhase, + perforationData.getNumPerforationsGlobal(), + wellControls, + 0.0, // initialization done at t = 0 + resCompFlowAccessors.get( extrinsicMeshData::flow::pressure{} ), + resCompFlowAccessors.get( extrinsicMeshData::flow::temperature{} ), + resCompFlowAccessors.get( extrinsicMeshData::flow::globalCompDensity{} ), + resCompFlowAccessors.get( extrinsicMeshData::flow::phaseVolumeFraction{} ), + resMultiFluidAccessors.get( extrinsicMeshData::multifluid::phaseMassDensity{} ), + resElementRegion, + resElementSubRegion, + resElementIndex, + wellElemGravCoef, + wellElemPressure, + wellElemTemp, + wellElemCompFrac ); // get well secondary variables on well elements string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); @@ -906,11 +912,13 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain ) { typename TYPEOFREF( castedFluid ) ::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); - compositionalMultiphaseBaseKernels::FluidUpdateKernel::launch< serialPolicy >( subRegion.size(), - fluidWrapper, - wellElemPressure, - wellElemTemp, - wellElemCompFrac ); + thermalCompositionalMultiphaseBaseKernels:: + FluidUpdateKernel:: + launch< serialPolicy >( subRegion.size(), + fluidWrapper, + wellElemPressure, + wellElemTemp, + wellElemCompFrac ); } ); CompDensInitializationKernel::launch( subRegion.size(), @@ -926,13 +934,15 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain ) // 6) Estimate the well rates // TODO: initialize rates using perforation rates - compositionalMultiphaseWellKernels::RateInitializationKernel::launch( subRegion.size(), - m_targetPhaseIndex, - wellControls, - 0.0, // initialization done at t = 0 - wellElemPhaseDens, - wellElemTotalDens, - connRate ); + compositionalMultiphaseWellKernels:: + RateInitializationKernel:: + launch( subRegion.size(), + m_targetPhaseIndex, + wellControls, + 0.0, // initialization done at t = 0 + wellElemPhaseDens, + wellElemTotalDens, + connRate ); } ); } ); } @@ -978,7 +988,7 @@ void CompositionalMultiphaseWell::assembleFluxTerms( real64 const GEOSX_UNUSED_P arrayView3d< real64 const, compflow::USD_COMP_DC > const & dWellElemCompFrac_dCompDens = subRegion.getExtrinsicData< extrinsicMeshData::well::dGlobalCompFraction_dGlobalCompDensity >(); - compositionalMultiphaseBaseKernels:: + isothermalCompositionalMultiphaseBaseKernels:: KernelLaunchSelector1< FluxKernel >( numFluidComponents(), subRegion.size(), dofManager.rankOffset(), @@ -1058,7 +1068,7 @@ void CompositionalMultiphaseWell::assembleAccumulationTerms( DomainPartition con arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const & wellElemPhaseCompFrac = fluid.phaseCompFraction(); arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > const & dWellElemPhaseCompFrac = fluid.dPhaseCompFraction(); - compositionalMultiphaseBaseKernels:: + isothermalCompositionalMultiphaseBaseKernels:: KernelLaunchSelector1< AccumulationKernel >( numFluidComponents(), subRegion.size(), numFluidPhases(), @@ -1121,7 +1131,7 @@ void CompositionalMultiphaseWell::assembleVolumeBalanceTerms( DomainPartition co arrayView1d< real64 const > const & wellElemVolume = subRegion.getReference< array1d< real64 > >( ElementSubRegionBase::viewKeyStruct::elementVolumeString() ); - compositionalMultiphaseBaseKernels:: + isothermalCompositionalMultiphaseBaseKernels:: KernelLaunchSelector1< VolumeBalanceKernel >( numFluidComponents(), subRegion.size(), numFluidPhases(), @@ -1174,25 +1184,24 @@ CompositionalMultiphaseWell::calculateResidualNorm( DomainPartition const & doma WellControls const & wellControls = getWellControls( subRegion ); - ResidualNormKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localRhs, - dofManager.rankOffset(), - subRegion.isLocallyOwned(), - subRegion.getTopWellElementIndex(), - m_numComponents, - numDofPerWellElement(), - m_targetPhaseIndex, - wellControls, - wellElemDofNumber, - wellElemGhostRank, - wellElemVolume, - wellElemPhaseDensOld, - wellElemTotalDensOld, - m_currentTime + m_currentDt, // residual normalized with rate of the end of the - // time - // interval - m_currentDt, - &localResidualNorm ); + ResidualNormKernel::launch< parallelDevicePolicy<> >( localRhs, + dofManager.rankOffset(), + subRegion.isLocallyOwned(), + subRegion.getTopWellElementIndex(), + m_numComponents, + numDofPerWellElement(), + m_targetPhaseIndex, + wellControls, + wellElemDofNumber, + wellElemGhostRank, + wellElemVolume, + wellElemPhaseDensOld, + wellElemTotalDensOld, + m_currentTime + m_currentDt, // residual normalized with rate of the end of the + // time + // interval + m_currentDt, + &localResidualNorm ); } ); } ); return sqrt( MpiWrapper::sum( localResidualNorm, MPI_COMM_GEOSX ) ); @@ -1243,18 +1252,17 @@ CompositionalMultiphaseWell::scalingForSystemSolution( DomainPartition const & d real64 const subRegionScalingFactor = - SolutionScalingKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localSolution, - dofManager.rankOffset(), - numFluidComponents(), - wellElemDofNumber, - wellElemGhostRank, - wellElemPressure, - dWellElemPressure, - wellElemCompDens, - dWellElemCompDens, - m_maxRelativePresChange, - m_maxCompFracChange ); + SolutionScalingKernel::launch< parallelDevicePolicy<> >( localSolution, + dofManager.rankOffset(), + numFluidComponents(), + wellElemDofNumber, + wellElemGhostRank, + wellElemPressure, + dWellElemPressure, + wellElemCompDens, + dWellElemCompDens, + m_maxRelativePresChange, + m_maxCompFracChange ); if( subRegionScalingFactor < scalingFactor ) @@ -1306,18 +1314,17 @@ CompositionalMultiphaseWell::checkSystemSolution( DomainPartition const & domain subRegion.getExtrinsicData< extrinsicMeshData::well::deltaGlobalCompDensity >(); localIndex const subRegionSolutionCheck = - SolutionCheckKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localSolution, - dofManager.rankOffset(), - numFluidComponents(), - wellElemDofNumber, - wellElemGhostRank, - wellElemPressure, - dWellElemPressure, - wellElemCompDens, - dWellElemCompDens, - m_allowCompDensChopping, - scalingFactor ); + SolutionCheckKernel::launch< parallelDevicePolicy<> >( localSolution, + dofManager.rankOffset(), + numFluidComponents(), + wellElemDofNumber, + wellElemGhostRank, + wellElemPressure, + dWellElemPressure, + wellElemCompDens, + dWellElemCompDens, + m_allowCompDensChopping, + scalingFactor ); if( subRegionSolutionCheck == 0 ) { @@ -1400,7 +1407,7 @@ void CompositionalMultiphaseWell::computePerforationRates( MeshLevel const & mes PerforationKernel::MultiFluidAccessors resMultiFluidAccessors( meshLevel.getElemManager(), flowSolver.getName() ); PerforationKernel::RelPermAccessors resRelPermAccessors( meshLevel.getElemManager(), flowSolver.getName() ); - compositionalMultiphaseBaseKernels:: + isothermalCompositionalMultiphaseBaseKernels:: KernelLaunchSelector2< PerforationKernel >( numFluidComponents(), numFluidPhases(), perforationData->size(), @@ -1438,7 +1445,6 @@ void CompositionalMultiphaseWell::computePerforationRates( MeshLevel const & mes compPerfRate, dCompPerfRate_dPres, dCompPerfRate_dComp ); - } @@ -1673,7 +1679,7 @@ void CompositionalMultiphaseWell::assemblePressureRelations( DomainPartition con bool controlHasSwitched = false; - compositionalMultiphaseBaseKernels:: + isothermalCompositionalMultiphaseBaseKernels:: KernelLaunchSelector1< PressureRelationKernel >( numFluidComponents(), subRegion.size(), dofManager.rankOffset(), @@ -1799,4 +1805,4 @@ void CompositionalMultiphaseWell::implicitStepComplete( real64 const & GEOSX_UNU } REGISTER_CATALOG_ENTRY( SolverBase, CompositionalMultiphaseWell, string const &, Group * const ) -}// namespace geosx +} // namespace geosx diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellExtrinsicData.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellExtrinsicData.hpp index af313c1dda3..a7a8f8ae9fc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellExtrinsicData.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellExtrinsicData.hpp @@ -48,6 +48,14 @@ EXTRINSIC_MESH_DATA_TRAIT( temperature, WRITE_AND_READ, "Temperature" ); +EXTRINSIC_MESH_DATA_TRAIT( deltaTemperature, + "deltaWellTemperature", + array1d< real64 >, + 0, + NOPLOT, + NO_WRITE, + "Accumulated temperature updates" ); + EXTRINSIC_MESH_DATA_TRAIT( globalCompDensity, "globalCompDensity", array2dLayoutComp, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp index aca6bc358fa..9ab30d56511 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.cpp @@ -482,9 +482,9 @@ FluxKernel:: // Apply equation/variable change transformation(s) real64 work[NC+1]{}; - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, 1, 2, localFluxJacobian_dRate, work ); - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC + 1, 2, localFluxJacobian_dPresCompUp, work ); - shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( NC, 2, localFlux ); + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC, 1, 2, localFluxJacobian_dRate, work ); + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC, NC + 1, 2, localFluxJacobian_dPresCompUp, work ); + shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( NC, NC, 2, localFlux ); for( integer i = 0; i < 2*NC; ++i ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp index f49bf82101e..361938c0c68 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellKernels.hpp @@ -27,9 +27,9 @@ #include "constitutive/relativePermeability/RelativePermeabilityExtrinsicData.hpp" #include "mesh/ElementRegionManager.hpp" #include "mesh/ObjectManagerBase.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseExtrinsicData.hpp" +#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" #include "physicsSolvers/fluidFlow/StencilAccessors.hpp" #include "physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWellExtrinsicData.hpp" #include "physicsSolvers/fluidFlow/wells/WellControls.hpp" @@ -509,12 +509,11 @@ struct RateInitializationKernel * @brief Define the interface for the property kernel in charge of computing the total mass density */ template< integer NUM_COMP, integer NUM_PHASE > -class TotalMassDensityKernel : public compositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > +class TotalMassDensityKernel : public isothermalCompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP > { public: - using Base = compositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; - + using Base = isothermalCompositionalMultiphaseBaseKernels::PropertyKernelBase< NUM_COMP >; using Base::numComp; /// Compile time value for the number of phases @@ -545,10 +544,10 @@ class TotalMassDensityKernel : public compositionalMultiphaseBaseKernels::Proper * @param[in] ei the element index * @param[in] totalMassDensityKernelOp the function used to customize the kernel */ - template< typename FUNC = compositionalMultiphaseBaseKernels::NoOpFunc > + template< typename FUNC = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc > GEOSX_HOST_DEVICE void compute( localIndex const ei, - FUNC && totalMassDensityKernelOp = compositionalMultiphaseBaseKernels::NoOpFunc{} ) const + FUNC && totalMassDensityKernelOp = isothermalCompositionalMultiphaseBaseKernels::NoOpFunc{} ) const { using Deriv = multifluid::DerivativeOffset; @@ -634,7 +633,7 @@ class TotalMassDensityKernelFactory { if( numPhase == 2 ) { - compositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { integer constexpr NUM_COMP = NC(); TotalMassDensityKernel< NUM_COMP, 2 > kernel( subRegion, fluid ); @@ -643,7 +642,7 @@ class TotalMassDensityKernelFactory } else if( numPhase == 3 ) { - compositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) + isothermalCompositionalMultiphaseBaseKernels::internal::kernelLaunchSelectorCompSwitch( numComp, [&] ( auto NC ) { integer constexpr NUM_COMP = NC(); TotalMassDensityKernel< NUM_COMP, 3 > kernel( subRegion, fluid ); @@ -659,9 +658,9 @@ class TotalMassDensityKernelFactory struct ResidualNormKernel { - template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > + template< typename POLICY > static void - launch( LOCAL_VECTOR const localResidual, + launch( arrayView1d< real64 const > const & localResidual, globalIndex const rankOffset, bool const isLocallyOwned, localIndex const iwelemControl, @@ -688,7 +687,7 @@ struct ResidualNormKernel real64 const absTargetTotalRate = LvArray::math::abs( targetTotalRate ); real64 const absTargetPhaseRate = LvArray::math::abs( targetPhaseRate ); - RAJA::ReduceSum< REDUCE_POLICY, real64 > sumScaled( 0.0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > sumScaled( 0.0 ); forAll< POLICY >( wellElemDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) { @@ -775,9 +774,9 @@ struct ResidualNormKernel struct SolutionScalingKernel { - template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > + template< typename POLICY > static real64 - launch( LOCAL_VECTOR const localSolution, + launch( arrayView1d< real64 const > const & localSolution, globalIndex const rankOffset, integer const numComponents, arrayView1d< globalIndex const > const & wellElemDofNumber, @@ -791,7 +790,7 @@ struct SolutionScalingKernel { real64 constexpr eps = minDensForDivision; - RAJA::ReduceMin< REDUCE_POLICY, real64 > minVal( 1.0 ); + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > minVal( 1.0 ); forAll< POLICY >( wellElemDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) { @@ -848,9 +847,9 @@ struct SolutionScalingKernel struct SolutionCheckKernel { - template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > + template< typename POLICY > static localIndex - launch( LOCAL_VECTOR const localSolution, + launch( arrayView1d< real64 const > const & localSolution, globalIndex const rankOffset, integer const numComponents, arrayView1d< globalIndex const > const & wellElemDofNumber, @@ -866,7 +865,7 @@ struct SolutionCheckKernel real64 constexpr eps = minDensForDivision; - RAJA::ReduceMin< REDUCE_POLICY, localIndex > minVal( 1 ); + RAJA::ReduceMin< ReducePolicy< POLICY >, localIndex > minVal( 1 ); forAll< POLICY >( wellElemDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 2d98c8d903d..a9148c91aa3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -707,20 +707,19 @@ SinglePhaseWell::calculateResidualNorm( DomainPartition const & domain, WellControls const & wellControls = getWellControls( subRegion ); - ResidualNormKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localRhs, - dofManager.rankOffset(), - subRegion.isLocallyOwned(), - subRegion.getTopWellElementIndex(), - wellControls, - wellElemDofNumber, - wellElemGhostRank, - wellElemVolume, - wellElemDensityOld, - m_currentTime + m_currentDt, // residual normalized with rate of the end of the - // time interval - m_currentDt, - &localResidualNorm ); + ResidualNormKernel::launch< parallelDevicePolicy<> >( localRhs, + dofManager.rankOffset(), + subRegion.isLocallyOwned(), + subRegion.getTopWellElementIndex(), + wellControls, + wellElemDofNumber, + wellElemGhostRank, + wellElemVolume, + wellElemDensityOld, + m_currentTime + m_currentDt, // residual normalized with rate of the end of the + // time interval + m_currentDt, + &localResidualNorm ); } ); } ); @@ -763,14 +762,14 @@ bool SinglePhaseWell::checkSystemSolution( DomainPartition const & domain, // here we can reuse the flow solver kernel checking that pressures are positive localIndex const subRegionSolutionCheck = - singlePhaseWellKernels::SolutionCheckKernel::launch< parallelDevicePolicy<>, - parallelDeviceReduce >( localSolution, - dofManager.rankOffset(), - wellElemDofNumber, - wellElemGhostRank, - wellElemPressure, - dWellElemPressure, - scalingFactor ); + singlePhaseWellKernels:: + SolutionCheckKernel::launch< parallelDevicePolicy<> >( localSolution, + dofManager.rankOffset(), + wellElemDofNumber, + wellElemGhostRank, + wellElemPressure, + dWellElemPressure, + scalingFactor ); if( subRegionSolutionCheck == 0 ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWellKernels.hpp index 3b0e1cfc980..4c502b1e477 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWellKernels.hpp @@ -308,9 +308,9 @@ struct RateInitializationKernel struct ResidualNormKernel { - template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > + template< typename POLICY > static void - launch( LOCAL_VECTOR const localResidual, + launch( arrayView1d< real64 const > const & localResidual, globalIndex const rankOffset, bool const isLocallyOwned, localIndex const iwelemControl, @@ -328,7 +328,7 @@ struct ResidualNormKernel real64 const targetRate = wellControls.getTargetTotalRate( timeAtEndOfStep ); real64 const absTargetRate = LvArray::math::abs( targetRate ); - RAJA::ReduceSum< REDUCE_POLICY, real64 > sumScaled( 0.0 ); + RAJA::ReduceSum< ReducePolicy< POLICY >, real64 > sumScaled( 0.0 ); forAll< POLICY >( wellElemDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const iwelem ) { @@ -379,9 +379,9 @@ struct ResidualNormKernel struct SolutionCheckKernel { - template< typename POLICY, typename REDUCE_POLICY, typename LOCAL_VECTOR > + template< typename POLICY > static localIndex - launch( LOCAL_VECTOR const localSolution, + launch( arrayView1d< real64 const > const & localSolution, globalIndex const rankOffset, arrayView1d< globalIndex const > const & presDofNumber, arrayView1d< integer const > const & ghostRank, @@ -389,7 +389,7 @@ struct SolutionCheckKernel arrayView1d< real64 const > const & dPres, real64 const scalingFactor ) { - RAJA::ReduceMin< REDUCE_POLICY, localIndex > minVal( 1 ); + RAJA::ReduceMin< ReducePolicy< POLICY >, localIndex > minVal( 1 ); forAll< POLICY >( presDofNumber.size(), [=] GEOSX_HOST_DEVICE ( localIndex const ei ) { diff --git a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoir.cpp b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoir.cpp index 2addaae1b63..78f2a2d70ac 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoir.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoir.cpp @@ -302,8 +302,8 @@ void CompositionalMultiphaseReservoir::assembleCouplingTerms( real64 const time_ // Apply equation/variable change transformation(s) stackArray1d< real64, 2 * MAX_NUM_DOF > work( 2 * resNumDofs ); - shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComps, resNumDofs*2, 2, localPerfJacobian, work ); - shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComps, 2, localPerf ); + shiftBlockRowsAheadByOneAndReplaceFirstRowWithColumnSum( numComps, numComps, resNumDofs*2, 2, localPerfJacobian, work ); + shiftBlockElementsAheadByOneAndReplaceFirstElementWithSum( numComps, numComps, 2, localPerf ); for( localIndex i = 0; i < localPerf.size(); ++i ) { diff --git a/src/coreComponents/schema/docs/AcousticSEM_other.rst b/src/coreComponents/schema/docs/AcousticSEM_other.rst index 71642d4c7ae..edd823f6468 100644 --- a/src/coreComponents/schema/docs/AcousticSEM_other.rst +++ b/src/coreComponents/schema/docs/AcousticSEM_other.rst @@ -1,20 +1,20 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================= -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================= -indexSeismoTrace integer Count for output pressure at receivers -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -pressureNp1AtReceivers real64_array2d Pressure value at each receiver for each timestep -receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank -receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point -sourceConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds -sourceIsLocal integer_array Flag that indicates whether the source is local to this MPI rank -sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point -sourceValue real64_array2d Source Value of the sources -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================= +========================= =================================================================================================================================================== ======================================================================= +Name Type Description +========================= =================================================================================================================================================== ======================================================================= +indexSeismoTrace integer Count for output pressure at receivers +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +pressureNp1AtReceivers real64_array2d Pressure value at each receiver for each timestep +receiverIsLocal integer_array Flag that indicates whether the receiver is local to this MPI rank +receiverNodeIds integer_array2d Indices of the nodes (in the right order) for each receiver point +sourceConstants real64_array2d Constant part of the receiver for the nodes listed in m_receiverNodeIds +sourceIsLocal integer_array Flag that indicates whether the source is local to this MPI rank +sourceNodeIds integer_array2d Indices of the nodes (in the right order) for each source point +sourceValue real64_array2d Source Value of the sources +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ======================================================================= diff --git a/src/coreComponents/schema/docs/BlackOilFluid_other.rst b/src/coreComponents/schema/docs/BlackOilFluid_other.rst index e6d3382d61a..b06aacaa1ac 100644 --- a/src/coreComponents/schema/docs/BlackOilFluid_other.rst +++ b/src/coreComponents/schema/docs/BlackOilFluid_other.rst @@ -8,7 +8,7 @@ dPhaseCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l dPhaseDensity real64_array4d Derivative of phase density with respect to pressure, temperature, and global component fractions dPhaseEnthalpy real64_array4d Derivative of phase enthalpy with respect to pressure, temperature, and global component fractions dPhaseFraction real64_array4d Derivative of phase fraction with respect to pressure, temperature, and global component fractions -dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fraction +dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fractions dPhaseMassDensity real64_array4d Derivative of phase mass density with respect to pressure, temperature, and global component fractions dPhaseViscosity real64_array4d Derivative of phase viscosity with respect to pressure, temperature, and global component fractions dTotalDensity real64_array3d Derivative of total density with respect to pressure, temperature, and global component fractions diff --git a/src/coreComponents/schema/docs/CO2BrineEzrokhiFluid_other.rst b/src/coreComponents/schema/docs/CO2BrineEzrokhiFluid_other.rst index b7df809a8fc..54d5845f0dc 100644 --- a/src/coreComponents/schema/docs/CO2BrineEzrokhiFluid_other.rst +++ b/src/coreComponents/schema/docs/CO2BrineEzrokhiFluid_other.rst @@ -7,7 +7,7 @@ dPhaseCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l dPhaseDensity real64_array4d Derivative of phase density with respect to pressure, temperature, and global component fractions dPhaseEnthalpy real64_array4d Derivative of phase enthalpy with respect to pressure, temperature, and global component fractions dPhaseFraction real64_array4d Derivative of phase fraction with respect to pressure, temperature, and global component fractions -dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fraction +dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fractions dPhaseMassDensity real64_array4d Derivative of phase mass density with respect to pressure, temperature, and global component fractions dPhaseViscosity real64_array4d Derivative of phase viscosity with respect to pressure, temperature, and global component fractions dTotalDensity real64_array3d Derivative of total density with respect to pressure, temperature, and global component fractions diff --git a/src/coreComponents/schema/docs/CO2BrineEzrokhiThermalFluid_other.rst b/src/coreComponents/schema/docs/CO2BrineEzrokhiThermalFluid_other.rst index b7df809a8fc..54d5845f0dc 100644 --- a/src/coreComponents/schema/docs/CO2BrineEzrokhiThermalFluid_other.rst +++ b/src/coreComponents/schema/docs/CO2BrineEzrokhiThermalFluid_other.rst @@ -7,7 +7,7 @@ dPhaseCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l dPhaseDensity real64_array4d Derivative of phase density with respect to pressure, temperature, and global component fractions dPhaseEnthalpy real64_array4d Derivative of phase enthalpy with respect to pressure, temperature, and global component fractions dPhaseFraction real64_array4d Derivative of phase fraction with respect to pressure, temperature, and global component fractions -dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fraction +dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fractions dPhaseMassDensity real64_array4d Derivative of phase mass density with respect to pressure, temperature, and global component fractions dPhaseViscosity real64_array4d Derivative of phase viscosity with respect to pressure, temperature, and global component fractions dTotalDensity real64_array3d Derivative of total density with respect to pressure, temperature, and global component fractions diff --git a/src/coreComponents/schema/docs/CO2BrinePhillipsFluid_other.rst b/src/coreComponents/schema/docs/CO2BrinePhillipsFluid_other.rst index b7df809a8fc..54d5845f0dc 100644 --- a/src/coreComponents/schema/docs/CO2BrinePhillipsFluid_other.rst +++ b/src/coreComponents/schema/docs/CO2BrinePhillipsFluid_other.rst @@ -7,7 +7,7 @@ dPhaseCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l dPhaseDensity real64_array4d Derivative of phase density with respect to pressure, temperature, and global component fractions dPhaseEnthalpy real64_array4d Derivative of phase enthalpy with respect to pressure, temperature, and global component fractions dPhaseFraction real64_array4d Derivative of phase fraction with respect to pressure, temperature, and global component fractions -dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fraction +dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fractions dPhaseMassDensity real64_array4d Derivative of phase mass density with respect to pressure, temperature, and global component fractions dPhaseViscosity real64_array4d Derivative of phase viscosity with respect to pressure, temperature, and global component fractions dTotalDensity real64_array3d Derivative of total density with respect to pressure, temperature, and global component fractions diff --git a/src/coreComponents/schema/docs/CO2BrinePhillipsThermalFluid_other.rst b/src/coreComponents/schema/docs/CO2BrinePhillipsThermalFluid_other.rst index b7df809a8fc..54d5845f0dc 100644 --- a/src/coreComponents/schema/docs/CO2BrinePhillipsThermalFluid_other.rst +++ b/src/coreComponents/schema/docs/CO2BrinePhillipsThermalFluid_other.rst @@ -7,7 +7,7 @@ dPhaseCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l dPhaseDensity real64_array4d Derivative of phase density with respect to pressure, temperature, and global component fractions dPhaseEnthalpy real64_array4d Derivative of phase enthalpy with respect to pressure, temperature, and global component fractions dPhaseFraction real64_array4d Derivative of phase fraction with respect to pressure, temperature, and global component fractions -dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fraction +dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fractions dPhaseMassDensity real64_array4d Derivative of phase mass density with respect to pressure, temperature, and global component fractions dPhaseViscosity real64_array4d Derivative of phase viscosity with respect to pressure, temperature, and global component fractions dTotalDensity real64_array3d Derivative of total density with respect to pressure, temperature, and global component fractions diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseFVM.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseFVM.rst index 902c77f2e61..ddbc407ff2b 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseFVM.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseFVM.rst @@ -9,6 +9,7 @@ computeCFLNumbers integer 0 Flag indicating whether CFL discretization string required Name of discretization object to use for this solver. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. inputFluxEstimate real64 1 Initial estimate of the input flux used only for residual scaling. This should be essentially equivalent to the input flux * dt. +isThermal integer 0 Flag indicating whether the problem is thermal or not. logLevel integer 0 Log level maxCompFractionChange real64 1 Maximum (absolute) change in a component fraction between two Newton iterations name string required A name is required for any non-unique nodes diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseFVM_other.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseFVM_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseFVM_other.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseFVM_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseFluid_other.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseFluid_other.rst index b7df809a8fc..54d5845f0dc 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseFluid_other.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseFluid_other.rst @@ -7,7 +7,7 @@ dPhaseCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l, 1l, 2l dPhaseDensity real64_array4d Derivative of phase density with respect to pressure, temperature, and global component fractions dPhaseEnthalpy real64_array4d Derivative of phase enthalpy with respect to pressure, temperature, and global component fractions dPhaseFraction real64_array4d Derivative of phase fraction with respect to pressure, temperature, and global component fractions -dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fraction +dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fractions dPhaseMassDensity real64_array4d Derivative of phase mass density with respect to pressure, temperature, and global component fractions dPhaseViscosity real64_array4d Derivative of phase viscosity with respect to pressure, temperature, and global component fractions dTotalDensity real64_array3d Derivative of total density with respect to pressure, temperature, and global component fractions diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM.rst index 96ee301c63e..9ac6d8bda01 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM.rst @@ -9,6 +9,7 @@ computeCFLNumbers integer 0 Flag indicating whether CFL discretization string required Name of discretization object to use for this solver. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. inputFluxEstimate real64 1 Initial estimate of the input flux used only for residual scaling. This should be essentially equivalent to the input flux * dt. +isThermal integer 0 Flag indicating whether the problem is thermal or not. logLevel integer 0 Log level maxCompFractionChange real64 1 Maximum (absolute) change in a component fraction between two Newton iterations maxRelativePressureChange real64 1 Maximum (relative) change in (face) pressure between two Newton iterations diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM_other.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM_other.rst index bf32e9bac59..8e2f2eb6a85 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM_other.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseHybridFVM_other.rst @@ -1,14 +1,14 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================ ================================================================ -Name Type Registered On Description -========================= =============================================================================================================================================================================================================================================================================================== ================================ ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -deltaFacePressure real64_array :ref:`DATASTRUCTURE_FaceManager` Accumulated face pressure updates -mimGravityCoefficient real64_array :ref:`DATASTRUCTURE_FaceManager` Mimetic gravity coefficient -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================ ================================================================ +========================= =================================================================================================================================================== ================================ ================================================================ +Name Type Registered On Description +========================= =================================================================================================================================================== ================================ ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +deltaFacePressure real64_array :ref:`DATASTRUCTURE_FaceManager` Accumulated face pressure updates +mimGravityCoefficient real64_array :ref:`DATASTRUCTURE_FaceManager` Mimetic gravity coefficient +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================ ================================================================ diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseReservoir_other.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseReservoir_other.rst index aa08bcefbb0..7acea5a6dbf 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseReservoir_other.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseReservoir_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/CompositionalMultiphaseWell_other.rst b/src/coreComponents/schema/docs/CompositionalMultiphaseWell_other.rst index b3804961a57..48a09333642 100644 --- a/src/coreComponents/schema/docs/CompositionalMultiphaseWell_other.rst +++ b/src/coreComponents/schema/docs/CompositionalMultiphaseWell_other.rst @@ -1,14 +1,14 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -WellControls node :ref:`DATASTRUCTURE_WellControls` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +WellControls node :ref:`DATASTRUCTURE_WellControls` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/ConstantThermalConductivity.rst b/src/coreComponents/schema/docs/ConstantThermalConductivity.rst index f3c7c40593e..a3d9cd2129f 100644 --- a/src/coreComponents/schema/docs/ConstantThermalConductivity.rst +++ b/src/coreComponents/schema/docs/ConstantThermalConductivity.rst @@ -1,11 +1,11 @@ -============================= ============ ======== ============================================================================= -Name Type Default Description -============================= ============ ======== ============================================================================= -name string required A name is required for any non-unique nodes -phaseNames string_array required List of fluid phases -thermalConductivityComponents R1Tensor required xx, yy, and zz components of a diagonal thermal conductivity tensor [W/(m.K)] -============================= ============ ======== ============================================================================= +============================= ============ ======== =============================================================================== +Name Type Default Description +============================= ============ ======== =============================================================================== +name string required A name is required for any non-unique nodes +phaseNames string_array required List of fluid phases +thermalConductivityComponents R1Tensor required xx, yy, and zz components of a diagonal thermal conductivity tensor [J/(s.m.K)] +============================= ============ ======== =============================================================================== diff --git a/src/coreComponents/schema/docs/Constitutive.rst b/src/coreComponents/schema/docs/Constitutive.rst index 6fa3cfd19db..cee410140ac 100644 --- a/src/coreComponents/schema/docs/Constitutive.rst +++ b/src/coreComponents/schema/docs/Constitutive.rst @@ -51,6 +51,7 @@ ProppantPorosity node :ref:`XML_ProppantPoros ProppantSlurryFluid node :ref:`XML_ProppantSlurryFluid` ProppantSolidProppantPermeability node :ref:`XML_ProppantSolidProppantPermeability` SlipDependentPermeability node :ref:`XML_SlipDependentPermeability` +SolidInternalEnergy node :ref:`XML_SolidInternalEnergy` TableCapillaryPressure node :ref:`XML_TableCapillaryPressure` TableRelativePermeability node :ref:`XML_TableRelativePermeability` TableRelativePermeabilityHysteresis node :ref:`XML_TableRelativePermeabilityHysteresis` diff --git a/src/coreComponents/schema/docs/Constitutive_other.rst b/src/coreComponents/schema/docs/Constitutive_other.rst index 2ca2f51060a..68a9ddde2d0 100644 --- a/src/coreComponents/schema/docs/Constitutive_other.rst +++ b/src/coreComponents/schema/docs/Constitutive_other.rst @@ -51,6 +51,7 @@ ProppantPorosity node :ref:`DATASTRUCTURE_ProppantPor ProppantSlurryFluid node :ref:`DATASTRUCTURE_ProppantSlurryFluid` ProppantSolidProppantPermeability node :ref:`DATASTRUCTURE_ProppantSolidProppantPermeability` SlipDependentPermeability node :ref:`DATASTRUCTURE_SlipDependentPermeability` +SolidInternalEnergy node :ref:`DATASTRUCTURE_SolidInternalEnergy` TableCapillaryPressure node :ref:`DATASTRUCTURE_TableCapillaryPressure` TableRelativePermeability node :ref:`DATASTRUCTURE_TableRelativePermeability` TableRelativePermeabilityHysteresis node :ref:`DATASTRUCTURE_TableRelativePermeabilityHysteresis` diff --git a/src/coreComponents/schema/docs/DeadOilFluid_other.rst b/src/coreComponents/schema/docs/DeadOilFluid_other.rst index 2b44e2d2747..083e50ec178 100644 --- a/src/coreComponents/schema/docs/DeadOilFluid_other.rst +++ b/src/coreComponents/schema/docs/DeadOilFluid_other.rst @@ -7,7 +7,7 @@ dPhaseCompFraction LvArray_Array< double, 5, camp_int_seq< long, 0l dPhaseDensity real64_array4d Derivative of phase density with respect to pressure, temperature, and global component fractions dPhaseEnthalpy real64_array4d Derivative of phase enthalpy with respect to pressure, temperature, and global component fractions dPhaseFraction real64_array4d Derivative of phase fraction with respect to pressure, temperature, and global component fractions -dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fraction +dPhaseInternalEnergy real64_array4d Derivative of phase internal energy with respect to pressure, temperature, and global component fractions dPhaseMassDensity real64_array4d Derivative of phase mass density with respect to pressure, temperature, and global component fractions dPhaseViscosity real64_array4d Derivative of phase viscosity with respect to pressure, temperature, and global component fractions dTotalDensity real64_array3d Derivative of total density with respect to pressure, temperature, and global component fractions diff --git a/src/coreComponents/schema/docs/EmbeddedSurfaceGenerator_other.rst b/src/coreComponents/schema/docs/EmbeddedSurfaceGenerator_other.rst index 8bd71bf6b41..b53cb9d3cce 100644 --- a/src/coreComponents/schema/docs/EmbeddedSurfaceGenerator_other.rst +++ b/src/coreComponents/schema/docs/EmbeddedSurfaceGenerator_other.rst @@ -1,14 +1,14 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================ ======================================================================================================================================================================================================================================================================================================================== -Name Type Registered On Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================ ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -parentEdgeIndex integer_array :ref:`DATASTRUCTURE_embeddedSurfacesNodeManager` Index of parent edge within the mesh object it is registered on. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================ ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ================================================ ======================================================================================================================================================================================================================================================================================================================== +Name Type Registered On Description +========================= =================================================================================================================================================== ================================================ ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +parentEdgeIndex integer_array :ref:`DATASTRUCTURE_embeddedSurfacesNodeManager` Index of parent edge within the mesh object it is registered on. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================ ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/FlowProppantTransport_other.rst b/src/coreComponents/schema/docs/FlowProppantTransport_other.rst index aa08bcefbb0..7acea5a6dbf 100644 --- a/src/coreComponents/schema/docs/FlowProppantTransport_other.rst +++ b/src/coreComponents/schema/docs/FlowProppantTransport_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/Hydrofracture_other.rst b/src/coreComponents/schema/docs/Hydrofracture_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/Hydrofracture_other.rst +++ b/src/coreComponents/schema/docs/Hydrofracture_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/LagrangianContact_other.rst b/src/coreComponents/schema/docs/LagrangianContact_other.rst index aa08bcefbb0..7acea5a6dbf 100644 --- a/src/coreComponents/schema/docs/LagrangianContact_other.rst +++ b/src/coreComponents/schema/docs/LagrangianContact_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/LaplaceFEM_other.rst b/src/coreComponents/schema/docs/LaplaceFEM_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/LaplaceFEM_other.rst +++ b/src/coreComponents/schema/docs/LaplaceFEM_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/MultiphasePoromechanics_other.rst b/src/coreComponents/schema/docs/MultiphasePoromechanics_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/MultiphasePoromechanics_other.rst +++ b/src/coreComponents/schema/docs/MultiphasePoromechanics_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/PhaseFieldDamageFEM_other.rst b/src/coreComponents/schema/docs/PhaseFieldDamageFEM_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/PhaseFieldDamageFEM_other.rst +++ b/src/coreComponents/schema/docs/PhaseFieldDamageFEM_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/PhaseFieldFracture_other.rst b/src/coreComponents/schema/docs/PhaseFieldFracture_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/PhaseFieldFracture_other.rst +++ b/src/coreComponents/schema/docs/PhaseFieldFracture_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/ProppantTransport_other.rst b/src/coreComponents/schema/docs/ProppantTransport_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/ProppantTransport_other.rst +++ b/src/coreComponents/schema/docs/ProppantTransport_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/SinglePhaseFVM_other.rst b/src/coreComponents/schema/docs/SinglePhaseFVM_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/SinglePhaseFVM_other.rst +++ b/src/coreComponents/schema/docs/SinglePhaseFVM_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/SinglePhaseHybridFVM_other.rst b/src/coreComponents/schema/docs/SinglePhaseHybridFVM_other.rst index 7cd3eeef08a..56309599fa0 100644 --- a/src/coreComponents/schema/docs/SinglePhaseHybridFVM_other.rst +++ b/src/coreComponents/schema/docs/SinglePhaseHybridFVM_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================ ================================================================ -Name Type Registered On Description -========================= =============================================================================================================================================================================================================================================================================================== ================================ ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -deltaFacePressure real64_array :ref:`DATASTRUCTURE_FaceManager` Accumulated face pressure updates -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================ ================================================================ +========================= =================================================================================================================================================== ================================ ================================================================ +Name Type Registered On Description +========================= =================================================================================================================================================== ================================ ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +deltaFacePressure real64_array :ref:`DATASTRUCTURE_FaceManager` Accumulated face pressure updates +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================ ================================================================ diff --git a/src/coreComponents/schema/docs/SinglePhasePoromechanicsEmbeddedFractures_other.rst b/src/coreComponents/schema/docs/SinglePhasePoromechanicsEmbeddedFractures_other.rst index aa08bcefbb0..7acea5a6dbf 100644 --- a/src/coreComponents/schema/docs/SinglePhasePoromechanicsEmbeddedFractures_other.rst +++ b/src/coreComponents/schema/docs/SinglePhasePoromechanicsEmbeddedFractures_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/SinglePhasePoromechanics_other.rst b/src/coreComponents/schema/docs/SinglePhasePoromechanics_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/SinglePhasePoromechanics_other.rst +++ b/src/coreComponents/schema/docs/SinglePhasePoromechanics_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/SinglePhaseProppantFVM_other.rst b/src/coreComponents/schema/docs/SinglePhaseProppantFVM_other.rst index c874c188464..5b4f62fa9a3 100644 --- a/src/coreComponents/schema/docs/SinglePhaseProppantFVM_other.rst +++ b/src/coreComponents/schema/docs/SinglePhaseProppantFVM_other.rst @@ -1,12 +1,12 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/SinglePhaseReservoir_other.rst b/src/coreComponents/schema/docs/SinglePhaseReservoir_other.rst index aa08bcefbb0..7acea5a6dbf 100644 --- a/src/coreComponents/schema/docs/SinglePhaseReservoir_other.rst +++ b/src/coreComponents/schema/docs/SinglePhaseReservoir_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/SinglePhaseWell_other.rst b/src/coreComponents/schema/docs/SinglePhaseWell_other.rst index b3804961a57..48a09333642 100644 --- a/src/coreComponents/schema/docs/SinglePhaseWell_other.rst +++ b/src/coreComponents/schema/docs/SinglePhaseWell_other.rst @@ -1,14 +1,14 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -WellControls node :ref:`DATASTRUCTURE_WellControls` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +WellControls node :ref:`DATASTRUCTURE_WellControls` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/SolidInternalEnergy.rst b/src/coreComponents/schema/docs/SolidInternalEnergy.rst new file mode 100644 index 00000000000..32cad1e4e49 --- /dev/null +++ b/src/coreComponents/schema/docs/SolidInternalEnergy.rst @@ -0,0 +1,12 @@ + + +======================= ====== ======== =================================================== +Name Type Default Description +======================= ====== ======== =================================================== +name string required A name is required for any non-unique nodes +referenceInternalEnergy real64 required Internal energy at the reference temperature [J/kg] +referenceTemperature real64 required Reference temperature [K] +volumetricHeatCapacity real64 required Solid volumetric heat capacity [J/(kg.K)] +======================= ====== ======== =================================================== + + diff --git a/src/coreComponents/schema/docs/SolidInternalEnergy_other.rst b/src/coreComponents/schema/docs/SolidInternalEnergy_other.rst new file mode 100644 index 00000000000..e245ee67a7a --- /dev/null +++ b/src/coreComponents/schema/docs/SolidInternalEnergy_other.rst @@ -0,0 +1,11 @@ + + +============================ ============== ============================================================================== +Name Type Description +============================ ============== ============================================================================== +dInternalEnergy_dTemperature real64_array2d Derivative of the solid internal energy w.r.t. temperature [J/(m^3.K)] +internalEnergy real64_array2d Internal energy of the solid per unit volume [J/m^3] +oldInternalEnergy real64_array2d Internal energy of the solid per unit volume at the previous time-step [J/m^3] +============================ ============== ============================================================================== + + diff --git a/src/coreComponents/schema/docs/SolidMechanicsEmbeddedFractures_other.rst b/src/coreComponents/schema/docs/SolidMechanicsEmbeddedFractures_other.rst index aa08bcefbb0..7acea5a6dbf 100644 --- a/src/coreComponents/schema/docs/SolidMechanicsEmbeddedFractures_other.rst +++ b/src/coreComponents/schema/docs/SolidMechanicsEmbeddedFractures_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/SolidMechanicsLagrangianSSLE_other.rst b/src/coreComponents/schema/docs/SolidMechanicsLagrangianSSLE_other.rst index d6506d6fedf..8ab18aa4617 100644 --- a/src/coreComponents/schema/docs/SolidMechanicsLagrangianSSLE_other.rst +++ b/src/coreComponents/schema/docs/SolidMechanicsLagrangianSSLE_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxForce real64 The maximum force contribution in the problem domain. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxForce real64 The maximum force contribution in the problem domain. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/SolidMechanics_LagrangianFEM_other.rst b/src/coreComponents/schema/docs/SolidMechanics_LagrangianFEM_other.rst index d6506d6fedf..8ab18aa4617 100644 --- a/src/coreComponents/schema/docs/SolidMechanics_LagrangianFEM_other.rst +++ b/src/coreComponents/schema/docs/SolidMechanics_LagrangianFEM_other.rst @@ -1,13 +1,13 @@ -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ -maxForce real64 The maximum force contribution in the problem domain. -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ================================================================ +========================= =================================================================================================================================================== ================================================================ +Name Type Description +========================= =================================================================================================================================================== ================================================================ +maxForce real64 The maximum force contribution in the problem domain. +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ================================================================ diff --git a/src/coreComponents/schema/docs/SurfaceGenerator_other.rst b/src/coreComponents/schema/docs/SurfaceGenerator_other.rst index 37b9b97583d..96cfc6ecec7 100644 --- a/src/coreComponents/schema/docs/SurfaceGenerator_other.rst +++ b/src/coreComponents/schema/docs/SurfaceGenerator_other.rst @@ -1,18 +1,18 @@ -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -Name Type Description -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== -discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. -failCriterion integer (no description available) -maxStableDt real64 Value of the Maximum Stable Timestep for this solver. -meshTargets geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. -tipEdges LvArray_SortedArray< int, int, LvArray_ChaiBuffer > Set containing all the tip edges -tipFaces LvArray_SortedArray< int, int, LvArray_ChaiBuffer > Set containing all the tip faces -tipNodes LvArray_SortedArray< int, int, LvArray_ChaiBuffer > Set containing all the nodes at the fracture tip -trailingFaces LvArray_SortedArray< int, int, LvArray_ChaiBuffer > Set containing all the trailing faces -LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` -NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` -========================= =============================================================================================================================================================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +Name Type Description +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== +discretization string Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. +failCriterion integer (no description available) +maxStableDt real64 Value of the Maximum Stable Timestep for this solver. +meshTargets geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > MeshBody/Region combinations that the solver will be applied to. +tipEdges LvArray_SortedArray< int, int, LvArray_ChaiBuffer > Set containing all the tip edges +tipFaces LvArray_SortedArray< int, int, LvArray_ChaiBuffer > Set containing all the tip faces +tipNodes LvArray_SortedArray< int, int, LvArray_ChaiBuffer > Set containing all the nodes at the fracture tip +trailingFaces LvArray_SortedArray< int, int, LvArray_ChaiBuffer > Set containing all the trailing faces +LinearSolverParameters node :ref:`DATASTRUCTURE_LinearSolverParameters` +NonlinearSolverParameters node :ref:`DATASTRUCTURE_NonlinearSolverParameters` +========================= =================================================================================================================================================== ======================================================================================================================================================================================================================================================================================================================== diff --git a/src/coreComponents/schema/docs/TwoPointFluxApproximation_other.rst b/src/coreComponents/schema/docs/TwoPointFluxApproximation_other.rst index b1b2c45ed80..825f9d230d6 100644 --- a/src/coreComponents/schema/docs/TwoPointFluxApproximation_other.rst +++ b/src/coreComponents/schema/docs/TwoPointFluxApproximation_other.rst @@ -1,15 +1,15 @@ -======================== =============================================================================================================================================================================================================================================================================================== ======================================== -Name Type Description -======================== =============================================================================================================================================================================================================================================================================================== ======================================== -cellStencil geosx_CellElementStencilTPFA (no description available) -coefficientName string Name of coefficient field -edfmStencil geosx_EmbeddedSurfaceToCellStencil (no description available) -faceElementToCellStencil geosx_FaceElementToCellStencil (no description available) -fieldName string Name of primary solution field -fractureStencil geosx_SurfaceElementStencil (no description available) -targetRegions geosx_mapBase< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, LvArray_Array< std___cxx11_basic_string< char, std_char_traits< char >, std_allocator< char > >, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > List of regions to build the stencil for -======================== =============================================================================================================================================================================================================================================================================================== ======================================== +======================== =================================================================================================================================================== ======================================== +Name Type Description +======================== =================================================================================================================================================== ======================================== +cellStencil geosx_CellElementStencilTPFA (no description available) +coefficientName string Name of coefficient field +edfmStencil geosx_EmbeddedSurfaceToCellStencil (no description available) +faceElementToCellStencil geosx_FaceElementToCellStencil (no description available) +fieldName string Name of primary solution field +fractureStencil geosx_SurfaceElementStencil (no description available) +targetRegions geosx_mapBase< std_string, LvArray_Array< std_string, 1, camp_int_seq< long, 0l >, int, LvArray_ChaiBuffer >, std_integral_constant< bool, true > > List of regions to build the stencil for +======================== =================================================================================================================================================== ======================================== diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index bef1c12123f..e620ba20ce5 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -1157,6 +1157,8 @@ the relative residual norm satisfies: + + @@ -1187,6 +1189,8 @@ the relative residual norm satisfies: + + @@ -1929,6 +1933,7 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> + @@ -2197,7 +2202,7 @@ The expected format is "{ waterMax, oilMax }", in that order--> - + @@ -2736,6 +2741,16 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia + + + + + + + + + + - + @@ -480,7 +480,7 @@ - + @@ -504,7 +504,7 @@ - + @@ -514,7 +514,7 @@ - + @@ -526,7 +526,7 @@ - + @@ -539,7 +539,7 @@ - + @@ -552,7 +552,7 @@ - + @@ -564,7 +564,7 @@ - + @@ -574,7 +574,7 @@ - + @@ -586,7 +586,7 @@ - + @@ -596,7 +596,7 @@ - + @@ -606,7 +606,7 @@ - + @@ -616,7 +616,7 @@ - + @@ -626,7 +626,7 @@ - + @@ -636,7 +636,7 @@ - + @@ -646,7 +646,7 @@ - + @@ -656,7 +656,7 @@ - + @@ -666,7 +666,7 @@ - + @@ -678,7 +678,7 @@ - + @@ -688,7 +688,7 @@ - + @@ -700,7 +700,7 @@ - + @@ -713,7 +713,7 @@ - + @@ -725,7 +725,7 @@ - + @@ -737,7 +737,7 @@ - + @@ -749,7 +749,7 @@ - + @@ -763,7 +763,7 @@ - + @@ -869,6 +869,7 @@ + @@ -902,7 +903,7 @@ - + @@ -996,7 +997,7 @@ - + @@ -1044,7 +1045,7 @@ - + @@ -1092,7 +1093,7 @@ - + @@ -1140,7 +1141,7 @@ - + @@ -1196,7 +1197,7 @@ - + @@ -1324,7 +1325,7 @@ - + @@ -1675,6 +1676,14 @@ + + + + + + + + diff --git a/src/coreComponents/unitTests/fluidFlowTests/CMakeLists.txt b/src/coreComponents/unitTests/fluidFlowTests/CMakeLists.txt index 6b53aa2b005..c7a67721aab 100644 --- a/src/coreComponents/unitTests/fluidFlowTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/fluidFlowTests/CMakeLists.txt @@ -6,6 +6,7 @@ set( gtest_geosx_tests testSinglePhaseBaseKernels.cpp testSinglePhaseFVMKernels.cpp testSinglePhaseHybridFVMKernels.cpp + testThermalCompMultiphaseFlow.cpp ) set( dependencyList gtest ) diff --git a/src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp b/src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp index 52dcf0665f0..6b48b5cf2db 100644 --- a/src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp +++ b/src/coreComponents/unitTests/fluidFlowTests/testCompFlowUtils.hpp @@ -780,7 +780,6 @@ void fillCellCenteredNumericalJacobian( COMPOSITIONAL_SOLVER & solver, elemDofNumber[ei] + numComp + 1, dT, jacobianFD.toViewConstSizes() ); - } } diff --git a/src/coreComponents/unitTests/fluidFlowTests/testThermalCompMultiphaseFlow.cpp b/src/coreComponents/unitTests/fluidFlowTests/testThermalCompMultiphaseFlow.cpp new file mode 100644 index 00000000000..4d2f2a6b86e --- /dev/null +++ b/src/coreComponents/unitTests/fluidFlowTests/testThermalCompMultiphaseFlow.cpp @@ -0,0 +1,421 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "constitutive/fluid/MultiFluidBase.hpp" +#include "finiteVolume/FiniteVolumeManager.hpp" +#include "finiteVolume/FluxApproximationBase.hpp" +#include "mainInterface/initialization.hpp" +#include "mainInterface/GeosxState.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseExtrinsicData.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseExtrinsicData.hpp" +#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp" + +using namespace geosx; +using namespace geosx::dataRepository; +using namespace geosx::constitutive; +using namespace geosx::constitutive::multifluid; +using namespace geosx::testing; + +CommandLineOptions g_commandLineOptions; + +// Sphinx start after input XML + +char const * pvtLiquid = "DensityFun PhillipsBrineDensity 1e6 7.5e7 5e5 295.15 370.15 25 0\n" + "ViscosityFun PhillipsBrineViscosity 0\n" + "EnthalpyFun BrineEnthalpy 1e6 7.5e7 5e5 295.15 370.15 25 0\n"; + +char const * pvtGas = "DensityFun SpanWagnerCO2Density 1e6 7.5e7 5e5 295.15 370.15 25\n" + "ViscosityFun FenghourCO2Viscosity 1e6 7.5e7 5e5 295.15 370.15 25\n" + "EnthalpyFun CO2Enthalpy 1e6 7.5e7 5e5 295.15 370.15 25\n"; + +char const * co2flash = "FlashModel CO2Solubility 1e6 7.5e7 5e5 295.15 370.15 25 0"; + +char const * xmlInput = + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n" + "\n"; + +// Sphinx end before input XML + +template< typename LAMBDA > +void testNumericalJacobian( CompositionalMultiphaseFVM & solver, + DomainPartition & domain, + real64 const perturbParameter, + real64 const relTol, + LAMBDA assembleFunction ) +{ + CRSMatrix< real64, globalIndex > const & jacobian = solver.getLocalMatrix(); + array1d< real64 > residual( jacobian.numRows() ); + + // assemble the analytical residual + solver.resetStateToBeginningOfStep( domain ); + + residual.zero(); + jacobian.zero(); + + assembleFunction( jacobian.toViewConstSizes(), residual.toView() ); + residual.move( LvArray::MemorySpace::host, false ); + + // copy the analytical residual + array1d< real64 > residualOrig( residual ); + + // create the numerical jacobian + jacobian.move( LvArray::MemorySpace::host ); + CRSMatrix< real64, globalIndex > jacobianFD( jacobian ); + jacobianFD.zero(); + + // fill jacobian FD + fillCellCenteredNumericalJacobian( solver, + domain, + true, + perturbParameter, + residual.toView(), + residualOrig.toView(), + jacobian.toView(), + jacobianFD.toView(), + assembleFunction ); + + // assemble the analytical jacobian + solver.resetStateToBeginningOfStep( domain ); + + residual.zero(); + jacobian.zero(); + assembleFunction( jacobian.toViewConstSizes(), residual.toView() ); + + compareLocalMatrices( jacobian.toViewConst(), jacobianFD.toViewConst(), relTol, 1e-6 ); +} + +void writeTableToFile( string const & filename, char const * str ) +{ + std::ofstream os( filename ); + ASSERT_TRUE( os.is_open() ); + os << str; + os.close(); +} + +void removeFile( string const & filename ) +{ + int const ret = std::remove( filename.c_str() ); + ASSERT_TRUE( ret == 0 ); +} + +class ThermalCompositionalMultiphaseFlowTest : public ::testing::Test +{ +public: + + ThermalCompositionalMultiphaseFlowTest(): + state( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ) + { + writeTableToFile( pvtLiquidFilename, pvtLiquid ); + writeTableToFile( pvtGasFilename, pvtGas ); + writeTableToFile( co2flashFilename, co2flash ); + } + + ~ThermalCompositionalMultiphaseFlowTest() override + { + removeFile( pvtLiquidFilename ); + removeFile( pvtGasFilename ); + removeFile( co2flashFilename ); + } + +protected: + + void SetUp() override + { + + setupProblemFromXML( state.getProblemManager(), xmlInput ); + solver = &state.getProblemManager().getPhysicsSolverManager().getGroup< CompositionalMultiphaseFVM >( "compflow" ); + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + solver->setupSystem( domain, + solver->getDofManager(), + solver->getLocalMatrix(), + solver->getSystemRhs(), + solver->getSystemSolution() ); + + solver->implicitStepSetup( time, dt, domain ); + } + + static real64 constexpr time = 0.0; + static real64 constexpr dt = 1e4; + static real64 constexpr eps = std::numeric_limits< real64 >::epsilon(); + + GeosxState state; + CompositionalMultiphaseFVM * solver; + + string const pvtLiquidFilename = "pvtliquid.txt"; + string const pvtGasFilename = "pvtgas.txt"; + string const co2flashFilename = "co2flash.txt"; +}; + +real64 constexpr ThermalCompositionalMultiphaseFlowTest::time; +real64 constexpr ThermalCompositionalMultiphaseFlowTest::dt; +real64 constexpr ThermalCompositionalMultiphaseFlowTest::eps; + +TEST_F( ThermalCompositionalMultiphaseFlowTest, derivativeNumericalCheck_composition ) +{ + real64 const perturb = std::sqrt( eps ); + real64 const tol = 1e-4; + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + testCompositionNumericalDerivatives( *solver, domain, perturb, tol ); +} + +TEST_F( ThermalCompositionalMultiphaseFlowTest, derivativeNumericalCheck_phaseVolumeFraction ) +{ + real64 const perturb = std::sqrt( eps ); + real64 const tol = 5e-2; // 5% error margin + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + testPhaseVolumeFractionNumericalDerivatives( *solver, domain, true, perturb, tol ); +} + +TEST_F( ThermalCompositionalMultiphaseFlowTest, derivativeNumericalCheck_phaseMobility ) +{ + real64 const perturb = std::sqrt( eps ); + real64 const tol = 5e-2; // 5% error margin + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + testPhaseMobilityNumericalDerivatives( *solver, domain, true, perturb, tol ); +} + +TEST_F( ThermalCompositionalMultiphaseFlowTest, jacobianNumericalCheck_flux ) +{ + real64 const perturb = std::sqrt( eps ); + real64 const tol = 1e-1; // 10% error margin + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + testNumericalJacobian( *solver, domain, perturb, tol, + [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + solver->assembleFluxTerms( dt, domain, solver->getDofManager(), localMatrix, localRhs ); + } ); +} + +#if 0 +TEST_F( ThermalCompositionalMultiphaseFlowTest, jacobianNumericalCheck_accumulationVolumeBalance ) +{ + real64 const perturb = sqrt( eps ); + real64 const tol = 1e-1; // 10% error margin + + DomainPartition & domain = state.getProblemManager().getDomainPartition(); + + testNumericalJacobian( *solver, domain, perturb, tol, + [&] ( CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) + { + solver->assembleAccumulationAndVolumeBalanceTerms( domain, solver->getDofManager(), localMatrix, localRhs ); + } ); +} +#endif + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geosx::basicSetup( argc, argv ); + int const result = RUN_ALL_TESTS(); + geosx::basicCleanup(); + return result; +} diff --git a/src/docs/sphinx/CompleteXMLSchema.rst b/src/docs/sphinx/CompleteXMLSchema.rst index 130b9440a03..d991f080a00 100644 --- a/src/docs/sphinx/CompleteXMLSchema.rst +++ b/src/docs/sphinx/CompleteXMLSchema.rst @@ -813,6 +813,13 @@ Element: SlipDependentPermeability .. include:: ../../coreComponents/schema/docs/SlipDependentPermeability.rst +.. _XML_SolidInternalEnergy: + +Element: SolidInternalEnergy +============================ +.. include:: ../../coreComponents/schema/docs/SolidInternalEnergy.rst + + .. _XML_SolidMechanicsEmbeddedFractures: Element: SolidMechanicsEmbeddedFractures @@ -1846,6 +1853,13 @@ Datastructure: SlipDependentPermeability .. include:: ../../coreComponents/schema/docs/SlipDependentPermeability_other.rst +.. _DATASTRUCTURE_SolidInternalEnergy: + +Datastructure: SolidInternalEnergy +================================== +.. include:: ../../coreComponents/schema/docs/SolidInternalEnergy_other.rst + + .. _DATASTRUCTURE_SolidMechanicsEmbeddedFractures: Datastructure: SolidMechanicsEmbeddedFractures