Add energy balance equation to Compositional solver.

inputFiles/thermalMultiphaseFlow/co2_thermal_2d.xml

@@ -0,0 +1,223 @@
+<?xml version="1.0" ?>
+
+<Problem>
+  <Solvers>
+    <CompositionalMultiphaseFVM
+      name="compflow"
+      logLevel="1"
+      discretization="fluidTPFA"
+      temperature="368.15"
+      useMass="1"
+      isThermal="1"
+      initialDt="1e3"
+      maxCompFractionChange="0.5"
+      targetRegions="{ region }">
+      <NonlinearSolverParameters
+        newtonTol="1.0e-6"
+        newtonMaxIter="100"
+        lineSearchAction="None"
+        maxTimeStepCuts="5"/>
+      <LinearSolverParameters
+        directParallel="0"/>
+    </CompositionalMultiphaseFVM>
+  </Solvers>
+
+  <Mesh>
+    <InternalMesh
+      name="mesh1"
+      elementTypes="{ C3D8 }"
+      xCoords="{ 0, 100 }"
+      yCoords="{ 0, 100 }"
+      zCoords="{ 0, 1 }"
+      nx="{ 10 }"
+      ny="{ 10 }"
+      nz="{ 1 }"
+      cellBlockNames="{ cb }"/>
+  </Mesh>
+
+  <Geometry>
+    <Box
+      name="sink"
+      xMin="{ 89.99, 89.99, -0.01 }"
+      xMax="{ 101.01, 101.01, 1.01 }"/>
+
+    <Box
+      name="source"
+      xMin="{ -0.01, -0.01, -0.01 }"
+      xMax="{ 10.01, 10.01, 1.01 }"/>
+  </Geometry>
+
+  <Events
+    maxTime="1.5e5">
+    <PeriodicEvent
+      name="outputs"
+      timeFrequency="2.5e4"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="solverApplications"
+      maxEventDt="2.5e4"
+      target="/Solvers/compflow"/>
+
+    <PeriodicEvent
+      name="restarts"
+      timeFrequency="7.5e4"
+      target="/Outputs/sidreRestart"/>
+  </Events>
+
+  <NumericalMethods>
+    <FiniteVolume>
+      <TwoPointFluxApproximation
+        name="fluidTPFA"/>
+    </FiniteVolume>
+  </NumericalMethods>
+
+  <ElementRegions>
+    <CellElementRegion
+      name="region"
+      cellBlocks="{ cb }"
+      materialList="{ fluid, rock, relperm, thermalCond }"/>
+  </ElementRegions>
+
+  <Constitutive>
+
+    <CompressibleSolidConstantPermeability
+      name="rock"
+      solidModelName="nullSolid"
+      porosityModelName="rockPorosity"
+      permeabilityModelName="rockPerm"
+      solidInternalEnergyModelName="rockInternalEnergy"/>
+    <NullModel
+      name="nullSolid"/>
+    <PressurePorosity
+      name="rockPorosity"
+      defaultReferencePorosity="0.2"
+      referencePressure="0.0"
+      compressibility="1.0e-9"/>
+    <SolidInternalEnergy
+      name="rockInternalEnergy"
+      volumetricHeatCapacity="1.95e6"
+      referenceTemperature="368.15"
+      referenceInternalEnergy="0"/>
+    <ConstantPermeability
+      name="rockPerm"
+      permeabilityComponents="{ 1.0e-13, 1.0e-13, 1.0e-13 }"/>
+
+    <CO2BrinePhillipsThermalFluid
+      name="fluid"
+      phaseNames="{ gas, water }"
+      componentNames="{ co2, water }"
+      componentMolarWeight="{ 44e-3, 18e-3 }"
+      phasePVTParaFiles="{ pvtgas.txt, pvtliquid.txt }"
+      flashModelParaFile="co2flash.txt"/>
+
+    <BrooksCoreyRelativePermeability
+      name="relperm"
+      phaseNames="{ gas, water }"
+      phaseMinVolumeFraction="{ 0.0, 0.0 }"
+      phaseRelPermExponent="{ 1.5, 1.5 }"
+      phaseRelPermMaxValue="{ 0.9, 0.9 }"/>
+
+    <ConstantThermalConductivity
+      name="thermalCond"
+      phaseNames="{ gas, water }"
+      thermalConductivityComponents="{ 0.6, 0.6, 0.6 }"/>
+
+  </Constitutive>
+
+  <FieldSpecifications>
+
+    <FieldSpecification
+      name="initialPressure"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/region/cb"
+      fieldName="pressure"
+      scale="9e6"/>
+    <FieldSpecification
+      name="initialTemperature"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/region/cb"
+      fieldName="temperature"
+      scale="368.15"/>
+    <FieldSpecification
+      name="initialComposition_co2"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/region/cb"
+      fieldName="globalCompFraction"
+      component="0"
+      scale="0.005"/>
+    <FieldSpecification
+      name="initialComposition_water"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/region/cb"
+      fieldName="globalCompFraction"
+      component="1"
+      scale="0.995"/>
+
+    <FieldSpecification
+      name="sinkPressure"
+      setNames="{ sink }"       
+      objectPath="ElementRegions/region/cb"
+      fieldName="pressure"
+      scale="7e6"/>
+    <FieldSpecification
+      name="sinkTemperature"
+      setNames="{ sink }"       
+      objectPath="ElementRegions/region/cb"
+      fieldName="temperature"
+      scale="368.15"/>
+     <FieldSpecification
+      name="sinkTermComposition_co2"
+      setNames="{ sink }"
+      objectPath="ElementRegions/region/cb"
+      fieldName="globalCompFraction"
+      component="0"
+      scale="0.005"/>
+    <FieldSpecification
+      name="sinkTermComposition_water"
+      setNames="{ sink }"
+      objectPath="ElementRegions/region/cb"
+      fieldName="globalCompFraction"
+      component="1"
+      scale="0.995"/>
+
+    <FieldSpecification
+      name="sourcePressure"
+      setNames="{ source }"     
+      objectPath="ElementRegions/region/cb"
+      fieldName="pressure"
+      scale="1.45e7"/>
+    <FieldSpecification
+      name="sourceTemperature"
+      setNames="{ source }"     
+      objectPath="ElementRegions/region/cb"
+      fieldName="temperature"
+      scale="300.15"/>
+    <FieldSpecification
+      name="sourceTermComposition_co2"
+      setNames="{ source }"
+      objectPath="ElementRegions/region/cb"
+      fieldName="globalCompFraction"
+      component="0"
+      scale="0.995"/>
+    <FieldSpecification
+      name="sourceTermComposition_water"
+      setNames="{ source }"
+      objectPath="ElementRegions/region/cb"
+      fieldName="globalCompFraction"
+      component="1"
+      scale="0.005"/>
+  </FieldSpecifications>
+
+  <Outputs>
+    <VTK
+      name="vtkOutput"/>
+
+    <Restart
+      name="sidreRestart"/>
+  </Outputs>
+</Problem>

