Break up multi fluid test (#3140)

GEOS-DEV · Jun 11, 2024 · 4de1337 · 4de1337
1 parent 010fc34
commit 4de1337
Show file tree

Hide file tree

Showing 16 changed files with 1,746 additions and 1,347 deletions.
diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3080-5511-8210a1a
+  baseline: integratedTests/baseline_integratedTests-pr3140-5514-fe26f80
 
 allow_fail:
   all: ''

diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md
@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3140 (2024-06-11)
+======================
+Fixed derivative in EzrokhiBrineDensity
+
 PR #3080 (2024-06-07)
 =====================
 Rebaseline after adding viscoelastic wave propagator.

diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp
@@ -194,11 +194,13 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::checkTablesParameters( real64 const
 template< typename PHASE1, typename PHASE2, typename FLASH >
 void CO2BrineFluid< PHASE1, PHASE2, FLASH >::initializePreSubGroups()
 {
+#if defined(GEOS_DEVICE_COMPILE)
   GEOS_THROW_IF( this->getCatalogName() == CO2BrineEzrokhiThermalFluid::catalogName(),
                  GEOS_FMT( "The `{}` model is disabled for now. Please use the other thermal CO2-brine model instead: `{}`",
                            CO2BrineEzrokhiThermalFluid::catalogName(),
                            CO2BrinePhillipsThermalFluid::catalogName() ),
                  InputError );
+#endif
 }
 
 template< typename PHASE1, typename PHASE2, typename FLASH >

diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/EzrokhiBrineDensity.hpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/EzrokhiBrineDensity.hpp
@@ -171,6 +171,7 @@ void EzrokhiBrineDensityUpdate::compute( real64 const & pressure,
                                          arraySlice1d< real64, USD3 > const & dValue,
                                          bool useMass ) const
 {
+  constexpr integer numDof = 4;
   using Deriv = multifluid::DerivativeOffset;
 
   real64 waterSatDensity_dTemperature = 0.0;
@@ -188,32 +189,39 @@ void EzrokhiBrineDensityUpdate::compute( real64 const & pressure,
 
   real64 const coefPhaseComposition = m_coef0 + temperature * ( m_coef1 + m_coef2 * temperature );
 
+  real64 const mw_co2 = m_componentMolarWeight[m_CO2Index];
+  real64 const mw_h2o = m_componentMolarWeight[m_waterIndex];
+
   // we have to convert molar component phase fraction (phaseComposition[m_CO2Index]) to mass fraction
-  real64 const waterMWInv = 1.0 / (phaseComposition[m_CO2Index] * m_componentMolarWeight[m_CO2Index]
-                                   + phaseComposition[m_waterIndex] * m_componentMolarWeight[m_waterIndex]);
+  real64 const waterMWInv = 1.0 / (phaseComposition[m_CO2Index] * mw_co2 + phaseComposition[m_waterIndex] * mw_h2o);
+  real64 dWaterMWInv[numDof]{};
+  for( integer dof = 0; dof < numDof; ++dof )
+  {
+    dWaterMWInv[dof] = -(dPhaseComposition[m_CO2Index][dof] * mw_co2 + dPhaseComposition[m_waterIndex][dof] * mw_h2o)*waterMWInv*waterMWInv;
+  }
 
-  real64 const massPhaseCompositionCO2 = phaseComposition[m_CO2Index] * m_componentMolarWeight[m_CO2Index] * waterMWInv;
+  real64 const massPhaseCompositionCO2 = phaseComposition[m_CO2Index] * mw_co2 * waterMWInv;
 
   real64 const exponent = coefPhaseComposition * massPhaseCompositionCO2;
+  real64 dExponent[numDof]{};
+  for( integer dof = 0; dof < numDof; ++dof )
+  {
+    dExponent[dof] = coefPhaseComposition * mw_co2 *
+                     (dPhaseComposition[m_CO2Index][dof] * waterMWInv + phaseComposition[m_CO2Index] * dWaterMWInv[dof]);
+  }
+  dExponent[Deriv::dT] += ( m_coef1 + 2.0 * m_coef2 * temperature) * massPhaseCompositionCO2;
 
-  real64 const exponent_dPressure = coefPhaseComposition * dPhaseComposition[m_CO2Index][Deriv::dP];
-  real64 const exponent_dTemperature = coefPhaseComposition * dPhaseComposition[m_CO2Index][Deriv::dT] +
-                                       ( m_coef1 + 2 * m_coef2 * temperature) * massPhaseCompositionCO2;
-  // compute only common part of derivatives w.r.t. CO2 and water phase compositions
-  // later to be multiplied by (phaseComposition[m_waterIndex]) and ( -phaseComposition[m_CO2Index] ) respectively
-  real64 const exponent_dPhaseComp = coefPhaseComposition * m_componentMolarWeight[m_CO2Index] * m_componentMolarWeight[m_waterIndex] * waterMWInv * waterMWInv;
-  real64 exponentPowered = pow( 10, exponent );
+  real64 exponentPowered = pow( 10.0, exponent );
 
   value = waterDensity * exponentPowered;
 
-  real64 const dValueCoef = LvArray::math::log( 10 ) * value;
-  dValue[Deriv::dP] = dValueCoef * exponent_dPressure + waterDensity_dPressure * exponentPowered;
-  dValue[Deriv::dT] = dValueCoef * exponent_dTemperature + waterDensity_dTemperature * exponentPowered;
+  real64 const dValueCoef = LvArray::math::log( 10.0 ) * value;
+
+  dValue[Deriv::dP] = dValueCoef * dExponent[Deriv::dP] + waterDensity_dPressure * exponentPowered;
+  dValue[Deriv::dT] = dValueCoef * dExponent[Deriv::dT] + waterDensity_dTemperature * exponentPowered;
 
-  // here, we multiply common part of derivatives by specific coefficients
-  real64 const dValue_dPhaseComp = dValueCoef * exponent_dPhaseComp;
-  dValue[Deriv::dC+m_CO2Index] = dValue_dPhaseComp * phaseComposition[m_waterIndex] * dPhaseComposition[m_CO2Index][Deriv::dC+m_CO2Index];
-  dValue[Deriv::dC+m_waterIndex] = dValue_dPhaseComp * ( -phaseComposition[m_CO2Index] ) * dPhaseComposition[m_waterIndex][Deriv::dC+m_waterIndex];
+  dValue[Deriv::dC+m_CO2Index] = dValueCoef * dExponent[Deriv::dC+m_CO2Index];
+  dValue[Deriv::dC+m_waterIndex] = dValueCoef * dExponent[Deriv::dC+m_waterIndex];
 
   if( !useMass )
   {

diff --git a/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp b/src/coreComponents/constitutive/fluid/multifluid/MultiFluidSelector.hpp
@@ -47,6 +47,8 @@ void constitutiveUpdatePassThru( MultiFluidBase const & fluid,
                                CO2BrinePhillipsFluid,
                                CO2BrineEzrokhiFluid,
                                CO2BrinePhillipsThermalFluid,
+// Including these in a CUDA build will lead to compiler segfault.
+// Need to split compilation units for all the options
 #if !defined(GEOS_DEVICE_COMPILE)
                                CO2BrineEzrokhiThermalFluid,
                                CompositionalTwoPhasePengRobinsonLBCViscosity,
@@ -69,6 +71,8 @@ void constitutiveUpdatePassThru( MultiFluidBase & fluid,
                                CO2BrinePhillipsFluid,
                                CO2BrineEzrokhiFluid,
                                CO2BrinePhillipsThermalFluid,
+// Including these in a CUDA build will lead to compiler segfault.
+// Need to split compilation units for all the options"
 #if !defined(GEOS_DEVICE_COMPILE)
                                CO2BrineEzrokhiThermalFluid,
                                CompositionalTwoPhasePengRobinsonLBCViscosity,

diff --git a/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt b/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt
@@ -4,6 +4,9 @@ set( gtest_geosx_tests
      testCO2BrinePVTModels.cpp
      testCO2SpycherPruessModels.cpp
      testDamage.cpp
+     testMultiFluidCO2Brine.cpp
+     testMultiFluidDeadOil.cpp
+     testMultiFluidLiveOil.cpp
      testRelPerm.cpp
      testRelPermHysteresis.cpp )
 
@@ -39,8 +42,7 @@ endif()
 
 if( ENABLE_PVTPackage )
     list( APPEND gtest_geosx_tests
-          testMultiFluid.cpp )
-
+          testMultiFluidCompositionalMultiphasePVTPackage.cpp )
     list( APPEND dependencyList PVTPackage )
 endif()