diff --git a/src/coreComponents/constitutive/fluid/multifluid/Layouts.hpp b/src/coreComponents/constitutive/fluid/multifluid/Layouts.hpp index 626f988104e..998eb89803e 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/Layouts.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/Layouts.hpp @@ -31,17 +31,6 @@ namespace geos namespace constitutive { -namespace singlefluid -{ -struct DerivativeOffset -{ - /// index of derivative wrt pressure - static integer constexpr dP = 0; - /// index of derivative wrt temperature - static integer constexpr dT = 1; - -}; -} namespace multifluid { diff --git a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.cpp index 28550afbfc6..9ed421a7942 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.cpp @@ -89,8 +89,8 @@ void CompressibleSinglePhaseFluid::allocateConstitutiveData( dataRepository::Gro getField< fields::singlefluid::density >().setApplyDefaultValue( m_defaultDensity ); getField< fields::singlefluid::viscosity >().setApplyDefaultValue( m_defaultViscosity ); - m_density.setValues< serialPolicy >( m_referenceDensity ); - m_viscosity.setValues< serialPolicy >( m_referenceViscosity ); + m_density.value.setValues< serialPolicy >( m_referenceDensity ); + m_viscosity.value.setValues< serialPolicy >( m_referenceViscosity ); } void CompressibleSinglePhaseFluid::postInputInitialization() @@ -131,6 +131,11 @@ void CompressibleSinglePhaseFluid::postInputInitialization() real64 dRho_dP; real64 dVisc_dP; createKernelWrapper().compute( m_referencePressure, m_referenceDensity, dRho_dP, m_referenceViscosity, dVisc_dP ); + + for( integer i=0; i().setDefaultValue( dRho_dP ); getField< fields::singlefluid::dViscosity_dPressure >().setDefaultValue( dVisc_dP ); } @@ -140,9 +145,11 @@ CompressibleSinglePhaseFluid::createKernelWrapper() { return KernelWrapper( KernelWrapper::DensRelationType( m_referencePressure, m_referenceDensity, m_compressibility ), KernelWrapper::ViscRelationType( m_referencePressure, m_referenceViscosity, m_viscosibility ), - m_density, + m_density.value, + m_density.derivs, m_dDensity_dPressure, - m_viscosity, + m_viscosity.value, + m_viscosity.derivs, m_dViscosity_dPressure ); } diff --git a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp index d4c545b002e..d4b6eda60ba 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp @@ -39,19 +39,23 @@ template< ExponentApproximationType DENS_EAT, ExponentApproximationType VISC_EAT class CompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate { public: - + using SingleFluidProp = SingleFluidVar< real64, 2, singlefluid::LAYOUT_FLUID, singlefluid::LAYOUT_FLUID_DC >; using DensRelationType = ExponentialRelation< real64, DENS_EAT >; using ViscRelationType = ExponentialRelation< real64, VISC_EAT >; - + using DerivOffset = singlefluid::DerivativeOffset; CompressibleSinglePhaseUpdate( DensRelationType const & densRelation, ViscRelationType const & viscRelation, arrayView2d< real64 > const & density, + arrayView3d< real64 > const & dDensity, arrayView2d< real64 > const & dDens_dPres, arrayView2d< real64 > const & viscosity, + arrayView3d< real64 > const & dViscosity, arrayView2d< real64 > const & dVisc_dPres ) : SingleFluidBaseUpdate( density, + dDensity, dDens_dPres, viscosity, + dViscosity, dVisc_dPres ), m_densRelation( densRelation ), m_viscRelation( viscRelation ) @@ -114,6 +118,14 @@ class CompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate m_dDens_dPres[k][q], m_viscosity[k][q], m_dVisc_dPres[k][q] ); + compute( pressure, + m_density[k][q], + m_dDensity[k][q][DerivOffset::dP], + m_viscosity[k][q], + m_dVisc_dPres[k][q] ); + // tjb std::cout << m_dDens_dPres[k][q]<< " " << m_density_c.derivs[k][q][DerivOffset::dP] << std::endl; + // std::cout.flush(); + // assert(fabs(m_dDens_dPres[k][q]-m_density_c.derivs[k][q][DerivOffset::dP])( m_referenceDensity ); - m_viscosity.setValues< serialPolicy >( m_referenceViscosity ); + m_density.value.setValues< serialPolicy >( m_referenceDensity ); + m_viscosity.value.setValues< serialPolicy >( m_referenceViscosity ); } @@ -109,20 +109,22 @@ ProppantSlurryFluid::createKernelWrapper() m_nIndices, m_Ks, m_isNewtonianFluid, - m_density, + m_density.value, + m_density.derivs, m_dDensity_dPressure, m_dDensity_dProppantConc, m_dDensity_dCompConc, m_componentDensity, m_dCompDens_dPres, m_dCompDens_dCompConc, - m_fluidDensity, + m_fluidDensity.value, m_dFluidDens_dPres, m_dFluidDens_dCompConc, m_fluidViscosity, m_dFluidVisc_dPres, m_dFluidVisc_dCompConc, - m_viscosity, + m_viscosity.value, + m_viscosity.derivs, m_dViscosity_dPressure, m_dViscosity_dProppantConc, m_dViscosity_dCompConc ); diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ProppantSlurryFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ProppantSlurryFluid.hpp index 0a9c1610635..97de80517bc 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/ProppantSlurryFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/ProppantSlurryFluid.hpp @@ -34,6 +34,8 @@ namespace constitutive class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate { public: + using SingleFluidProp = SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >; + using DerivOffset = constitutive::singlefluid::DerivativeOffsetC<0>; /** * @brief @@ -80,6 +82,7 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate arrayView1d< real64 const > const & Ks, bool const isNewtonianFluid, arrayView2d< real64 > const & density, + arrayView3d< real64 > const & dDensity, arrayView2d< real64 > const & dDens_dPres, arrayView2d< real64 > const & dDens_dProppantConc, arrayView3d< real64 > const & dDens_dCompConc, @@ -93,6 +96,7 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate arrayView2d< real64 > const & dFluidVisc_dPres, arrayView3d< real64 > const & dFluidVisc_dCompConc, arrayView2d< real64 > const & viscosity, + arrayView3d< real64 > const & dViscosity, arrayView2d< real64 > const & dVisc_dPres, arrayView2d< real64 > const & dVisc_dProppantConc, arrayView3d< real64 > const & dVisc_dCompConc ) @@ -103,6 +107,7 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate Ks, isNewtonianFluid, density, + dDensity, dDens_dPres, dDens_dProppantConc, dDens_dCompConc, @@ -116,6 +121,7 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate dFluidVisc_dPres, dFluidVisc_dCompConc, viscosity, + dViscosity, dVisc_dPres, dVisc_dProppantConc, dVisc_dCompConc ), @@ -164,10 +170,12 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate m_dFluidVisc_dCompConc[k][q], isProppantBoundary, m_density[k][q], + m_dDensity[k][q][DerivOffset::dP], // tjb add deriv:dp m_dDensity_dPressure[k][q], m_dDensity_dProppantConc[k][q], m_dDensity_dCompConc[k][q], m_viscosity[k][q], + m_dViscosity[k][q][DerivOffset::dP],// tjb add deriv:dp m_dViscosity_dPressure[k][q], m_dViscosity_dProppantConc[k][q], m_dViscosity_dCompConc[k][q] ); @@ -260,10 +268,12 @@ class ProppantSlurryFluidUpdate final : public SlurryFluidBaseUpdate arraySlice1d< real64 const > const & dFluidViscosity_dComponentConcentration, integer const & isProppantBoundary, real64 & density, + real64 & dDensity_dp, // tjb real64 & dDensity_dPressure, real64 & dDensity_dProppantConcentration, arraySlice1d< real64 > const & dDensity_dComponentConcentration, real64 & viscosity, + real64 & dViscosity_dp, // tjb real64 & dViscosity_dPressure, real64 & dViscosity_dProppantConcentration, arraySlice1d< real64 > const & dViscosity_dComponentConcentration ) const; @@ -477,10 +487,12 @@ ProppantSlurryFluidUpdate:: arraySlice1d< real64 const > const & GEOS_UNUSED_PARAM( dFluidViscosity_dComponentConcentration ), integer const & isProppantBoundary, real64 & density, + real64 & dDensity_dp, //tjb real64 & dDensity_dPressure, real64 & dDensity_dProppantConcentration, arraySlice1d< real64 > const & dDensity_dComponentConcentration, real64 & viscosity, + real64 & dViscosity_dp, //tjb real64 & dViscosity_dPressure, real64 & dViscosity_dProppantConcentration, arraySlice1d< real64 > const & dViscosity_dComponentConcentration ) const @@ -495,6 +507,7 @@ ProppantSlurryFluidUpdate:: density = (1.0 - effectiveConcentration) * fluidDensity + effectiveConcentration * m_referenceProppantDensity; dDensity_dPressure = (1.0 - effectiveConcentration) * dFluidDensity_dPressure; + dDensity_dp = (1.0 - effectiveConcentration) * dFluidDensity_dPressure; dDensity_dProppantConcentration = 0.0; for( localIndex c = 0; c < NC; ++c ) @@ -505,7 +518,7 @@ ProppantSlurryFluidUpdate:: real64 const coef = pow( 1.0 + 1.25 * effectiveConcentration / (1.0 - effectiveConcentration / m_maxProppantConcentration), 2.0 ); viscosity = fluidViscosity * coef; dViscosity_dPressure = dFluidViscosity_dPressure * coef; - + dViscosity_dp = dFluidViscosity_dPressure * coef; dViscosity_dProppantConcentration = 0.0; for( localIndex c = 0; c < NC; ++c ) { diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.cpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.cpp index e9870a759a4..2e70650d141 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.cpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.cpp @@ -30,25 +30,36 @@ namespace constitutive { SingleFluidBase::SingleFluidBase( string const & name, Group * const parent ) - : ConstitutiveBase( name, parent ) + : ConstitutiveBase( name, parent ), + m_numDOF( 1 ) { - registerField( fields::singlefluid::density{}, &m_density ); + //registerField( fields::singlefluid::density{}, &m_density.value ); registerField( fields::singlefluid::dDensity_dPressure{}, &m_dDensity_dPressure ); registerField( fields::singlefluid::dDensity_dTemperature{}, &m_dDensity_dTemperature ); registerField( fields::singlefluid::density_n{}, &m_density_n ); - registerField( fields::singlefluid::viscosity{}, &m_viscosity ); + registerField( fields::singlefluid::density{}, &m_density.value ); + registerField( fields::singlefluid::dDensity{}, &m_density.derivs ); + + registerField( fields::singlefluid::viscosity{}, &m_viscosity.value ); + registerField( fields::singlefluid::dViscosity{}, &m_viscosity.derivs ); + registerField( fields::singlefluid::dViscosity_dPressure{}, &m_dViscosity_dPressure ); registerField( fields::singlefluid::dViscosity_dTemperature{}, &m_dViscosity_dTemperature ); - registerField( fields::singlefluid::internalEnergy{}, &m_internalEnergy ); + registerField( fields::singlefluid::internalEnergy{}, &m_internalEnergy.value ); + registerField( fields::singlefluid::dInternalEnergy{}, &m_internalEnergy.derivs ); + registerField( fields::singlefluid::internalEnergy_n{}, &m_internalEnergy_n ); registerField( fields::singlefluid::dInternalEnergy_dPressure{}, &m_dInternalEnergy_dPressure ); registerField( fields::singlefluid::dInternalEnergy_dTemperature{}, &m_dInternalEnergy_dTemperature ); - registerField( fields::singlefluid::enthalpy{}, &m_enthalpy ); + registerField( fields::singlefluid::enthalpy{}, &m_enthalpy.value ); + registerField( fields::singlefluid::dEnthalpy{}, &m_enthalpy.derivs ); +#if 1 registerField( fields::singlefluid::dEnthalpy_dPressure{}, &m_dEnthalpy_dPressure ); registerField( fields::singlefluid::dEnthalpy_dTemperature{}, &m_dEnthalpy_dTemperature ); +#endif } @@ -68,8 +79,8 @@ void SingleFluidBase::initializeState() const void SingleFluidBase::saveConvergedState() const { - m_density_n.setValues< parallelDevicePolicy<> >( m_density.toViewConst() ); - m_internalEnergy_n.setValues< parallelDevicePolicy<> >( m_internalEnergy.toViewConst() ); + m_density_n.setValues< parallelDevicePolicy<> >( m_density.value.toViewConst() ); + m_internalEnergy_n.setValues< parallelDevicePolicy<> >( m_internalEnergy.value.toViewConst() ); } //START_SPHINX_INCLUDE_00 @@ -80,21 +91,33 @@ void SingleFluidBase::allocateConstitutiveData( Group & parent, resize( parent.size() ); - m_density.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dDensity_dPressure.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dDensity_dTemperature.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_density_n.resize( parent.size(), numConstitutivePointsPerParentIndex ); - m_viscosity.resize( parent.size(), numConstitutivePointsPerParentIndex ); + // tjb new density + m_density.value.resize( parent.size(), numConstitutivePointsPerParentIndex ); + m_density.derivs.resize( parent.size(), numConstitutivePointsPerParentIndex, m_numDOF ); + + // tjb new viscosity + m_viscosity.value.resize( parent.size(), numConstitutivePointsPerParentIndex ); + m_viscosity.derivs.resize( parent.size(), numConstitutivePointsPerParentIndex, m_numDOF ); + m_dViscosity_dPressure.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dViscosity_dTemperature.resize( parent.size(), numConstitutivePointsPerParentIndex ); - m_internalEnergy.resize( parent.size(), numConstitutivePointsPerParentIndex ); + // tjb new internal energy + m_internalEnergy.value.resize( parent.size(), numConstitutivePointsPerParentIndex ); + m_internalEnergy.derivs.resize( parent.size(), numConstitutivePointsPerParentIndex, m_numDOF ); + m_internalEnergy_n.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dInternalEnergy_dPressure.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dInternalEnergy_dTemperature.resize( parent.size(), numConstitutivePointsPerParentIndex ); - m_enthalpy.resize( parent.size(), numConstitutivePointsPerParentIndex ); + // tjb new enthalpy + m_enthalpy.value.resize( parent.size(), numConstitutivePointsPerParentIndex ); + m_enthalpy.derivs.resize( parent.size(), numConstitutivePointsPerParentIndex, m_numDOF ); + m_dEnthalpy_dPressure.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dEnthalpy_dTemperature.resize( parent.size(), numConstitutivePointsPerParentIndex ); } diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp index d5f63c8dedb..98d85ff9b0a 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidBase.hpp @@ -21,7 +21,8 @@ #define GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_SINGLEFLUIDBASE_HPP #include "constitutive/ConstitutiveBase.hpp" - +#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidUtils.hpp" namespace geos { @@ -34,6 +35,7 @@ namespace constitutive class SingleFluidBaseUpdate { public: + using SingleFluidProp = SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >; /** * @brief Get number of elements in this wrapper. @@ -58,13 +60,18 @@ class SingleFluidBaseUpdate * @param viscosity fluid viscosity * @param dVisc_dPres derivative of viscosity w.r.t. pressure */ - SingleFluidBaseUpdate( arrayView2d< real64 > const & density, + SingleFluidBaseUpdate( arrayView2d< real64 > const & density, + arrayView3d< real64 > const & dDensity, arrayView2d< real64 > const & dDens_dPres, arrayView2d< real64 > const & viscosity, + arrayView3d< real64 > const & dViscosity, arrayView2d< real64 > const & dVisc_dPres ) - : m_density( density ), + : + m_density( density ), + m_dDensity( dDensity ), m_dDens_dPres( dDens_dPres ), m_viscosity( viscosity ), + m_dViscosity( dViscosity ), m_dVisc_dPres( dVisc_dPres ) {} @@ -92,13 +99,15 @@ class SingleFluidBaseUpdate /// Fluid density - arrayView2d< real64 > m_density; + arrayView2d< real64 > m_density; + arrayView3d< real64 > m_dDensity; /// Derivative of density w.r.t. pressure arrayView2d< real64 > m_dDens_dPres; /// Fluid viscosity arrayView2d< real64 > m_viscosity; + arrayView3d< real64 > m_dViscosity; /// Derivative of viscosity w.r.t. pressure arrayView2d< real64 > m_dVisc_dPres; @@ -215,6 +224,9 @@ class SingleFluidBase : public ConstitutiveBase { public: + using SingleFluidProp = SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >; + //using SingleFluidPropConst = SingleFluidVar< real64 const, 2, constitutive::singlefluid::LAYOUT_FLUID, + // constitutive::singlefluid::LAYOUT_FLUID_DC >; /** * @brief Constructor. * @param name name of the group @@ -236,7 +248,11 @@ class SingleFluidBase : public ConstitutiveBase // *** SingleFluid-specific interface - arrayView2d< real64 const > density() const { return m_density; } + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > density() const { return m_density.value; } + arrayView2d< real64, constitutive::singlefluid::USD_FLUID > density() { return m_density.value; } + + arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DC > dDensity() const + { return m_density.derivs; } arrayView2d< real64 const > dDensity_dPressure() const { return m_dDensity_dPressure; } @@ -245,15 +261,22 @@ class SingleFluidBase : public ConstitutiveBase arrayView2d< real64 const > density_n() const { return m_density_n; } - arrayView2d< real64 const > viscosity() const { return m_viscosity; } + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > viscosity() const { return m_viscosity.value; } + arrayView2d< real64, constitutive::singlefluid::USD_FLUID > viscosity() { return m_viscosity.value; } + + arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DC > dViscosity() const + { return m_viscosity.derivs; } arrayView2d< real64 const > dViscosity_dPressure() const { return m_dViscosity_dPressure; } arrayView2d< real64 > dViscosity_dTemperature() { return m_dViscosity_dTemperature; } arrayView2d< real64 const > dViscosity_dTemperature() const { return m_dViscosity_dTemperature; } - arrayView2d< real64 > internalEnergy() { return m_internalEnergy; } - arrayView2d< real64 const > internalEnergy() const { return m_internalEnergy; } + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > internalEnergy() const { return m_internalEnergy.value; } + arrayView2d< real64, constitutive::singlefluid::USD_FLUID > internalEnergy() { return m_internalEnergy.value; } + + arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DC > dInternalEnergy() const + { return m_internalEnergy.derivs; } arrayView2d< real64 > internalEnergy_n() { return m_internalEnergy_n; } arrayView2d< real64 const > internalEnergy_n() const { return m_internalEnergy_n; } @@ -264,8 +287,11 @@ class SingleFluidBase : public ConstitutiveBase arrayView2d< real64 > dInternalEnergy_dTemperature() { return m_dInternalEnergy_dTemperature; } arrayView2d< real64 const > dInternalEnergy_dTemperature() const { return m_dInternalEnergy_dTemperature; } - arrayView2d< real64 > enthalpy() { return m_enthalpy; } - arrayView2d< real64 const > enthalpy() const { return m_enthalpy; } + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > enthalpy() const { return m_enthalpy.value; } + arrayView2d< real64, constitutive::singlefluid::USD_FLUID > enthalpy() { return m_enthalpy.value; } + + arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DC > dEnthalpy() const + { return m_enthalpy.derivs; } arrayView2d< real64 > dEnthalpy_dPressure() { return m_dEnthalpy_dPressure; } arrayView2d< real64 const > dEnthalpy_dPressure() const { return m_dEnthalpy_dPressure; } @@ -288,23 +314,30 @@ class SingleFluidBase : public ConstitutiveBase virtual void postInputInitialization() override; + // Degrees of freedom in fluid characterization + integer m_numDOF; //START_SPHINX_INCLUDE_00 - array2d< real64 > m_density; + SingleFluidProp m_density; + SingleFluidProp m_viscosity; + SingleFluidProp m_internalEnergy; + SingleFluidProp m_enthalpy; + + //array2d< real64 > m_density; array2d< real64 > m_dDensity_dPressure; array2d< real64 > m_dDensity_dTemperature; array2d< real64 > m_density_n; - array2d< real64 > m_viscosity; + //array2d< real64 > m_viscosity; array2d< real64 > m_dViscosity_dPressure; array2d< real64 > m_dViscosity_dTemperature; - array2d< real64 > m_internalEnergy; + //array2d< real64 > m_internalEnergy; array2d< real64 > m_internalEnergy_n; array2d< real64 > m_dInternalEnergy_dPressure; array2d< real64 > m_dInternalEnergy_dTemperature; - array2d< real64 > m_enthalpy; + //array2d< real64 > m_enthalpy; array2d< real64 > m_dEnthalpy_dPressure; array2d< real64 > m_dEnthalpy_dTemperature; //END_SPHINX_INCLUDE_00 diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidFields.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidFields.hpp index 61949828116..9c51c78ae38 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidFields.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidFields.hpp @@ -21,7 +21,8 @@ #define GEOS_CONSTITUTIVE_FLUID_SINGLEFLUIDFIELDS_HPP_ #include "mesh/MeshFields.hpp" - +#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidUtils.hpp" namespace geos { @@ -31,14 +32,25 @@ namespace fields namespace singlefluid { +using array2dLayoutFluid = array2d< real64, constitutive::singlefluid::LAYOUT_FLUID >; +using array3dLayoutFluid_dC = array3d< real64, constitutive::singlefluid::LAYOUT_FLUID_DC >; + DECLARE_FIELD( density, "density", - array2d< real64 >, + array2dLayoutFluid, 0, LEVEL_0, WRITE_AND_READ, "Density" ); +DECLARE_FIELD( dDensity, + "dDensity", + array3dLayoutFluid_dC, + 0, + LEVEL_0, + WRITE_AND_READ, + "dDensity" ); + DECLARE_FIELD( dDensity_dPressure, "dDensity_dPressure", array2d< real64 >, @@ -65,12 +77,20 @@ DECLARE_FIELD( density_n, DECLARE_FIELD( viscosity, "viscosity", - array2d< real64 >, + array2dLayoutFluid, 0, - NOPLOT, + LEVEL_0, WRITE_AND_READ, "Viscosity" ); +DECLARE_FIELD( dViscosity, + "dViscosity", + array3dLayoutFluid_dC, + 0, + LEVEL_0, + WRITE_AND_READ, + "dViscosity" ); + DECLARE_FIELD( dViscosity_dPressure, "dViscosity_dPressure", array2d< real64 >, @@ -87,13 +107,22 @@ DECLARE_FIELD( dViscosity_dTemperature, WRITE_AND_READ, "Derivative of viscosity with respect to temperature" ); + DECLARE_FIELD( internalEnergy, "internalEnergy", - array2d< real64 >, + array2dLayoutFluid, 0, - NOPLOT, + LEVEL_0, + WRITE_AND_READ, + "InternalEnergy" ); + +DECLARE_FIELD( dInternalEnergy, + "dInternalEnergy", + array3dLayoutFluid_dC, + 0, + LEVEL_0, WRITE_AND_READ, - "Internal energy" ); + "dInternalEnergy" ); DECLARE_FIELD( internalEnergy_n, "internalEnergy_n", @@ -121,12 +150,20 @@ DECLARE_FIELD( dInternalEnergy_dTemperature, DECLARE_FIELD( enthalpy, "enthalpy", - array2d< real64 >, + array2dLayoutFluid, 0, - NOPLOT, + LEVEL_0, WRITE_AND_READ, "Enthalpy" ); +DECLARE_FIELD( dEnthalpy, + "dEnthalpy", + array3dLayoutFluid_dC, + 0, + LEVEL_0, + WRITE_AND_READ, + "dEnthalpy" ); +#if 1 DECLARE_FIELD( dEnthalpy_dPressure, "dEnthalpy_dPressure", array2d< real64 >, @@ -142,7 +179,7 @@ DECLARE_FIELD( dEnthalpy_dTemperature, NOPLOT, WRITE_AND_READ, "Derivative of enthalpy with respect to temperature" ); - +#endif } } diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidLayouts.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidLayouts.hpp new file mode 100644 index 00000000000..360b9a786b3 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidLayouts.hpp @@ -0,0 +1,131 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 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 Layouts.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_LAYOUTS_HPP +#define GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_LAYOUTS_HPP + +#include "common/DataTypes.hpp" +#include "common/GeosxConfig.hpp" + +#include "LvArray/src/typeManipulation.hpp" +#include "RAJA/RAJA.hpp" + +namespace geos +{ +namespace constitutive +{ + +namespace singlefluid +{ +struct DerivativeOffset +{ + /// index of derivative wrt pressure + static integer constexpr dP = 0; + /// index of derivative wrt temperature + static integer constexpr dT = 1; + +}; + + +/// indices of pressure, temperature +template< integer IS_THERMAL > +struct DerivativeOffsetC {}; + +template<> +struct DerivativeOffsetC< 1 > +{ + /// index of derivative wrt pressure + static integer constexpr dP = 0; + /// index of derivative wrt temperature + static integer constexpr dT = dP + 1; + /// number of derivatives + static integer constexpr nDer = 2; +}; +template<> +struct DerivativeOffsetC< 0 > +{ + /// index of derivative wrt pressure + static integer constexpr dP = 0; + /// number of derivatives + static integer constexpr nDer = 1; +}; + +#if defined( GEOS_USE_DEVICE ) + +/// Constitutive model phase property array layout +using LAYOUT_PHASE = RAJA::PERM_JKI; +/// Constitutive model phase property compositional derivative array layout +using LAYOUT_PHASE_DC = RAJA::PERM_JKLI; + +/// Constitutive model phase composition array layout +using LAYOUT_PHASE_COMP = RAJA::PERM_JKLI; +/// Constitutive model phase composition compositional derivative array layout +using LAYOUT_PHASE_COMP_DC = RAJA::PERM_JKLMI; + +/// Constitutive model fluid property array layout +using LAYOUT_FLUID = RAJA::PERM_JI; +/// Constitutive model fluid property compositional derivative array layout +using LAYOUT_FLUID_DC = RAJA::PERM_JKI; + +///Constitutive model singe fluid property derivative array layout +using LAYOUT_SINGLEFLUID_DC = RAJA::PERMJI; + +#else + +/// Constitutive model single phase property array layout with derivatives +/// Constitutive model phase property array layout +using LAYOUT_PHASE = RAJA::PERM_IJK; +/// Constitutive model phase property compositional derivative array layout +using LAYOUT_PHASE_DC = RAJA::PERM_IJKL; + +/// Constitutive model phase composition array layout +using LAYOUT_PHASE_COMP = RAJA::PERM_IJKL; +/// Constitutive model phase composition compositional derivative array layout +using LAYOUT_PHASE_COMP_DC = RAJA::PERM_IJKLM; + +/// Constitutive model fluid property array layout +using LAYOUT_FLUID = RAJA::PERM_IJ; +/// Constitutive model fluid property compositional derivative array layout +using LAYOUT_FLUID_DC = RAJA::PERM_IJK; + +/// Constitutive model singe fluid property derivative array layout +using LAYOUT_SINGLEFLUID_DC = RAJA::PERM_IJ; + +#endif + +/// Constitutive model phase property unit stride dimension +static constexpr int USD_PHASE = LvArray::typeManipulation::getStrideOneDimension( LAYOUT_PHASE{} ); +/// Constitutive model phase property compositional derivative unit stride dimension +static constexpr int USD_PHASE_DC = LvArray::typeManipulation::getStrideOneDimension( LAYOUT_PHASE_DC{} ); + +/// Constitutive model phase composition unit stride dimension +static constexpr int USD_PHASE_COMP = LvArray::typeManipulation::getStrideOneDimension( LAYOUT_PHASE_COMP{} ); +/// Constitutive model phase composition compositional derivative unit stride dimension +static constexpr int USD_PHASE_COMP_DC = LvArray::typeManipulation::getStrideOneDimension( LAYOUT_PHASE_COMP_DC{} ); + +/// Constitutive model fluid property unit stride dimension +static constexpr int USD_FLUID = LvArray::typeManipulation::getStrideOneDimension( LAYOUT_FLUID{} ); +/// Constitutive model fluid property compositional derivative unit stride dimension +static constexpr int USD_FLUID_DC = LvArray::typeManipulation::getStrideOneDimension( LAYOUT_FLUID_DC{} ); + +} // namespace singefluid +} // namespace constitutive +} // namespace geos + +#endif //GEOS_CONSTITUTIVE_SINGLEFLUID_LAYOUTS_HPP diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidUtils.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidUtils.hpp new file mode 100644 index 00000000000..98fdb2a733c --- /dev/null +++ b/src/coreComponents/constitutive/fluid/singlefluid/SingleFluidUtils.hpp @@ -0,0 +1,131 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 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 SingleFluidUtils.hpp + */ +#ifndef GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_SINGLEFLUIDUTILS_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_SINGLEFLUIDUTILS_HPP_ + +#include "common/DataTypes.hpp" +#include "common/DataLayouts.hpp" +#include "constitutive/fluid/multifluid/MultiFluidUtils.hpp" + +namespace geos +{ + +namespace constitutive +{ + + +/** + * @brief Helper struct used to represent a variable and its compositional derivatives + * @tparam DIM number of dimensions + */ +template< typename T, int DIM, int USD, int USD_DC > +struct SingleFluidVarSlice +{ + GEOS_HOST_DEVICE + SingleFluidVarSlice( internal::ArraySliceOrRef< T, DIM, USD > inputValue, + internal::ArraySliceOrRef< T, DIM+1, USD_DC > inputDerivs ): + value( inputValue ), + derivs( inputDerivs ) + {} + + internal::ArraySliceOrRef< T, DIM, USD > value; /// variable value + internal::ArraySliceOrRef< T, DIM + 1, USD_DC > derivs; /// derivative w.r.t. pressure, temperature, compositions + + using ValueType = internal::ArraySliceOrRef< T, DIM, USD >; +}; + +/** + * @brief Struct holding views into fluid data, used to simplify parameter passing in kernel wrapper constructors. + * @tparam NDIM number of dimensions + * @tparam USD unit-stride-dim of primary property + * @tparam USD_DC unit-stride-dim of derivatives + */ +template< typename T, int NDIM, int USD, int USD_DC > +struct SingleFluidVarView +{ + SingleFluidVarView() = default; + + GEOS_HOST_DEVICE + SingleFluidVarView ( SingleFluidVarView const & src ): + value( src.value ), + derivs( src.derivs ) + {} + + GEOS_HOST_DEVICE + SingleFluidVarView ( ArrayView< T, NDIM, USD > const & valueSrc, + ArrayView< T, NDIM + 1, USD_DC > const & derivsSrc ): + value( valueSrc ), + derivs( derivsSrc ) + {} + + ArrayView< T, NDIM, USD > value; ///< View into property values + ArrayView< T, NDIM + 1, USD_DC > derivs; ///< View into property derivatives w.r.t. pressure, temperature, compositions + + using SliceType = SingleFluidVarSlice< T, NDIM - 2, USD - 2, USD_DC - 2 >; + + using ValueType = ArrayView< T, NDIM, USD >; + + GEOS_HOST_DEVICE + SliceType operator()( localIndex const k, localIndex const q ) const + { + return { value[k][q], derivs[k][q] }; + } +}; + +/** + * @brief Struct holding views into fluid data, used to simplify parameter passing in kernel wrapper constructors. + * @tparam NDIM number of dimensions + * @tparam PERM unit-stride-dim of primary property + * @tparam PERM_DC unit-stride-dim of derivatives + */ +template< typename T, int NDIM, typename PERM, typename PERM_DC > +struct SingleFluidVar +{ + Array< real64, NDIM, PERM > value; ///< Property values + Array< real64, NDIM + 1, PERM_DC > derivs; ///< Property derivatives w.r.t. pressure, temperature, compositions + + using ViewType = SingleFluidVarView< T, NDIM, getUSD< PERM >, getUSD< PERM_DC > >; + using ViewTypeConst = SingleFluidVarView< T const, NDIM, getUSD< PERM >, getUSD< PERM_DC > >; + + using SliceType = typename ViewType::SliceType; + using SliceTypeConst = typename ViewTypeConst::SliceType; + + using ValueType = Array< real64, NDIM, PERM >; + template< int MAXSIZE > + using StackValueType = StackArray< real64, NDIM, MAXSIZE, PERM >; + using ViewValueType = typename ViewType::ValueType; + using SliceValueType = typename SliceType::ValueType; + + ViewType toView() + { + return { value.toView(), derivs.toView() }; + } + + ViewTypeConst toViewConst() const + { + return { value.toViewConst(), derivs.toViewConst() }; + } +}; + + +} // namespace constitutive + +} // namespace geos + +#endif //GEOS_CONSTITUTIVE_FLUID_SINGLEFLUIDUTILS_HPP_ diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.cpp b/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.cpp index c51906552a5..25e9b2b8811 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.cpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.cpp @@ -61,10 +61,12 @@ SlurryFluidBase::SlurryFluidBase( string const & name, Group * const parent ): setInputFlag( InputFlags::OPTIONAL ). setDescription( "Flow consistency index" ); + // these would be in dDensity registerField( fields::slurryfluid::dDensity_dProppantConcentration{}, &m_dDensity_dProppantConc ); registerField( fields::slurryfluid::dDensity_dComponentConcentration{}, &m_dDensity_dCompConc ); - registerField( fields::slurryfluid::fluidDensity{}, &m_fluidDensity ); + + registerField( fields::slurryfluid::fluidDensity{}, &m_fluidDensity.value ); registerField( fields::slurryfluid::dFluidDensity_dPressure{}, &m_dFluidDens_dPres ); registerField( fields::slurryfluid::dFluidDensity_dComponentConcentration{}, &m_dFluidDens_dCompConc ); @@ -107,12 +109,15 @@ localIndex SlurryFluidBase::numFluidComponents() const void SlurryFluidBase::allocateConstitutiveData( Group & parent, localIndex const numConstitutivePointsPerParentIndex ) { + localIndex const NC = numFluidComponents(); + m_numDOF = 2 + NC; // pressure,proppantconc, NC compconc + SingleFluidBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex ); this->resize( parent.size() ); - localIndex const NC = numFluidComponents(); + // tjb these are also sized in m_dDenisty in base class m_dDensity_dProppantConc.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dDensity_dCompConc.resize( parent.size(), numConstitutivePointsPerParentIndex, NC ); @@ -120,7 +125,7 @@ void SlurryFluidBase::allocateConstitutiveData( Group & parent, m_dCompDens_dPres.resize( parent.size(), numConstitutivePointsPerParentIndex, NC ); m_dCompDens_dCompConc.resize( parent.size(), numConstitutivePointsPerParentIndex, NC, NC ); - m_fluidDensity.resize( parent.size(), numConstitutivePointsPerParentIndex ); + m_fluidDensity.value.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dFluidDens_dPres.resize( parent.size(), numConstitutivePointsPerParentIndex ); m_dFluidDens_dCompConc.resize( parent.size(), numConstitutivePointsPerParentIndex, NC ); diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.hpp index 2ab967932c9..ab8363ab45d 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidBase.hpp @@ -34,6 +34,7 @@ namespace constitutive class SlurryFluidBaseUpdate { public: + using SingleFluidProp = SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >; /** * @brief Get number of elements in this wrapper. @@ -94,6 +95,7 @@ class SlurryFluidBaseUpdate arrayView1d< real64 const > const & Ks, bool const isNewtonianFluid, arrayView2d< real64 > const & density, + arrayView3d< real64 > const & dDensity, arrayView2d< real64 > const & dDens_dPres, arrayView2d< real64 > const & dDens_dProppantConc, arrayView3d< real64 > const & dDens_dCompConc, @@ -107,6 +109,7 @@ class SlurryFluidBaseUpdate arrayView2d< real64 > const & dFluidVisc_dPres, arrayView3d< real64 > const & dFluidVisc_dCompConc, arrayView2d< real64 > const & viscosity, + arrayView3d< real64 > const & dViscosity, arrayView2d< real64 > const & dVisc_dPres, arrayView2d< real64 > const & dVisc_dProppantConc, arrayView3d< real64 > const & dVisc_dCompConc ) @@ -117,6 +120,7 @@ class SlurryFluidBaseUpdate m_Ks( Ks ), m_isNewtonianFluid( isNewtonianFluid ), m_density( density ), + m_dDensity( dDensity ), m_dDensity_dPressure( dDens_dPres ), m_dDensity_dProppantConc( dDens_dProppantConc ), m_dDensity_dCompConc( dDens_dCompConc ), @@ -130,6 +134,7 @@ class SlurryFluidBaseUpdate m_dFluidVisc_dPres( dFluidVisc_dPres ), m_dFluidVisc_dCompConc( dFluidVisc_dCompConc ), m_viscosity( viscosity ), + m_dViscosity( dViscosity ), m_dViscosity_dPressure( dVisc_dPres ), m_dViscosity_dProppantConc( dVisc_dProppantConc ), m_dViscosity_dCompConc( dVisc_dCompConc ) @@ -169,6 +174,8 @@ class SlurryFluidBaseUpdate bool m_isNewtonianFluid; arrayView2d< real64 > m_density; + arrayView3d< real64 > m_dDensity; + arrayView2d< real64 > m_dDensity_dPressure; arrayView2d< real64 > m_dDensity_dProppantConc; arrayView3d< real64 > m_dDensity_dCompConc; @@ -177,7 +184,7 @@ class SlurryFluidBaseUpdate arrayView3d< real64 > m_dCompDens_dPres; arrayView4d< real64 > m_dCompDens_dCompConc; - arrayView2d< real64 > m_fluidDensity; + arrayView2d< real64, constitutive::singlefluid::USD_FLUID > m_fluidDensity; arrayView2d< real64 > m_dFluidDens_dPres; arrayView3d< real64 > m_dFluidDens_dCompConc; @@ -186,6 +193,8 @@ class SlurryFluidBaseUpdate arrayView3d< real64 > m_dFluidVisc_dCompConc; arrayView2d< real64 > m_viscosity; + arrayView3d< real64 > m_dViscosity; + arrayView2d< real64 > m_dViscosity_dPressure; arrayView2d< real64 > m_dViscosity_dProppantConc; arrayView3d< real64 > m_dViscosity_dCompConc; @@ -247,6 +256,7 @@ class SlurryFluidBaseUpdate class SlurryFluidBase : public SingleFluidBase { public: + using SingleFluidProp = SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >; SlurryFluidBase( string const & name, Group * const parent ); @@ -269,8 +279,8 @@ class SlurryFluidBase : public SingleFluidBase arrayView3d< real64 > dDensity_dComponentConcentration() { return m_dDensity_dCompConc; } arrayView3d< real64 const > dDensity_dComponentConcentration() const { return m_dDensity_dCompConc; } - arrayView2d< real64 > fluidDensity() { return m_fluidDensity; } - arrayView2d< real64 const > fluidDensity() const { return m_fluidDensity; } + arrayView2d< real64, constitutive::singlefluid::USD_FLUID > fluidDensity() { return m_fluidDensity.value; } + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > fluidDensity() const { return m_fluidDensity.value; } arrayView2d< real64 > dFluidDensity_dPressure() { return m_dFluidDens_dPres; } arrayView2d< real64 const > dFluidDensity_dPressure() const { return m_dFluidDens_dPres; } @@ -317,6 +327,7 @@ class SlurryFluidBase : public SingleFluidBase array1d< real64 > m_nIndices; array1d< real64 > m_Ks; + // these can be folded into base but specializations will need input into DOFS for sizing array2d< real64 > m_dDensity_dProppantConc; array3d< real64 > m_dDensity_dCompConc; @@ -324,7 +335,7 @@ class SlurryFluidBase : public SingleFluidBase array3d< real64 > m_dCompDens_dPres; array4d< real64 > m_dCompDens_dCompConc; - array2d< real64 > m_fluidDensity; + SingleFluidProp m_fluidDensity; array2d< real64 > m_dFluidDens_dPres; array3d< real64 > m_dFluidDens_dCompConc; diff --git a/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidFields.hpp b/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidFields.hpp index fe45d82982f..85af1661f61 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidFields.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/SlurryFluidFields.hpp @@ -21,7 +21,8 @@ #define GEOS_CONSTITUTIVE_FLUID_SINGLEFLUID_SLURRYFLUIDFIELDS_HPP_ #include "mesh/MeshFields.hpp" - +#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidUtils.hpp" namespace geos { @@ -30,6 +31,8 @@ namespace fields namespace slurryfluid { +using array2dLayoutFluid = array2d< real64, constitutive::singlefluid::LAYOUT_FLUID >; +using array3dLayoutFluid_dC = array3d< real64, constitutive::singlefluid::LAYOUT_FLUID_DC >; DECLARE_FIELD( dDensity_dProppantConcentration, "dDens_dProppantConc", @@ -73,7 +76,15 @@ DECLARE_FIELD( dComponentDensity_dComponentConcentration, DECLARE_FIELD( fluidDensity, "FluidDensity", - array2d< real64 >, + array2dLayoutFluid, + 0, + LEVEL_0, + WRITE_AND_READ, + "Fluid density" ); + +DECLARE_FIELD( dFluidDensity, + "dFluidDensity", + array3dLayoutFluid_dC, 0, LEVEL_0, WRITE_AND_READ, diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp index 181af8bcba2..595d117f7a9 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.cpp @@ -34,7 +34,7 @@ ThermalCompressibleSinglePhaseFluid::ThermalCompressibleSinglePhaseFluid( string m_internalEnergyModelType( ExponentApproximationType::Linear ) { m_densityModelType = ExponentApproximationType::Full; - + m_numDOF=2; registerWrapper( viewKeyStruct::thermalExpansionCoeffString(), &m_thermalExpansionCoeff ). setApplyDefaultValue( 0.0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -69,7 +69,7 @@ void ThermalCompressibleSinglePhaseFluid::allocateConstitutiveData( dataReposito { CompressibleSinglePhaseFluid::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex ); - m_internalEnergy.setValues< serialPolicy >( m_referenceInternalEnergy ); + m_internalEnergy.value.setValues< serialPolicy >( m_referenceInternalEnergy ); } void ThermalCompressibleSinglePhaseFluid::postInputInitialization() @@ -103,16 +103,20 @@ ThermalCompressibleSinglePhaseFluid::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, + m_density.value, + m_density.derivs, m_dDensity_dPressure, m_dDensity_dTemperature, - m_viscosity, + m_viscosity.value, + m_viscosity.derivs, m_dViscosity_dPressure, m_dViscosity_dTemperature, - m_internalEnergy, + m_internalEnergy.value, + m_internalEnergy.derivs, m_dInternalEnergy_dPressure, m_dInternalEnergy_dTemperature, - m_enthalpy, + m_enthalpy.value, + m_enthalpy.derivs, m_dEnthalpy_dPressure, m_dEnthalpy_dTemperature, m_referenceInternalEnergy ); diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp index 01d7a802e6d..f3415db180f 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp @@ -50,28 +50,36 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate ViscRelationType const & viscRelation, IntEnergyRelationType const & intEnergyRelation, arrayView2d< real64 > const & density, + arrayView3d< real64 > const & dDensity, arrayView2d< real64 > const & dDens_dPres, arrayView2d< real64 > const & dDens_dTemp, arrayView2d< real64 > const & viscosity, + arrayView3d< real64 > const & dViscosity, arrayView2d< real64 > const & dVisc_dPres, arrayView2d< real64 > const & dVisc_dTemp, arrayView2d< real64 > const & internalEnergy, + arrayView3d< real64 > const & dInternalEnergy, arrayView2d< real64 > const & dIntEnergy_dPres, arrayView2d< real64 > const & dIntEnergy_dTemp, arrayView2d< real64 > const & enthalpy, + arrayView3d< real64 > const & dEnthalpy, arrayView2d< real64 > const & dEnthalpy_dPres, arrayView2d< real64 > const & dEnthalpy_dTemp, real64 const & refIntEnergy ) : SingleFluidBaseUpdate( density, + dDensity, dDens_dPres, viscosity, + dViscosity, dVisc_dPres ), m_dDens_dTemp( dDens_dTemp ), m_dVisc_dTemp( dVisc_dTemp ), m_internalEnergy( internalEnergy ), + m_dInternalEnergy( dInternalEnergy ), m_dIntEnergy_dPres( dIntEnergy_dPres ), m_dIntEnergy_dTemp( dIntEnergy_dTemp ), m_enthalpy( enthalpy ), + m_dEnthalpy( dEnthalpy ), m_dEnthalpy_dPres( dEnthalpy_dPres ), m_dEnthalpy_dTemp( dEnthalpy_dTemp ), m_densRelation( densRelation ), @@ -147,6 +155,9 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate m_dDens_dPres[k][q], m_viscosity[k][q], m_dVisc_dPres[k][q] ); + // tjb + m_dDensity[k][q][0] = m_dDens_dPres[k][q]; + m_dViscosity[k][q][0] = m_dVisc_dPres[k][q]; } GEOS_HOST_DEVICE @@ -170,6 +181,16 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate m_enthalpy[k][q], m_dEnthalpy_dPres[k][q], m_dEnthalpy_dTemp[k][q] ); + // tjb + m_dDensity[k][q][0] = m_dDens_dPres[k][q]; + m_dDensity[k][q][1] = m_dDens_dTemp[k][q]; + m_dViscosity[k][q][0] = m_dVisc_dPres[k][q]; + m_dViscosity[k][q][1] = m_dVisc_dTemp[k][q]; + m_dInternalEnergy[k][q][0] = m_dIntEnergy_dPres[k][q]; + m_dInternalEnergy[k][q][1] = m_dIntEnergy_dTemp[k][q]; + m_dEnthalpy[k][q][0] = m_dEnthalpy_dPres[k][q]; + m_dEnthalpy[k][q][1] = m_dEnthalpy_dTemp[k][q]; + } private: @@ -182,6 +203,7 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate /// Fluid internal energy arrayView2d< real64 > m_internalEnergy; + arrayView3d< real64 > m_dInternalEnergy; /// Derivative of internal energy w.r.t. pressure arrayView2d< real64 > m_dIntEnergy_dPres; @@ -191,7 +213,7 @@ class ThermalCompressibleSinglePhaseUpdate : public SingleFluidBaseUpdate /// Fluid enthalpy arrayView2d< real64 > m_enthalpy; - + arrayView3d< real64 > m_dEnthalpy; /// Derivative of enthalpy w.r.t. pressure arrayView2d< real64 > m_dEnthalpy_dPres; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 429ea336ab7..68671bd2b20 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -38,6 +38,7 @@ #include "mainInterface/ProblemManager.hpp" #include "mesh/DomainPartition.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" +#include "physicsSolvers/KernelLaunchSelectors.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp" @@ -97,7 +98,12 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies ) { subRegion.registerField< deltaVolume >( getName() ); + //registerField( fields::flow::mobility{}, &m_mobility.value ); + //registerField( fields::flow::dMobility{}, &m_mobility.derivs ); subRegion.registerField< mobility >( getName() ); + // tjb new mobility + subRegion.registerField< dMobility >( getName()).reference().resizeDimension< 1 >( m_numDofPerCell ); + subRegion.registerField< dMobility_dPressure >( getName() ); subRegion.registerField< fields::flow::mass >( getName() ); @@ -108,6 +114,8 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies ) { subRegion.registerField< dMobility_dTemperature >( getName() ); } + + } ); elemManager.forElementSubRegions< SurfaceElementSubRegion >( regionNames, @@ -210,8 +218,10 @@ SinglePhaseBase::FluidPropViews SinglePhaseBase::getFluidProperties( Constitutiv { SingleFluidBase const & singleFluid = dynamicCast< SingleFluidBase const & >( fluid ); return { singleFluid.density(), + singleFluid.dDensity(), singleFluid.dDensity_dPressure(), singleFluid.viscosity(), + singleFluid.dViscosity(), singleFluid.dViscosity_dPressure(), singleFluid.defaultDensity(), singleFluid.defaultViscosity() }; @@ -360,6 +370,7 @@ void SinglePhaseBase::updateMobility( ObjectManagerBase & dataGroup ) const // output arrayView1d< real64 > const mob = dataGroup.getField< fields::flow::mobility >(); + arrayView2d< real64 > const dMobility = dataGroup.getField< fields::flow::dMobility >(); // tjb arrayView1d< real64 > const dMob_dPres = dataGroup.getField< fields::flow::dMobility_dPressure >(); // input @@ -368,6 +379,17 @@ void SinglePhaseBase::updateMobility( ObjectManagerBase & dataGroup ) const getConstitutiveModel< SingleFluidBase >( dataGroup, dataGroup.getReference< string >( viewKeyStruct::fluidNamesString() ) ); FluidPropViews fluidProps = getFluidProperties( fluid ); + geos::internal::kernelLaunchSelectorThermalSwitch( m_isThermal, [&] ( auto ISTHERMAL ) { + integer constexpr NUMDOF = ISTHERMAL() + 1; + singlePhaseBaseKernels::MobilityKernel::compute_value_and_derivatives< parallelDevicePolicy<>, NUMDOF >( dataGroup.size(), + fluidProps.dens, + fluidProps.dDens, + fluidProps.visc, + fluidProps.dVisc, + mob, + dMobility ); + } ); + // tjb delete these if( m_isThermal ) { arrayView1d< real64 > const dMob_dTemp = @@ -377,25 +399,43 @@ void SinglePhaseBase::updateMobility( ObjectManagerBase & dataGroup ) const singlePhaseBaseKernels::MobilityKernel::launch< parallelDevicePolicy<> >( dataGroup.size(), fluidProps.dens, + fluidProps.dDens, fluidProps.dDens_dPres, thermalFluidProps.dDens_dTemp, fluidProps.visc, + fluidProps.dVisc, fluidProps.dVisc_dPres, thermalFluidProps.dVisc_dTemp, mob, dMob_dPres, dMob_dTemp ); + for( localIndex i=0; i >( dataGroup.size(), fluidProps.dens, + fluidProps.dDens, fluidProps.dDens_dPres, fluidProps.visc, + fluidProps.dVisc, fluidProps.dVisc_dPres, mob, dMob_dPres ); + + for( localIndex i=0; i; + /** * @brief main constructor for Group Objects * @param name the name of this instantiation of Group in the repository @@ -393,8 +397,10 @@ class SinglePhaseBase : public FlowSolverBase struct FluidPropViews { arrayView2d< real64 const > const dens; ///< density + arrayView3d< real64 const > const dDens; ///< density derivatives arrayView2d< real64 const > const dDens_dPres; ///< derivative of density w.r.t. pressure arrayView2d< real64 const > const visc; ///< viscosity + arrayView3d< real64 const > const dVisc; ///< viscosity derivatives arrayView2d< real64 const > const dVisc_dPres; ///< derivative of viscosity w.r.t. pressure real64 const defaultDensity; ///< default density to use for new elements real64 const defaultViscosity; ///< default vi to use for new elements diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp index 8b3f80b95c9..b5f1b58d265 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp @@ -21,6 +21,8 @@ #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASEBASEFIELDS_HPP_ #include "mesh/MeshFields.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidUtils.hpp" namespace geos { @@ -33,6 +35,8 @@ namespace fields namespace flow { +using array2dLayoutFluid = array2d< real64, constitutive::singlefluid::LAYOUT_FLUID >; + DECLARE_FIELD( mobility, "mobility", array1d< real64 >, @@ -41,6 +45,14 @@ DECLARE_FIELD( mobility, WRITE_AND_READ, "Mobility" ); +DECLARE_FIELD( dMobility, + "dMobility", + array2dLayoutFluid, + 0, + NOPLOT, + WRITE_AND_READ, + "dMobility" ); + DECLARE_FIELD( dMobility_dPressure, "dMobility_dPressure", array1d< real64 >, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index b87ce6bcf75..e7069a22926 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -429,8 +429,10 @@ void SinglePhaseFVM< SinglePhaseProppantBase >::assembleFluxTerms( real64 const flowAccessors.get< fields::flow::pressure >(), flowAccessors.get< fields::flow::gravityCoefficient >(), fluidAccessors.get< fields::singlefluid::density >(), + fluidAccessors.get< fields::singlefluid::dDensity >(), fluidAccessors.get< fields::singlefluid::dDensity_dPressure >(), flowAccessors.get< fields::flow::mobility >(), + flowAccessors.get< fields::flow::dMobility >(), flowAccessors.get< fields::flow::dMobility_dPressure >(), permAccessors.get< fields::permeability::permeability >(), permAccessors.get< fields::permeability::dPerm_dPressure >(), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.cpp index 0a5ccd5d41e..5d0d520f636 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.cpp @@ -96,8 +96,10 @@ SinglePhaseBase::FluidPropViews SinglePhaseProppantBase::getFluidProperties( con { SlurryFluidBase const & slurryFluid = dynamicCast< SlurryFluidBase const & >( fluid ); return { slurryFluid.density(), + slurryFluid.dDensity(), slurryFluid.dDensity_dPressure(), slurryFluid.viscosity(), + slurryFluid.dViscosity(), slurryFluid.dViscosity_dPressure(), slurryFluid.getField< fields::singlefluid::density >().getDefaultValue(), slurryFluid.getField< fields::singlefluid::viscosity >().getDefaultValue() }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp index d7bcb6068f1..59a50afbe84 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/AccumulationKernels.hpp @@ -23,6 +23,8 @@ #include "common/DataTypes.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidUtils.hpp" +///#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "codingUtilities/Utilities.hpp" @@ -44,6 +46,7 @@ class AccumulationKernel { public: + using SingleFluidProp = constitutive::SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >; /// Compute time value for the number of degrees of freedom static constexpr integer numDof = NUM_DOF; @@ -51,6 +54,11 @@ class AccumulationKernel /// Compute time value for the number of equations static constexpr integer numEqn = NUM_DOF; + + /// Note: Derivative lineup only supports dP & dT, not component terms + static constexpr integer isThermal = NUM_DOF-1; + using DerivOffset = constitutive::singlefluid::DerivativeOffsetC; + /** * @brief Constructor * @param[in] rankOffset the offset of my MPI rank @@ -77,6 +85,7 @@ class AccumulationKernel m_porosity( solid.getPorosity() ), m_dPoro_dPres( solid.getDporosity_dPressure() ), m_density( fluid.density() ), + m_dDensity ( fluid.dDensity() ), m_dDensity_dPres( fluid.dDensity_dPressure() ), m_mass_n( subRegion.template getField< fields::flow::mass_n >() ), m_localMatrix( localMatrix ), @@ -161,10 +170,15 @@ class AccumulationKernel { // Residual contribution is mass conservation in the cell stack.localResidual[0] = stack.poreVolume * m_density[ei][0] - m_mass_n[ei]; - // Derivative of residual wrt to pressure in the cell - stack.localJacobian[0][0] = stack.dPoreVolume_dPres * m_density[ei][0] + m_dDensity_dPres[ei][0] * stack.poreVolume; + //std::cout << m_dDensity_dPres[ei][0]<< " " << m_dDensity[ei][0][DerivOffset::dP] << std::endl; + //tjb use DerivOffset::dP + // std::cout.flush(); + assert( fabs( m_dDensity_dPres[ei][0]-m_dDensity[ei][0][DerivOffset::dP] ) const m_dPoro_dPres; /// Views on density - arrayView2d< real64 const > const m_density; + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const m_density; + arrayView3d< real64 const, constitutive::singlefluid::USD_FLUID_DC > const m_dDensity; arrayView2d< real64 const > const m_dDensity_dPres; /// View on mass diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp index 8083ad10746..59f15d5fe71 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernel.hpp @@ -223,8 +223,10 @@ class FluxComputeKernel : public FluxComputeKernelBase m_pres, m_gravCoef, m_dens, + m_dDens, m_dDens_dPres, m_mob, + m_dMob, m_dMob_dPres, alpha, mobility, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernelBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernelBase.hpp index 5a3896fdd79..f118d3698d9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernelBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxComputeKernelBase.hpp @@ -62,22 +62,28 @@ class FluxComputeKernelBase using DofNumberAccessor = ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > >; + using SingleFluidProp = constitutive::SingleFluidVar< real64, 2, constitutive::singlefluid::LAYOUT_FLUID, constitutive::singlefluid::LAYOUT_FLUID_DC >; + using DerivOffset = constitutive::singlefluid::DerivativeOffset; + using SinglePhaseFlowAccessors = StencilAccessors< fields::ghostRank, fields::flow::pressure, fields::flow::pressure_n, fields::flow::gravityCoefficient, fields::flow::mobility, + fields::flow::dMobility, fields::flow::dMobility_dPressure >; using SinglePhaseFluidAccessors = StencilMaterialAccessors< constitutive::SingleFluidBase, fields::singlefluid::density, + fields::singlefluid::dDensity, fields::singlefluid::dDensity_dPressure >; using SlurryFluidAccessors = StencilMaterialAccessors< constitutive::SlurryFluidBase, fields::singlefluid::density, + fields::singlefluid::dDensity, fields::singlefluid::dDensity_dPressure >; using PermeabilityAccessors = @@ -120,8 +126,10 @@ class FluxComputeKernelBase m_gravCoef( singlePhaseFlowAccessors.get( fields::flow::gravityCoefficient {} ) ), m_pres( singlePhaseFlowAccessors.get( fields::flow::pressure {} ) ), m_mob( singlePhaseFlowAccessors.get( fields::flow::mobility {} ) ), + m_dMob( singlePhaseFlowAccessors.get( fields::flow::dMobility {} ) ), m_dMob_dPres( singlePhaseFlowAccessors.get( fields::flow::dMobility_dPressure {} ) ), m_dens( singlePhaseFluidAccessors.get( fields::singlefluid::density {} ) ), + m_dDens( singlePhaseFluidAccessors.get( fields::singlefluid::dDensity {} ) ), m_dDens_dPres( singlePhaseFluidAccessors.get( fields::singlefluid::dDensity_dPressure {} ) ), m_localMatrix( localMatrix ), m_localRhs( localRhs ) @@ -152,10 +160,12 @@ class FluxComputeKernelBase /// Views on fluid mobility ElementViewConst< arrayView1d< real64 const > > const m_mob; + ElementViewConst< arrayView2d< real64 const > > const m_dMob; ElementViewConst< arrayView1d< real64 const > > const m_dMob_dPres; /// Views on fluid density ElementViewConst< arrayView2d< real64 const > > const m_dens; + ElementViewConst< arrayView3d< real64 const > > const m_dDens; ElementViewConst< arrayView2d< real64 const > > const m_dDens_dPres; // Residual and jacobian diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp index 07f81dfd2f0..18e6620ed7f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/FluxKernelsHelper.hpp @@ -22,6 +22,7 @@ #include "common/DataTypes.hpp" #include "mesh/ElementRegionManager.hpp" +#include "constitutive/fluid/singlefluid/SingleFluidLayouts.hpp" namespace geos { @@ -38,6 +39,7 @@ namespace singlePhaseFluxKernelsHelper template< typename VIEWTYPE > using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + GEOS_HOST_DEVICE inline void computeSinglePhaseFlux( localIndex const ( &seri )[2], @@ -48,8 +50,10 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2], ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, ElementViewConst< arrayView2d< real64 const > > const & dens, + ElementViewConst< arrayView3d< real64 const > > const & dDens, ElementViewConst< arrayView2d< real64 const > > const & dDens_dPres, ElementViewConst< arrayView1d< real64 const > > const & mob, + ElementViewConst< arrayView2d< real64 const > > const & dMob, ElementViewConst< arrayView1d< real64 const > > const & dMob_dPres, real64 & alpha, real64 & mobility, @@ -58,6 +62,7 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2], real64 ( & dFlux_dP )[2], real64 & dFlux_dTrans ) { + using DerivOffset = constitutive::singlefluid::DerivativeOffsetC<0>; // average density real64 densMean = 0.0; real64 dDensMean_dP[2]; @@ -66,6 +71,9 @@ void computeSinglePhaseFlux( localIndex const ( &seri )[2], { densMean += 0.5 * dens[seri[ke]][sesri[ke]][sei[ke]][0]; dDensMean_dP[ke] = 0.5 * dDens_dPres[seri[ke]][sesri[ke]][sei[ke]][0]; + // tjb - add derivid + dDensMean_dP[ke] = 0.5 * dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dP]; + assert( fabs( dDens_dPres[seri[ke]][sesri[ke]][sei[ke]][0]-dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dP] ) static void launch( localIndex const size, arrayView2d< real64 const > const & dens, + arrayView3d< real64 const > const & dDens, arrayView2d< real64 const > const & dDens_dPres, arrayView2d< real64 const > const & visc, + arrayView3d< real64 const > const & dVisc, arrayView2d< real64 const > const & dVisc_dPres, arrayView1d< real64 > const & mob, arrayView1d< real64 > const & dMob_dPres ) @@ -91,21 +108,46 @@ struct MobilityKernel forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a ) { compute( dens[a][0], + dDens[a][0][0], // tjb use deriv::dp dDens_dPres[a][0], visc[a][0], + dVisc[a][0][0], // tjb use deriv::dp dVisc_dPres[a][0], mob[a], dMob_dPres[a] ); } ); } + // Generic version + template< typename POLICY, integer NUMDOF > + static void compute_value_and_derivatives( localIndex const size, + arrayView2d< real64 const > const & density, + arrayView3d< real64 const > const & dDensity, + arrayView2d< real64 const > const & viscosity, + arrayView3d< real64 const > const & dViscosity, + arrayView1d< real64 > const & mobility, + arrayView2d< real64 > const & dMobility ) + { + forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + mobility[a] = density[a][0] / viscosity[a][0]; + for( int i=0; i static void launch( localIndex const size, arrayView2d< real64 const > const & dens, + arrayView3d< real64 const > const & dDens, // tjb arrayView2d< real64 const > const & dDens_dPres, arrayView2d< real64 const > const & dDens_dTemp, arrayView2d< real64 const > const & visc, + arrayView3d< real64 const > const & dVisc, // tjb arrayView2d< real64 const > const & dVisc_dPres, arrayView2d< real64 const > const & dVisc_dTemp, arrayView1d< real64 > const & mob, @@ -114,15 +156,28 @@ struct MobilityKernel { forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const a ) { - compute( dens[a][0], - dDens_dPres[a][0], - dDens_dTemp[a][0], - visc[a][0], - dVisc_dPres[a][0], - dVisc_dTemp[a][0], - mob[a], - dMob_dPres[a], - dMob_dTemp[a] ); + old_compute( dens[a][0], + dDens[a][0][0], // tjb use deriv::dp + dDens[a][0][1], // tjb use deriv::dt + dDens_dPres[a][0], + dDens_dTemp[a][0], + visc[a][0], + dVisc[a][0][0], // tjb use deriv::dp + dVisc[a][0][1], // tjb use deriv::dt + dVisc_dPres[a][0], + dVisc_dTemp[a][0], + mob[a], + dMob_dPres[a], + dMob_dTemp[a] ); + /* + compute( dens[a][0], + dDens[a][0], + visc[a][0], + dVisc[a][0], + mob[a], + dMob_dPres[a], + dMob_dTemp[a] ); + */ } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp index 792558b5b2e..128b64f3576 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalAccumulationKernels.hpp @@ -51,10 +51,14 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU using Base::m_porosity; using Base::m_dPoro_dPres; using Base::m_density; + using Base::m_dDensity; using Base::m_dDensity_dPres; using Base::m_localMatrix; using Base::m_localRhs; + /// Note: Derivative lineup only supports dP & dT, not component terms + static constexpr integer isThermal = NUM_DOF-1; + using DerivOffset = constitutive::singlefluid::DerivativeOffsetC; /** * @brief Constructor * @param[in] rankOffset the offset of my MPI rank @@ -76,6 +80,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU m_dDensity_dTemp( fluid.dDensity_dTemperature() ), m_dPoro_dTemp( solid.getDporosity_dTemperature() ), m_internalEnergy( fluid.internalEnergy() ), + m_dInternalEnergy( fluid.dInternalEnergy() ), m_dInternalEnergy_dPres( fluid.dInternalEnergy_dPressure() ), m_dInternalEnergy_dTemp( fluid.dInternalEnergy_dTemperature() ), m_rockInternalEnergy( solid.getInternalEnergy() ), @@ -163,15 +168,30 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU // Step 2: assemble the fluid part of the accumulation term of the energy equation real64 const fluidEnergy = stack.poreVolume * m_density[ei][0] * m_internalEnergy[ei][0]; - + // tjb + assert( fabs( m_dDensity_dPres[ei][0]-m_dDensity[ei][0][DerivOffset::dP] ) const m_internalEnergy; + arrayView3d< real64 const > const m_dInternalEnergy; + + //tjb remove arrayView2d< real64 const > const m_dInternalEnergy_dPres; arrayView2d< real64 const > const m_dInternalEnergy_dTemp; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp index 8e6bc1ba8fe..9416c579a97 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp @@ -89,6 +89,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo StencilMaterialAccessors< constitutive::SingleFluidBase, fields::singlefluid::dDensity_dTemperature, fields::singlefluid::enthalpy, + fields::singlefluid::dEnthalpy, fields::singlefluid::dEnthalpy_dPressure, fields::singlefluid::dEnthalpy_dTemperature >; @@ -145,6 +146,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo m_dMob_dTemp( thermalSinglePhaseFlowAccessors.get( fields::flow::dMobility_dTemperature {} ) ), m_dDens_dTemp( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dDensity_dTemperature {} ) ), m_enthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::enthalpy {} ) ), + m_dEnthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) ), m_dEnthalpy_dPres( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy_dPressure {} ) ), m_dEnthalpy_dTemp( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy_dTemperature {} ) ), m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) ), @@ -298,6 +300,8 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo /// Views on enthalpies ElementViewConst< arrayView2d< real64 const > > const m_enthalpy; + ElementViewConst< arrayView3d< real64 const > > const m_dEnthalpy; + // tjb remove ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dPres; ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dTemp; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp index cba7c440e09..58f75f6decd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalFluxComputeKernel.hpp @@ -63,7 +63,9 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E 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 = singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, STENCILWRAPPER >; using Base::numDof; @@ -84,6 +86,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E StencilMaterialAccessors< constitutive::SingleFluidBase, fields::singlefluid::dDensity_dTemperature, fields::singlefluid::enthalpy, + fields::singlefluid::dEnthalpy, fields::singlefluid::dEnthalpy_dPressure, fields::singlefluid::dEnthalpy_dTemperature >; @@ -133,6 +136,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E m_dMob_dTemp( thermalSinglePhaseFlowAccessors.get( fields::flow::dMobility_dTemperature {} ) ), m_dDens_dTemp( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dDensity_dTemperature {} ) ), m_enthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::enthalpy {} ) ), + m_dEnthalpy( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy {} ) ), m_dEnthalpy_dPres( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy_dPressure {} ) ), m_dEnthalpy_dTemp( thermalSinglePhaseFluidAccessors.get( fields::singlefluid::dEnthalpy_dTemperature {} ) ), m_thermalConductivity( thermalConductivityAccessors.get( fields::thermalconductivity::effectiveConductivity {} ) ), @@ -185,6 +189,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E void computeFlux( localIndex const iconn, StackVariables & stack ) const { + using DerivOffset = constitutive::singlefluid::DerivativeOffsetC<1>; // *********************************************** // First, we call the base computeFlux to compute: // 1) compFlux and its derivatives (including derivatives wrt temperature), @@ -212,7 +217,10 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E for( integer ke = 0; ke < 2; ++ke ) { - real64 const dDens_dT = m_dDens_dTemp[seri[ke]][sesri[ke]][sei[ke]][0]; + //real64 const dDens_dT = m_dDens_dTemp[seri[ke]][sesri[ke]][sei[ke]][0]; + //tjb derivid + real64 const dDens_dT = m_dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT]; + assert( fabs( m_dDens_dTemp[seri[ke]][sesri[ke]][sei[ke]][0]-m_dDens[seri[ke]][sesri[ke]][sei[ke]][0][DerivOffset::dT] )= 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]; + // tjb remove + assert( fabs( m_dMob_dTemp[seri[k_up]][sesri[k_up]][sei[k_up]]-m_dMob[seri[k_up]][sesri[k_up]][sei[k_up]][DerivOffset::dT] ) > const m_enthalpy; + ElementViewConst< arrayView3d< real64 const > > const m_dEnthalpy; + // tjb remove ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dPres; ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dTemp; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/proppant/ProppantFluxKernels.cpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/proppant/ProppantFluxKernels.cpp index 22dec62d10f..59dc46b5572 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/proppant/ProppantFluxKernels.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/proppant/ProppantFluxKernels.cpp @@ -39,8 +39,10 @@ void FaceElementFluxKernel:: ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, ElementViewConst< arrayView2d< real64 const > > const & dens, + ElementViewConst< arrayView3d< real64 const > > const & dDens, ElementViewConst< arrayView2d< real64 const > > const & dDens_dPres, ElementViewConst< arrayView1d< real64 const > > const & mob, + ElementViewConst< arrayView2d< real64 const > > const & dMob, ElementViewConst< arrayView1d< real64 const > > const & dMob_dPres, ElementViewConst< arrayView3d< real64 const > > const & permeability, ElementViewConst< arrayView3d< real64 const > > const & dPerm_dPres, @@ -103,8 +105,10 @@ void FaceElementFluxKernel:: pres, gravCoef, dens, + dDens, dDens_dPres, mob, + dMob, dMob_dPres, dt, localFlux, @@ -154,8 +158,10 @@ FaceElementFluxKernel::compute( localIndex const numFluxElems, ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, ElementViewConst< arrayView2d< real64 const > > const & dens, + ElementViewConst< arrayView3d< real64 const > > const & dDens, ElementViewConst< arrayView2d< real64 const > > const & dDens_dPres, ElementViewConst< arrayView1d< real64 const > > const & mob, + ElementViewConst< arrayView2d< real64 const > > const & dMob, ElementViewConst< arrayView1d< real64 const > > const & dMob_dPres, real64 const dt, arraySlice1d< real64 > const & flux, @@ -187,8 +193,10 @@ FaceElementFluxKernel::compute( localIndex const numFluxElems, pres, gravCoef, dens, + dDens, dDens_dPres, mob, + dMob, dMob_dPres, alpha, mobility, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/proppant/ProppantFluxKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/proppant/ProppantFluxKernels.hpp index 69b50bb2f82..48235e9c6b3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/proppant/ProppantFluxKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/proppant/ProppantFluxKernels.hpp @@ -73,8 +73,10 @@ struct FaceElementFluxKernel ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, ElementViewConst< arrayView2d< real64 const > > const & dens, + ElementViewConst< arrayView3d< real64 const > > const & dDens, ElementViewConst< arrayView2d< real64 const > > const & dDens_dPres, ElementViewConst< arrayView1d< real64 const > > const & mob, + ElementViewConst< arrayView2d< real64 const > > const & dMob, ElementViewConst< arrayView1d< real64 const > > const & dMob_dPres, ElementViewConst< arrayView3d< real64 const > > const & permeability, ElementViewConst< arrayView3d< real64 const > > const & dPerm_dPres, @@ -103,8 +105,10 @@ struct FaceElementFluxKernel ElementViewConst< arrayView1d< real64 const > > const & pres, ElementViewConst< arrayView1d< real64 const > > const & gravCoef, ElementViewConst< arrayView2d< real64 const > > const & dens, + ElementViewConst< arrayView3d< real64 const > > const & dDens, ElementViewConst< arrayView2d< real64 const > > const & dDens_dPres, ElementViewConst< arrayView1d< real64 const > > const & mob, + ElementViewConst< arrayView2d< real64 const > > const & dMob, ElementViewConst< arrayView1d< real64 const > > const & dMob_dPres, real64 const dt, arraySlice1d< real64 > const & flux, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp index d2a2b413126..35d59489f6f 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsConformingFractures.hpp @@ -60,8 +60,10 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe using AbstractBase::m_gravCoef; using AbstractBase::m_pres; using AbstractBase::m_mob; + using AbstractBase::m_dMob; using AbstractBase::m_dMob_dPres; using AbstractBase::m_dens; + using AbstractBase::m_dDens; using AbstractBase::m_dDens_dPres; using Base = singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper >; @@ -197,8 +199,10 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe m_pres, m_gravCoef, m_dens, + m_dDens, m_dDens_dPres, m_mob, + m_dMob, m_dMob_dPres, alpha, mobility, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEmbeddedFractures.hpp index 2cb1dec4ecc..c1485f398bf 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEmbeddedFractures.hpp @@ -59,8 +59,10 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe using AbstractBase::m_gravCoef; using AbstractBase::m_pres; using AbstractBase::m_mob; + using AbstractBase::m_dMob; using AbstractBase::m_dMob_dPres; using AbstractBase::m_dens; + using AbstractBase::m_dDens; using AbstractBase::m_dDens_dPres; using Base = singlePhaseFVMKernels::FluxComputeKernel< NUM_EQN, NUM_DOF, SurfaceElementStencilWrapper >; @@ -187,8 +189,10 @@ class ConnectorBasedAssemblyKernel : public singlePhaseFVMKernels::FluxComputeKe m_pres, m_gravCoef, m_dens, + m_dDens, m_dDens_dPres, m_mob, + m_dMob, m_dMob_dPres, alpha, mobility, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp index 03690cbf27d..4997b2c8d20 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEmbeddedFractures.hpp @@ -332,6 +332,7 @@ class ConnectorBasedAssemblyKernel : public singlePhasePoromechanicsEmbeddedFrac /// Views on enthalpies ElementViewConst< arrayView2d< real64 const > > const m_enthalpy; + ElementViewConst< arrayView3d< real64 const > > const m_dEnthalpy; ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dPres; ElementViewConst< arrayView2d< real64 const > > const m_dEnthalpy_dTemp; diff --git a/src/coreComponents/unitTests/fluidFlowTests/testSinglePhaseMobilityKernel.cpp b/src/coreComponents/unitTests/fluidFlowTests/testSinglePhaseMobilityKernel.cpp index b48cb9d4cf8..043a6720dac 100644 --- a/src/coreComponents/unitTests/fluidFlowTests/testSinglePhaseMobilityKernel.cpp +++ b/src/coreComponents/unitTests/fluidFlowTests/testSinglePhaseMobilityKernel.cpp @@ -30,8 +30,10 @@ TEST( SinglePhaseBaseKernels, mobility ) int constexpr NTEST = 3; real64 const dens[NTEST] = { 800.0, 1000.0, 1500.0 }; + real64 const dDens_dPres[NTEST] = { 1e-5, 1e-10, 0.0 }; real64 const visc[NTEST] = { 5.0, 2.0, 1.0 }; + real64 const dVisc_dPres[NTEST] = { 1e-7, 0.0, 0.0 }; for( int i = 0; i < NTEST; ++i ) @@ -41,7 +43,7 @@ TEST( SinglePhaseBaseKernels, mobility ) real64 mob; real64 dMob_dPres; - MobilityKernel::compute( dens[i], dDens_dPres[i], visc[i], dVisc_dPres[i], mob, dMob_dPres ); + MobilityKernel::compute( dens[i], dDens_dPres[i], dDens_dPres[i], visc[i], dVisc_dPres[i], dVisc_dPres[i], mob, dMob_dPres ); // compute etalon real64 const mob_et = dens[i] / visc[i];