inputFiles/thermalMultiphaseFlow/co2flash.txt

@@ -0,0 +1 @@
+FlashModel CO2Solubility  1e6 7.5e7 5e5 299.15 369.15 10 0

inputFiles/thermalMultiphaseFlow/pvtgas.txt

@@ -0,0 +1,3 @@
+DensityFun SpanWagnerCO2Density 1e6 7.5e7 5e5 299.15 369.15 10
+ViscosityFun FenghourCO2Viscosity 1e6 7.5e7 5e5 299.15 369.15 10
+EnthalpyFun CO2Enthalpy 1e6 7.5e7 5e5 299.15 369.15 10

inputFiles/thermalMultiphaseFlow/pvtliquid.txt

@@ -0,0 +1,3 @@
+DensityFun PhillipsBrineDensity 1e6 7.5e7 5e5 299.15 369.15 10 0
+ViscosityFun PhillipsBrineViscosity 0
+EnthalpyFun BrineEnthalpy 1e6 7.5e7 5e5 299.15 369.15 10 0

src/coreComponents/constitutive/fluid/CO2BrineFluid.cpp

@@ -185,7 +185,7 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::createPVTModels()
     {
       string_array const strs = stringutilities::tokenize( str, " " );
 
-      if( strs[0] == toString( SubModelInputNames::DENSITY ) )
+      if( strs[0] == "DensityFun" )
       {
         if( strs[1] == PHASE1::Density::catalogName() )
         {
@@ -196,7 +196,7 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::createPVTModels()
           phase2InputParams[PHASE2::InputParamOrder::DENSITY] = strs;
         }
       }
-      else if( strs[0] == toString( SubModelInputNames::VISCOSITY ) )
+      else if( strs[0] == "ViscosityFun" )
       {
         if( strs[1] == PHASE1::Viscosity::catalogName() )
         {
@@ -207,7 +207,7 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::createPVTModels()
           phase2InputParams[PHASE2::InputParamOrder::VISCOSITY] = strs;
         }
       }
-      else if( strs[0] == toString( SubModelInputNames::ENTHALPY ) )
+      else if( strs[0] == "EnthalpyFun" )
       {
         if( strs[1] == PHASE1::Enthalpy::catalogName() )
         {

src/coreComponents/constitutive/fluid/CO2BrineFluid.hpp

@@ -569,25 +569,6 @@ CO2BrineFluid< PHASE1, PHASE2, FLASH >::KernelWrapper::
            m_totalDensity( k, q ) );
 }
 
-/// Declare strings associated with enumeration values
-/// Needed for now, because we don't use the catalogNames for input (yet)
-ENUM_STRINGS( CO2BrinePhillipsFluid::SubModelInputNames,
-              "DensityFun",
-              "ViscosityFun",
-              "EnthalpyFun" );
-ENUM_STRINGS( CO2BrinePhillipsThermalFluid::SubModelInputNames,
-              "DensityFun",
-              "ViscosityFun",
-              "EnthalpyFun" );
-ENUM_STRINGS( CO2BrineEzrokhiFluid::SubModelInputNames,
-              "DensityFun",
-              "ViscosityFun",
-              "EnthalpyFun" );
-ENUM_STRINGS( CO2BrineEzrokhiThermalFluid::SubModelInputNames,
-              "DensityFun",
-              "ViscosityFun",
-              "EnthalpyFun" );
-
 
 } // namespace constitutive
 

