diff --git a/inputFiles/singlePhaseFlow/ReactiveCompressible_1d.xml b/inputFiles/singlePhaseFlow/ReactiveCompressible_1d.xml
new file mode 100644
index 00000000000..a9fb09731d4
--- /dev/null
+++ b/inputFiles/singlePhaseFlow/ReactiveCompressible_1d.xml
@@ -0,0 +1,283 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt
index a65eb5eec19..49157dc8f8c 100644
--- a/src/coreComponents/constitutive/CMakeLists.txt
+++ b/src/coreComponents/constitutive/CMakeLists.txt
@@ -135,6 +135,9 @@ set( constitutive_headers
fluid/singlefluid/SingleFluidSelector.hpp
fluid/singlefluid/SlurryFluidSelector.hpp
fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
+ fluid/singlefluid/reactive/ReactiveCompressibleSinglePhaseFluid.hpp
+ fluid/singlefluid/reactive/ReactiveSingleFluid.hpp
+ fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.hpp
permeability/CarmanKozenyPermeability.hpp
permeability/ConstantPermeability.hpp
permeability/DamagePermeability.hpp
@@ -279,6 +282,9 @@ set( constitutive_sources
fluid/singlefluid/SingleFluidBase.cpp
fluid/singlefluid/SlurryFluidBase.cpp
fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp
+ fluid/singlefluid/reactive/ReactiveCompressibleSinglePhaseFluid.cpp
+ fluid/singlefluid/reactive/ReactiveSingleFluid.cpp
+ fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.cpp
permeability/CarmanKozenyPermeability.cpp
permeability/ConstantPermeability.cpp
permeability/DamagePermeability.cpp
diff --git a/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveFluidSelector.hpp b/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveFluidSelector.hpp
index 03439e9629b..71b42c912f0 100644
--- a/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveFluidSelector.hpp
+++ b/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveFluidSelector.hpp
@@ -21,6 +21,8 @@
#include "constitutive/ConstitutivePassThruHandler.hpp"
#include "constitutive/fluid/multifluid/reactive/ReactiveBrineFluid.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ReactiveCompressibleSinglePhaseFluid.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.hpp"
#include "common/GeosxConfig.hpp"
@@ -46,6 +48,22 @@ void constitutiveUpdatePassThru( ReactiveMultiFluid & fluid,
ReactiveBrineThermal >::execute( fluid, std::forward< LAMBDA >( lambda ) );
}
+template< typename LAMBDA >
+void constitutiveUpdatePassThru( ReactiveSingleFluid const & fluid,
+ LAMBDA && lambda )
+{
+ ConstitutivePassThruHandler< ReactiveCompressibleSinglePhase,
+ ThermalReactiveCompressibleSinglePhase >::execute( fluid, std::forward< LAMBDA >( lambda ) );
+}
+
+template< typename LAMBDA >
+void constitutiveUpdatePassThru( ReactiveSingleFluid & fluid,
+ LAMBDA && lambda )
+{
+ ConstitutivePassThruHandler< ReactiveCompressibleSinglePhase,
+ ThermalReactiveCompressibleSinglePhase >::execute( fluid, std::forward< LAMBDA >( lambda ) );
+}
+
} // namespace constitutive
} // namespace geos
diff --git a/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveMultiFluid.cpp b/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveMultiFluid.cpp
index 09ac848a16a..a4d9deb8885 100644
--- a/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveMultiFluid.cpp
+++ b/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveMultiFluid.cpp
@@ -39,7 +39,7 @@ ReactiveMultiFluid::
registerField( fields::reactivefluid::primarySpeciesConcentration{}, &m_primarySpeciesConcentration );
registerField( fields::reactivefluid::secondarySpeciesConcentration{}, &m_secondarySpeciesConcentration );
- registerField( fields::reactivefluid::primarySpeciesTotalConcentration{}, &m_primarySpeciesTotalConcentration );
+ registerField( fields::reactivefluid::primarySpeciesAggregateConcentration{}, &m_primarySpeciesTotalConcentration );
registerField( fields::reactivefluid::kineticReactionRates{}, &m_kineticReactionRates );
}
diff --git a/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveMultiFluidFields.hpp b/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveMultiFluidFields.hpp
index d3a3a77f967..e10a004bdde 100644
--- a/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveMultiFluidFields.hpp
+++ b/src/coreComponents/constitutive/fluid/multifluid/reactive/ReactiveMultiFluidFields.hpp
@@ -33,6 +33,7 @@ namespace reactivefluid
{
using array2dLayoutComp = array2d< real64, compflow::LAYOUT_COMP >;
+using array3dLayoutComp_dC = array3d< real64, compflow::LAYOUT_COMP_DC >;
DECLARE_FIELD( primarySpeciesConcentration,
"primarySpeciesConcentration",
@@ -42,13 +43,29 @@ DECLARE_FIELD( primarySpeciesConcentration,
WRITE_AND_READ,
"primarySpeciesConcentration" );
-DECLARE_FIELD( primarySpeciesTotalConcentration,
- "primarySpeciesTotalConcentration",
+DECLARE_FIELD( primarySpeciesAggregateConcentration,
+ "primarySpeciesAggregateConcentration",
array2dLayoutComp,
0,
LEVEL_0,
WRITE_AND_READ,
- "primarySpeciesTotalConcentration" );
+ "primarySpeciesAggregateConcentration" );
+
+DECLARE_FIELD( primarySpeciesAggregateConcentration_n,
+ "primarySpeciesAggregateConcentration_n",
+ array2dLayoutComp,
+ 0,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "primarySpeciesAggregateConcentration at the previous timestep" );
+
+DECLARE_FIELD( dPrimarySpeciesAggregateConcentration_dLogPrimaryConc,
+ "dPrimarySpeciesAggregateConcentration_dLogPrimaryConc",
+ array3dLayoutComp_dC,
+ 0,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "Deivatives of primarySpeciesAggregateConcentration w.r.t log primary species concentration" );
DECLARE_FIELD( secondarySpeciesConcentration,
"secondarySpeciesConcentration",
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp
index c86f9c25b19..ce243afd080 100644
--- a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp
+++ b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp
@@ -130,6 +130,21 @@ class CompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate
m_dViscosity[k][q][DerivOffset::dP] );
}
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const GEOS_UNUSED_PARAM( temperature ),
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & GEOS_UNUSED_PARAM( composition ) ) const override
+ {
+ compute( pressure,
+ m_density[k][q],
+ m_dDensity[k][q][DerivOffset::dP],
+ m_viscosity[k][q],
+ m_dViscosity[k][q][DerivOffset::dP] );
+ }
+
private:
/// Relationship between the fluid density and pressure
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp
index 9d8d3cbfb2d..af8c4e3fddf 100644
--- a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp
+++ b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp
@@ -20,6 +20,7 @@
#ifndef GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_SINGLEFLUIDBASE_HPP
#define GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_SINGLEFLUIDBASE_HPP
+#include "common/DataLayouts.hpp"
#include "constitutive/ConstitutiveBase.hpp"
#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp"
#include "constitutive/fluid/singlefluid/SingleFluidUtils.hpp"
@@ -205,6 +206,21 @@ class SingleFluidBaseUpdate
real64 const pressure,
real64 const temperature ) const = 0;
+ /**
+ * @brief Update fluid state at a single point.
+ * @param[in] k element index
+ * @param[in] q gauss point index
+ * @param[in] pressure the target pressure value
+ * @param[in] temperature the target temperature value
+ * @param[in] logPrimaryConc the target logPrimaryConc value
+ */
+ GEOS_HOST_DEVICE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const temperature,
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & logPrimaryConc ) const = 0;
+
};
//END_SPHINX_INCLUDE_02
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidSelector.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidSelector.hpp
index be6560c6c21..ba378a397e1 100644
--- a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidSelector.hpp
+++ b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidSelector.hpp
@@ -22,6 +22,8 @@
#include "constitutive/ConstitutivePassThruHandler.hpp"
#include "constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp"
#include "constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ReactiveCompressibleSinglePhaseFluid.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.hpp"
namespace geos
{
@@ -33,7 +35,9 @@ template< typename LAMBDA >
void constitutiveUpdatePassThru( SingleFluidBase const & fluid,
LAMBDA && lambda )
{
- ConstitutivePassThruHandler< ThermalCompressibleSinglePhaseFluid,
+ ConstitutivePassThruHandler< ThermalReactiveCompressibleSinglePhase,
+ ReactiveCompressibleSinglePhase,
+ ThermalCompressibleSinglePhaseFluid,
CompressibleSinglePhaseFluid >::execute( fluid, std::forward< LAMBDA >( lambda ) );
}
@@ -41,7 +45,9 @@ template< typename LAMBDA >
void constitutiveUpdatePassThru( SingleFluidBase & fluid,
LAMBDA && lambda )
{
- ConstitutivePassThruHandler< ThermalCompressibleSinglePhaseFluid,
+ ConstitutivePassThruHandler< ThermalReactiveCompressibleSinglePhase,
+ ReactiveCompressibleSinglePhase,
+ ThermalCompressibleSinglePhaseFluid,
CompressibleSinglePhaseFluid >::execute( fluid, std::forward< LAMBDA >( lambda ) );
}
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
index 2a97afae9e5..59ef2713a0b 100644
--- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
+++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp
@@ -165,6 +165,21 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate
m_dEnthalpy[k][q][DerivOffset::dT] );
}
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const GEOS_UNUSED_PARAM( temperature ),
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & GEOS_UNUSED_PARAM( composition ) ) const override
+ {
+ compute( pressure,
+ m_density[k][q],
+ m_dDensity[k][q][DerivOffset::dP],
+ m_viscosity[k][q],
+ m_dViscosity[k][q][DerivOffset::dP] );
+ }
+
private:
/// Fluid internal energy and derivatives
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveCompressibleSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveCompressibleSinglePhaseFluid.cpp
new file mode 100644
index 00000000000..c5af8c0fe0d
--- /dev/null
+++ b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveCompressibleSinglePhaseFluid.cpp
@@ -0,0 +1,163 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ReactiveCompressibleSinglePhaseFluid.cpp
+ */
+
+#include "ReactiveCompressibleSinglePhaseFluid.hpp"
+
+#include "constitutive/fluid/singlefluid/SingleFluidFields.hpp"
+
+namespace geos
+{
+
+using namespace dataRepository;
+
+namespace constitutive
+{
+
+ReactiveCompressibleSinglePhase::ReactiveCompressibleSinglePhase( string const & name, Group * const parent ):
+ ReactiveSingleFluid( name, parent ),
+ m_densityModelType( ExponentApproximationType::Linear ),
+ m_viscosityModelType( ExponentApproximationType::Linear )
+{
+ registerWrapper( viewKeyStruct::defaultDensityString(), &m_defaultDensity ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Default value for density." );
+
+ registerWrapper( viewKeyStruct::defaultViscosityString(), &m_defaultViscosity ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Default value for viscosity." );
+
+ registerWrapper( viewKeyStruct::compressibilityString(), &m_compressibility ).
+ setApplyDefaultValue( 0.0 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Fluid compressibility" );
+
+ registerWrapper( viewKeyStruct::viscosibilityString(), &m_viscosibility ).
+ setApplyDefaultValue( 0.0 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Fluid viscosity exponential coefficient" );
+
+ registerWrapper( viewKeyStruct::referencePressureString(), &m_referencePressure ).
+ setApplyDefaultValue( 0.0 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Reference pressure" );
+
+ registerWrapper( viewKeyStruct::referenceDensityString(), &m_referenceDensity ).
+ setApplyDefaultValue( 1000.0 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Reference fluid density" );
+
+ registerWrapper( viewKeyStruct::referenceViscosityString(), &m_referenceViscosity ).
+ setApplyDefaultValue( 0.001 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Reference fluid viscosity" );
+
+ registerWrapper( viewKeyStruct::densityModelTypeString(), &m_densityModelType ).
+ setApplyDefaultValue( m_densityModelType ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Type of density model. Valid options:\n* " + EnumStrings< ExponentApproximationType >::concat( "\n* " ) );
+
+ registerWrapper( viewKeyStruct::viscosityModelTypeString(), &m_viscosityModelType ).
+ setApplyDefaultValue( m_viscosityModelType ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Type of viscosity model. Valid options:\n* " + EnumStrings< ExponentApproximationType >::concat( "\n* " ) );
+
+}
+
+ReactiveCompressibleSinglePhase::~ReactiveCompressibleSinglePhase() = default;
+
+void ReactiveCompressibleSinglePhase::allocateConstitutiveData( dataRepository::Group & parent,
+ localIndex const numConstitutivePointsPerParentIndex )
+{
+ ReactiveSingleFluid::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex );
+
+ getField< fields::singlefluid::density >().setApplyDefaultValue( m_defaultDensity );
+ getField< fields::singlefluid::viscosity >().setApplyDefaultValue( m_defaultViscosity );
+
+ m_density.value.setValues< serialPolicy >( m_referenceDensity );
+ m_viscosity.value.setValues< serialPolicy >( m_referenceViscosity );
+}
+
+void ReactiveCompressibleSinglePhase::postInputInitialization()
+{
+ ReactiveSingleFluid::postInputInitialization();
+
+ auto const checkNonnegative = [&]( real64 const value, auto const & attribute )
+ {
+ GEOS_THROW_IF_LT_MSG( value, 0.0,
+ GEOS_FMT( "{}: invalid value of attribute '{}'", getFullName(), attribute ),
+ InputError );
+ };
+ checkNonnegative( m_compressibility, viewKeyStruct::compressibilityString() );
+ checkNonnegative( m_viscosibility, viewKeyStruct::viscosibilityString() );
+
+ auto const checkPositive = [&]( real64 const value, auto const & attribute )
+ {
+ GEOS_THROW_IF_LE_MSG( value, 0.0,
+ GEOS_FMT( "{}: invalid value of attribute '{}'", getFullName(), attribute ),
+ InputError );
+ };
+ checkPositive( m_referenceDensity, viewKeyStruct::referenceDensityString() );
+ checkPositive( m_referenceViscosity, viewKeyStruct::referenceViscosityString() );
+
+ // Due to the way update wrapper is currently implemented, we can only support one model type
+ auto const checkModelType = [&]( ExponentApproximationType const value, auto const & attribute )
+ {
+ GEOS_THROW_IF_NE_MSG( value, ExponentApproximationType::Linear,
+ GEOS_FMT( "{}: invalid model type in attribute '{}' (only linear currently supported)", getFullName(), attribute ),
+ InputError );
+ };
+ checkModelType( m_densityModelType, viewKeyStruct::densityModelTypeString() );
+ checkModelType( m_viscosityModelType, viewKeyStruct::viscosityModelTypeString() );
+
+ // Set default values for derivatives (cannot be done in base class)
+ // TODO: reconsider the necessity of this
+
+ real64 dRho_dP;
+ real64 dVisc_dP;
+ createKernelWrapper().compute( m_referencePressure, m_referenceDensity, dRho_dP, m_referenceViscosity, dVisc_dP );
+ for( integer i=0; i
+
+namespace geos
+{
+
+namespace constitutive
+{
+
+/**
+ * @brief Update class for the model suitable for lambda capture.
+ * @tparam DENS_EAT type of density exponent approximation
+ * @tparam VISC_EAT type of viscosity exponent approximation
+ */
+template< ExponentApproximationType DENS_EAT, ExponentApproximationType VISC_EAT >
+class ReactiveCompressibleSinglePhaseUpdate : public ReactiveSingleFluidUpdate
+{
+public:
+
+ using DensRelationType = ExponentialRelation< real64, DENS_EAT >;
+ using ViscRelationType = ExponentialRelation< real64, VISC_EAT >;
+ using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
+
+ ReactiveCompressibleSinglePhaseUpdate( DensRelationType const & densRelation,
+ ViscRelationType const & viscRelation,
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const & density,
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > const & dDensity,
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const & viscosity,
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > const & dViscosity,
+ integer const numPrimarySpecies,
+ // chemicalReactions::EquilibriumReactions const & equilibriumReactions,
+ // chemicalReactions::KineticReactions const & kineticReactions,
+ arrayView2d< real64, compflow::USD_COMP > const & primarySpeciesConcentration,
+ arrayView2d< real64, compflow::USD_COMP > const & secondarySpeciesConcentration,
+ arrayView2d< real64, compflow::USD_COMP > const & primarySpeciesAggregateConcentration,
+ arrayView3d< real64, compflow::USD_COMP_DC > const & dPrimarySpeciesAggregateConcentration_dLogPrimaryConc,
+ arrayView2d< real64, compflow::USD_COMP > const & kineticReactionRates )
+ : ReactiveSingleFluidUpdate( density, dDensity, viscosity, dViscosity, numPrimarySpecies,
+ // equilibriumReactions, kineticReactions,
+ primarySpeciesConcentration, secondarySpeciesConcentration, primarySpeciesAggregateConcentration,
+ dPrimarySpeciesAggregateConcentration_dLogPrimaryConc, kineticReactionRates ),
+ m_densRelation( densRelation ),
+ m_viscRelation( viscRelation )
+ {}
+
+ /// Default copy constructor
+ ReactiveCompressibleSinglePhaseUpdate( ReactiveCompressibleSinglePhaseUpdate const & ) = default;
+
+ /// Default move constructor
+ ReactiveCompressibleSinglePhaseUpdate( ReactiveCompressibleSinglePhaseUpdate && ) = default;
+
+ /// Deleted copy assignment operator
+ ReactiveCompressibleSinglePhaseUpdate & operator=( ReactiveCompressibleSinglePhaseUpdate const & ) = delete;
+
+ /// Deleted move assignment operator
+ ReactiveCompressibleSinglePhaseUpdate & operator=( ReactiveCompressibleSinglePhaseUpdate && ) = delete;
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void compute( real64 const pressure,
+ real64 & density,
+ real64 & dDensity_dPressure,
+ real64 & viscosity,
+ real64 & dViscosity_dPressure ) const override
+ {
+ m_densRelation.compute( pressure, density, dDensity_dPressure );
+ m_viscRelation.compute( pressure, viscosity, dViscosity_dPressure );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void compute( real64 const pressure,
+ real64 const GEOS_UNUSED_PARAM( temperature ),
+ real64 & density,
+ real64 & dDensity_dPressure,
+ real64 & GEOS_UNUSED_PARAM( dDensity_dTemperature ),
+ real64 & viscosity,
+ real64 & dViscosity_dPressure,
+ real64 & GEOS_UNUSED_PARAM( dViscosity_dTemperature ),
+ real64 & GEOS_UNUSED_PARAM( internalEnergy ),
+ real64 & GEOS_UNUSED_PARAM( dInternalEnergy_dPressure ),
+ real64 & GEOS_UNUSED_PARAM( dInternalEnergy_dTemperature ),
+ real64 & GEOS_UNUSED_PARAM( enthalpy ),
+ real64 & GEOS_UNUSED_PARAM( dEnthalpy_dPressure ),
+ real64 & GEOS_UNUSED_PARAM( dEnthalpy_dTemperature ) ) const override
+ {
+ m_densRelation.compute( pressure, density, dDensity_dPressure );
+ m_viscRelation.compute( pressure, viscosity, dViscosity_dPressure );
+ }
+
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure ) const override
+ {
+ compute( pressure,
+ m_density[k][q],
+ m_dDensity[k][q][DerivOffset::dP],
+ m_viscosity[k][q],
+ m_dViscosity[k][q][DerivOffset::dP] );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const GEOS_UNUSED_PARAM( temperature ) ) const override
+ {
+ compute( pressure,
+ m_density[k][q],
+ m_dDensity[k][q][DerivOffset::dP],
+ m_viscosity[k][q],
+ m_dViscosity[k][q][DerivOffset::dP] );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const GEOS_UNUSED_PARAM( temperature ),
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & GEOS_UNUSED_PARAM( logPrimaryConc ) ) const override
+ {
+ compute( pressure,
+ m_density[k][q],
+ m_dDensity[k][q][DerivOffset::dP],
+ m_viscosity[k][q],
+ m_dViscosity[k][q][DerivOffset::dP] );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void updateChemistry( localIndex const k,
+ localIndex const GEOS_UNUSED_PARAM( q ),
+ real64 const pressure,
+ real64 const temperature,
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & primaryConc ) const override
+ {
+ for( int i=0; i < m_numPrimarySpecies; i++ )
+ {
+ m_primarySpeciesAggregateConcentration[k][i] = primaryConc[i];
+ }
+
+ computeChemistry( pressure,
+ temperature,
+ m_primarySpeciesAggregateConcentration[k],
+ m_primarySpeciesConcentration[k],
+ m_secondarySpeciesConcentration[k],
+ m_kineticReactionRates[k] );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void updateChemistryLogConc( localIndex const k,
+ localIndex const GEOS_UNUSED_PARAM( q ),
+ real64 const GEOS_UNUSED_PARAM( pressure ),
+ real64 const GEOS_UNUSED_PARAM( temperature ),
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & logPrimaryConc ) const override
+ {
+ for( int i=0; i < m_numPrimarySpecies; i++ )
+ {
+ m_primarySpeciesConcentration[k][i] = std::exp( logPrimaryConc[i] );
+
+ m_primarySpeciesAggregateConcentration[k][i] = m_primarySpeciesConcentration[k][i];
+
+ m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc[k][i][i] = m_primarySpeciesConcentration[k][i];
+ }
+
+ // computeChemistry( pressure,
+ // temperature,
+ // m_primarySpeciesAggregateConcentration[k],
+ // m_primarySpeciesConcentration[k],
+ // m_secondarySpeciesConcentration[k],
+ // m_kineticReactionRates[k] );
+ }
+
+private:
+
+ /// Relationship between the fluid density and pressure
+ DensRelationType m_densRelation;
+
+ /// Relationship between the fluid viscosity and pressure
+ ViscRelationType m_viscRelation;
+
+};
+
+class ReactiveCompressibleSinglePhase : public ReactiveSingleFluid
+{
+public:
+ using DerivOffset = singlefluid::DerivativeOffset;
+ ReactiveCompressibleSinglePhase( string const & name, Group * const parent );
+
+ virtual ~ReactiveCompressibleSinglePhase() override;
+
+ static string catalogName() { return "ReactiveCompressibleSinglePhase"; }
+
+ virtual string getCatalogName() const override { return catalogName(); }
+
+ virtual void allocateConstitutiveData( dataRepository::Group & parent,
+ localIndex const numConstitutivePointsPerParentIndex ) override;
+
+ /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT, not just linear)
+ using KernelWrapper = ReactiveCompressibleSinglePhaseUpdate< ExponentApproximationType::Linear, ExponentApproximationType::Linear >;
+
+ /**
+ * @brief Create an update kernel wrapper.
+ * @return the wrapper
+ */
+ KernelWrapper createKernelWrapper();
+
+ struct viewKeyStruct
+ {
+ static constexpr char const * defaultDensityString() { return "defaultDensity"; }
+ static constexpr char const * defaultViscosityString() { return "defaultViscosity"; }
+ static constexpr char const * compressibilityString() { return "compressibility"; }
+ static constexpr char const * viscosibilityString() { return "viscosibility"; }
+ static constexpr char const * referencePressureString() { return "referencePressure"; }
+ static constexpr char const * referenceDensityString() { return "referenceDensity"; }
+ static constexpr char const * referenceViscosityString() { return "referenceViscosity"; }
+ static constexpr char const * densityModelTypeString() { return "densityModelType"; }
+ static constexpr char const * viscosityModelTypeString() { return "viscosityModelType"; }
+ };
+
+ real64 defaultDensity() const override final { return m_defaultDensity; }
+ real64 defaultViscosity() const override final { return m_defaultViscosity; }
+
+protected:
+
+ virtual void postInputInitialization() override;
+
+ /// default density value
+ real64 m_defaultDensity;
+
+ /// default viscosity value
+ real64 m_defaultViscosity;
+
+ /// scalar fluid bulk modulus parameter
+ real64 m_compressibility;
+
+ /// scalar fluid viscosity exponential coefficient
+ real64 m_viscosibility;
+
+ /// reference pressure parameter
+ real64 m_referencePressure;
+
+ /// reference density parameter
+ real64 m_referenceDensity;
+
+ /// reference viscosity parameter
+ real64 m_referenceViscosity;
+
+ /// type of density model in terms of pressure
+ ExponentApproximationType m_densityModelType;
+
+ /// type of viscosity model
+ ExponentApproximationType m_viscosityModelType;
+
+};
+
+} /* namespace constitutive */
+
+} /* namespace geos */
+
+#endif /* GEOS_CONSTITUTIVE_FLUID_REACTIVECOMPRESSIBLESINGLEPHASEFLUID_HPP_ */
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.cpp b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.cpp
new file mode 100644
index 00000000000..f406cf3e3a4
--- /dev/null
+++ b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.cpp
@@ -0,0 +1,105 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ReactiveSingleFluid.cpp
+ */
+#include "ReactiveSingleFluid.hpp"
+
+#include "constitutive/fluid/singlefluid/SingleFluidFields.hpp"
+#include "constitutive/fluid/multifluid/reactive/ReactiveMultiFluidFields.hpp"
+
+
+namespace geos
+{
+
+using namespace dataRepository;
+
+namespace constitutive
+{
+
+ReactiveSingleFluid::
+ ReactiveSingleFluid( string const & name, Group * const parent ):
+ SingleFluidBase( name, parent )
+{
+ registerWrapper( viewKeyStruct::primarySpeciesNamesString(), &m_primarySpeciesNames ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "List of primary species names" );
+
+ // For now this is being hardcoded. We will see where this should come from.
+ m_numPrimarySpecies = 3;
+ m_numSecondarySpecies = 11;
+ m_numKineticReactions = 2;
+
+ registerField( fields::reactivefluid::primarySpeciesConcentration{}, &m_primarySpeciesConcentration );
+ registerField( fields::reactivefluid::secondarySpeciesConcentration{}, &m_secondarySpeciesConcentration );
+ registerField( fields::reactivefluid::primarySpeciesAggregateConcentration{}, &m_primarySpeciesAggregateConcentration );
+ registerField( fields::reactivefluid::primarySpeciesAggregateConcentration_n{}, &m_primarySpeciesAggregateConcentration_n );
+ registerField( fields::reactivefluid::dPrimarySpeciesAggregateConcentration_dLogPrimaryConc{}, &m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc );
+ registerField( fields::reactivefluid::kineticReactionRates{}, &m_kineticReactionRates );
+}
+
+void ReactiveSingleFluid::postInputInitialization()
+{
+ SingleFluidBase::postInputInitialization();
+
+ // // call to correctly set member array tertiary sizes on the 'main' material object
+ // resizeFields( 0, 0 );
+
+ // createChemicalReactions();
+}
+
+void ReactiveSingleFluid::saveConvergedState() const
+{
+ SingleFluidBase::saveConvergedState();
+
+ m_primarySpeciesAggregateConcentration_n.setValues< parallelDevicePolicy<> >( m_primarySpeciesAggregateConcentration.toViewConst() );
+}
+
+void ReactiveSingleFluid::resizeFields( localIndex const size, localIndex const numPts )
+{
+ GEOS_UNUSED_VAR( numPts );
+
+ integer const numPrimarySpecies = this->numPrimarySpecies();
+ integer const numSecondarySpecies = this->numSecondarySpecies();
+ integer const numKineticReactions = this->numKineticReactions();
+
+ m_primarySpeciesConcentration.resize( size, numPrimarySpecies );
+ m_secondarySpeciesConcentration.resize( size, numSecondarySpecies );
+ m_primarySpeciesAggregateConcentration.resize( size, numPrimarySpecies );
+ m_primarySpeciesAggregateConcentration_n.resize( size, numPrimarySpecies );
+ m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc.resize( size, numPrimarySpecies, numPrimarySpecies );
+ m_kineticReactionRates.resize( size, numKineticReactions );
+}
+
+void ReactiveSingleFluid::allocateConstitutiveData( dataRepository::Group & parent,
+ localIndex const numConstitutivePointsPerParentIndex )
+{
+ SingleFluidBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex );
+ resizeFields( parent.size(), numConstitutivePointsPerParentIndex );
+}
+
+// void ReactiveSingleFluid::createChemicalReactions()
+// {
+// // instantiate reactions objects
+// m_equilibriumReactions = std::make_unique< chemicalReactions::EquilibriumReactions >( getName() + "_equilibriumReactions",
+// m_numPrimarySpecies, m_numSecondarySpecies );
+// m_kineticReactions = std::make_unique< chemicalReactions::KineticReactions >( getName() + "_kineticReactions", m_numPrimarySpecies,
+// m_numSecondarySpecies, m_numKineticReactions );
+// }
+
+} //namespace constitutive
+
+} //namespace geos
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.hpp
new file mode 100644
index 00000000000..48a5457faba
--- /dev/null
+++ b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.hpp
@@ -0,0 +1,238 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ReactiveSingleFluid.hpp
+ */
+
+#ifndef GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_REACTIVE_REACTIVESINGLEFLUID_HPP_
+#define GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_REACTIVE_REACTIVESINGLEFLUID_HPP_
+
+
+#include "common/format/EnumStrings.hpp"
+#include "constitutive/fluid/singlefluid/SingleFluidBase.hpp"
+#include "constitutive/fluid/multifluid/reactive/chemicalReactions/EquilibriumReactions.hpp"
+#include "constitutive/fluid/multifluid/reactive/chemicalReactions/KineticReactions.hpp"
+
+#include
+
+namespace geos
+{
+
+namespace constitutive
+{
+
+class ReactiveSingleFluidUpdate : public SingleFluidBaseUpdate
+{
+public:
+
+ GEOS_HOST_DEVICE
+ void computeChemistry( real64 const pressure,
+ real64 const temperature,
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & primarySpeciesAggregateConcentration,
+ arraySlice1d< real64, compflow::USD_COMP - 1 > const & primarySpeciesConcentration,
+ arraySlice1d< real64, compflow::USD_COMP - 1 > const & secondarySpeciesConcentration,
+ arraySlice1d< real64, compflow::USD_COMP - 1 > const & kineticReactionRates ) const;
+
+ GEOS_HOST_DEVICE
+ virtual void updateChemistry( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const temperature,
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & primaryConc ) const = 0;
+
+ GEOS_HOST_DEVICE
+ virtual void updateChemistryLogConc( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const temperature,
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & logPrimaryConc ) const = 0;
+
+protected:
+
+ ReactiveSingleFluidUpdate( arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const & density,
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > const & dDensity,
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const & viscosity,
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > const & dViscosity,
+ integer const numPrimarySpecies,
+ // chemicalReactions::EquilibriumReactions const & equilibriumReactions,
+ // chemicalReactions::KineticReactions const & kineticReactions,
+ arrayView2d< real64, compflow::USD_COMP > const & primarySpeciesConcentration,
+ arrayView2d< real64, compflow::USD_COMP > const & secondarySpeciesConcentration,
+ arrayView2d< real64, compflow::USD_COMP > const & primarySpeciesAggregateConcentration,
+ arrayView3d< real64, compflow::USD_COMP_DC > const & dPrimarySpeciesAggregateConcentration_dLogPrimaryConc,
+ arrayView2d< real64, compflow::USD_COMP > const & kineticReactionRates )
+ : SingleFluidBaseUpdate( density,
+ dDensity,
+ viscosity,
+ dViscosity ),
+ m_numPrimarySpecies( numPrimarySpecies ),
+ // m_equilibriumReactions( equilibriumReactions.createKernelWrapper() ),
+ // m_kineticReactions( kineticReactions.createKernelWrapper() ),
+ m_primarySpeciesConcentration( primarySpeciesConcentration ),
+ m_secondarySpeciesConcentration( secondarySpeciesConcentration ),
+ m_primarySpeciesAggregateConcentration( primarySpeciesAggregateConcentration ),
+ m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc( dPrimarySpeciesAggregateConcentration_dLogPrimaryConc ),
+ m_kineticReactionRates( kineticReactionRates )
+ {}
+
+ /**
+ * @brief Copy constructor.
+ */
+ ReactiveSingleFluidUpdate( ReactiveSingleFluidUpdate const & ) = default;
+
+ /**
+ * @brief Move constructor.
+ */
+ ReactiveSingleFluidUpdate( ReactiveSingleFluidUpdate && ) = default;
+
+ /**
+ * @brief Deleted copy assignment operator
+ * @return reference to this object
+ */
+ ReactiveSingleFluidUpdate & operator=( ReactiveSingleFluidUpdate const & ) = delete;
+
+ /**
+ * @brief Deleted move assignment operator
+ * @return reference to this object
+ */
+ ReactiveSingleFluidUpdate & operator=( ReactiveSingleFluidUpdate && ) = delete;
+
+ /// Reaction related terms
+ integer m_numPrimarySpecies;
+
+ // chemicalReactions::EquilibriumReactions::KernelWrapper m_equilibriumReactions;
+
+ // chemicalReactions::KineticReactions::KernelWrapper m_kineticReactions;
+
+ arrayView2d< real64, compflow::USD_COMP > m_primarySpeciesConcentration;
+
+ arrayView2d< real64, compflow::USD_COMP > m_secondarySpeciesConcentration;
+
+ arrayView2d< real64, compflow::USD_COMP > m_primarySpeciesAggregateConcentration;
+
+ arrayView3d< real64, compflow::USD_COMP_DC > m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc;
+
+ arrayView2d< real64, compflow::USD_COMP > m_kineticReactionRates;
+};
+
+class ReactiveSingleFluid : public SingleFluidBase
+{
+public:
+
+ using exec_policy = serialPolicy;
+
+ ReactiveSingleFluid( string const & name,
+ Group * const parent );
+
+ virtual void saveConvergedState() const override;
+
+ virtual void allocateConstitutiveData( dataRepository::Group & parent,
+ localIndex const numConstitutivePointsPerParentIndex ) override;
+
+ arrayView2d< real64 const, compflow::USD_COMP > primarySpeciesConcentration() const
+ { return m_primarySpeciesConcentration; }
+
+ arrayView2d< real64 const, compflow::USD_COMP > primarySpeciesAggregateConcentration() const
+ { return m_primarySpeciesAggregateConcentration; }
+
+ arrayView2d< real64 const, compflow::USD_COMP > primarySpeciesAggregateConcentration_n() const
+ { return m_primarySpeciesAggregateConcentration_n; }
+
+ arrayView3d< real64 const, compflow::USD_COMP_DC > dPrimarySpeciesAggregateConcentration_dLogPrimaryConc() const
+ { return m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc; }
+
+ arrayView2d< real64 const, compflow::USD_COMP > secondarySpeciesConcentration() const
+ { return m_secondarySpeciesConcentration; }
+
+ arrayView2d< real64 const, compflow::USD_COMP > kineticReactionRates() const
+ { return m_kineticReactionRates; }
+
+ integer numPrimarySpecies() const { return m_numPrimarySpecies; }
+
+ integer numSecondarySpecies() const { return m_numSecondarySpecies; }
+
+ integer numKineticReactions() const { return m_numKineticReactions; }
+
+
+ struct viewKeyStruct : ConstitutiveBase::viewKeyStruct
+ {
+ static constexpr char const * primarySpeciesNamesString() { return "primarySpeciesNames"; }
+ };
+
+protected:
+
+ virtual void postInputInitialization() override;
+
+ void createChemicalReactions();
+
+ virtual void resizeFields( localIndex const size, localIndex const numPts );
+
+ /// Reaction related terms
+ array1d< string > m_primarySpeciesNames;
+
+ integer m_numPrimarySpecies;
+
+ integer m_numSecondarySpecies;
+
+ integer m_numKineticReactions;
+
+ // std::unique_ptr< chemicalReactions::EquilibriumReactions > m_equilibriumReactions;
+
+ // std::unique_ptr< chemicalReactions::KineticReactions > m_kineticReactions;
+
+ array2d< real64, constitutive::multifluid::LAYOUT_FLUID > m_primarySpeciesConcentration;
+
+ array2d< real64, constitutive::multifluid::LAYOUT_FLUID > m_secondarySpeciesConcentration;
+
+ array2d< real64, constitutive::multifluid::LAYOUT_FLUID > m_primarySpeciesAggregateConcentration;
+
+ array2d< real64, constitutive::multifluid::LAYOUT_FLUID > m_primarySpeciesAggregateConcentration_n;
+
+ array3d< real64, constitutive::multifluid::LAYOUT_FLUID_DC > m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc;
+
+ array2d< real64, constitutive::multifluid::LAYOUT_FLUID > m_kineticReactionRates;
+};
+
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+void ReactiveSingleFluidUpdate::
+ computeChemistry( real64 const pressure,
+ real64 const temperature,
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & primarySpeciesAggregateConcentration,
+ arraySlice1d< real64, compflow::USD_COMP - 1 > const & primarySpeciesConcentration,
+ arraySlice1d< real64, compflow::USD_COMP - 1 > const & secondarySpeciesConcentration,
+ arraySlice1d< real64, compflow::USD_COMP - 1 > const & kineticReactionRates ) const
+{
+ GEOS_UNUSED_VAR( pressure, temperature, primarySpeciesAggregateConcentration, primarySpeciesConcentration, secondarySpeciesConcentration, kineticReactionRates );
+
+ // // 2. solve for equilibrium
+ // m_equilibriumReactions.updateConcentrations( temperature,
+ // primarySpeciesAggregateConcentration,
+ // primarySpeciesConcentration,
+ // secondarySpeciesConcentration );
+
+ // // 3. compute kinetic reaction rates
+ // m_kineticReactions.computeReactionRates( temperature,
+ // primarySpeciesConcentration,
+ // secondarySpeciesConcentration,
+ // kineticReactionRates );
+}
+
+} // namespace constitutive
+
+} // namespace geos
+
+#endif //GEOS_CONSTITUTIVE_FLUID_REACTIVEMULTIFLUID_HPP
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.cpp
new file mode 100644
index 00000000000..3f7a746cf72
--- /dev/null
+++ b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.cpp
@@ -0,0 +1,128 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ThermalReactiveCompressibleSinglePhaseFluid.cpp
+ */
+
+#include "ThermalReactiveCompressibleSinglePhaseFluid.hpp"
+
+#include "constitutive/fluid/singlefluid/SingleFluidFields.hpp"
+
+namespace geos
+{
+
+using namespace dataRepository;
+
+namespace constitutive
+{
+
+ThermalReactiveCompressibleSinglePhase::ThermalReactiveCompressibleSinglePhase( string const & name, Group * const parent ):
+ ReactiveCompressibleSinglePhase( name, parent ),
+ m_internalEnergyModelType( ExponentApproximationType::Linear )
+{
+ m_densityModelType = ExponentApproximationType::Full;
+ registerWrapper( viewKeyStruct::thermalExpansionCoeffString(), &m_thermalExpansionCoeff ).
+ setApplyDefaultValue( 0.0 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Fluid thermal expansion coefficient. Unit: 1/K" );
+
+ registerWrapper( viewKeyStruct::specificHeatCapacityString(), &m_specificHeatCapacity ).
+ setApplyDefaultValue( 0.0 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Fluid heat capacity. Unit: J/kg/K" );
+
+ registerWrapper( viewKeyStruct::referenceTemperatureString(), &m_referenceTemperature ).
+ setApplyDefaultValue( 0.0 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Reference temperature" );
+
+ registerWrapper( viewKeyStruct::referenceInternalEnergyString(), &m_referenceInternalEnergy ).
+ setApplyDefaultValue( 0.001 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Reference fluid internal energy" );
+
+ registerWrapper( viewKeyStruct::internalEnergyModelTypeString(), &m_internalEnergyModelType ).
+ setApplyDefaultValue( m_internalEnergyModelType ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Type of internal energy model. Valid options:\n* " + EnumStrings< ExponentApproximationType >::concat( "\n* " ) );
+
+}
+
+ThermalReactiveCompressibleSinglePhase::~ThermalReactiveCompressibleSinglePhase() = default;
+
+void ThermalReactiveCompressibleSinglePhase::allocateConstitutiveData( dataRepository::Group & parent,
+ localIndex const numConstitutivePointsPerParentIndex )
+{
+ ReactiveCompressibleSinglePhase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex );
+
+ m_internalEnergy.value.setValues< serialPolicy >( m_referenceInternalEnergy );
+}
+
+void ThermalReactiveCompressibleSinglePhase::postInputInitialization()
+{
+ ReactiveCompressibleSinglePhase::postInputInitialization();
+
+ auto const checkNonnegative = [&]( real64 const value, auto const & attribute )
+ {
+ GEOS_THROW_IF_LT_MSG( value, 0.0,
+ GEOS_FMT( "{}: invalid value of attribute '{}'", getFullName(), attribute ),
+ InputError );
+ };
+
+ checkNonnegative( m_thermalExpansionCoeff, viewKeyStruct::thermalExpansionCoeffString() );
+ checkNonnegative( m_specificHeatCapacity, viewKeyStruct::specificHeatCapacityString() );
+ checkNonnegative( m_referenceInternalEnergy, viewKeyStruct::referenceInternalEnergyString() );
+
+ // Due to the way update wrapper is currently implemented, we can only support one model type
+ auto const checkModelType = [&]( ExponentApproximationType const value, auto const & attribute )
+ {
+ GEOS_THROW_IF( value != ExponentApproximationType::Linear && value != ExponentApproximationType::Full,
+ GEOS_FMT( "{}: invalid model type in attribute '{}' (only linear or fully exponential currently supported)", getFullName(), attribute ),
+ InputError );
+ };
+ checkModelType( m_internalEnergyModelType, viewKeyStruct::internalEnergyModelTypeString() );
+}
+
+ThermalReactiveCompressibleSinglePhase::KernelWrapper
+ThermalReactiveCompressibleSinglePhase::createKernelWrapper()
+{
+ return KernelWrapper( KernelWrapper::DensRelationType( m_referencePressure, m_referenceTemperature, m_referenceDensity, m_compressibility, -m_thermalExpansionCoeff ),
+ KernelWrapper::ViscRelationType( m_referencePressure, m_referenceViscosity, m_viscosibility ),
+ KernelWrapper::IntEnergyRelationType( m_referenceTemperature, m_referenceInternalEnergy, m_specificHeatCapacity/m_referenceInternalEnergy ),
+ m_density.value,
+ m_density.derivs,
+ m_viscosity.value,
+ m_viscosity.derivs,
+ m_numPrimarySpecies,
+ // *m_equilibriumReactions,
+ // *m_kineticReactions,
+ m_primarySpeciesConcentration.toView(),
+ m_secondarySpeciesConcentration.toView(),
+ m_primarySpeciesAggregateConcentration.toView(),
+ m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc.toView(),
+ m_kineticReactionRates.toView(),
+ m_internalEnergy.value,
+ m_internalEnergy.derivs,
+ m_enthalpy.value,
+ m_enthalpy.derivs,
+ m_referenceInternalEnergy );
+}
+
+REGISTER_CATALOG_ENTRY( ConstitutiveBase, ThermalReactiveCompressibleSinglePhase, string const &, Group * const )
+
+} /* namespace constitutive */
+
+} /* namespace geos */
diff --git a/src/coreComponents/constitutive/fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.hpp
new file mode 100644
index 00000000000..8e0cdbdcfb1
--- /dev/null
+++ b/src/coreComponents/constitutive/fluid/singlefluid/reactive/ThermalReactiveCompressibleSinglePhaseFluid.hpp
@@ -0,0 +1,333 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ThermalReactiveCompressibleSinglePhaseFluid.hpp
+ */
+
+#ifndef GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_THERMALREACTIVECOMPRESSIBLESINGLEPHASEFLUID_HPP_
+#define GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_THERMALREACTIVECOMPRESSIBLESINGLEPHASEFLUID_HPP_
+
+#include "constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ReactiveCompressibleSinglePhaseFluid.hpp"
+
+#include "constitutive/ExponentialRelation.hpp"
+
+namespace geos
+{
+
+namespace constitutive
+{
+
+/**
+ * @brief Update class for the model suitable for lambda capture.
+ * @tparam DENS_EAT type of density exponent approximation for the pressure part
+ * @tparam VISC_EAT type of viscosity exponent approximation
+ * @tparam INTENERGY_EAT type of internal energy exponent approximation
+ */
+template< ExponentApproximationType DENS_EAT, ExponentApproximationType VISC_EAT, ExponentApproximationType INTENERGY_EAT >
+class ThermalReactiveCompressibleSinglePhaseUpdate : public ReactiveSingleFluidUpdate
+{
+public:
+
+ using DensRelationType = ExponentialRelation< real64, DENS_EAT, 3 >;
+ using ViscRelationType = ExponentialRelation< real64, VISC_EAT >;
+ using IntEnergyRelationType = ExponentialRelation< real64, INTENERGY_EAT >;
+ using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
+
+ ThermalReactiveCompressibleSinglePhaseUpdate( DensRelationType const & densRelation,
+ ViscRelationType const & viscRelation,
+ IntEnergyRelationType const & intEnergyRelation,
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const & density,
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > const & dDensity,
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const & viscosity,
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > const & dViscosity,
+ integer const numPrimarySpecies,
+ // chemicalReactions::EquilibriumReactions const & equilibriumReactions,
+ // chemicalReactions::KineticReactions const & kineticReactions,
+ arrayView2d< real64, compflow::USD_COMP > const & primarySpeciesConcentration,
+ arrayView2d< real64, compflow::USD_COMP > const & secondarySpeciesConcentration,
+ arrayView2d< real64, compflow::USD_COMP > const & primarySpeciesAggregateConcentration,
+ arrayView3d< real64, compflow::USD_COMP_DC > const & dPrimarySpeciesAggregateConcentration_dLogPrimaryConc,
+ arrayView2d< real64, compflow::USD_COMP > const & kineticReactionRates,
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const & internalEnergy,
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > const & dInternalEnergy,
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > const & enthalpy,
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > const & dEnthalpy,
+ real64 const & refIntEnergy )
+ : ReactiveSingleFluidUpdate( density, dDensity, viscosity, dViscosity, numPrimarySpecies,
+ // equilibriumReactions, kineticReactions,
+ primarySpeciesConcentration, secondarySpeciesConcentration, primarySpeciesAggregateConcentration,
+ dPrimarySpeciesAggregateConcentration_dLogPrimaryConc, kineticReactionRates ),
+ m_internalEnergy( internalEnergy ),
+ m_dInternalEnergy( dInternalEnergy ),
+ m_enthalpy( enthalpy ),
+ m_dEnthalpy( dEnthalpy ),
+ m_densRelation( densRelation ),
+ m_viscRelation( viscRelation ),
+ m_intEnergyRelation( intEnergyRelation ),
+ m_refIntEnergy( refIntEnergy )
+ {}
+
+ /// Default copy constructor
+ ThermalReactiveCompressibleSinglePhaseUpdate( ThermalReactiveCompressibleSinglePhaseUpdate const & ) = default;
+
+ /// Default move constructor
+ ThermalReactiveCompressibleSinglePhaseUpdate( ThermalReactiveCompressibleSinglePhaseUpdate && ) = default;
+
+ /// Deleted copy assignment operator
+ ThermalReactiveCompressibleSinglePhaseUpdate & operator=( ThermalReactiveCompressibleSinglePhaseUpdate const & ) = delete;
+
+ /// Deleted move assignment operator
+ ThermalReactiveCompressibleSinglePhaseUpdate & operator=( ThermalReactiveCompressibleSinglePhaseUpdate && ) = delete;
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void compute( real64 const pressure,
+ real64 & density,
+ real64 & dDensity_dPressure,
+ real64 & viscosity,
+ real64 & dViscosity_dPressure ) const override
+ {
+ m_densRelation.compute( pressure, density, dDensity_dPressure );
+ m_viscRelation.compute( pressure, viscosity, dViscosity_dPressure );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void compute( real64 const pressure,
+ real64 const temperature,
+ real64 & density,
+ real64 & dDensity_dPressure,
+ real64 & dDensity_dTemperature,
+ real64 & viscosity,
+ real64 & dViscosity_dPressure,
+ real64 & dViscosity_dTemperature,
+ real64 & internalEnergy,
+ real64 & dInternalEnergy_dPressure,
+ real64 & dInternalEnergy_dTemperature,
+ real64 & enthalpy,
+ real64 & dEnthalpy_dPressure,
+ real64 & dEnthalpy_dTemperature ) const override
+ {
+ m_viscRelation.compute( pressure, viscosity, dViscosity_dPressure );
+ dViscosity_dTemperature = 0.0;
+
+ m_densRelation.compute( pressure, temperature, density, dDensity_dPressure, dDensity_dTemperature );
+
+ /// Compute the internal energy (only sensitive to temperature)
+ m_intEnergyRelation.compute( temperature, internalEnergy, dInternalEnergy_dTemperature );
+ dInternalEnergy_dPressure = 0.0;
+
+ enthalpy = internalEnergy - m_refIntEnergy;
+ dEnthalpy_dPressure = 0.0;
+ dEnthalpy_dTemperature = dInternalEnergy_dTemperature;
+
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure ) const override
+ {
+ compute( pressure,
+ m_density[k][q],
+ m_dDensity[k][q][DerivOffset::dP],
+ m_viscosity[k][q],
+ m_dViscosity[k][q][DerivOffset::dP] );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const temperature ) const override
+ {
+ compute( pressure,
+ temperature,
+ m_density[k][q],
+ m_dDensity[k][q][DerivOffset::dP],
+ m_dDensity[k][q][DerivOffset::dT],
+ m_viscosity[k][q],
+ m_dViscosity[k][q][DerivOffset::dP],
+ m_dViscosity[k][q][DerivOffset::dT],
+ m_internalEnergy[k][q],
+ m_dInternalEnergy[k][q][DerivOffset::dP],
+ m_dInternalEnergy[k][q][DerivOffset::dT],
+ m_enthalpy[k][q],
+ m_dEnthalpy[k][q][DerivOffset::dP],
+ m_dEnthalpy[k][q][DerivOffset::dT] );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void update( localIndex const k,
+ localIndex const q,
+ real64 const pressure,
+ real64 const temperature,
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & GEOS_UNUSED_PARAM( composition ) ) const override
+ {
+ compute( pressure,
+ temperature,
+ m_density[k][q],
+ m_dDensity[k][q][DerivOffset::dP],
+ m_dDensity[k][q][DerivOffset::dT],
+ m_viscosity[k][q],
+ m_dViscosity[k][q][DerivOffset::dP],
+ m_dViscosity[k][q][DerivOffset::dT],
+ m_internalEnergy[k][q],
+ m_dInternalEnergy[k][q][DerivOffset::dP],
+ m_dInternalEnergy[k][q][DerivOffset::dT],
+ m_enthalpy[k][q],
+ m_dEnthalpy[k][q][DerivOffset::dP],
+ m_dEnthalpy[k][q][DerivOffset::dT] );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void updateChemistry( localIndex const k,
+ localIndex const GEOS_UNUSED_PARAM( q ),
+ real64 const pressure,
+ real64 const temperature,
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & primaryConc ) const override
+ {
+ for( int i=0; i < m_numPrimarySpecies; i++ )
+ {
+ m_primarySpeciesAggregateConcentration[k][i] = primaryConc[i];
+ }
+
+ computeChemistry( pressure,
+ temperature,
+ m_primarySpeciesAggregateConcentration[k],
+ m_primarySpeciesConcentration[k],
+ m_secondarySpeciesConcentration[k],
+ m_kineticReactionRates[k] );
+ }
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ virtual void updateChemistryLogConc( localIndex const k,
+ localIndex const GEOS_UNUSED_PARAM( q ),
+ real64 const GEOS_UNUSED_PARAM( pressure ),
+ real64 const GEOS_UNUSED_PARAM( temperature ),
+ arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & logPrimaryConc ) const override
+ {
+ for( int i=0; i < m_numPrimarySpecies; i++ )
+ {
+ m_primarySpeciesConcentration[k][i] = std::exp( logPrimaryConc[i] );
+
+ m_primarySpeciesAggregateConcentration[k][i] = m_primarySpeciesConcentration[k][i];
+
+ m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc[k][i][i] = m_primarySpeciesConcentration[k][i];
+ }
+
+ // computeChemistry( pressure,
+ // temperature,
+ // m_primarySpeciesAggregateConcentration[k],
+ // m_primarySpeciesConcentration[k],
+ // m_secondarySpeciesConcentration[k],
+ // m_kineticReactionRates[k] );
+ }
+
+private:
+
+ /// Fluid internal energy and derivatives
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > m_internalEnergy;
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > m_dInternalEnergy;
+
+ /// Fluid enthalpy and derivatives
+ arrayView2d< real64, constitutive::singlefluid::USD_FLUID > m_enthalpy;
+ arrayView3d< real64, constitutive::singlefluid::USD_FLUID_DER > m_dEnthalpy;
+
+ /// Relationship between the fluid density and pressure & temperature
+ DensRelationType m_densRelation;
+
+ /// Relationship between the fluid viscosity and pressure
+ ViscRelationType m_viscRelation;
+
+ /// Relationship between the fluid internal energy and temperature
+ IntEnergyRelationType m_intEnergyRelation;
+
+ /// Reference internal energy of the fluid
+ real64 const m_refIntEnergy;
+
+};
+
+class ThermalReactiveCompressibleSinglePhase : public ReactiveCompressibleSinglePhase
+{
+public:
+
+ ThermalReactiveCompressibleSinglePhase( string const & name, Group * const parent );
+
+ virtual ~ThermalReactiveCompressibleSinglePhase() override;
+
+ static string catalogName() { return "ThermalReactiveCompressibleSinglePhase"; }
+
+ virtual string getCatalogName() const override { return catalogName(); }
+
+ virtual void allocateConstitutiveData( dataRepository::Group & parent,
+ localIndex const numConstitutivePointsPerParentIndex ) override;
+
+ using ReactiveCompressibleSinglePhase::m_densityModelType;
+
+ /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT, not just linear)
+ using KernelWrapper = ThermalReactiveCompressibleSinglePhaseUpdate< ExponentApproximationType::Full, ExponentApproximationType::Linear, ExponentApproximationType::Linear >;
+
+ /**
+ * @brief Create an update kernel wrapper.
+ * @return the wrapper
+ */
+ KernelWrapper createKernelWrapper();
+
+ struct viewKeyStruct : public ReactiveCompressibleSinglePhase::viewKeyStruct
+ {
+ static constexpr char const * thermalExpansionCoeffString() { return "thermalExpansionCoeff"; }
+ static constexpr char const * specificHeatCapacityString() { return "specificHeatCapacity"; }
+ static constexpr char const * referenceTemperatureString() { return "referenceTemperature"; }
+ static constexpr char const * referenceInternalEnergyString() { return "referenceInternalEnergy"; }
+ static constexpr char const * internalEnergyModelTypeString() { return "internalEnergyModelType"; }
+ };
+
+ virtual bool isThermal() const override { return true; }
+
+protected:
+
+ virtual void postInputInitialization() override;
+
+private:
+
+ /// scalar fluid thermal expansion coefficient
+ real64 m_thermalExpansionCoeff;
+
+ /// scalar fluid volumetric heat capacity coefficient
+ real64 m_specificHeatCapacity;
+
+ /// reference temperature parameter
+ real64 m_referenceTemperature;
+
+ /// reference internal energy parameter
+ real64 m_referenceInternalEnergy;
+
+ /// type of internal energy model
+ ExponentApproximationType m_internalEnergyModelType;
+};
+
+} /* namespace constitutive */
+
+} /* namespace geos */
+
+#endif /* GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_THERMALREACTIVECOMPRESSIBLESINGLEPHASEFLUID_HPP_ */
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt
index ce0829dd415..f7f6af1aca4 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt
@@ -63,6 +63,13 @@ set( fluidFlowSolvers_headers
kernels/singlePhase/ThermalFluxComputeKernel.hpp
kernels/singlePhase/proppant/ProppantBaseKernels.hpp
kernels/singlePhase/proppant/ProppantFluxKernels.hpp
+ kernels/singlePhase/reactive/AccumulationKernels.hpp
+ kernels/singlePhase/reactive/FluidUpdateKernel.hpp
+ kernels/singlePhase/reactive/FluxComputeKernel.hpp
+ kernels/singlePhase/reactive/KernelLaunchSelectors.hpp
+ kernels/singlePhase/reactive/ResidualNormKernel.hpp
+ kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp
+ kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp
kernels/compositional/AccumulationKernel.hpp
kernels/compositional/AquiferBCKernel.hpp
kernels/compositional/PPUPhaseFlux.hpp
@@ -134,6 +141,7 @@ set( fluidFlowSolvers_sources
SinglePhaseFVM.cpp
SinglePhaseHybridFVM.cpp
SinglePhaseProppantBase.cpp
+ SinglePhaseReactiveTransport.cpp
SourceFluxStatistics.cpp
StencilDataCollection.cpp
kernels/singlePhase/proppant/ProppantFluxKernels.cpp
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index 6a17a1ef7d4..3a4f750000d 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -195,11 +195,11 @@ void SinglePhaseBase::validateConstitutiveModels( DomainPartition & domain ) con
constitutiveUpdatePassThru( fluid, [&] ( auto & castedFluid )
{
string const fluidModelName = castedFluid.getCatalogName();
- GEOS_THROW_IF( m_isThermal && (fluidModelName != "ThermalCompressibleSinglePhaseFluid"),
+ GEOS_THROW_IF( m_isThermal && ((fluidModelName != "ThermalCompressibleSinglePhaseFluid") && (fluidModelName != "ThermalReactiveCompressibleSinglePhase")),
GEOS_FMT( "SingleFluidBase {}: the thermal option is enabled in the solver, but the fluid model {} is not for thermal fluid",
getDataContext(), fluid.getDataContext() ),
InputError );
- GEOS_THROW_IF( !m_isThermal && (fluidModelName == "ThermalCompressibleSinglePhaseFluid"),
+ GEOS_THROW_IF( !m_isThermal && ((fluidModelName == "ThermalCompressibleSinglePhaseFluid") || (fluidModelName == "ThermalReactiveCompressibleSinglePhase")),
GEOS_FMT( "SingleFluidBase {}: the fluid model is for thermal fluid {}, but the solver option is incompatible with the fluid model",
getDataContext(), fluid.getDataContext() ),
InputError );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
index 245aedaf584..056e5564c21 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
@@ -222,7 +222,7 @@ class SinglePhaseBase : public FlowSolverBase
* @param localMatrix local system matrix
* @param localRhs local system right-hand side vector
*/
- void
+ virtual void
applyDirichletBC( real64 const time_n,
real64 const dt,
DomainPartition & domain,
@@ -265,7 +265,7 @@ class SinglePhaseBase : public FlowSolverBase
arrayView1d< real64 > const & localRhs ) const = 0;
virtual void
- updateState ( DomainPartition & domain ) override final;
+ updateState ( DomainPartition & domain ) override;
/**
* @brief Function to update all constitutive state and dependent variables
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp
new file mode 100644
index 00000000000..705ae361db5
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp
@@ -0,0 +1,1030 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file SinglePhaseReactiveTransport.cpp
+ */
+
+#include "SinglePhaseReactiveTransport.hpp"
+
+#include "constitutive/ConstitutiveManager.hpp"
+#include "constitutive/ConstitutivePassThru.hpp"
+#include "constitutive/diffusion/DiffusionFields.hpp"
+#include "constitutive/diffusion/DiffusionSelector.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.hpp"
+#include "constitutive/fluid/multifluid/reactive/ReactiveFluidSelector.hpp"
+#include "finiteVolume/FluxApproximationBase.hpp"
+#include "mesh/DomainPartition.hpp"
+#include "physicsSolvers/LogLevelsInfo.hpp"
+
+/**
+ * @namespace the geos namespace that encapsulates the majority of the code
+ */
+namespace geos
+{
+
+using namespace dataRepository;
+using namespace constitutive;
+
+SinglePhaseReactiveTransport::SinglePhaseReactiveTransport( const string & name,
+ Group * const parent ):
+ SinglePhaseBase( name, parent ),
+ m_numPrimarySpecies( 0 ),
+ m_hasDiffusion( 0 )
+{
+ // To add modeling parameters we want to add here
+
+ addLogLevel< logInfo::Convergence >();
+}
+
+// TODO: we need to update the class of ReactiveSingleFluid to be consistent with the chemistry module!!!
+void SinglePhaseReactiveTransport::registerDataOnMesh( Group & meshBodies )
+{
+ using namespace fields::flow;
+
+ SinglePhaseBase::registerDataOnMesh( meshBodies );
+
+ DomainPartition const & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+ ConstitutiveManager const & cm = domain.getConstitutiveManager();
+
+ // 0. Find a reactive fluid model name (at this point, models are already attached to subregions)
+ forDiscretizationOnMeshTargets( meshBodies, [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ mesh.getElemManager().forElementSubRegions( regionNames,
+ [&]( localIndex const,
+ ElementSubRegionBase & subRegion )
+ {
+ if( m_reactiveFluidModelName.empty() )
+ {
+ m_reactiveFluidModelName = getConstitutiveName< ReactiveSingleFluid >( subRegion );
+ }
+
+ // If at least one region has a diffusion model, consider it enabled for all
+ string const diffusionName = getConstitutiveName< DiffusionBase >( subRegion );
+ if( !diffusionName.empty() )
+ {
+ m_hasDiffusion = true;
+ }
+ } );
+ } );
+
+ // 1. Set key dimensions of the problem
+ // Check needed to avoid errors when running in schema generation mode.
+ if( !m_reactiveFluidModelName.empty() )
+ {
+ ReactiveSingleFluid const & reactiveFluid = cm.getConstitutiveRelation< ReactiveSingleFluid >( m_reactiveFluidModelName );
+ m_numPrimarySpecies = reactiveFluid.numPrimarySpecies();
+ m_isThermal = reactiveFluid.isThermal();
+ }
+
+ // n_c components + one pressure ( + one temperature if needed )
+ m_numDofPerCell = m_isThermal ? m_numPrimarySpecies + 2 : m_numPrimarySpecies + 1;
+
+ // 2. Register and resize all fields as necessary
+ forDiscretizationOnMeshTargets( meshBodies, [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ ElementRegionManager & elemManager = mesh.getElemManager();
+
+ elemManager.forElementSubRegions< ElementSubRegionBase >( regionNames,
+ [&]( localIndex const,
+ ElementSubRegionBase & subRegion )
+ {
+ if( m_hasDiffusion )
+ {
+ subRegion.registerWrapper< string >( viewKeyStruct::diffusionNamesString() ).
+ setPlotLevel( PlotLevel::NOPLOT ).
+ setRestartFlags( RestartFlags::NO_WRITE ).
+ setSizedFromParent( 0 ).
+ setDescription( "Name of the diffusion constitutive model to use" );
+
+ string & diffusionName = subRegion.getReference< string >( viewKeyStruct::diffusionNamesString() );
+ diffusionName = getConstitutiveName< DiffusionBase >( subRegion );
+ GEOS_THROW_IF( diffusionName.empty(),
+ GEOS_FMT( "Diffusion model not found on subregion {}", subRegion.getName() ),
+ InputError );
+ }
+
+ subRegion.registerField< logPrimarySpeciesConcentration >( getName() ).
+ reference().resizeDimension< 1 >( m_numPrimarySpecies );
+
+ subRegion.registerField< logPrimarySpeciesConcentration_n >( getName() ).
+ reference().resizeDimension< 1 >( m_numPrimarySpecies );
+
+ subRegion.registerField< primarySpeciesAggregateMole >( getName() ).
+ reference().resizeDimension< 1 >( m_numPrimarySpecies );
+
+ subRegion.registerField< primarySpeciesAggregateMole_n >( getName() ).
+ reference().resizeDimension< 1 >( m_numPrimarySpecies );
+
+ subRegion.registerField< bcLogPrimarySpeciesConcentration >( getName() ).
+ reference().resizeDimension< 1 >( m_numPrimarySpecies );
+ } );
+ } );
+}
+
+void SinglePhaseReactiveTransport::setupDofs( DomainPartition const & domain,
+ DofManager & dofManager ) const
+{
+ // add a field for the cell-centered degrees of freedom
+ dofManager.addField( viewKeyStruct::elemDofFieldString(),
+ FieldLocation::Elem,
+ m_numDofPerCell,
+ getMeshTargets() );
+
+ NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
+ FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
+ FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName );
+
+ dofManager.addCoupling( viewKeyStruct::elemDofFieldString(), fluxApprox );
+}
+
+void SinglePhaseReactiveTransport::resetStateToBeginningOfStep( DomainPartition & domain )
+{
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const,
+ auto & subRegion )
+ {
+ arrayView1d< real64 > const pres = subRegion.template getField< fields::flow::pressure >();
+ arrayView1d< real64 const > const pres_n = subRegion.template getField< fields::flow::pressure_n >();
+ pres.setValues< parallelDevicePolicy<> >( pres_n );
+
+ arrayView2d< real64, compflow::USD_COMP > const logPrimarySpeciesConc = subRegion.template getField< fields::flow::logPrimarySpeciesConcentration >();
+ arrayView2d< real64 const, compflow::USD_COMP > const logPrimarySpeciesConc_n = subRegion.template getField< fields::flow::logPrimarySpeciesConcentration_n >();
+ logPrimarySpeciesConc.setValues< parallelDevicePolicy<> >( logPrimarySpeciesConc_n );
+
+ if( m_isThermal )
+ {
+ arrayView1d< real64 > const temp = subRegion.template getField< fields::flow::temperature >();
+ arrayView1d< real64 const > const temp_n = subRegion.template getField< fields::flow::temperature_n >();
+ temp.setValues< parallelDevicePolicy<> >( temp_n );
+ }
+
+ updatePorosityAndPermeability( subRegion );
+ updateFluidState( subRegion );
+
+ if( m_isThermal )
+ {
+ updateSolidInternalEnergyModel( subRegion );
+ updateEnergy( subRegion );
+ }
+ } );
+ } );
+}
+
+void SinglePhaseReactiveTransport::implicitStepComplete( real64 const & time,
+ real64 const & dt,
+ DomainPartition & domain )
+{
+ GEOS_MARK_FUNCTION;
+
+ SinglePhaseBase::implicitStepComplete( time, dt, domain );
+
+ if( m_hasDiffusion )
+ {
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ mesh.getElemManager().forElementSubRegions( regionNames,
+ [&]( localIndex const,
+ ElementSubRegionBase & subRegion )
+ {
+ string const & diffusionName = subRegion.getReference< string >( viewKeyStruct::diffusionNamesString() );
+ DiffusionBase const & diffusionMaterial = getConstitutiveModel< DiffusionBase >( subRegion, diffusionName );
+ arrayView1d< real64 const > const temperature = subRegion.template getField< fields::flow::temperature >();
+ diffusionMaterial.saveConvergedTemperatureState( temperature );
+ } );
+ } );
+ }
+}
+
+void SinglePhaseReactiveTransport::assembleSystem( real64 const GEOS_UNUSED_PARAM( time_n ),
+ real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+{
+ GEOS_MARK_FUNCTION;
+
+ assembleAccumulationTermsInMassBalanceAndSpeciesAmountEqs( dt,
+ domain,
+ dofManager,
+ localMatrix,
+ localRhs );
+ assembleFluxTerms( dt,
+ domain,
+ dofManager,
+ localMatrix,
+ localRhs );
+}
+
+void SinglePhaseReactiveTransport::assembleAccumulationTermsInMassBalanceAndSpeciesAmountEqs( real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) const
+{
+ GEOS_MARK_FUNCTION;
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ mesh.getElemManager().forElementSubRegions( regionNames,
+ [&]( localIndex const,
+ ElementSubRegionBase const & subRegion )
+ {
+ geos::constitutive::ReactiveSingleFluid const & fluid =
+ getConstitutiveModel< geos::constitutive::ReactiveSingleFluid >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) );
+ geos::constitutive::CoupledSolidBase const & solid =
+ getConstitutiveModel< geos::constitutive::CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
+
+ string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
+
+ if( m_isThermal )
+ {
+ thermalSinglePhaseReactiveBaseKernels::
+ AccumulationKernelFactory::
+ createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies,
+ dt,
+ dofManager.rankOffset(),
+ dofKey,
+ subRegion,
+ fluid,
+ solid,
+ localMatrix,
+ localRhs );
+ }
+ else
+ {
+ singlePhaseReactiveBaseKernels::
+ AccumulationKernelFactory::
+ createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies,
+ dt,
+ dofManager.rankOffset(),
+ dofKey,
+ subRegion,
+ fluid,
+ solid,
+ localMatrix,
+ localRhs );
+ }
+ } );
+ } );
+}
+
+void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt,
+ DomainPartition const & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+{
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel const & mesh,
+ string_array const & )
+ {
+ NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
+ FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
+ FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName );
+
+ string const & dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
+
+ fluxApprox.forAllStencils( mesh, [&] ( auto & stencil )
+ {
+ typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper();
+
+ if( m_isThermal )
+ {
+ thermalSinglePhaseReactiveFVMKernels::
+ FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies,
+ m_hasDiffusion,
+ dofManager.rankOffset(),
+ dofKey,
+ getName(),
+ mesh.getElemManager(),
+ stencilWrapper,
+ dt,
+ localMatrix.toViewConstSizes(),
+ localRhs.toView() );
+ }
+ else
+ {
+ singlePhaseReactiveFVMKernels::
+ FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies,
+ m_hasDiffusion,
+ dofManager.rankOffset(),
+ dofKey,
+ getName(),
+ mesh.getElemManager(),
+ stencilWrapper,
+ dt,
+ localMatrix.toViewConstSizes(),
+ localRhs.toView() );
+ }
+
+ // To add diffusion
+ } );
+ } );
+}
+
+SinglePhaseBase::FluidPropViews SinglePhaseReactiveTransport::getFluidProperties( constitutive::ConstitutiveBase const & fluid ) const
+{
+ ReactiveSingleFluid const & reactiveFluid = dynamicCast< ReactiveSingleFluid const & >( fluid );
+ return { reactiveFluid.density(),
+ reactiveFluid.dDensity(),
+ reactiveFluid.viscosity(),
+ reactiveFluid.dViscosity(),
+ reactiveFluid.defaultDensity(),
+ reactiveFluid.defaultViscosity() };
+}
+
+void SinglePhaseReactiveTransport::updateState( DomainPartition & domain )
+{
+ GEOS_MARK_FUNCTION;
+
+ SinglePhaseBase::updateState( domain );
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const,
+ auto & subRegion )
+ {
+ updateSpeciesAmount( subRegion );
+ } );
+ } );
+}
+
+void SinglePhaseReactiveTransport::updateSpeciesAmount( ElementSubRegionBase & subRegion ) const
+{
+ GEOS_MARK_FUNCTION;
+
+ arrayView2d< real64, compflow::USD_COMP > const primarySpeciesAggregateMole = subRegion.getField< fields::flow::primarySpeciesAggregateMole >();
+ arrayView2d< real64, compflow::USD_COMP > const primarySpeciesAggregateMole_n = subRegion.getField< fields::flow::primarySpeciesAggregateMole_n >();
+
+ CoupledSolidBase const & porousSolid =
+ getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
+ arrayView2d< real64 const > const porosity = porousSolid.getPorosity();
+ arrayView2d< real64 const > const porosity_n = porousSolid.getPorosity_n();
+
+ arrayView1d< real64 const > const volume = subRegion.getElementVolume();
+ arrayView1d< real64 > const deltaVolume = subRegion.getField< fields::flow::deltaVolume >();
+
+ ReactiveSingleFluid & fluid =
+ getConstitutiveModel< ReactiveSingleFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) );
+ arrayView2d< real64 const, compflow::USD_COMP > const primarySpeciesAggregateConcentration = fluid.primarySpeciesAggregateConcentration();
+ arrayView2d< real64 const, compflow::USD_COMP > const primarySpeciesAggregateConcentration_n = fluid.primarySpeciesAggregateConcentration_n();
+
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+ {
+ for( integer is = 0; is < m_numPrimarySpecies; ++is )
+ {
+ primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][is];
+
+ if( isZero( primarySpeciesAggregateMole_n[ei][is] ) )
+ primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][is];
+ }
+ } );
+}
+
+void SinglePhaseReactiveTransport::updateFluidModel( ObjectManagerBase & dataGroup ) const
+{
+ GEOS_MARK_FUNCTION;
+
+ arrayView1d< real64 const > const pres = dataGroup.getField< fields::flow::pressure >();
+ arrayView1d< real64 const > const temp = dataGroup.getField< fields::flow::temperature >();
+ arrayView2d< real64 const, compflow::USD_COMP > const logPrimaryConc = dataGroup.getField< fields::flow::logPrimarySpeciesConcentration >();
+
+ ReactiveSingleFluid & fluid =
+ getConstitutiveModel< ReactiveSingleFluid >( dataGroup, dataGroup.getReference< string >( viewKeyStruct::fluidNamesString() ) );
+
+ constitutive::constitutiveUpdatePassThru( fluid, [&]( auto & castedFluid )
+ {
+ typename TYPEOFREF( castedFluid ) ::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper();
+ singlePhaseReactiveBaseKernels::FluidUpdateKernel::launch( fluidWrapper, pres, temp, logPrimaryConc );
+ } );
+}
+
+void SinglePhaseReactiveTransport::initializeFluidState( MeshLevel & mesh, string_array const & regionNames )
+{
+ mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const,
+ auto & subRegion )
+ {
+ ReactiveSingleFluid const & fluid =
+ getConstitutiveModel< ReactiveSingleFluid >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString()));
+ updateFluidState( subRegion );
+
+ // 2. save the initial density (for use in the single-phase poromechanics solver to compute the deltaBodyForce)
+ fluid.initializeState();
+
+ SinglePhaseBase::updateMass( subRegion );
+ updateSpeciesAmount( subRegion );
+
+ // If the diffusion is supported, initialize the model
+ if( m_hasDiffusion )
+ {
+ string const & diffusionName = subRegion.template getReference< string >( viewKeyStruct::diffusionNamesString() );
+ DiffusionBase const & diffusionMaterial = getConstitutiveModel< DiffusionBase >( subRegion, diffusionName );
+ arrayView1d< real64 const > const temperature = subRegion.template getField< fields::flow::temperature >();
+ diffusionMaterial.initializeTemperatureState( temperature );
+ }
+ } );
+}
+
+void SinglePhaseReactiveTransport::initializePostInitialConditionsPreSubGroups()
+{
+ GEOS_MARK_FUNCTION;
+
+ FlowSolverBase::initializePostInitialConditionsPreSubGroups();
+
+ DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ FieldIdentifiers fieldsToBeSync;
+ fieldsToBeSync.addElementFields( { fields::flow::logPrimarySpeciesConcentration::key() },
+ regionNames );
+
+ CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false );
+ } );
+
+ FlowSolverBase::initializeState( domain );
+}
+
+void SinglePhaseReactiveTransport::applyBoundaryConditions( real64 const time_n,
+ real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+{
+ GEOS_MARK_FUNCTION;
+
+ // apply pressure boundary conditions.
+ applyDirichletBC( time_n, dt, domain, dofManager, localMatrix.toViewConstSizes(), localRhs.toView() );
+
+ // // apply flux boundary conditions (To finish)
+ // applySourceFluxBC( time_n, dt, dofManager, domain, localMatrix.toViewConstSizes(), localRhs.toView() );
+
+ // // apply aquifer boundary conditions (To finish)
+ // applyAquiferBC( time_n, dt, dofManager, domain, localMatrix.toViewConstSizes(), localRhs.toView() );
+}
+
+// // To finish
+// void SinglePhaseReactiveTransport::applySourceFluxBC( real64 const time,
+// real64 const dt,
+// DofManager const & dofManager,
+// DomainPartition & domain,
+// CRSMatrixView< real64, globalIndex const > const & localMatrix,
+// arrayView1d< real64 > const & localRhs ) const
+// {
+// GEOS_MARK_FUNCTION;
+
+// FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance();
+
+// string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
+
+// // Step 1: count individual source flux boundary conditions
+
+// std::map< string, localIndex > bcNameToBcId;
+// localIndex bcCounter = 0;
+
+// fsManager.forSubGroups< SourceFluxBoundaryCondition >( [&] ( SourceFluxBoundaryCondition const & bc )
+// {
+// // collect all the bc names to idx
+// bcNameToBcId[bc.getName()] = bcCounter;
+// bcCounter++;
+// } );
+
+// if( bcCounter == 0 )
+// {
+// return;
+// }
+
+// // Step 2: count the set size for each source flux (each source flux may have multiple target sets)
+
+// array1d< globalIndex > bcAllSetsSize( bcNameToBcId.size() );
+
+// computeSourceFluxSizeScalingFactor( time_n,
+// dt,
+// domain,
+// bcNameToBcId,
+// bcAllSetsSize.toView() );
+
+// // Step 3: we are ready to impose the boundary condition, normalized by the set size
+
+// forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+// MeshLevel & mesh,
+// arrayView1d< string const > const & )
+// {
+// integer const isThermal = m_isThermal;
+
+// fsManager.apply< ElementSubRegionBase,
+// SourceFluxBoundaryCondition >( time + dt,
+// mesh,
+// SourceFluxBoundaryCondition::catalogName(),
+// [&, isThermal]( SourceFluxBoundaryCondition const & fs,
+// string const & setName,
+// SortedArrayView< localIndex const > const & targetSet,
+// ElementSubRegionBase & subRegion,
+// string const & )
+// {
+// if( fs.getLogLevel() >= 1 && m_nonlinearSolverParameters.m_numNewtonIterations == 0 )
+// {
+// globalIndex const numTargetElems = MpiWrapper::sum< globalIndex >( targetSet.size() );
+// GEOS_LOG_RANK_0( GEOS_FMT( bcLogMessage,
+// getName(), time+dt, fs.getCatalogName(), fs.getName(),
+// setName, subRegion.getName(), fs.getScale(), numTargetElems ) );
+// }
+
+// if( targetSet.size() == 0 )
+// {
+// return;
+// }
+// if( !subRegion.hasWrapper( dofKey ) )
+// {
+// if( fs.getLogLevel() >= 1 )
+// {
+// GEOS_LOG_RANK( GEOS_FMT( "{}: trying to apply SourceFlux, but its targetSet named '{}' intersects with non-simulated region
+// named '{}'.",
+// getDataContext(), setName, subRegion.getName() ) );
+// }
+// return;
+// }
+
+// arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
+// arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+
+// // Step 3.1: get the values of the source boundary condition that need to be added to the rhs
+
+// array1d< globalIndex > dofArray( targetSet.size() );
+// array1d< real64 > rhsContributionArray( targetSet.size() );
+// arrayView1d< real64 > rhsContributionArrayView = rhsContributionArray.toView();
+// localIndex const rankOffset = dofManager.rankOffset();
+
+// RAJA::ReduceSum< parallelDeviceReduce, real64 > massProd( 0.0 );
+
+// // note that the dofArray will not be used after this step (simpler to use dofNumber instead)
+// fs.computeRhsContribution< FieldSpecificationAdd,
+// parallelDevicePolicy<> >( targetSet.toViewConst(),
+// time + dt,
+// dt,
+// subRegion,
+// dofNumber,
+// rankOffset,
+// localMatrix,
+// dofArray.toView(),
+// rhsContributionArrayView,
+// [] GEOS_HOST_DEVICE ( localIndex const )
+// {
+// return 0.0;
+// } );
+
+// // Step 3.2: we are ready to add the right-hand side contributions, taking into account our equation layout
+
+// // get the normalizer
+// real64 const sizeScalingFactor = bcAllSetsSize[bcNameToBcId.at( fs.getName())];
+
+// if( isThermal )
+// {
+
+// }
+// else
+// {
+// integer const fluidComponentId = fs.getComponent();
+// integer const numFluidSpecies = m_numPrimarySpecies;
+// forAll< parallelDevicePolicy<> >( targetSet.size(), [sizeScalingFactor,
+// targetSet,
+// rankOffset,
+// ghostRank,
+// fluidComponentId,
+// numFluidSpecies,
+// dofNumber,
+// rhsContributionArrayView,
+// localRhs,
+// massProd] GEOS_HOST_DEVICE ( localIndex const a )
+// {
+// // we need to filter out ghosts here, because targetSet may contain them
+// localIndex const ei = targetSet[a];
+// if( ghostRank[ei] >= 0 )
+// {
+// return;
+// }
+
+// real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor
+// here!
+// massProd += rhsValue;
+
+// globalIndex const totalMassBalanceRow = dofNumber[ei] - rankOffset;
+// globalIndex const speciesMassBalanceRow = dofNumber[ei] - rankOffset + fluidComponentId + 1;
+// localRhs[totalMassBalanceRow] += rhsValue;
+// } );
+// }
+// } );
+// } );
+// }
+
+namespace
+{
+char const bcLogMessage[] =
+ "SinglePhaseReactiveTransport {}: at time {}s, "
+ "the <{}> boundary condition '{}' is applied to the element set '{}' in subRegion '{}'. "
+ "\nThe scale of this boundary condition is {} and multiplies the value of the provided function (if any). "
+ "\nThe total number of target elements (including ghost elements) is {}. "
+ "\nNote that if this number is equal to zero for all subRegions, the boundary condition will not be applied on this element set.";
+}
+
+void SinglePhaseReactiveTransport::applyDirichletBC( real64 const time_n,
+ real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) const
+{
+ FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance();
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & )
+ {
+ // 1. Apply pressure Dirichlet BCs, store in a separate field
+ applyFieldValue< ElementSubRegionBase >( time_n, dt, mesh, bcLogMessage,
+ fields::flow::pressure::key(), fields::flow::bcPressure::key() );
+ // 2. Apply primary species BC (log promary species concentration) and store them for constitutive call
+ applyFieldValue< ElementSubRegionBase >( time_n, dt, mesh, bcLogMessage,
+ fields::flow::logPrimarySpeciesConcentration::key(), fields::flow::bcLogPrimarySpeciesConcentration::key() );
+ // 3. Apply temperature Dirichlet BCs, store in a separate field
+ if( m_isThermal )
+ {
+ applyFieldValue< ElementSubRegionBase >( time_n, dt, mesh, bcLogMessage,
+ fields::flow::temperature::key(), fields::flow::bcTemperature::key() );
+ }
+
+ globalIndex const rankOffset = dofManager.rankOffset();
+ string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() );
+
+ // 4. Apply pressure and log primary species concentration to the system
+ fsManager.apply< ElementSubRegionBase >( time_n + dt,
+ mesh,
+ fields::flow::pressure::key(),
+ [&] ( FieldSpecificationBase const &,
+ string const &,
+ SortedArrayView< localIndex const > const & targetSet,
+ ElementSubRegionBase & subRegion,
+ string const & )
+ {
+ arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+ arrayView1d< globalIndex const > const dofNumber =
+ subRegion.getReference< array1d< globalIndex > >( dofKey );
+
+ // in the isothermal case, we use the reservoir temperature to enforce the boundary condition
+ // in the thermal case, the validation function guarantees that temperature has been provided
+ arrayView1d< real64 const > const bcPres =
+ subRegion.getReference< array1d< real64 > >( fields::flow::bcPressure::key() );
+ arrayView2d< real64 const, compflow::USD_COMP > const bcLogPrimaryConc =
+ subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( fields::flow::bcLogPrimarySpeciesConcentration::key() );
+
+ arrayView1d< real64 const > const pres =
+ subRegion.getReference< array1d< real64 > >( fields::flow::pressure::key() );
+ arrayView2d< real64 const, compflow::USD_COMP > const logPrimaryConc =
+ subRegion.getReference< array2d< real64, compflow::LAYOUT_COMP > >( fields::flow::logPrimarySpeciesConcentration::key() );
+
+ integer const numPrimarySpecies = m_numPrimarySpecies;
+ forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
+ {
+ localIndex const ei = targetSet[a];
+ if( ghostRank[ei] >= 0 )
+ {
+ return;
+ }
+
+ globalIndex const dofIndex = dofNumber[ei];
+ localIndex const localRow = dofIndex - rankOffset;
+ real64 rhsValue;
+
+ // 4.1. Apply pressure value to the matrix/rhs
+ FieldSpecificationEqual::SpecifyFieldValue( dofIndex,
+ rankOffset,
+ localMatrix,
+ rhsValue,
+ bcPres[ei],
+ pres[ei] );
+ localRhs[localRow] = rhsValue;
+
+ integer const speciesDofBeginIndex = m_isThermal? 2:1;
+
+ // 4.2. For each component, apply target global density value
+ for( integer is = 0; is < numPrimarySpecies; ++is )
+ {
+ FieldSpecificationEqual::SpecifyFieldValue( dofIndex + is + speciesDofBeginIndex,
+ rankOffset,
+ localMatrix,
+ rhsValue,
+ bcLogPrimaryConc[ei][is],
+ logPrimaryConc[ei][is] );
+ localRhs[localRow + is + speciesDofBeginIndex] = rhsValue;
+ }
+ } );
+ } );
+
+ // 5. Apply temperature to the system
+ if( m_isThermal )
+ {
+ fsManager.apply< ElementSubRegionBase >( time_n + dt,
+ mesh,
+ fields::flow::temperature::key(),
+ [&] ( FieldSpecificationBase const &,
+ string const &,
+ SortedArrayView< localIndex const > const & targetSet,
+ ElementSubRegionBase & subRegion,
+ string const & )
+ {
+ arrayView1d< integer const > const ghostRank =
+ subRegion.getReference< array1d< integer > >( ObjectManagerBase::viewKeyStruct::ghostRankString() );
+ arrayView1d< globalIndex const > const dofNumber =
+ subRegion.getReference< array1d< globalIndex > >( dofKey );
+ arrayView1d< real64 const > const bcTemp =
+ subRegion.getReference< array1d< real64 > >( fields::flow::bcTemperature::key() );
+ arrayView1d< real64 const > const temp =
+ subRegion.getReference< array1d< real64 > >( fields::flow::temperature::key() );
+
+ forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
+ {
+ localIndex const ei = targetSet[a];
+ if( ghostRank[ei] >= 0 )
+ {
+ return;
+ }
+
+ globalIndex const dofIndex = dofNumber[ei];
+ localIndex const localRow = dofIndex - rankOffset;
+ real64 rhsValue;
+
+ // 4.2. Apply temperature value to the matrix/rhs
+ FieldSpecificationEqual::SpecifyFieldValue( dofIndex + 1,
+ rankOffset,
+ localMatrix,
+ rhsValue,
+ bcTemp[ei],
+ temp[ei] );
+ localRhs[localRow + 1] = rhsValue;
+ } );
+ } );
+ }
+ } );
+}
+
+real64 SinglePhaseReactiveTransport::calculateResidualNorm( real64 const & GEOS_UNUSED_PARAM( time_n ),
+ real64 const & GEOS_UNUSED_PARAM( dt ),
+ DomainPartition const & domain,
+ DofManager const & dofManager,
+ arrayView1d< real64 const > const & localRhs )
+{
+ GEOS_MARK_FUNCTION;
+
+ integer constexpr numNorm = 3; // total mass balance, energy balance, and species mass balance
+ array1d< real64 > localResidualNorm;
+ array1d< real64 > localResidualNormalizer;
+ localResidualNorm.resize( numNorm );
+ localResidualNormalizer.resize( numNorm );
+
+ physicsSolverBaseKernels::NormType const normType = SinglePhaseBase::getNonlinearSolverParameters().normType();
+
+ globalIndex const rankOffset = dofManager.rankOffset();
+ string const dofKey = dofManager.getKey( SinglePhaseBase::viewKeyStruct::elemDofFieldString() );
+
+ this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
+ MeshLevel const & mesh,
+ string_array const & regionNames )
+ {
+ mesh.getElemManager().forElementSubRegions( regionNames,
+ [&]( localIndex const,
+ ElementSubRegionBase const & subRegion )
+ {
+ real64 subRegionResidualNorm[numNorm]{};
+ real64 subRegionResidualNormalizer[numNorm]{};
+
+ // step 1: compute the norm in the subRegion
+
+ if( m_isThermal )
+ {
+ singlePhaseReactiveBaseKernels::
+ ResidualNormKernelFactory::
+ createAndLaunch< parallelDevicePolicy<> >( normType,
+ m_numPrimarySpecies,
+ rankOffset,
+ dofKey,
+ localRhs,
+ subRegion,
+ m_nonlinearSolverParameters.m_minNormalizer,
+ subRegionResidualNorm,
+ subRegionResidualNormalizer );
+ }
+ else
+ {
+ real64 subRegionFlowResidualNorm[2]{};
+ real64 subRegionFlowResidualNormalizer[2]{};
+
+ singlePhaseReactiveBaseKernels::
+ ResidualNormKernelFactory::
+ createAndLaunch< parallelDevicePolicy<> >( normType,
+ m_numPrimarySpecies,
+ rankOffset,
+ dofKey,
+ localRhs,
+ subRegion,
+ m_nonlinearSolverParameters.m_minNormalizer,
+ subRegionFlowResidualNorm,
+ subRegionFlowResidualNormalizer );
+ subRegionResidualNorm[0] = subRegionFlowResidualNorm[0];
+ subRegionResidualNorm[1] = subRegionFlowResidualNorm[1];
+ subRegionResidualNormalizer[0] = subRegionFlowResidualNormalizer[0];
+ subRegionResidualNormalizer[1] = subRegionFlowResidualNormalizer[1];
+ }
+
+ // step 2: first reduction across meshBodies/regions/subRegions
+
+ if( normType == physicsSolverBaseKernels::NormType::Linf )
+ {
+ physicsSolverBaseKernels::LinfResidualNormHelper::
+ updateLocalNorm< numNorm >( subRegionResidualNorm, localResidualNorm );
+ }
+ else
+ {
+ physicsSolverBaseKernels::L2ResidualNormHelper::
+ updateLocalNorm< numNorm >( subRegionResidualNorm, subRegionResidualNormalizer, localResidualNorm, localResidualNormalizer );
+ }
+ } );
+ } );
+
+ // step 3: second reduction across MPI ranks
+
+ real64 residualNorm = 0.0;
+ array1d< real64 > globalResidualNorm;
+ globalResidualNorm.resize( numNorm );
+ if( m_isThermal )
+ {
+ if( normType == physicsSolverBaseKernels::NormType::Linf )
+ {
+ physicsSolverBaseKernels::LinfResidualNormHelper::
+ computeGlobalNorm( localResidualNorm, globalResidualNorm );
+ }
+ else
+ {
+ physicsSolverBaseKernels::L2ResidualNormHelper::
+ computeGlobalNorm( localResidualNorm, localResidualNormalizer, globalResidualNorm );
+ }
+ residualNorm = sqrt( globalResidualNorm[0] * globalResidualNorm[0] + globalResidualNorm[1] * globalResidualNorm[1] + globalResidualNorm[2] * globalResidualNorm[2] );
+
+ GEOS_LOG_LEVEL_RANK_0_NLR( logInfo::Convergence, GEOS_FMT( " ( RtotalMass RspeciesAmount ) = ( {:4.2e} {:4.2e} ) ( Renergy ) = ( {:4.2e} )",
+ globalResidualNorm[0], globalResidualNorm[2], globalResidualNorm[1] ));
+ }
+ else
+ {
+ if( normType == physicsSolverBaseKernels::NormType::Linf )
+ {
+ physicsSolverBaseKernels::LinfResidualNormHelper::
+ computeGlobalNorm( localResidualNorm, globalResidualNorm );
+ }
+ else
+ {
+ physicsSolverBaseKernels::L2ResidualNormHelper::
+ computeGlobalNorm( localResidualNorm, localResidualNormalizer, globalResidualNorm );
+ }
+ residualNorm = sqrt( globalResidualNorm[0] * globalResidualNorm[0] + globalResidualNorm[1] * globalResidualNorm[1] );
+
+ GEOS_LOG_LEVEL_RANK_0_NLR( logInfo::Convergence, GEOS_FMT( " ( RtotalMass RspeciesAmount ) = ( {:4.2e} {:4.2e} )",
+ globalResidualNorm[0], globalResidualNorm[1] ) );
+ }
+ return residualNorm;
+}
+
+void SinglePhaseReactiveTransport::applySystemSolution( DofManager const & dofManager,
+ arrayView1d< real64 const > const & localSolution,
+ real64 const scalingFactor,
+ real64 const dt,
+ DomainPartition & domain )
+{
+ GEOS_UNUSED_VAR( dt );
+
+ if( m_isThermal )
+ {
+ DofManager::CompMask pressureMask( m_numDofPerCell, 0, 1 );
+ DofManager::CompMask temperatureMask( m_numDofPerCell, 1, 2 );
+ DofManager::CompMask speciesMask( m_numDofPerCell, 2, m_numPrimarySpecies+2 );
+
+ dofManager.addVectorToField( localSolution,
+ viewKeyStruct::elemDofFieldString(),
+ fields::flow::pressure::key(),
+ scalingFactor,
+ pressureMask );
+
+ dofManager.addVectorToField( localSolution,
+ viewKeyStruct::elemDofFieldString(),
+ fields::flow::temperature::key(),
+ scalingFactor,
+ temperatureMask );
+
+ dofManager.addVectorToField( localSolution,
+ viewKeyStruct::elemDofFieldString(),
+ fields::flow::logPrimarySpeciesConcentration::key(),
+ scalingFactor,
+ speciesMask );
+ }
+ else
+ {
+ DofManager::CompMask pressureMask( m_numDofPerCell, 0, 1 );
+ DofManager::CompMask speciesMask( m_numDofPerCell, 1, m_numPrimarySpecies+1 );
+
+ dofManager.addVectorToField( localSolution,
+ viewKeyStruct::elemDofFieldString(),
+ fields::flow::pressure::key(),
+ scalingFactor,
+ pressureMask );
+
+ dofManager.addVectorToField( localSolution,
+ viewKeyStruct::elemDofFieldString(),
+ fields::flow::logPrimarySpeciesConcentration::key(),
+ scalingFactor,
+ speciesMask );
+ }
+
+ this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ std::vector< string > fields{ fields::flow::pressure::key() };
+
+ if( m_isThermal )
+ {
+ fields.emplace_back( fields::flow::temperature::key() );
+ }
+
+ fields.emplace_back( fields::flow::logPrimarySpeciesConcentration::key() );
+
+ FieldIdentifiers fieldsToBeSync;
+ fieldsToBeSync.addElementFields( fields, regionNames );
+
+ CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), true );
+ } );
+}
+
+void SinglePhaseReactiveTransport::saveConvergedState( ElementSubRegionBase & subRegion ) const
+{
+ SinglePhaseBase::saveConvergedState( subRegion );
+
+ arrayView2d< real64 const, compflow::USD_COMP > const primarySpeciesAggregateMole = subRegion.template getField< fields::flow::primarySpeciesAggregateMole >();
+ arrayView2d< real64, compflow::USD_COMP > const primarySpeciesAggregateMole_n = subRegion.template getField< fields::flow::primarySpeciesAggregateMole_n >();
+ primarySpeciesAggregateMole_n.setValues< parallelDevicePolicy<> >( primarySpeciesAggregateMole );
+}
+
+void SinglePhaseReactiveTransport::assembleEDFMFluxTerms( real64 const GEOS_UNUSED_PARAM( time_n ),
+ real64 const GEOS_UNUSED_PARAM( dt ),
+ DomainPartition const & GEOS_UNUSED_PARAM( domain ),
+ DofManager const & GEOS_UNUSED_PARAM( dofManager ),
+ CRSMatrixView< real64, globalIndex const > const & GEOS_UNUSED_PARAM( localMatrix ),
+ arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ),
+ string const & GEOS_UNUSED_PARAM( jumpDofKey ) )
+{}
+
+void SinglePhaseReactiveTransport::applyAquiferBC( real64 const GEOS_UNUSED_PARAM( time ),
+ real64 const GEOS_UNUSED_PARAM( dt ),
+ DomainPartition & GEOS_UNUSED_PARAM( domain ),
+ DofManager const & GEOS_UNUSED_PARAM( dofManager ),
+ CRSMatrixView< real64, globalIndex const > const & GEOS_UNUSED_PARAM( localMatrix ),
+ arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) ) const
+{}
+
+void SinglePhaseReactiveTransport::assembleStabilizedFluxTerms( real64 const GEOS_UNUSED_PARAM( dt ),
+ DomainPartition const & GEOS_UNUSED_PARAM( domain ),
+ DofManager const & GEOS_UNUSED_PARAM( dofManager ),
+ CRSMatrixView< real64, globalIndex const > const & GEOS_UNUSED_PARAM( localMatrix ),
+ arrayView1d< real64 > const & GEOS_UNUSED_PARAM( localRhs ) )
+{}
+
+REGISTER_CATALOG_ENTRY( PhysicsSolverBase, SinglePhaseReactiveTransport, string const &, Group * const )
+
+} /* namespace geos */
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp
new file mode 100644
index 00000000000..61aa4732fff
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp
@@ -0,0 +1,291 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file SinglePhaseReactiveTransport.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVETRANSPORT_HPP_
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVETRANSPORT_HPP_
+
+#include "fieldSpecification/FieldSpecificationManager.hpp"
+#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp"
+#include "physicsSolvers/fluidFlow/SinglePhaseFVM.hpp"
+#include "physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluidUpdateKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.hpp"
+
+
+namespace geos
+{
+
+/**
+ * @class SinglePhaseReactiveTransport
+ *
+ * A solver for single phase reactive transport
+ */
+class SinglePhaseReactiveTransport : public SinglePhaseBase
+{
+
+public:
+
+ /**
+ * @brief main constructor for Group Objects
+ * @param name the name of this instantiation of Group in the repository
+ * @param parent the parent group of this instantiation of Group
+ */
+ SinglePhaseReactiveTransport( const string & name,
+ dataRepository::Group * const parent );
+
+ SinglePhaseReactiveTransport() = delete;
+
+ /// deleted copy constructor
+ SinglePhaseReactiveTransport( SinglePhaseReactiveTransport const & ) = delete;
+
+ /// default move constructor
+ SinglePhaseReactiveTransport( SinglePhaseReactiveTransport && ) = default;
+
+ /// deleted assignment operator
+ SinglePhaseReactiveTransport & operator=( SinglePhaseReactiveTransport const & ) = delete;
+
+ /// deleted move operator
+ SinglePhaseReactiveTransport & operator=( SinglePhaseReactiveTransport && ) = delete;
+
+ /**
+ * @brief default destructor
+ */
+ virtual ~SinglePhaseReactiveTransport() override = default;
+
+ /**
+ * @brief name of the node manager in the object catalog
+ * @return string that contains the catalog name to generate a new NodeManager object through the object catalog.
+ */
+ static string catalogName()
+ {
+ return "SinglePhaseReactiveTransport";
+ }
+
+ /**
+ * @copydoc PhysicsSolverBase::getCatalogName()
+ */
+ string getCatalogName() const override { return catalogName(); }
+
+ virtual void registerDataOnMesh( dataRepository::Group & meshBodies ) override;
+
+ virtual void
+ setupDofs( DomainPartition const & domain,
+ DofManager & dofManager ) const override;
+
+ virtual void
+ assembleSystem( real64 const time_n,
+ real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) override;
+
+ virtual void
+ applyBoundaryConditions( real64 const time_n,
+ real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) override;
+
+ virtual real64
+ calculateResidualNorm( real64 const & time_n,
+ real64 const & dt,
+ DomainPartition const & domain,
+ DofManager const & dofManager,
+ arrayView1d< real64 const > const & localRhs ) override;
+
+ virtual void
+ applySystemSolution( DofManager const & dofManager,
+ arrayView1d< real64 const > const & localSolution,
+ real64 const scalingFactor,
+ real64 const dt,
+ DomainPartition & domain ) override;
+
+ virtual void
+ resetStateToBeginningOfStep( DomainPartition & domain ) override;
+
+ virtual void
+ implicitStepComplete( real64 const & time,
+ real64 const & dt,
+ DomainPartition & domain ) override final;
+
+ virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override final;
+
+ virtual void
+ updateState ( DomainPartition & domain ) override final;
+
+ void updateSpeciesAmount( ElementSubRegionBase & subRegion ) const;
+
+ virtual void updateFluidModel( ObjectManagerBase & dataGroup ) const override;
+
+ virtual void initializePostInitialConditionsPreSubGroups() override;
+
+ virtual void initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) override;
+
+ /**
+ * @brief assembles the accumulation terms in total mass balance and primary species amount equation for all cells
+ * @param dt time step
+ * @param domain the physical domain object
+ * @param dofManager degree-of-freedom manager associated with the linear system
+ * @param localMatrix the system matrix
+ * @param localRhs the system right-hand side vector
+ */
+ void assembleAccumulationTermsInMassBalanceAndSpeciesAmountEqs( real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) const;
+
+ /**
+ * @brief assembles the flux terms for all cells
+ * @param dt time step
+ * @param domain the physical domain object
+ * @param dofManager degree-of-freedom manager associated with the linear system
+ * @param matrix the system matrix
+ * @param rhs the system right-hand side vector
+ */
+ virtual void
+ assembleFluxTerms( real64 const dt,
+ DomainPartition const & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) override;
+
+ /**
+ * @brief Function to perform the Application of Dirichlet type BC's
+ * @param time current time
+ * @param dt time step
+ * @param dofManager degree-of-freedom manager associated with the linear system
+ * @param domain the domain
+ * @param localMatrix local system matrix
+ * @param localRhs local system right-hand side vector
+ */
+ virtual void
+ applyDirichletBC( real64 const time_n,
+ real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) const override;
+
+ /**
+ * @brief Utility function that encapsulates the call to FieldSpecificationBase::applyFieldValue in BC application
+ * @param[in] time_n the time at the beginning of the step
+ * @param[in] dt the time step
+ * @param[in] mesh the mesh level object
+ * @param[in] logMessage the log message issued by the solver if the bc is called
+ * @param[in] fieldKey the key of the field specified in the xml file
+ * @param[in] boundaryFieldKey the key of the boundary field
+ */
+ template< typename OBJECT_TYPE >
+ void applyFieldValue( real64 const & time_n,
+ real64 const & dt,
+ MeshLevel & mesh,
+ char const logMessage[],
+ string const fieldKey,
+ string const boundaryFieldKey ) const;
+
+ virtual void
+ applyAquiferBC( real64 const time,
+ real64 const dt,
+ DomainPartition & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) const override;
+
+ virtual void
+ assembleEDFMFluxTerms( real64 const time_n,
+ real64 const dt,
+ DomainPartition const & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs,
+ string const & jumpDofKey ) override final;
+
+ virtual void
+ assembleStabilizedFluxTerms( real64 const dt,
+ DomainPartition const & domain,
+ DofManager const & dofManager,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs ) override;
+
+ struct viewKeyStruct : SinglePhaseBase::viewKeyStruct
+ {
+ static constexpr char const * diffusionNamesString() { return "diffusionNames"; }
+ };
+
+protected:
+
+ virtual FluidPropViews getFluidProperties( constitutive::ConstitutiveBase const & fluid ) const override;
+
+ /// the number of primary species in the fluid
+ integer m_numPrimarySpecies;
+
+ /// name of the reactive fluid constitutive model
+ string m_reactiveFluidModelName;
+
+ /// flag to determine whether or not to apply diffusion
+ integer m_hasDiffusion;
+};
+
+template< typename OBJECT_TYPE >
+void SinglePhaseReactiveTransport::applyFieldValue( real64 const & time_n,
+ real64 const & dt,
+ MeshLevel & mesh,
+ char const logMessage[],
+ string const fieldKey,
+ string const boundaryFieldKey ) const
+{
+ FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance();
+
+ fsManager.apply< OBJECT_TYPE >( time_n + dt,
+ mesh,
+ fieldKey,
+ [&]( FieldSpecificationBase const & fs,
+ string const & setName,
+ SortedArrayView< localIndex const > const & lset,
+ OBJECT_TYPE & targetGroup,
+ string const & )
+ {
+ if( fs.getLogLevel() >= 1 && m_nonlinearSolverParameters.m_numNewtonIterations == 0 )
+ {
+ globalIndex const numTargetElems = MpiWrapper::sum< globalIndex >( lset.size() );
+ GEOS_LOG_RANK_0( GEOS_FMT( logMessage,
+ getName(), time_n+dt, fs.getCatalogName(), fs.getName(),
+ setName, targetGroup.getName(), fs.getScale(), numTargetElems ) );
+ }
+
+ // Specify the bc value of the field
+ fs.applyFieldValue< FieldSpecificationEqual,
+ parallelDevicePolicy<> >( lset,
+ time_n + dt,
+ targetGroup,
+ boundaryFieldKey );
+ } );
+}
+
+} /* namespace geos */
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVETRANSPORT_HPP_
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp
new file mode 100644
index 00000000000..3aa6c8f8cd4
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp
@@ -0,0 +1,92 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file SinglePhaseReactiveTransportFields.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVETRANSPORTFIELDS_HPP_
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVETRANSPORTFIELDS_HPP_
+
+#include "mesh/MeshFields.hpp"
+
+namespace geos
+{
+/**
+ * A scope for field traits.
+ */
+namespace fields
+{
+
+namespace flow
+{
+using array2dLayoutComp = array2d< real64, compflow::LAYOUT_COMP >;
+using array2dLayoutFluid_dC = array2d< real64, compflow::LAYOUT_FLUID_DC >;
+
+DECLARE_FIELD( logPrimarySpeciesConcentration,
+ "logPrimarySpeciesConcentration",
+ array2dLayoutComp,
+ 0,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "Natural log of primary species concentration (molarity)" );
+
+DECLARE_FIELD( logPrimarySpeciesConcentration_n,
+ "logPrimarySpeciesConcentration_n",
+ array2dLayoutComp,
+ 0,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "Natural log of primary species concentration (molarity) at the previous converged time step" );
+
+DECLARE_FIELD( bcLogPrimarySpeciesConcentration,
+ "bcLogPrimarySpeciesConcentration",
+ array2dLayoutComp,
+ 0,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "Boundary condition for natural log of primary species concentration (molarity)" );
+
+DECLARE_FIELD( primarySpeciesAggregateMole,
+ "primarySpeciesAggregateMole",
+ array2dLayoutComp,
+ 0,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "Aggregate amount of primary species in mole" );
+
+DECLARE_FIELD( primarySpeciesAggregateMole_n,
+ "primarySpeciesAggregateMole_n",
+ array2dLayoutComp,
+ 0,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "Aggregate amount of primary species in mole at the previous converged time step" );
+
+DECLARE_FIELD( dMobility_dLogPrimaryConc,
+ "dMobility_dLogPrimaryConc",
+ array2dLayoutFluid_dC,
+ 0,
+ NOPLOT,
+ NO_WRITE,
+ "Derivative of fluid mobility with respect to log of primary species concentration" );
+
+}
+
+}
+
+}
+
+#endif // GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVETRANSPORTFIELDS_HPP_
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp
index b80166eb1b8..bf4800761b3 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp
@@ -56,8 +56,8 @@ class AccumulationKernel
/// Note: Derivative lineup only supports dP & dT, not component terms
- static constexpr integer isThermal = NUM_DOF-1;
- using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< isThermal >;
+ // static constexpr integer isThermal = NUM_DOF-1;
+ using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 0 >; // Why not just input 0
/**
* @brief Constructor
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp
index e2e92b14792..c470b9790a1 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp
@@ -51,8 +51,8 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU
using Base::m_dMass;
/// Note: Derivative lineup only supports dP & dT, not component terms
- static constexpr integer isThermal = NUM_DOF-1;
- using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< isThermal >;
+ // static constexpr integer isThermal = NUM_DOF-1;
+ using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
/**
* @brief Constructor
* @param[in] rankOffset the offset of my MPI rank
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp
new file mode 100644
index 00000000000..188939ac350
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp
@@ -0,0 +1,331 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file AccumulationKernels.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_ACCUMULATIONKERNELS_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_ACCUMULATIONKERNELS_HPP
+
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.hpp"
+#include "constitutive/solid/CoupledSolidBase.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp"
+
+namespace geos
+{
+
+namespace singlePhaseReactiveBaseKernels
+{
+
+/******************************** AccumulationKernel ********************************/
+
+/**
+ * @class AccumulationKernel
+ * @brief Define the interface for the assembly kernel in charge of accumulation
+ */
+template< typename SUBREGION_TYPE, integer NUM_DOF, integer NUM_SPECIES >
+class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SUBREGION_TYPE, NUM_DOF >
+{
+
+public:
+
+ using Base = singlePhaseBaseKernels::AccumulationKernel< SUBREGION_TYPE, NUM_DOF >;
+ using Base::numDof;
+ using Base::numEqn;
+ using Base::m_rankOffset;
+ using Base::m_dofNumber;
+ using Base::m_elemGhostRank;
+ using Base::m_localMatrix;
+ using Base::m_localRhs;
+ using Base::m_dMass;
+
+ /// Note: Derivative lineup only supports dP & dT, not component terms
+ using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 0 >;
+
+ /// Compile time value for the number of primary species
+ static constexpr integer numSpecies = NUM_SPECIES;
+
+ /**
+ * @brief Constructor
+ * @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
+ */
+ AccumulationKernel( globalIndex const rankOffset,
+ string const dofKey,
+ SUBREGION_TYPE const & subRegion,
+ constitutive::ReactiveSingleFluid const & fluid,
+ constitutive::CoupledSolidBase const & solid,
+ real64 const & dt,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+ : Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ),
+ m_dt( dt ),
+ m_volume( subRegion.getElementVolume() ),
+ m_deltaVolume( subRegion.template getField< fields::flow::deltaVolume >() ),
+ m_porosity( solid.getPorosity() ),
+ m_dPoro_dPres( solid.getDporosity_dPressure() ),
+ // m_dDensity_dLogPrimaryConc( fluid.dDensity_dLogPrimaryConc() ),
+ // m_dPoro_dLogPrimaryConc( solid.getDporosity_dLogPrimaryConc() ),
+ m_primarySpeciesAggregateConcentration( fluid.primarySpeciesAggregateConcentration() ),
+ // m_dPrimarySpeciesAggregateConcentration_dPres( fluid.dPrimarySpeciesAggregateConcentration_dPres() ),
+ m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc( fluid.dPrimarySpeciesAggregateConcentration_dLogPrimaryConc() ),
+ m_primarySpeciesAggregateKineticRate( fluid.kineticReactionRates() ),
+ // m_dPrimarySpeciesAggregateKineticRate_dPres( fluid.dPrimarySpeciesAggregateKineticRate_dPres() ),
+ // m_dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc( fluid.dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc() ),
+ m_primarySpeciesAggregateMole_n( subRegion.template getField< fields::flow::primarySpeciesAggregateMole_n >() )
+ {}
+
+ /**
+ * @struct StackVariables
+ * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack
+ */
+ struct StackVariables : public Base::StackVariables
+ {
+public:
+
+ GEOS_HOST_DEVICE
+ StackVariables()
+ : Base::StackVariables()
+ {}
+
+ using Base::StackVariables::localRow;
+ using Base::StackVariables::dofIndices;
+ using Base::StackVariables::localResidual;
+ using Base::StackVariables::localJacobian;
+
+ // Pore volume information
+
+ /// Pore volume at time n+1
+ real64 poreVolume = 0.0;
+
+ /// Derivative of pore volume with respect to pressure
+ real64 dPoreVolume_dPres = 0.0;
+
+ /// Derivative of pore volume with respect to each primary species concentration
+ real64 dPoreVolume_dLogPrimaryConc[numSpecies]{};
+
+ };
+
+ /**
+ * @brief Performs the setup phase for the kernel.
+ * @param[in] ei the element index
+ * @param[in] stack the stack variables
+ */
+ GEOS_HOST_DEVICE
+ void setup( localIndex const ei,
+ StackVariables & stack ) const
+ {
+ // initialize the pore volume
+ stack.poreVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * m_porosity[ei][0];
+ stack.dPoreVolume_dPres = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dPres[ei][0];
+
+ Base::setup( ei, stack );
+
+ // // is - index of the primary species
+ // for( integer is = 0; is < numSpecies; ++is )
+ // {
+ // stack.dPoreVolume_dLogPrimaryConc[is] = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dLogPrimaryConc[ei][is]
+ // }
+ }
+
+ /**
+ * @brief Compute the local accumulation contributions to the residual and Jacobian
+ * @param[in] ei the element index
+ * @param[inout] stack the stack variables
+ * @param[in] kernelOp the function used to customize the kernel
+ */
+ GEOS_HOST_DEVICE
+ void computeAccumulation( localIndex const ei,
+ StackVariables & stack ) const
+ {
+ // Residual[is] += (primarySpeciesAggregateConcentration[is] * stack.poreVolume - primarySpeciesAggregateMole_n[is])
+ // - dt * m_volume * primarySpeciesKineticRate[is] // To Check: what's the unit of the kinetic rate
+
+ Base::computeAccumulation( ei, stack );
+
+ // Step 1: assemble the derivatives of the total mass equation w.r.t log of primary species concentration
+ // for( integer is = 0; is < numSpecies; ++is )
+ // {
+ // stack.localJacobian[0][is+numDof-numSpecies] = stack.poreVolume * m_dDensity_dLogPrimaryConc[ei][is] +
+ // stack.dPoreVolume_dLogPrimaryConc[is] * m_density[ei][0];
+ // }
+
+ arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dPrimarySpeciesAggregateConcentration_dLogPrimaryConc = m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc[ei];
+
+ for( integer is = 0; is < numSpecies; ++is )
+ {
+ // Step 2: assemble the accumulation term of the species mass balance equation
+ // Step 2.1: residual
+ // Primary species mole amount in pore volume
+ stack.localResidual[is+numEqn-numSpecies] -= m_primarySpeciesAggregateMole_n[ei][is];
+ stack.localResidual[is+numEqn-numSpecies] += m_primarySpeciesAggregateConcentration[ei][is] * stack.poreVolume;
+
+ // // Reaction term
+ // stack.localResidual[is+numEqn-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) *
+ // m_primarySpeciesAggregateKineticRate[is];
+
+ // Step 2.1: jacobian
+ // Drivative of primary species amount in pore volume wrt pressure
+ stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * m_primarySpeciesAggregateConcentration[ei][is]
+ /* + stack.poreVolume * m_dTotalPrimarySpeciesConcentration_dPres[ei][is] */;
+ // // Derivative of reaction term wrt pressure
+ // stack.localJacobian[is+numEqn-numSpecies][0] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) *
+ // m_dPrimarySpeciesTotalKineticRate_dPres[is];
+
+ // Derivative wrt log of primary species concentration
+ // arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dPrimarySpeciesTotalKineticRate_dLogPrimaryConc =
+ // m_dPrimarySpeciesTotalKineticRate_dLogPrimaryConc[ei];
+
+ for( integer js = 0; js < numSpecies; ++js )
+ {
+ stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] += /* stack.dPoreVolume_dLogPrimaryConc[js] *
+ m_primarySpeciesAggregateConcentration[ei][is]
+ + */stack.poreVolume * dPrimarySpeciesAggregateConcentration_dLogPrimaryConc[is][js]; // To
+ // check
+ // if
+ // the
+ // permutation
+ // is
+ // consistent
+
+ // stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) *
+ // dPrimarySpeciesTotalKineticRate_dLogPrimaryConc[is][js];
+ }
+ }
+ }
+
+ /**
+ * @brief Performs the complete phase for the kernel.
+ * @param[in] ei the element index
+ * @param[inout] stack the stack variables
+ */
+ GEOS_HOST_DEVICE
+ void complete( localIndex const ei,
+ StackVariables & stack ) const
+ {
+ // Step 1: assemble the total mass balance equation
+ // - the total mass balance equations (i = 0)
+ Base::complete( ei, stack );
+
+ // Step 2: assemble the primary species mole amount balance equation
+ // - the species mole amount balance equations (i = numEqn-numSpecies to i = numEqn-1)
+ integer const beginRowSpecies = numEqn-numSpecies;
+ for( integer i = 0; i < numSpecies; ++i )
+ {
+ m_localRhs[stack.localRow + beginRowSpecies + i] += stack.localResidual[beginRowSpecies+i];
+ m_localMatrix.template addToRow< serialAtomic >( stack.localRow + beginRowSpecies + i,
+ stack.dofIndices,
+ stack.localJacobian[beginRowSpecies + i],
+ numDof );
+ }
+ }
+
+protected:
+
+ /// Time step size
+ real64 const m_dt;
+
+ /// View on the element volumes
+ arrayView1d< real64 const > const m_volume;
+ arrayView1d< real64 const > const m_deltaVolume;
+
+ /// Views on the porosity
+ arrayView2d< real64 const > const m_porosity;
+ arrayView2d< real64 const > const m_dPoro_dPres;
+
+ // // View on the derivatives of fluid density wrt log of primary species concentration
+ // arrayView2d< real64 const, compflow::USD_COMP > m_dDensity_dLogPrimaryConc;
+
+ // // View on the derivatives of porosity wrt log of primary species concentration
+ // arrayView2d< real64 const, compflow::USD_COMP > m_dPoro_dLogPrimaryConc;
+
+ // View on the total concentration of ions that contain the primary species
+ arrayView2d< real64 const, compflow::USD_COMP > m_primarySpeciesAggregateConcentration;
+
+ // // View on the derivatives of aggregate concentration for the primary species wrt pressure
+ // arrayView2d< real64 const, compflow::USD_COMP > m_dPrimarySpeciesAggregateConcentration_dPres;
+
+ // View on the derivatives of total ion concentration for the primary species wrt log of primary species concentration
+ arrayView3d< real64 const, compflow::USD_COMP_DC > m_dPrimarySpeciesAggregateConcentration_dLogPrimaryConc;
+
+ // View on the aggregate kinetic rate of primary species from all reactions
+ arrayView2d< real64 const, compflow::USD_COMP > m_primarySpeciesAggregateKineticRate;
+
+ // // View on the derivatives of aggregate kinetic rate of primary species wrt pressure
+ // arrayView2d< real64 const, compflow::USD_COMP > m_dPrimarySpeciesAggregateKineticRate_dPres;
+
+ // // View on the derivatives of aggregate kinetic rate of primary species wrt log of primary species concentration
+ // arrayView3d< real64 const, compflow::USD_COMP_DC > m_dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc;
+
+ // View on primary species mole amount from previous time step
+ arrayView2d< real64 const, compflow::USD_COMP > m_primarySpeciesAggregateMole_n;
+};
+
+/**
+ * @class AccumulationKernelFactory
+ */
+class AccumulationKernelFactory
+{
+public:
+
+ /**
+ * @brief Create a new kernel and launch
+ * @tparam POLICY the policy used in the RAJA kernel
+ * @param[in] numSpecies the number of primary species
+ * @param[in] dt time step
+ * @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, typename SUBREGION_TYPE >
+ static void
+ createAndLaunch( integer const numSpecies,
+ real64 const dt,
+ globalIndex const rankOffset,
+ string const dofKey,
+ SUBREGION_TYPE const & subRegion,
+ constitutive::ReactiveSingleFluid const & fluid,
+ constitutive::CoupledSolidBase const & solid,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+ {
+ internal::kernelLaunchSelectorCompSwitch( numSpecies, [&] ( auto NS )
+ {
+ integer constexpr NUM_SPECIES = NS();
+ integer constexpr NUM_DOF = 1+NS();
+ AccumulationKernel< SUBREGION_TYPE, NUM_DOF, NUM_SPECIES > kernel( rankOffset, dofKey, subRegion, fluid, solid, dt, localMatrix, localRhs );
+ AccumulationKernel< SUBREGION_TYPE, NUM_DOF, NUM_SPECIES >::template launch< POLICY >( subRegion.size(), kernel );
+ } );
+ }
+};
+
+} // namespace singlePhaseReactiveBaseKernels
+
+} // namespace geos
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_ACCUMULATIONKERNELS_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluidUpdateKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluidUpdateKernel.hpp
new file mode 100644
index 00000000000..e2561c89c8e
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluidUpdateKernel.hpp
@@ -0,0 +1,58 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file FluidUpdateKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_FLUIDUPDATEKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_FLUIDUPDATEKERNEL_HPP
+
+#include "common/DataTypes.hpp"
+#include "common/GEOS_RAJA_Interface.hpp"
+
+namespace geos
+{
+
+namespace singlePhaseReactiveBaseKernels
+{
+
+/******************************** FluidUpdateKernel ********************************/
+
+struct FluidUpdateKernel
+{
+ template< typename FLUID_WRAPPER >
+ static void launch( FLUID_WRAPPER const & fluidWrapper,
+ arrayView1d< real64 const > const & pres,
+ arrayView1d< real64 const > const & temp,
+ arrayView2d< real64 const, compflow::USD_COMP > const logPrimaryConc )
+ {
+ forAll< parallelDevicePolicy<> >( fluidWrapper.numElems(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ for( localIndex q = 0; q < fluidWrapper.numGauss(); ++q )
+ {
+ fluidWrapper.update( k, q, pres[k], temp[k], logPrimaryConc[k] );
+
+ fluidWrapper.updateChemistryLogConc( k, q, pres[k], temp[k], logPrimaryConc[k] );
+ }
+ } );
+ }
+};
+
+} // namespace singlePhaseReactiveBaseKernels
+
+} // namespace geos
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_FLUIDUPDATEKERNEL_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp
new file mode 100644
index 00000000000..9902addd22f
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp
@@ -0,0 +1,564 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file FluxComputeKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_FLUXCOMPUTEKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_FLUXCOMPUTEKERNEL_HPP
+
+#include "constitutive/diffusion/DiffusionFields.hpp"
+#include "constitutive/diffusion/DiffusionBase.hpp"
+#include "constitutive/solid/porosity/PorosityBase.hpp"
+#include "constitutive/solid/porosity/PorosityFields.hpp"
+#include "constitutive/fluid/singlefluid/reactive/ReactiveSingleFluid.hpp"
+#include "constitutive/fluid/multifluid/reactive/ReactiveMultiFluidFields.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp"
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp"
+
+
+namespace geos
+{
+
+namespace singlePhaseReactiveFVMKernels
+{
+
+/**
+ * @class FluxComputeKernel
+ * @tparam NUM_SPECIES number of fluid primary species
+ * @tparam NUM_EQN number of equations
+ * @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_SPECIES, integer NUM_EQN, integer NUM_DOF, typename STENCILWRAPPER >
+class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, STENCILWRAPPER >
+{
+public:
+
+ /// Compile time value for the number of primary species
+ static constexpr integer numSpecies = NUM_SPECIES;
+
+ /// Number of flux support points (hard-coded for TFPA)
+ static constexpr integer numFluxSupportPoints = 2;
+
+ /**
+ * @brief The type for element-based data. Consists entirely of ArrayView's.
+ *
+ * Can be converted from ElementRegionManager::ElementViewConstAccessor
+ * by calling .toView() or .toViewConst() on an accessor instance
+ */
+ template< typename VIEWTYPE >
+ using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
+
+ using AbstractBase = singlePhaseFVMKernels::FluxComputeKernelBase;
+ using DofNumberAccessor = AbstractBase::DofNumberAccessor;
+ using SinglePhaseFlowAccessors = AbstractBase::SinglePhaseFlowAccessors;
+ using SinglePhaseFluidAccessors = AbstractBase::SinglePhaseFluidAccessors;
+ using PermeabilityAccessors = AbstractBase::PermeabilityAccessors;
+
+ using AbstractBase::m_dt;
+ using AbstractBase::m_rankOffset;
+ using AbstractBase::m_dofNumber;
+ using AbstractBase::m_gravCoef;
+ using AbstractBase::m_mob;
+ using AbstractBase::m_dens;
+ using AbstractBase::m_dDens;
+
+ using Base = singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, STENCILWRAPPER >;
+ 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 ReactiveSinglePhaseFlowAccessors =
+ StencilAccessors< fields::flow::logPrimarySpeciesConcentration,
+ fields::flow::dMobility_dLogPrimaryConc >;
+
+ using ReactiveSinglePhaseFluidAccessors =
+ StencilMaterialAccessors< constitutive::ReactiveSingleFluid,
+ fields::reactivefluid::primarySpeciesAggregateConcentration,
+ fields::reactivefluid::dPrimarySpeciesAggregateConcentration_dLogPrimaryConc >;
+
+ using DiffusionAccessors =
+ StencilMaterialAccessors< constitutive::DiffusionBase,
+ fields::diffusion::diffusivity,
+ fields::diffusion::dDiffusivity_dTemperature >;
+
+ using PorosityAccessors =
+ StencilMaterialAccessors< constitutive::PorosityBase,
+ fields::porosity::referencePorosity >;
+
+ /**
+ * @brief Constructor for the kernel interface
+ * @param[in] rankOffset the offset of my MPI rank
+ * @param[in] stencilWrapper reference to the stencil wrapper
+ * @param[in] dofNumberAccessor
+ * @param[in] singlePhaseFlowAccessors
+ * @param[in] reactiveSinglePhaseFlowAccessors
+ * @param[in] singlePhaseFluidAccessors
+ * @param[in] reactiveSinglePhaseFluidAccessors
+ * @param[in] permeabilityAccessors
+ * @param[in] diffusionAccessors
+ * @param[in] porosityAccessors
+ * @param[in] hasDiffusion the flag to turn on diffusion calculation
+ * @param[in] dt time step size
+ * @param[inout] localMatrix the local CRS matrix
+ * @param[inout] localRhs the local right-hand side vector
+ */
+ FluxComputeKernel( globalIndex const rankOffset,
+ STENCILWRAPPER const & stencilWrapper,
+ DofNumberAccessor const & dofNumberAccessor,
+ SinglePhaseFlowAccessors const & singlePhaseFlowAccessors,
+ ReactiveSinglePhaseFlowAccessors const & reactiveSinglePhaseFlowAccessors,
+ SinglePhaseFluidAccessors const & singlePhaseFluidAccessors,
+ ReactiveSinglePhaseFluidAccessors const & reactiveSinglePhaseFluidAccessors,
+ PermeabilityAccessors const & permeabilityAccessors,
+ DiffusionAccessors const & diffusionAccessors,
+ PorosityAccessors const & porosityAccessors,
+ integer const & hasDiffusion,
+ real64 const & dt,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+ : Base( rankOffset,
+ stencilWrapper,
+ dofNumberAccessor,
+ singlePhaseFlowAccessors,
+ singlePhaseFluidAccessors,
+ permeabilityAccessors,
+ dt,
+ localMatrix,
+ localRhs ),
+ m_logPrimarySpeciesConc( reactiveSinglePhaseFlowAccessors.get( fields::flow::logPrimarySpeciesConcentration {} ) ),
+ m_dMob_dLogPrimaryConc( reactiveSinglePhaseFlowAccessors.get( fields::flow::dMobility_dLogPrimaryConc {} ) ),
+ m_primarySpeciesAggregateConc( reactiveSinglePhaseFluidAccessors.get( fields::reactivefluid::primarySpeciesAggregateConcentration {} ) ),
+ m_dPrimarySpeciesAggregateConc_dLogPrimaryConc( reactiveSinglePhaseFluidAccessors.get( fields::reactivefluid::dPrimarySpeciesAggregateConcentration_dLogPrimaryConc {} ) ),
+ m_diffusivity( diffusionAccessors.get( fields::diffusion::diffusivity {} ) ),
+ m_dDiffusivity_dTemp( diffusionAccessors.get( fields::diffusion::dDiffusivity_dTemperature {} ) ),
+ m_referencePorosity( porosityAccessors.get( fields::porosity::referencePorosity {} ) ),
+ m_hasDiffusion( hasDiffusion )
+ {}
+
+ /**
+ * @struct StackVariables
+ * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack
+ */
+ struct StackVariables : public Base::StackVariables
+ {
+public:
+
+ /**
+ * @brief Constructor for the stack variables
+ * @param[in] size size of the stencil for this connection
+ * @param[in] numElems number of elements for this connection
+ */
+ GEOS_HOST_DEVICE
+ StackVariables( localIndex const size, localIndex numElems )
+ : Base::StackVariables( size, numElems )
+ {}
+
+ 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;
+
+ /// Diffusion transmissibility
+ real64 diffusionTransmissibility[maxNumConns][numFluxSupportPoints]{};
+ /// Derivatives of diffusion transmissibility with respect to temperature
+ real64 dDiffusionTrans_dT[maxNumConns][numFluxSupportPoints]{};
+
+ };
+
+ /**
+ * @brief Compute the local flux contributions to the residual and Jacobian
+ * @tparam FUNC the type of the function that can be used to customize the computation of the flux
+ * @param[in] iconn the connection index
+ * @param[inout] stack the stack variables
+ * @param[in] NoOpFunc the function used to customize the computation of the flux
+ */
+ template< typename FUNC = NoOpFunc >
+ GEOS_HOST_DEVICE
+ void computeFlux( localIndex const iconn,
+ StackVariables & stack,
+ FUNC && kernelOp = NoOpFunc{} ) const
+ {
+ using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
+ // ***********************************************
+ // First, we call the base computeFlux to compute:
+ // 1) massFlux and its derivatives,
+ // 2) speciesFlux and its derivatives
+ Base::computeFlux( iconn, stack, [&] ( localIndex const (&k)[2],
+ localIndex const (&seri)[2],
+ localIndex const (&sesri)[2],
+ localIndex const (&sei)[2],
+ localIndex const connectionIndex,
+ real64 const alpha,
+ real64 const mobility,
+ real64 const & potGrad,
+ real64 const & fluxVal,
+ real64 const (&dFlux_dP)[2] )
+ {
+ // Step 1: compute the derivatives of the fluid density, potential difference,
+ // and the massFlux wrt log of primary species concentration (to complete)
+ real64 dFlux_dLogConc[numFluxSupportPoints][numSpecies]{};
+
+ GEOS_UNUSED_VAR( dFlux_dLogConc ); // Todo: to add the massFlux derivatives wrt speciesConc
+
+ // Step 2: compute the speciesFlux
+ real64 speciesFlux[numSpecies]{};
+ real64 dSpeciesFlux_dP[numFluxSupportPoints][numSpecies]{};
+ real64 dSpeciesFlux_dLogConc[numFluxSupportPoints][numSpecies][numSpecies]{};
+ // real64 dSpeciesFlux_dTrans[numSpecies]{};
+
+ // choose upstream cell
+ localIndex const k_up = (potGrad >= 0) ? 0 : 1;
+
+ localIndex const er_up = seri[k_up];
+ localIndex const esr_up = sesri[k_up];
+ localIndex const ei_up = sei[k_up];
+
+ real64 const fluidDens_up = m_dens[er_up][esr_up][ei_up][0];
+ real64 const dDens_dPres = m_dDens[er_up][esr_up][ei_up][0][DerivOffset::dP];
+
+ // compute species fluxes and derivatives using upstream cell concentration
+ for( integer is = 0; is < numSpecies; ++is )
+ {
+ real64 const aggregateConc_i = m_primarySpeciesAggregateConc[er_up][esr_up][ei_up][is];
+ speciesFlux[is] = aggregateConc_i / fluidDens_up * fluxVal;
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ dSpeciesFlux_dP[ke][is] += aggregateConc_i / fluidDens_up * dFlux_dP[ke];
+ }
+
+ dSpeciesFlux_dP[k_up][is] += -aggregateConc_i * fluxVal * dDens_dPres / (fluidDens_up * fluidDens_up);
+
+ for( integer js = 0; js < numSpecies; ++js )
+ {
+ real64 const dAggregateConc_i_dLogConc_j = m_dPrimarySpeciesAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][is][js];
+ dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConc_i_dLogConc_j / fluidDens_up * fluxVal;
+ }
+ }
+
+ /// populate local flux vector and derivatives
+ for( integer is = 0; is < numSpecies; ++is )
+ {
+ integer const eqIndex0 = k[0] * numEqn + numEqn - numSpecies + is;
+ integer const eqIndex1 = k[1] * numEqn + numEqn - numSpecies + is;
+
+ stack.localFlux[eqIndex0] += m_dt * speciesFlux[is];
+ stack.localFlux[eqIndex1] -= m_dt * speciesFlux[is];
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const localDofIndexPres = k[ke] * numDof;
+ stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dSpeciesFlux_dP[ke][is];
+ stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dSpeciesFlux_dP[ke][is];
+
+ for( integer js = 0; js < numSpecies; ++js )
+ {
+ localIndex const localDofIndexSpecies = localDofIndexPres + js + numDof - numSpecies;
+ stack.localFluxJacobian[eqIndex0][localDofIndexSpecies] += m_dt * dSpeciesFlux_dLogConc[ke][is][js];
+ stack.localFluxJacobian[eqIndex1][localDofIndexSpecies] -= m_dt * dSpeciesFlux_dLogConc[ke][is][js];
+ }
+ }
+ }
+
+ // Customize the kernel with this lambda
+ kernelOp( k, seri, sesri, sei, connectionIndex, alpha, mobility, potGrad, fluxVal, dFlux_dP, fluidDens_up );
+ } );
+ }
+
+ /**
+ * @brief Compute the local diffusion contributions to the residual and Jacobian
+ * @tparam FUNC the type of the function that can be used to customize the computation of the flux
+ * @param[in] iconn the connection index
+ * @param[inout] stack the stack variables
+ * @param[in] NoOpFunc the function used to customize the computation of the flux
+ */
+ template< typename FUNC = NoOpFunc >
+ GEOS_HOST_DEVICE
+ void computeDiffusion( localIndex const iconn,
+ StackVariables & stack,
+ FUNC && kernelOp = NoOpFunc{} ) const
+ {
+ if( m_hasDiffusion )
+ {
+ // *****************************************************
+ // Computation of the diffusion term in the species flux
+
+ // Step 1: compute the diffusion transmissibilities at this face
+ m_stencilWrapper.computeWeights( iconn,
+ m_diffusivity,
+ m_dDiffusivity_dTemp,
+ stack.diffusionTransmissibility,
+ stack.dDiffusionTrans_dT );
+
+ localIndex k[numFluxSupportPoints]{};
+ localIndex connectionIndex = 0;
+
+ for( k[0] = 0; k[0] < stack.numFluxElems; ++k[0] )
+ {
+ for( k[1] = k[0] + 1; k[1] < stack.numFluxElems; ++k[1] )
+ {
+ localIndex const seri[numFluxSupportPoints] = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )};
+ localIndex const sesri[numFluxSupportPoints] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )};
+ localIndex const sei[numFluxSupportPoints] = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )};
+
+ // clear working arrays
+ real64 diffusionFlux[numSpecies]{};
+ real64 speciesGrad[numSpecies]{};
+ // real64 dDiffusionFlux_dP[numFluxSupportPoints][numSpecies]{}; // Turn on if diffusionFlux is pressure-dependent
+ real64 dDiffusionFlux_dLogConc[numFluxSupportPoints][numSpecies][numSpecies]{};
+
+ real64 const diffusionTrans[numFluxSupportPoints] = { stack.diffusionTransmissibility[connectionIndex][0],
+ stack.diffusionTransmissibility[connectionIndex][1] };
+
+ //***** calculation of flux *****
+ // loop over primary species
+ for( integer is = 0; is < numSpecies; ++is )
+ {
+ // real64 dSpeciesGrad_i_dP[numFluxSupportPoints]{}; // Turn on if speciesGrad is pressure-dependent
+ real64 dSpeciesGrad_i_dLogConc[numFluxSupportPoints][numSpecies]{};
+
+ // Step 2: compute species gradient
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const er = seri[ke];
+ localIndex const esr = sesri[ke];
+ localIndex const ei = sei[ke];
+
+ real64 const aggregateConc_i = m_primarySpeciesAggregateConc[er][esr][ei][is];
+
+ speciesGrad[is] += diffusionTrans[ke] * aggregateConc_i;
+
+ for( integer js = 0; js < numSpecies; ++js )
+ {
+ real64 const dAggregateConc_i_dLogConc_j = m_dPrimarySpeciesAggregateConc_dLogPrimaryConc[er][esr][ei][is][js];
+
+ dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConc_i_dLogConc_j;
+ }
+ }
+
+ // choose upstream cell for species upwinding
+ localIndex const k_up = (speciesGrad[is] >= 0) ? 0 : 1;
+
+ localIndex const er_up = seri[k_up];
+ localIndex const esr_up = sesri[k_up];
+ localIndex const ei_up = sei[k_up];
+
+ // computation of the upwinded species flux
+ diffusionFlux[is] += m_referencePorosity[er_up][esr_up][ei_up] * speciesGrad[is];
+
+ // add contributions of the derivatives of component fractions wrt pressure/component fractions
+ for( integer ke = 0; ke < numFluxSupportPoints; ke++ )
+ {
+ for( integer js = 0; js < numSpecies; ++js )
+ {
+ dDiffusionFlux_dLogConc[ke][is][js] += m_referencePorosity[er_up][esr_up][ei_up] * dSpeciesGrad_i_dLogConc[ke][js];
+ }
+ }
+
+ // Add the local diffusion flux contribution to the residual and Jacobian
+ integer const eqIndex0 = k[0] * numEqn + numEqn - numSpecies + is;
+ integer const eqIndex1 = k[1] * numEqn + numEqn - numSpecies + is;
+
+ stack.localFlux[eqIndex0] += m_dt * diffusionFlux[is];
+ stack.localFlux[eqIndex1] -= m_dt * diffusionFlux[is];
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const localDofIndexPres = k[ke] * numDof;
+ // stack.localFluxJacobian[eqIndex0][localDofIndexPres] += m_dt * dDiffusionFlux_dP[ke][is];
+ // stack.localFluxJacobian[eqIndex1][localDofIndexPres] -= m_dt * dDiffusionFlux_dP[ke][is];
+
+ for( integer js = 0; js < numSpecies; ++js )
+ {
+ localIndex const localDofIndexSpecies = localDofIndexPres + js + numDof - numSpecies;
+ stack.localFluxJacobian[eqIndex0][localDofIndexSpecies] += m_dt * dDiffusionFlux_dLogConc[ke][is][js];
+ stack.localFluxJacobian[eqIndex1][localDofIndexSpecies] -= m_dt * dDiffusionFlux_dLogConc[ke][is][js];
+ }
+ }
+
+ // Customize the kernel with this lambda
+ kernelOp( is, k, seri, sesri, sei, connectionIndex, k_up );
+ } // loop over primary species
+ connectionIndex++;
+ }
+ }
+ }
+ }
+
+ /**
+ * @brief Performs the complete phase for the kernel.
+ * @param[in] iconn the connection index
+ * @param[inout] stack the stack variables
+ */
+ template< typename FUNC = NoOpFunc >
+ GEOS_HOST_DEVICE
+ void complete( localIndex const iconn,
+ StackVariables & stack,
+ FUNC && kernelOp = NoOpFunc{} ) const
+ {
+ // Call Base::complete to assemble the total mass balance equation
+ // In the lambda, add contribution to residual and jacobian into the species amount balance equation
+ Base::complete( iconn, stack, [&] ( integer const i,
+ localIndex const localRow )
+ {
+ // The no. of fluxes is equal to the no. of equations in m_localRhs and m_localMatrix
+ for( integer is = 0; is < numSpecies; ++is )
+ {
+ RAJA::atomicAdd( parallelDeviceAtomic{}, &AbstractBase::m_localRhs[localRow + numEqn - numSpecies + is],
+ stack.localFlux[i * numEqn + numEqn - numSpecies + is] );
+ AbstractBase::m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >
+ ( localRow + numEqn - numSpecies + is,
+ stack.dofColIndices.data(),
+ stack.localFluxJacobian[i * numEqn + numEqn - numSpecies + is].dataIfContiguous(),
+ stack.stencilSize * numDof );
+ }
+
+ // call the lambda to assemble additional terms, such as thermal terms
+ kernelOp( i, localRow );
+ } );
+ }
+
+ /**
+ * @brief Performs the kernel launch
+ * @tparam POLICY the policy used in the RAJA kernels
+ * @tparam KERNEL_TYPE the kernel type
+ * @param[in] numConnections the number of connections
+ * @param[inout] kernelComponent the kernel component providing access to setup/compute/complete functions and stack variables
+ */
+ template< typename POLICY, typename KERNEL_TYPE >
+ static void
+ launch( localIndex const numConnections,
+ KERNEL_TYPE const & kernelComponent )
+ {
+ GEOS_MARK_FUNCTION;
+
+ forAll< POLICY >( numConnections, [=] GEOS_HOST_DEVICE ( localIndex const iconn )
+ {
+ typename KERNEL_TYPE::StackVariables stack( kernelComponent.stencilSize( iconn ),
+ kernelComponent.numPointsInFlux( iconn ) );
+
+ kernelComponent.setup( iconn, stack );
+ kernelComponent.computeFlux( iconn, stack );
+ kernelComponent.computeDiffusion( iconn, stack );
+ kernelComponent.complete( iconn, stack );
+ } );
+ }
+
+protected:
+
+ /// Views on log of primary species concentration
+ ElementViewConst< arrayView2d< real64 const, compflow::USD_COMP > > const m_logPrimarySpeciesConc;
+
+ /// Views on derivatives of fluid mobilities
+ ElementViewConst< arrayView2d< real64 const, compflow::USD_FLUID_DC > > const m_dMob_dLogPrimaryConc;
+
+ /// Views on primary species aggregate concentration
+ ElementViewConst< arrayView2d< real64 const, compflow::USD_COMP > > const m_primarySpeciesAggregateConc;
+
+ /// Views on the derivative of primary species aggregate concentration wrt log of primary concentration
+ ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const m_dPrimarySpeciesAggregateConc_dLogPrimaryConc;
+
+ /// Views on diffusivity
+ ElementViewConst< arrayView3d< real64 const > > const m_diffusivity;
+ ElementViewConst< arrayView3d< real64 const > > const m_dDiffusivity_dTemp;
+
+ /// View on the reference porosity
+ ElementViewConst< arrayView1d< real64 const > > const m_referencePorosity;
+
+ /// Flag of adding the diffusion term
+ integer const m_hasDiffusion;
+};
+
+/**
+ * @class FluxComputeKernelFactory
+ */
+class FluxComputeKernelFactory
+{
+public:
+
+ /**
+ * @brief Create a new kernel and launch
+ * @tparam POLICY the policy used in the RAJA kernel
+ * @tparam STENCILWRAPPER the type of the stencil wrapper
+ * @param[in] numSpecies the number of primary species
+ * @param[in] hasDiffusion the flag of adding diffusion term
+ * @param[in] rankOffset the offset of my MPI rank
+ * @param[in] dofKey string to get the element degrees of freedom numbers
+ * @param[in] solverName name of the solver (to name accessors)
+ * @param[in] elemManager reference to the element region manager
+ * @param[in] stencilWrapper reference to the stencil wrapper
+ * @param[in] dt time step size
+ * @param[inout] localMatrix the local CRS matrix
+ * @param[inout] localRhs the local right-hand side vector
+ */
+ template< typename POLICY, typename STENCILWRAPPER >
+ static void
+ createAndLaunch( integer const numSpecies,
+ integer const hasDiffusion,
+ globalIndex const rankOffset,
+ string const & dofKey,
+ string const & solverName,
+ ElementRegionManager const & elemManager,
+ STENCILWRAPPER const & stencilWrapper,
+ real64 const & dt,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+ {
+ singlePhaseReactiveBaseKernels::internal::kernelLaunchSelectorCompSwitch( numSpecies, [&]( auto NS )
+ {
+ integer constexpr NUM_SPECIES = NS();
+ integer constexpr NUM_DOF = 1+NS();
+ integer constexpr NUM_EQN = 1+NS();
+
+ ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > dofNumberAccessor =
+ elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey );
+ dofNumberAccessor.setName( solverName + "/accessors/" + dofKey );
+
+ using KernelType = FluxComputeKernel< NUM_SPECIES, NUM_EQN, NUM_DOF, STENCILWRAPPER >;
+ typename KernelType::SinglePhaseFlowAccessors flowAccessors( elemManager, solverName );
+ typename KernelType::ReactiveSinglePhaseFlowAccessors reactiveFlowAccessors( elemManager, solverName );
+ typename KernelType::SinglePhaseFluidAccessors fluidAccessors( elemManager, solverName );
+ typename KernelType::ReactiveSinglePhaseFluidAccessors reactiveFluidAccessors( elemManager, solverName );
+ typename KernelType::PermeabilityAccessors permAccessors( elemManager, solverName );
+ typename KernelType::DiffusionAccessors diffusionAccessors( elemManager, solverName );
+ typename KernelType::PorosityAccessors porosityAccessors( elemManager, solverName );
+
+ KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor,
+ flowAccessors, reactiveFlowAccessors, fluidAccessors, reactiveFluidAccessors,
+ permAccessors, diffusionAccessors, porosityAccessors, hasDiffusion,
+ dt, localMatrix, localRhs );
+ KernelType::template launch< POLICY >( stencilWrapper.size(), kernel );
+ } );
+ }
+};
+
+} // namespace singlePhaseReactiveFVMKernels
+
+} // namespace geos
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_FLUXCOMPUTEKERNEL_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp
new file mode 100644
index 00000000000..a1eff26eb32
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp
@@ -0,0 +1,73 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file KernelLaunchSelector.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_KERNELLAUNCHSELECTOR_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_KERNELLAUNCHSELECTOR_HPP
+
+#include "physicsSolvers/KernelLaunchSelectors.hpp"
+#include "codingUtilities/Utilities.hpp"
+#include "common/DataLayouts.hpp"
+#include "common/DataTypes.hpp"
+#include "common/GEOS_RAJA_Interface.hpp"
+
+namespace geos
+{
+
+namespace singlePhaseReactiveBaseKernels
+{
+
+/******************************** Kernel launch machinery ********************************/
+
+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; }
+ case 6:
+ { lambda( std::integral_constant< T, 5 >() ); return; }
+ case 7:
+ { lambda( std::integral_constant< T, 5 >() ); return; }
+ default:
+ { GEOS_ERROR( "Unsupported number of primary species: " << value ); }
+ }
+}
+
+} // namespace internal
+
+} // namespace singlePhaseReactiveBaseKernels
+
+} // namespace geos
+
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_KERNELLAUNCHSELECTOR_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp
new file mode 100644
index 00000000000..c9fc35acce1
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ResidualNormKernel.hpp
@@ -0,0 +1,331 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ResidualNormKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_RESIDUALNORMKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_RESIDUALNORMKERNEL_HPP
+
+#include "physicsSolvers/PhysicsSolverBaseKernels.hpp"
+#include "physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp"
+
+namespace geos
+{
+
+namespace singlePhaseReactiveBaseKernels
+{
+
+/******************************** ResidualNormKernel ********************************/
+
+/**
+ * @class IsothermalResidualNormKernel
+ */
+class IsothermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBase< 2 >
+{
+public:
+
+ using Base = physicsSolverBaseKernels::ResidualNormKernelBase< 2 >;
+ using Base::m_minNormalizer;
+ using Base::m_rankOffset;
+ using Base::m_localResidual;
+ using Base::m_dofNumber;
+
+ IsothermalResidualNormKernel( globalIndex const rankOffset,
+ arrayView1d< real64 const > const & localResidual,
+ arrayView1d< globalIndex const > const & dofNumber,
+ arrayView1d< localIndex const > const & ghostRank,
+ integer const numPrimarySpecies,
+ ElementSubRegionBase const & subRegion,
+ real64 const minNormalizer )
+ : Base( rankOffset,
+ localResidual,
+ dofNumber,
+ ghostRank,
+ minNormalizer ),
+ m_numPrimarySpecies( numPrimarySpecies ),
+ m_mass_n( subRegion.template getField< fields::flow::mass_n >() ),
+ m_primarySpeciesAggregateMole_n( subRegion.getField< fields::flow::primarySpeciesAggregateMole_n >() )
+ {}
+
+ GEOS_HOST_DEVICE
+ virtual void computeLinf( localIndex const ei,
+ LinfStackVariables & stack ) const override
+ {
+ real64 const totalMassNormalizer = LvArray::math::max( m_minNormalizer, m_mass_n[ei] );
+
+ // step 1: total mass residuals
+ real64 const valMass = LvArray::math::abs( m_localResidual[stack.localRow] ) / totalMassNormalizer;
+ if( valMass > stack.localValue[0] )
+ {
+ stack.localValue[0] = valMass;
+ }
+
+ // step 2: species amount residuals
+ for( integer idof = 0; idof < m_numPrimarySpecies; ++idof )
+ {
+ real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] );
+ real64 const valAmount = LvArray::math::abs( m_localResidual[stack.localRow + idof + 1] ) / speciesAmountNormalizer;
+ if( valAmount > stack.localValue[1] )
+ {
+ stack.localValue[1] = valAmount;
+ }
+ }
+ }
+
+ GEOS_HOST_DEVICE
+ virtual void computeL2( localIndex const ei,
+ L2StackVariables & stack ) const override
+ {
+ real64 const totalMassNormalizer = LvArray::math::max( m_minNormalizer, m_mass_n[ei] );
+
+ // step 1: total mass residuals
+ stack.localValue[0] += m_localResidual[stack.localRow] * m_localResidual[stack.localRow];
+ stack.localNormalizer[0] += totalMassNormalizer;
+
+ // step 2: species amount residuals
+ for( integer idof = 0; idof < m_numPrimarySpecies; ++idof )
+ {
+ real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] );
+
+ stack.localValue[1] += m_localResidual[stack.localRow + idof + 1] * m_localResidual[stack.localRow + idof + 1];
+ stack.localNormalizer[1] += speciesAmountNormalizer;
+ }
+ }
+
+
+protected:
+
+ /// Number of primary species
+ integer const m_numPrimarySpecies;
+
+ /// View on mass at the previous converged time step
+ arrayView1d< real64 const > const m_mass_n;
+
+ // View on primary species aggregate amount (moles) from previous time step
+ arrayView2d< real64 const, compflow::USD_COMP > m_primarySpeciesAggregateMole_n;
+
+};
+
+/**
+ * @class ThermalResidualNormKernel
+ */
+class ThermalResidualNormKernel : public physicsSolverBaseKernels::ResidualNormKernelBase< 3 >
+{
+public:
+
+ using Base = physicsSolverBaseKernels::ResidualNormKernelBase< 3 >;
+ using Base::m_minNormalizer;
+ using Base::m_rankOffset;
+ using Base::m_localResidual;
+ using Base::m_dofNumber;
+
+ ThermalResidualNormKernel( globalIndex const rankOffset,
+ arrayView1d< real64 const > const & localResidual,
+ arrayView1d< globalIndex const > const & dofNumber,
+ arrayView1d< localIndex const > const & ghostRank,
+ integer const numPrimarySpecies,
+ ElementSubRegionBase const & subRegion,
+ real64 const minNormalizer )
+ : Base( rankOffset,
+ localResidual,
+ dofNumber,
+ ghostRank,
+ minNormalizer ),
+ m_numPrimarySpecies( numPrimarySpecies ),
+ m_mass_n( subRegion.template getField< fields::flow::mass_n >() ),
+ m_primarySpeciesAggregateMole_n( subRegion.getField< fields::flow::primarySpeciesAggregateMole_n >() ),
+ m_energy_n( subRegion.template getField< fields::flow::energy_n >() )
+ {}
+
+ GEOS_HOST_DEVICE
+ void computeMassEnergyNormalizers( localIndex const ei,
+ real64 & totalMassNormalizer,
+ real64 & energyNormalizer ) const
+ {
+ totalMassNormalizer = LvArray::math::max( m_minNormalizer, m_mass_n[ei] );
+ energyNormalizer = LvArray::math::max( m_minNormalizer, LvArray::math::abs( m_energy_n[ei] ) ); // energy can be negative
+ }
+
+ GEOS_HOST_DEVICE
+ virtual void computeLinf( localIndex const ei,
+ LinfStackVariables & stack ) const override
+ {
+ real64 totalMassNormalizer = 0.0, energyNormalizer = 0.0;
+ computeMassEnergyNormalizers( ei, totalMassNormalizer, energyNormalizer );
+
+ // step 1: mass residual
+
+ real64 const valMass = LvArray::math::abs( m_localResidual[stack.localRow] ) / totalMassNormalizer;
+ if( valMass > stack.localValue[0] )
+ {
+ stack.localValue[0] = valMass;
+ }
+
+ // step 2: energy residual
+ real64 const valEnergy = LvArray::math::abs( m_localResidual[stack.localRow + 1] ) / energyNormalizer;
+ if( valEnergy > stack.localValue[1] )
+ {
+ stack.localValue[1] = valEnergy;
+ }
+
+ // step 3: species amount residuals
+ for( integer idof = 0; idof < m_numPrimarySpecies; ++idof )
+ {
+ real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] );
+ real64 const valAmount = LvArray::math::abs( m_localResidual[stack.localRow + idof + 2] ) / speciesAmountNormalizer;
+ if( valAmount > stack.localValue[2] )
+ {
+ stack.localValue[2] = valAmount;
+ }
+ }
+ }
+
+ GEOS_HOST_DEVICE
+ virtual void computeL2( localIndex const ei,
+ L2StackVariables & stack ) const override
+ {
+ real64 totalMassNormalizer = 0.0, energyNormalizer = 0.0;
+ computeMassEnergyNormalizers( ei, totalMassNormalizer, energyNormalizer );
+
+ // step 1: mass residual
+
+ stack.localValue[0] += m_localResidual[stack.localRow] * m_localResidual[stack.localRow];
+ stack.localNormalizer[0] += totalMassNormalizer;
+
+ // step 2: energy residual
+
+ stack.localValue[1] += m_localResidual[stack.localRow + 1] * m_localResidual[stack.localRow + 1];
+ stack.localNormalizer[1] += energyNormalizer;
+
+ // step 3: species amount residuals
+ for( integer idof = 0; idof < m_numPrimarySpecies; ++idof )
+ {
+ real64 const speciesAmountNormalizer = LvArray::math::max( m_minNormalizer, m_primarySpeciesAggregateMole_n[ei][idof] );
+
+ stack.localValue[2] += m_localResidual[stack.localRow + idof + 2] * m_localResidual[stack.localRow + idof + 2];
+ stack.localNormalizer[2] += speciesAmountNormalizer;
+ }
+ }
+
+
+protected:
+
+ /// Number of primary species
+ integer const m_numPrimarySpecies;
+
+ /// View on mass at the previous converged time step
+ arrayView1d< real64 const > const m_mass_n;
+
+ // View on primary species aggregate amount (moles) from previous time step
+ arrayView2d< real64 const, compflow::USD_COMP > m_primarySpeciesAggregateMole_n;
+
+ /// View on energy at the previous converged time step
+ arrayView1d< real64 const > const m_energy_n;
+
+};
+
+/**
+ * @class ResidualNormKernelFactory
+ */
+class ResidualNormKernelFactory
+{
+public:
+
+ /**
+ * @brief Create a new kernel and launch
+ * @tparam POLICY the policy used in the RAJA kernel
+ * @param[in] normType the type of norm used (Linf or L2)
+ * @param[in] numPrimarySpecies the number of primary species
+ * @param[in] rankOffset the offset of my MPI rank
+ * @param[in] dofKey the string key to retrieve the degress of freedom numbers
+ * @param[in] localResidual the residual vector on my MPI rank
+ * @param[in] subRegion the element subregion
+ * @param[out] residualNorm the residual norm on the subRegion
+ * @param[out] residualNormalizer the residual normalizer on the subRegion
+ */
+ template< typename POLICY >
+ static void
+ createAndLaunch( physicsSolverBaseKernels::NormType const normType,
+ integer const numPrimarySpecies,
+ globalIndex const rankOffset,
+ string const dofKey,
+ arrayView1d< real64 const > const & localResidual,
+ ElementSubRegionBase const & subRegion,
+ real64 const minNormalizer,
+ real64 (& residualNorm)[2],
+ real64 (& residualNormalizer)[2] )
+ {
+ arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
+ arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+
+ IsothermalResidualNormKernel kernel( rankOffset, localResidual, dofNumber, ghostRank, numPrimarySpecies, subRegion, minNormalizer );
+ if( normType == physicsSolverBaseKernels::NormType::Linf )
+ {
+ IsothermalResidualNormKernel::launchLinf< POLICY >( subRegion.size(), kernel, residualNorm );
+ }
+ else // L2 norm
+ {
+ IsothermalResidualNormKernel::launchL2< POLICY >( subRegion.size(), kernel, residualNorm, residualNormalizer );
+ }
+ }
+
+ /**
+ * @brief Create a new kernel and launch (thermal version)
+ * @tparam POLICY the policy used in the RAJA kernel
+ * @param[in] normType the type of norm used (Linf or L2)
+ * @param[in] numPrimarySpecies the number of primary species
+ * @param[in] rankOffset the offset of my MPI rank
+ * @param[in] dofKey the string key to retrieve the degress of freedom numbers
+ * @param[in] localResidual the residual vector on my MPI rank
+ * @param[in] subRegion the element subregion
+ * @param[out] residualNorm the residual norm on the subRegion
+ * @param[out] residualNormalizer the residual normalizer on the subRegion
+ */
+ template< typename POLICY >
+ static void
+ createAndLaunch( physicsSolverBaseKernels::NormType const normType,
+ integer const numPrimarySpecies,
+ globalIndex const rankOffset,
+ string const & dofKey,
+ arrayView1d< real64 const > const & localResidual,
+ ElementSubRegionBase const & subRegion,
+ real64 const minNormalizer,
+ real64 (& residualNorm)[3],
+ real64 (& residualNormalizer)[3] )
+ {
+ arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey );
+ arrayView1d< integer const > const ghostRank = subRegion.ghostRank();
+
+ ThermalResidualNormKernel kernel( rankOffset, localResidual, dofNumber, ghostRank, numPrimarySpecies, subRegion, minNormalizer );
+ if( normType == physicsSolverBaseKernels::NormType::Linf )
+ {
+ ThermalResidualNormKernel::launchLinf< POLICY >( subRegion.size(), kernel, residualNorm );
+ }
+ else // L2 norm
+ {
+ ThermalResidualNormKernel::launchL2< POLICY >( subRegion.size(), kernel, residualNorm, residualNormalizer );
+ }
+ }
+
+};
+
+} // namespace singlePhaseReactiveBaseKernels
+
+} // namespace geos
+
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_REACTIVE_RESIDUALNORMKERNEL_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp
new file mode 100644
index 00000000000..8911ddde9b9
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp
@@ -0,0 +1,259 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ThermalAccumulationKernels.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_THERMALACCUMULATIONKERNELS_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_THERMALACCUMULATIONKERNELS_HPP
+
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp"
+
+namespace geos
+{
+
+namespace thermalSinglePhaseReactiveBaseKernels
+{
+
+/******************************** AccumulationKernel ********************************/
+
+/**
+ * @class AccumulationKernel
+ * @brief Define the interface for the assembly kernel in charge of accumulation
+ */
+template< typename SUBREGION_TYPE, integer NUM_DOF, integer NUM_SPECIES >
+class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKernel< SUBREGION_TYPE, NUM_DOF, NUM_SPECIES >
+{
+
+public:
+
+ using Base = singlePhaseReactiveBaseKernels::AccumulationKernel< SUBREGION_TYPE, NUM_DOF, NUM_SPECIES >;
+ using Base::numDof;
+ using Base::numEqn;
+ using Base::numSpecies;
+ using Base::m_rankOffset;
+ using Base::m_dofNumber;
+ using Base::m_elemGhostRank;
+ using Base::m_localMatrix;
+ using Base::m_localRhs;
+ using Base::m_dMass;
+ using Base::m_volume;
+ using Base::m_deltaVolume;
+ using Base::m_primarySpeciesAggregateConcentration;
+
+ /// Note: Derivative lineup only supports dP & dT, not component terms
+ using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
+ /**
+ * @brief Constructor
+ * @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
+ */
+ AccumulationKernel( globalIndex const rankOffset,
+ string const dofKey,
+ SUBREGION_TYPE const & subRegion,
+ constitutive::ReactiveSingleFluid const & fluid,
+ constitutive::CoupledSolidBase const & solid,
+ real64 const & dt,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+ : Base( rankOffset, dofKey, subRegion, fluid, solid, dt, localMatrix, localRhs ),
+ m_energy( subRegion.template getField< fields::flow::energy >() ),
+ m_energy_n( subRegion.template getField< fields::flow::energy_n >() ),
+ m_dEnergy( subRegion.template getField< fields::flow::dEnergy >() ),
+ m_dPoro_dTemp( solid.getDporosity_dTemperature() )
+ // m_dPrimarySpeciesAggregateConcentration_dTemp( fluid.dPrimarySpeciesAggregateConcentration_dTemp() ),
+ // m_dPrimarySpeciesTotalKineticRate_dTemp( fluid.dPrimarySpeciesTotalKineticRate_dTemp() ),
+ {}
+
+ /**
+ * @struct StackVariables
+ * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack
+ */
+ struct StackVariables : public Base::StackVariables
+ {
+public:
+
+ GEOS_HOST_DEVICE
+ StackVariables()
+ : Base::StackVariables()
+ {}
+
+ using Base::StackVariables::localRow;
+ using Base::StackVariables::dofIndices;
+ using Base::StackVariables::localResidual;
+ using Base::StackVariables::localJacobian;
+ using Base::StackVariables::poreVolume;
+ using Base::StackVariables::dPoreVolume_dLogPrimaryConc;
+
+ // Pore volume information
+
+ /// Derivative of pore volume with respect to temperature
+ real64 dPoreVolume_dTemp = 0.0;
+ };
+
+ /**
+ * @brief Performs the setup phase for the kernel.
+ * @param[in] ei the element index
+ * @param[in] stack the stack variables
+ */
+ GEOS_HOST_DEVICE
+ void setup( localIndex const ei,
+ StackVariables & stack ) const
+ {
+ Base::setup( ei, stack );
+
+ stack.dPoreVolume_dTemp = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dTemp[ei][0];
+ }
+
+ /**
+ * @brief Compute the local accumulation contributions to the residual and Jacobian
+ * @param[in] ei the element index
+ * @param[inout] stack the stack variables
+ * @param[in] kernelOp the function used to customize the kernel
+ */
+ GEOS_HOST_DEVICE
+ void computeAccumulation( localIndex const ei,
+ StackVariables & stack ) const
+ {
+ Base::computeAccumulation( ei, stack );
+
+ // Step 1: assemble the derivatives of the mass balance equation w.r.t temperature
+ stack.localJacobian[0][numDof-numSpecies-1] = m_dMass[ei][DerivOffset::dT];
+
+ // Step 2: assemble the accumulation term of the energy equation
+ // Step 2.1: assemble the residual and derivatives wrt pressure and temperature
+ stack.localResidual[numEqn-numSpecies-1] = m_energy[ei] - m_energy_n[ei];
+ stack.localJacobian[numEqn-numSpecies-1][0] += m_dEnergy[ei][DerivOffset::dP];
+ stack.localJacobian[numEqn-numSpecies-1][numDof-numSpecies-1] += m_dEnergy[ei][DerivOffset::dT];
+
+ // Step 2.2: assemble the derivatives of the energy equation w.r.t log primary species concentration
+ // for( integer is = 0; is < numSpecies; ++is )
+ // {
+ // stack.localJacobian[numEqn-numSpecies-1][is+numDof-numSpecies] += stack.dPoreVolume_dLogPrimaryConc[is] * m_density[ei][0] *
+ // m_fluidInternalEnergy[ei][0]
+ // - stack.dPoreVolume_dLogPrimaryConc[is] *
+ // m_rockInternalEnergy[ei][0]
+ // + stack.poreVolume * m_dDensity_dLogPrimaryConc[ei][is] *
+ // m_fluidInternalEnergy[ei][0]
+ // + stack.poreVolume * m_density[ei][0] *
+ // m_dFluidInternalEnergy_dLogPrimaryConc[ei][is];
+ // }
+
+ // Step 3: assemble the derivatives of the species amount balance equation w.r.t temperature
+ for( integer is = 0; is < numSpecies; ++is )
+ {
+ // Drivative of primary species amount in pore volume wrt temperature
+ stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += stack.dPoreVolume_dTemp * m_primarySpeciesAggregateConcentration[ei][is]
+ /* + stack.poreVolume *
+ m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] */;
+ // // Derivative of reaction term wrt temperature
+ // stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) *
+ // m_dPrimarySpeciesTotalKineticRate_dTemp[is];
+ }
+ }
+
+ /**
+ * @brief Performs the complete phase for the kernel.
+ * @param[in] ei the element index
+ * @param[inout] stack the stack variables
+ */
+ GEOS_HOST_DEVICE
+ void complete( localIndex const ei,
+ StackVariables & stack ) const
+ {
+ // Step 1: assemble the total mass balance equation (i = 0)
+ // and species amount balance equation (i = numEqn-numSpecies to i = numEqn-1)
+ Base::complete( ei, stack );
+
+ // Step 2: assemble the energy equation (i = numEqn-numSpecies-1)
+ m_localRhs[stack.localRow + numEqn-numSpecies-1] += stack.localResidual[numEqn-numSpecies-1];
+ m_localMatrix.template addToRow< serialAtomic >( stack.localRow + numEqn-numSpecies-1,
+ stack.dofIndices,
+ stack.localJacobian[numEqn-numSpecies-1],
+ numDof );
+ }
+
+protected:
+
+ /// View on energy
+ arrayView1d< real64 const > const m_energy;
+ arrayView1d< real64 const > const m_energy_n;
+ arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_dEnergy;
+
+ /// Views on the porosity derivative
+ arrayView2d< real64 const > const m_dPoro_dTemp;
+
+ // // View on the derivatives of aggregate concentration for the primary species wrt temperature
+ // arrayView2d< real64 const, compflow::USD_COMP > m_dPrimarySpeciesAggregateConcentration_dTemp;
+
+ // // View on the derivatives of total kinetic rate of primary species wrt temperature
+ // arrayView2d< real64 const, compflow::USD_COMP > m_dPrimarySpeciesTotalKineticRate_dTemp;
+
+};
+
+/**
+ * @class AccumulationKernelFactory
+ */
+class AccumulationKernelFactory
+{
+public:
+
+ /**
+ * @brief Create a new kernel and launch
+ * @tparam POLICY the policy used in the RAJA kernel
+ * @param[in] numSpecies the number of primary species
+ * @param[in] dt time step
+ * @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, typename SUBREGION_TYPE >
+ static void
+ createAndLaunch( integer const numSpecies,
+ real64 const dt,
+ globalIndex const rankOffset,
+ string const dofKey,
+ SUBREGION_TYPE const & subRegion,
+ constitutive::ReactiveSingleFluid const & fluid,
+ constitutive::CoupledSolidBase const & solid,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+ {
+ singlePhaseReactiveBaseKernels::
+ internal::kernelLaunchSelectorCompSwitch( numSpecies, [&] ( auto NS )
+ {
+ integer constexpr NUM_SPECIES = NS();
+ integer constexpr NUM_DOF = 2+NS();
+ AccumulationKernel< SUBREGION_TYPE, NUM_DOF, NUM_SPECIES > kernel( rankOffset, dofKey, subRegion, fluid, solid, dt, localMatrix, localRhs );
+ AccumulationKernel< SUBREGION_TYPE, NUM_DOF, NUM_SPECIES >::template launch< POLICY >( subRegion.size(), kernel );
+ } );
+ }
+};
+
+} // namespace thermalSinglePhaseReactiveBaseKernels
+
+} // namespace geos
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_THERMALACCUMULATIONKERNELS_HPP
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp
new file mode 100644
index 00000000000..a61571baea6
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp
@@ -0,0 +1,628 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 TotalEnergies
+ * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2023-2024 Chevron
+ * Copyright (c) 2019- GEOS/GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file ThermalFluxComputeKernel.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_THERMALFLUXCOMPUTEKERNEL_HPP
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_THERMALFLUXCOMPUTEKERNEL_HPP
+
+#include "physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp"
+
+#include "constitutive/thermalConductivity/SinglePhaseThermalConductivityBase.hpp"
+#include "constitutive/thermalConductivity/ThermalConductivityFields.hpp"
+
+namespace geos
+{
+
+namespace thermalSinglePhaseReactiveFVMKernels
+{
+/******************************** FluxComputeKernel ********************************/
+
+/**
+ * @class FluxComputeKernel
+ * @tparam NUM_SPECIES number of fluid primary species
+ * @tparam NUM_EQN number of equations
+ * @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_SPECIES, integer NUM_EQN, integer NUM_DOF, typename STENCILWRAPPER >
+class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKernel< NUM_SPECIES, NUM_EQN, 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 = singlePhaseFVMKernels::FluxComputeKernelBase;
+ using DofNumberAccessor = AbstractBase::DofNumberAccessor;
+ using SinglePhaseFlowAccessors = AbstractBase::SinglePhaseFlowAccessors;
+ using SinglePhaseFluidAccessors = AbstractBase::SinglePhaseFluidAccessors;
+ using PermeabilityAccessors = AbstractBase::PermeabilityAccessors;
+
+ using AbstractBase::m_dt;
+ using AbstractBase::m_rankOffset;
+ using AbstractBase::m_dofNumber;
+ using AbstractBase::m_gravCoef;
+ using AbstractBase::m_mob;
+ using AbstractBase::m_dMob;
+ using AbstractBase::m_dens;
+ using AbstractBase::m_dDens;
+
+ using Base = singlePhaseReactiveFVMKernels::FluxComputeKernel< NUM_SPECIES, NUM_EQN, NUM_DOF, STENCILWRAPPER >;
+ using ReactiveSinglePhaseFlowAccessors = typename Base::ReactiveSinglePhaseFlowAccessors;
+ using ReactiveSinglePhaseFluidAccessors = typename Base::ReactiveSinglePhaseFluidAccessors;
+ using DiffusionAccessors = typename Base::DiffusionAccessors;
+ using PorosityAccessors = typename Base::PorosityAccessors;
+ using Base::numSpecies;
+ using Base::numFluxSupportPoints;
+ using Base::numDof;
+ using Base::numEqn;
+ using Base::maxNumElems;
+ using Base::maxNumConns;
+ using Base::maxStencilSize;
+ using Base::m_stencilWrapper;
+ using Base::m_seri;
+ using Base::m_sesri;
+ using Base::m_sei;
+ using Base::m_primarySpeciesAggregateConc;
+ using Base::m_referencePorosity;
+
+ using ThermalSinglePhaseFlowAccessors =
+ StencilAccessors< fields::flow::temperature >;
+
+ using ThermalReactiveSinglePhaseFluidAccessors =
+ StencilMaterialAccessors< constitutive::ReactiveSingleFluid,
+ fields::singlefluid::enthalpy,
+ fields::singlefluid::dEnthalpy >;
+
+ using ThermalConductivityAccessors =
+ StencilMaterialAccessors< constitutive::SinglePhaseThermalConductivityBase,
+ fields::thermalconductivity::effectiveConductivity,
+ fields::thermalconductivity::dEffectiveConductivity_dT >;
+
+
+ /**
+ * @brief Constructor for the kernel interface
+ * @param[in] rankOffset the offset of my MPI rank
+ * @param[in] stencilWrapper reference to the stencil wrapper
+ * @param[in] dofNumberAccessor accessor for the dofs numbers
+ * @param[in] singlePhaseFlowAccessors accessor for wrappers registered by the solver
+ * @param[in] reactiveSinglePhaseFlowAccessors accessor for *reactive* wrappers registered by the solver
+ * @param[in] thermalSinglePhaseFlowAccessors accessor for *thermal* wrappers registered by the solver
+ * @param[in] singlePhaseFluidAccessors accessor for wrappers registered by the single fluid model
+ * @param[in] reactiveSinglePhaseFluidAccessors accessor for *reactive* wrappers registered by the single fluid model
+ * @param[in] thermalReactiveSinglePhaseFluidAccessors accessor for *thermal reactive* wrappers registered by the single fluid model
+ * @param[in] permeabilityAccessors accessor for wrappers registered by the permeability model
+ * @param[in] diffusionAccessors accessor for wrappers registered by the diffusion model
+ * @param[in] porosityAccessors accessor for wrappers registered by the porosity model
+ * @param[in] thermalConductivityAccessors accessor for wrappers registered by the thermal conductivity model
+ * @param[in] hasDiffusion the flag to turn on diffusion calculation
+ * @param[in] dt time step size
+ * @param[inout] localMatrix the local CRS matrix
+ * @param[inout] localRhs the local right-hand side vector
+ */
+ FluxComputeKernel( globalIndex const rankOffset,
+ STENCILWRAPPER const & stencilWrapper,
+ DofNumberAccessor const & dofNumberAccessor,
+ SinglePhaseFlowAccessors const & singlePhaseFlowAccessors,
+ ReactiveSinglePhaseFlowAccessors const & reactiveSinglePhaseFlowAccessors,
+ ThermalSinglePhaseFlowAccessors const & thermalSinglePhaseFlowAccessors,
+ SinglePhaseFluidAccessors const & singlePhaseFluidAccessors,
+ ReactiveSinglePhaseFluidAccessors const & reactiveSinglePhaseFluidAccessors,
+ ThermalReactiveSinglePhaseFluidAccessors const & thermalReactiveSinglePhaseFluidAccessors,
+ PermeabilityAccessors const & permeabilityAccessors,
+ DiffusionAccessors const & diffusionAccessors,
+ PorosityAccessors const & porosityAccessors,
+ ThermalConductivityAccessors const & thermalConductivityAccessors,
+ integer const & hasDiffusion,
+ real64 const & dt,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+ : Base( rankOffset,
+ stencilWrapper,
+ dofNumberAccessor,
+ singlePhaseFlowAccessors,
+ reactiveSinglePhaseFlowAccessors,
+ singlePhaseFluidAccessors,
+ reactiveSinglePhaseFluidAccessors,
+ permeabilityAccessors,
+ diffusionAccessors,
+ porosityAccessors,
+ hasDiffusion,
+ dt,
+ localMatrix,
+ localRhs ),
+ m_temp( thermalSinglePhaseFlowAccessors.get( fields::flow::temperature {} ) ),
+ m_enthalpy( thermalReactiveSinglePhaseFluidAccessors.get( fields::singlefluid::enthalpy {} ) ),
+ m_dEnthalpy( thermalReactiveSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) ),
+ // m_dPrimarySpeciesAggregateConcentration_dTemp( fluid.dPrimarySpeciesAggregateConcentration_dTemp() ),
+ m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) ),
+ m_dThermalCond_dT( thermalConductivityAccessors.get( fields::thermalconductivity::dEffectiveConductivity_dT {} ) )
+ {}
+
+ struct StackVariables : public Base::StackVariables
+ {
+public:
+
+ GEOS_HOST_DEVICE
+ StackVariables( localIndex const size, localIndex numElems )
+ : Base::StackVariables( size, numElems ),
+ energyFlux( 0.0 ),
+ dEnergyFlux_dP( size ),
+ dEnergyFlux_dT( size )
+ {}
+
+ 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;
+ using Base::StackVariables::diffusionTransmissibility;
+ using Base::StackVariables::dDiffusionTrans_dT;
+
+
+ // Thermal transmissibility
+ real64 thermalTransmissibility[maxNumConns][2]{};
+
+ /// Derivatives of thermal transmissibility with respect to temperature
+ real64 dThermalTrans_dT[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;
+
+ };
+
+ /**
+ * @brief Compute the local flux contributions to the residual and Jacobian
+ * @param[in] iconn the connection index
+ * @param[inout] stack the stack variables
+ */
+ GEOS_HOST_DEVICE
+ void computeFlux( localIndex const iconn,
+ StackVariables & stack ) const
+ {
+ using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >;
+ // ***********************************************
+ // First, we call the base computeFlux to compute:
+ // 1) massFlux and speciesFlux and their derivatives (including derivatives wrt temperature),
+ // 2) enthalpy part of energyFlux and its derivatives (including derivatives wrt temperature)
+ //
+ // Computing dFlux_dT and the enthalpy flux requires quantities already computed in the base computeFlux,
+ // such as potGrad, fluxVal, and the indices of the upwind cell
+ // We use the lambda below (called **inside** the phase loop of the base computeFlux) to access these variables
+ Base::computeFlux( iconn, stack, [&] ( localIndex const (&k)[2],
+ localIndex const (&seri)[2],
+ localIndex const (&sesri)[2],
+ localIndex const (&sei)[2],
+ localIndex const connectionIndex,
+ real64 const alpha,
+ real64 const mobility,
+ real64 const & potGrad,
+ real64 const & fluxVal,
+ real64 const (&dFlux_dP)[2],
+ real64 const fluidDens_up )
+ {
+ // Step 1: compute the derivatives of the (upwinded) massFlux wrt temperature
+ // --------------------------------------------------------------------------
+ // Step 1.1: compute the derivatives of the mean density at the interface wrt temperature
+ real64 dDensMean_dT[numFluxSupportPoints]{0.0, 0.0};
+
+ real64 const trans[numFluxSupportPoints] = { stack.transmissibility[connectionIndex][0], stack.transmissibility[connectionIndex][1] };
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ dDensMean_dT[ke] = 0.5 * m_dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT];
+ }
+
+ // Step 1.2: compute the derivatives of the potential difference wrt temperature
+ real64 dGravHead_dT[numFluxSupportPoints]{0.0, 0.0};
+
+ // compute derivative of gravity potential difference wrt temperature
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const er = seri[ke];
+ localIndex const esr = sesri[ke];
+ localIndex const ei = sei[ke];
+
+ real64 const gravD = trans[ke] * m_gravCoef[er][esr][ei];
+
+ for( integer i = 0; i < numFluxSupportPoints; ++i )
+ {
+ dGravHead_dT[i] += dDensMean_dT[i] * gravD;
+ }
+ }
+
+ real64 dFlux_dT[numFluxSupportPoints]{0.0, 0.0};
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ dFlux_dT[ke] -= dGravHead_dT[ke];
+ }
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ dFlux_dT[ke] *= mobility;
+ }
+
+ // compute the derivatives of the mobility wrt temperature
+ // *** upwinding ***
+ real64 dMob_dT[numFluxSupportPoints]{};
+
+ if( alpha <= 0.0 || alpha >= 1.0 )
+ {
+ localIndex const k_up = 1 - localIndex( fmax( fmin( alpha, 1.0 ), 0.0 ) );
+ dMob_dT[k_up] = m_dMob[seri[k_up]][sesri[k_up]][sei[k_up]][DerivOffset::dT];
+ }
+ else
+ {
+ real64 const mobWeights[numFluxSupportPoints] = { alpha, 1.0 - alpha };
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ dMob_dT[ke] = mobWeights[ke] * m_dMob[seri[ke]][sesri[ke]][sei[ke]][DerivOffset::dT];
+ }
+ }
+
+ // add contribution from upstream cell mobility derivatives
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ dFlux_dT[ke] += dMob_dT[ke] * potGrad;
+ }
+
+ // Step 1.3: populate local jacobian
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const localDofIndexTemp = k[ke] * numDof + numDof - numSpecies - 1;
+ stack.localFluxJacobian[k[0]*numEqn][localDofIndexTemp] += m_dt * dFlux_dT[ke];
+ stack.localFluxJacobian[k[1]*numEqn][localDofIndexTemp] -= m_dt * dFlux_dT[ke];
+ }
+
+ // Step 2: compute the derivatives of the speciesFlux wrt temperature
+ // -------------------------------------------------------------------
+ real64 dSpeciesFlux_dT[numFluxSupportPoints][numSpecies]{};
+
+ {
+ // Step 2.1: compute the derivatives of the upstream density wrt temperature
+ // choose upstream cell
+ localIndex const k_up = (potGrad >= 0) ? 0 : 1;
+
+ localIndex const er_up = seri[k_up];
+ localIndex const esr_up = sesri[k_up];
+ localIndex const ei_up = sei[k_up];
+
+ real64 const dDens_dTemp = m_dDens[er_up][esr_up][ei_up][0][DerivOffset::dT];
+
+ // Step 2.2: compute speciesFlux derivative wrt temperature
+ for( integer is = 0; is < numSpecies; ++is )
+ {
+ real64 const aggregateConc_i = m_primarySpeciesAggregateConc[er_up][esr_up][ei_up][is];
+
+ // real64 const dAggregateConc_i_dTemp = m_dPrimarySpeciesAggregateConcentration_dTemp[er_up][esr_up][ei_up][is];
+ // dSpeciesFlux_dT[k_up][is] += dAggregateConc_i_dTemp * fluxVal / fluidDens_up;
+ dSpeciesFlux_dT[k_up][is] += -aggregateConc_i * fluxVal * dDens_dTemp / (fluidDens_up * fluidDens_up);
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ dSpeciesFlux_dT[ke][is] += aggregateConc_i / fluidDens_up * dFlux_dT[ke];
+ }
+ }
+ }
+
+ // Step 2.3: populate local jacobian
+ for( integer is = 0; is < numSpecies; ++is )
+ {
+ integer const eqIndex0 = k[0] * numEqn + numEqn - numSpecies + is;
+ integer const eqIndex1 = k[1] * numEqn + numEqn - numSpecies + is;
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const localDofIndexTemp = k[ke] * numDof + numDof - numSpecies - 1;
+
+ stack.localFluxJacobian[eqIndex0][localDofIndexTemp] += m_dt * dSpeciesFlux_dT[ke][is];
+ stack.localFluxJacobian[eqIndex1][localDofIndexTemp] -= m_dt * dSpeciesFlux_dT[ke][is];
+ }
+ }
+
+ // Step 3: compute the enthalpy flux
+ // ----------------------------------
+ real64 enthalpy = 0.0;
+ real64 dEnthalpy_dP[numFluxSupportPoints]{0.0, 0.0};
+ real64 dEnthalpy_dT[numFluxSupportPoints]{0.0, 0.0};
+ // Todo: to add the enthalpy derivatives wrt speciesConc if needed
+ // real64 dEnthalpy_dLogConc[numFluxSupportPoints][numSpecies]{};
+
+ if( alpha <= 0.0 || alpha >= 1.0 )
+ {
+ localIndex const k_up = 1 - localIndex( fmax( fmin( alpha, 1.0 ), 0.0 ) );
+
+ enthalpy = m_enthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0];
+ dEnthalpy_dP[k_up] = m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][DerivOffset::dP];
+ dEnthalpy_dT[k_up] = m_dEnthalpy[seri[k_up]][sesri[k_up]][sei[k_up]][0][DerivOffset::dT];
+ }
+ else
+ {
+ real64 const mobWeights[numFluxSupportPoints] = { alpha, 1.0 - alpha };
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ enthalpy += mobWeights[ke] * m_enthalpy[seri[ke]][sesri[ke]][sei[ke]][0];
+ dEnthalpy_dP[ke] = mobWeights[ke] * m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dP];
+ dEnthalpy_dT[ke] = mobWeights[ke] * m_dEnthalpy[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT];
+ }
+ }
+
+ stack.energyFlux += fluxVal * enthalpy;
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ stack.dEnergyFlux_dP[ke] += dFlux_dP[ke] * enthalpy;
+ stack.dEnergyFlux_dT[ke] += dFlux_dT[ke] * enthalpy;
+ }
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ stack.dEnergyFlux_dP[ke] += fluxVal * dEnthalpy_dP[ke];
+ stack.dEnergyFlux_dT[ke] += fluxVal * dEnthalpy_dT[ke];
+ }
+
+ } );
+
+ // *****************************************************
+ // Computation of the conduction term in the energy flux
+ // Note that the enthalpy term in the energy was computed above
+ // Note that this term is computed using an explicit treatment of conductivity for now
+
+ // Step 1: compute the thermal transmissibilities at this face
+ // We follow how the thermal compositional multi-phase solver does to update the thermal transmissibility
+ m_stencilWrapper.computeWeights( iconn,
+ m_thermalConductivity,
+ m_dThermalCond_dT,
+ stack.thermalTransmissibility,
+ stack.dThermalTrans_dT );
+
+ localIndex k[numFluxSupportPoints];
+ localIndex connectionIndex = 0;
+
+ for( k[0] = 0; k[0] < stack.numFluxElems; ++k[0] )
+ {
+ for( k[1] = k[0] + 1; k[1] < stack.numFluxElems; ++k[1] )
+ {
+ real64 const thermalTrans[numFluxSupportPoints] = { stack.thermalTransmissibility[connectionIndex][0], stack.thermalTransmissibility[connectionIndex][1] };
+ real64 const dThermalTrans_dT[numFluxSupportPoints] = { stack.dThermalTrans_dT[connectionIndex][0], stack.dThermalTrans_dT[connectionIndex][1] };
+
+ localIndex const seri[numFluxSupportPoints] = {m_seri( iconn, k[0] ), m_seri( iconn, k[1] )};
+ localIndex const sesri[numFluxSupportPoints] = {m_sesri( iconn, k[0] ), m_sesri( iconn, k[1] )};
+ localIndex const sei[numFluxSupportPoints] = {m_sei( iconn, k[0] ), m_sei( iconn, k[1] )};
+
+ // Step 2: compute temperature difference at the interface
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const er = seri[ke];
+ localIndex const esr = sesri[ke];
+ localIndex const ei = sei[ke];
+
+ stack.energyFlux += thermalTrans[ke] * m_temp[er][esr][ei];
+ stack.dEnergyFlux_dT[ke] += thermalTrans[ke] + dThermalTrans_dT[ke] * m_temp[er][esr][ei];
+ }
+
+ integer const eqIndex0 = k[0] * numEqn + numEqn - numSpecies - 1;
+ integer const eqIndex1 = k[1] * numEqn + numEqn - numSpecies - 1;
+
+ // add energyFlux and its derivatives to localFlux and localFluxJacobian
+ stack.localFlux[eqIndex0] += m_dt * stack.energyFlux;
+ stack.localFlux[eqIndex1] -= m_dt * stack.energyFlux;
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ integer const localDofIndexPres = k[ke] * numDof;
+ stack.localFluxJacobian[eqIndex0][localDofIndexPres] = m_dt * stack.dEnergyFlux_dP[ke];
+ stack.localFluxJacobian[eqIndex1][localDofIndexPres] = -m_dt * stack.dEnergyFlux_dP[ke];
+ integer const localDofIndexTemp = localDofIndexPres + numDof - numSpecies - 1;
+ stack.localFluxJacobian[eqIndex0][localDofIndexTemp] = m_dt * stack.dEnergyFlux_dT[ke];
+ stack.localFluxJacobian[eqIndex1][localDofIndexTemp] = -m_dt * stack.dEnergyFlux_dT[ke];
+ }
+
+ connectionIndex++;
+ }
+ }
+ }
+
+ /**
+ * @brief Compute the local flux contributions to the residual and Jacobian
+ * @param[in] iconn the connection index
+ * @param[inout] stack the stack variables
+ */
+ GEOS_HOST_DEVICE
+ void computeDiffusion( localIndex const iconn,
+ StackVariables & stack ) const
+ {
+ Base::computeDiffusion( iconn, stack, [&] ( integer const is,
+ localIndex const (&k)[2],
+ localIndex const (&seri)[2],
+ localIndex const (&sesri)[2],
+ localIndex const (&sei)[2],
+ localIndex const connectionIndex,
+ localIndex const k_up )
+ {
+ real64 dDiffusionFlux_dT[numFluxSupportPoints]{};
+ real64 dSpeciesGrad_dT[numFluxSupportPoints]{};
+
+ // Calculate diffusion derivative wrt temperature
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const er = seri[ke];
+ localIndex const esr = sesri[ke];
+ localIndex const ei = sei[ke];
+
+ // dSpeciesGrad_dT[ke] += stack.diffusionTransmissibility[connectionIndex][ke]
+ // * m_dPrimarySpeciesAggregateConcentration_dTemp[er][esr][ei][is];
+
+ dSpeciesGrad_dT[ke] += stack.dDiffusionTrans_dT[connectionIndex][ke] * m_primarySpeciesAggregateConc[er][esr][ei][is];
+ }
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ke++ )
+ {
+ localIndex const er_up = seri[k_up];
+ localIndex const esr_up = sesri[k_up];
+ localIndex const ei_up = sei[k_up];
+
+ dDiffusionFlux_dT[ke] += m_referencePorosity[er_up][esr_up][ei_up] * dSpeciesGrad_dT[ke];
+ }
+
+ // populate local Jacobian
+ integer const eqIndex0 = k[0] * numEqn + numEqn - numSpecies + is;
+ integer const eqIndex1 = k[1] * numEqn + numEqn - numSpecies + is;
+
+ for( integer ke = 0; ke < numFluxSupportPoints; ++ke )
+ {
+ localIndex const localDofIndexTemp = k[ke] * numDof + numDof - numSpecies - 1;
+ stack.localFluxJacobian[eqIndex0][localDofIndexTemp] += m_dt * dDiffusionFlux_dT[ke];
+ stack.localFluxJacobian[eqIndex1][localDofIndexTemp] -= m_dt * dDiffusionFlux_dT[ke];
+ }
+ } );
+ }
+
+ /**
+ * @brief Performs the complete phase for the kernel.
+ * @param[in] iconn the connection index
+ * @param[inout] stack the stack variables
+ */
+ GEOS_HOST_DEVICE
+ void complete( localIndex const iconn,
+ StackVariables & stack ) const
+ {
+ // Call Case::complete to assemble the mass balance equations
+ // In the lambda, add contribution to residual and jacobian into the energy balance equation
+ Base::complete( iconn, stack, [&] ( integer const i,
+ localIndex const localRow )
+ {
+ // The no. of fluxes is equal to the no. of equations in m_localRhs and m_localMatrix
+ // Different from the one in compositional multi-phase flow, which has a volume balance eqn.
+ RAJA::atomicAdd( parallelDeviceAtomic{}, &AbstractBase::m_localRhs[localRow + numEqn - numSpecies - 1], stack.localFlux[i * numEqn + numEqn - numSpecies - 1] );
+
+ AbstractBase::m_localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow + numEqn - numSpecies - 1,
+ stack.dofColIndices.data(),
+ stack.localFluxJacobian[i * numEqn + numEqn - numSpecies - 1].dataIfContiguous(),
+ stack.stencilSize * numDof );
+
+ } );
+ }
+
+protected:
+
+ /// Views on temperature
+ ElementViewConst< arrayView1d< real64 const > > const m_temp;
+
+ /// Views on enthalpies
+ ElementViewConst< arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > > const m_enthalpy;
+ ElementViewConst< arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DER > > const m_dEnthalpy;
+
+ // /// Views on the derivative of primary species aggregate concentration wrt temperature
+ // ElementViewConst< arrayView2d< real64 const, compflow::USD_COMP > > const m_dPrimarySpeciesAggregateConc_dTemp;
+
+ /// View on thermal conductivity
+ ElementViewConst< arrayView3d< real64 const > > m_thermalConductivity;
+
+ /// View on derivatives of thermal conductivity w.r.t. temperature
+ ElementViewConst< arrayView3d< real64 const > > m_dThermalCond_dT;
+
+};
+
+/**
+ * @class FluxComputeKernelFactory
+ */
+class FluxComputeKernelFactory
+{
+public:
+
+ /**
+ * @brief Create a new kernel and launch
+ * @tparam POLICY the policy used in the RAJA kernel
+ * @tparam STENCILWRAPPER the type of the stencil wrapper
+ * @param[in] numSpecies the number of primary species
+ * @param[in] hasDiffusion the flag of adding diffusion term
+ * @param[in] rankOffset the offset of my MPI rank
+ * @param[in] dofKey string to get the element degrees of freedom numbers
+ * @param[in] solverName name of the solver (to name accessors)
+ * @param[in] elemManager reference to the element region manager
+ * @param[in] stencilWrapper reference to the stencil wrapper
+ * @param[in] dt time step size
+ * @param[inout] localMatrix the local CRS matrix
+ * @param[inout] localRhs the local right-hand side vector
+ */
+ template< typename POLICY, typename STENCILWRAPPER >
+ static void
+ createAndLaunch( integer const numSpecies,
+ integer const hasDiffusion,
+ globalIndex const rankOffset,
+ string const & dofKey,
+ string const & solverName,
+ ElementRegionManager const & elemManager,
+ STENCILWRAPPER const & stencilWrapper,
+ real64 const & dt,
+ CRSMatrixView< real64, globalIndex const > const & localMatrix,
+ arrayView1d< real64 > const & localRhs )
+ {
+ singlePhaseReactiveBaseKernels::internal::kernelLaunchSelectorCompSwitch( numSpecies, [&]( auto NS )
+ {
+ integer constexpr NUM_SPECIES = NS();
+ integer constexpr NUM_DOF = 2+NS();
+ integer constexpr NUM_EQN = 2+NS();
+
+ ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > dofNumberAccessor =
+ elemManager.constructArrayViewAccessor< globalIndex, 1 >( dofKey );
+ dofNumberAccessor.setName( solverName + "/accessors/" + dofKey );
+
+ using KernelType = FluxComputeKernel< NUM_SPECIES, NUM_EQN, NUM_DOF, STENCILWRAPPER >;
+ typename KernelType::SinglePhaseFlowAccessors flowAccessors( elemManager, solverName );
+ typename KernelType::ReactiveSinglePhaseFlowAccessors reactiveFlowAccessors( elemManager, solverName );
+ typename KernelType::ThermalSinglePhaseFlowAccessors thermalFlowAccessors( elemManager, solverName );
+ typename KernelType::SinglePhaseFluidAccessors fluidAccessors( elemManager, solverName );
+ typename KernelType::ReactiveSinglePhaseFluidAccessors reactiveFluidAccessors( elemManager, solverName );
+ typename KernelType::ThermalReactiveSinglePhaseFluidAccessors thermalFluidAccessors( elemManager, solverName );
+ typename KernelType::PermeabilityAccessors permAccessors( elemManager, solverName );
+ typename KernelType::DiffusionAccessors diffusionAccessors( elemManager, solverName );
+ typename KernelType::PorosityAccessors porosityAccessors( elemManager, solverName );
+ typename KernelType::ThermalConductivityAccessors thermalConductivityAccessors( elemManager, solverName );
+
+ KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor,
+ flowAccessors, reactiveFlowAccessors, thermalFlowAccessors, fluidAccessors, reactiveFluidAccessors, thermalFluidAccessors,
+ permAccessors, diffusionAccessors, porosityAccessors, thermalConductivityAccessors,
+ hasDiffusion, dt, localMatrix, localRhs );
+ KernelType::template launch< POLICY >( stencilWrapper.size(), kernel );
+ } );
+ }
+};
+
+} // namespace thermalSinglePhaseReactiveFVMKernels
+
+} // namespace geos
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEREACTIVE_THERMALFLUXCOMPUTEKERNEL_HPP