src/coreComponents/constitutive/fluid/MultiFluidExtrinsicData.hpp

@@ -156,7 +156,7 @@ EXTRINSIC_MESH_DATA_TRAIT( dPhaseInternalEnergy,
                            0,
                            NOPLOT,
                            NO_WRITE,
-                           "Derivative of phase internal energy with respect to pressure, temperature, and global component fraction" );
+                           "Derivative of phase internal energy with respect to pressure, temperature, and global component fractions" );
 
 EXTRINSIC_MESH_DATA_TRAIT( phaseCompFraction,
                            "phaseCompFraction",

src/coreComponents/constitutive/fluid/multiFluidSelector.hpp

@@ -44,7 +44,10 @@ void constitutiveUpdatePassThru( MultiFluidBase const & fluid,
                                CompositionalMultiphaseFluid,
 #endif
                                CO2BrinePhillipsFluid,
-                               CO2BrineEzrokhiFluid >::execute( fluid, std::forward< LAMBDA >( lambda ) );
+                               CO2BrineEzrokhiFluid,
+                               CO2BrinePhillipsThermalFluid /*, // if I uncomment the two models at the same time, the compiler segfaults on
+                                                               Lassen!
+                                                               CO2BrineEzrokhiThermalFluid*/>::execute( fluid, std::forward< LAMBDA >( lambda ) );
 }
 
 template< typename LAMBDA >
@@ -57,7 +60,10 @@ void constitutiveUpdatePassThru( MultiFluidBase & fluid,
                                CompositionalMultiphaseFluid,
 #endif
                                CO2BrinePhillipsFluid,
-                               CO2BrineEzrokhiFluid >::execute( fluid, std::forward< LAMBDA >( lambda ) );
+                               CO2BrineEzrokhiFluid,
+                               CO2BrinePhillipsThermalFluid /*, // if I uncomment the two models at the same time, the compiler segfaults on
+                                                               Lassen!
+                                                               CO2BrineEzrokhiThermalFluid*/>::execute( fluid, std::forward< LAMBDA >( lambda ) );
 }
 
 } // namespace constitutive

src/coreComponents/constitutive/solid/SolidInternalEnergy.cpp

@@ -30,34 +30,34 @@ SolidInternalEnergy::SolidInternalEnergy( string const & name, Group * const par
   ConstitutiveBase( name, parent ),
   m_internalEnergy(),
   m_dInternalEnergy_dTemperature(),
-  m_specificHeatCapacity(),
+  m_volumetricHeatCapacity(),
   m_referenceTemperature(),
   m_referenceInternalEnergy()
 {
   registerWrapper( viewKeyStruct::internalEnergyString(), &m_internalEnergy ).
     setPlotLevel( PlotLevel::LEVEL_0 ).
     setApplyDefaultValue( 0.0 ).
-    setDescription( "Internal energy of the solid" );
+    setDescription( "Internal energy of the solid per unit volume [J/m^3]" );
 
   registerWrapper( viewKeyStruct::oldInternalEnergyString(), &m_oldInternalEnergy ).
     setApplyDefaultValue( 0.0 ).
-    setDescription( "Internal energy of the solid at the previous time-step" );
+    setDescription( "Internal energy of the solid per unit volume at the previous time-step [J/m^3]" );
 
   registerWrapper( viewKeyStruct::dInternalEnergy_dTemperatureString(), &m_dInternalEnergy_dTemperature ).
     setApplyDefaultValue( 0.0 ).
-    setDescription( "Derivative of the solid internal energy w.r.t. temperature" );
+    setDescription( "Derivative of the solid internal energy w.r.t. temperature [J/(m^3.K)]" );
 
-  registerWrapper( viewKeyStruct::specificHeatCapacityString(), &m_specificHeatCapacity ).
-    setApplyDefaultValue( 0.0 ).
-    setDescription( "Solid specific heat capacity" );
+  registerWrapper( viewKeyStruct::volumetricHeatCapacityString(), &m_volumetricHeatCapacity ).
+    setInputFlag( InputFlags::REQUIRED ).
+    setDescription( "Solid volumetric heat capacity [J/(kg.K)]" );
 
   registerWrapper( viewKeyStruct::referenceTemperatureString(), &m_referenceTemperature ).
-    setApplyDefaultValue( 0.0 ).
-    setDescription( "Reference temperature" );
+    setInputFlag( InputFlags::REQUIRED ).
+    setDescription( "Reference temperature [K]" );
 
   registerWrapper( viewKeyStruct::referenceInternalEnergyString(), &m_referenceInternalEnergy ).
-    setApplyDefaultValue( 0.0 ).
-    setDescription( "Internal energy at the reference temperature" );
+    setInputFlag( InputFlags::REQUIRED ).
+    setDescription( "Internal energy at the reference temperature [J/kg]" );
 }
 
 void SolidInternalEnergy::allocateConstitutiveData( Group & parent,
@@ -81,6 +81,8 @@ void SolidInternalEnergy::saveConvergedState() const
   } );
 }
 
-}
+REGISTER_CATALOG_ENTRY( ConstitutiveBase, SolidInternalEnergy, string const &, Group * const )
 
-}
+} /* namespace constitutive */
+
+} /* namespace geosx */