From 65064b0ce78cdab2035d935c93a7317293da05fa Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 18 Apr 2024 17:16:41 -0500 Subject: [PATCH] Make sequential poromechanics to conserve mass (#2952) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Pavel Tomin <“paveltomin@users.noreply.github.com”> Co-authored-by: Matteo Cusini <49037133+CusiniM@users.noreply.github.com> Co-authored-by: Randolph R. Settgast Co-authored-by: Nicola Castelletto --- integratedTests | 2 +- .../fluidFlow/CompositionalMultiphaseBase.cpp | 86 ++++++++++++++++--- .../fluidFlow/CompositionalMultiphaseBase.hpp | 14 ++- .../CompositionalMultiphaseBaseFields.hpp | 25 ++++-- .../fluidFlow/FlowSolverBase.cpp | 12 +++ .../fluidFlow/FlowSolverBaseFields.hpp | 50 +++++++++-- ...rmalCompositionalMultiphaseBaseKernels.hpp | 37 +++----- .../fluidFlow/SinglePhaseBase.cpp | 83 +++++++++++++++++- .../fluidFlow/SinglePhaseBase.hpp | 23 ++++- .../fluidFlow/SinglePhaseBaseKernels.hpp | 46 ++++------ .../fluidFlow/SinglePhaseFVM.cpp | 8 +- .../fluidFlow/SinglePhaseHybridFVM.cpp | 5 -- ...rmalCompositionalMultiphaseBaseKernels.hpp | 32 +++---- .../ThermalSinglePhaseBaseKernels.hpp | 35 +++----- ...mpositionalMultiphaseReservoirAndWells.hpp | 2 +- .../multiphysics/CoupledSolver.hpp | 8 +- .../multiphysics/PoromechanicsSolver.hpp | 19 ++++ .../SinglePhaseReservoirAndWells.hpp | 2 +- 18 files changed, 341 insertions(+), 148 deletions(-) diff --git a/integratedTests b/integratedTests index 9e5112893b1..671feaaf28d 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit 9e5112893b1fc371c058a4da293771e8e18aee2b +Subproject commit 671feaaf28d51f101144145afa00a6c535669ff1 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 0513f370c70..f023c1a9ec3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -359,10 +359,15 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< dPhaseMobility >( getName() ). reference().resizeDimension< 1, 2 >( m_numPhases, m_numComponents + 2 ); // dP, dT, dC + // needed for time step selector subRegion.registerField< phaseVolumeFraction_n >( getName() ). reference().resizeDimension< 1 >( m_numPhases ); - subRegion.registerField< phaseMobility_n >( getName() ). - reference().resizeDimension< 1 >( m_numPhases ); + + subRegion.registerField< compAmount >( getName() ). + reference().resizeDimension< 1 >( m_numComponents ); + subRegion.registerField< compAmount_n >( getName() ). + reference().resizeDimension< 1 >( m_numComponents ); + } ); FaceManager & faceManager = mesh.getFaceManager(); @@ -697,6 +702,57 @@ void CompositionalMultiphaseBase::updateCapPressureModel( ObjectManagerBase & da } } +void CompositionalMultiphaseBase::updateCompAmount( ElementSubRegionBase & subRegion ) const +{ + GEOS_MARK_FUNCTION; + + string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ); + CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); + arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >(); + arrayView2d< real64, compflow::USD_COMP > const compAmount = subRegion.getField< fields::flow::compAmount >(); + + integer const numComp = m_numComponents; + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + for( integer ic = 0; ic < numComp; ++ic ) + { + compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compDens[ei][ic]; + } + } ); +} + +void CompositionalMultiphaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const +{ + GEOS_MARK_FUNCTION; + + string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ); + CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); + arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); + arrayView2d< real64 const > rockInternalEnergy = porousMaterial.getInternalEnergy(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = subRegion.getField< fields::flow::phaseVolumeFraction >(); + string const & fluidName = getConstitutiveName< MultiFluidBase >( subRegion ); + MultiFluidBase & fluid = subRegion.getConstitutiveModel< MultiFluidBase >( fluidName ); + arrayView3d< real64 const, multifluid::USD_PHASE > const phaseDens = fluid.phaseDensity(); + arrayView3d< real64 const, multifluid::USD_PHASE > const phaseInternalEnergy = fluid.phaseInternalEnergy(); + + arrayView1d< real64 > const energy = subRegion.getField< fields::flow::energy >(); + + integer const numPhases = m_numPhases; + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + energy[ei] = volume[ei] * (1.0 - porosity[ei][0]) * rockInternalEnergy[ei][0]; + for( integer ip = 0; ip < numPhases; ++ip ) + { + energy[ei] += volume[ei] * porosity[ei][0] * phaseVolFrac[ei][ip] * phaseDens[ei][0][ip] * phaseInternalEnergy[ei][0][ip]; + } + } ); +} + void CompositionalMultiphaseBase::updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const { arrayView1d< real64 const > const temp = dataGroup.getField< fields::flow::temperature >(); @@ -715,12 +771,13 @@ void CompositionalMultiphaseBase::updateSolidInternalEnergyModel( ObjectManagerB temp ); } -real64 CompositionalMultiphaseBase::updateFluidState( ObjectManagerBase & subRegion ) const +real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & subRegion ) const { GEOS_MARK_FUNCTION; updateGlobalComponentFraction( subRegion ); updateFluidModel( subRegion ); + updateCompAmount( subRegion ); real64 const maxDeltaPhaseVolFrac = updatePhaseVolumeFraction( subRegion ); updateRelPermModel( subRegion ); updatePhaseMobility( subRegion ); @@ -768,7 +825,6 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, compDens[ei][ic] = totalDens[ei][0] * compFrac[ei][ic]; } } ); - } ); // with initial component densities defined - check if they need to be corrected to avoid zero diags etc @@ -796,6 +852,7 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, porousSolid.scaleReferencePorosity( netToGross ); saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes updatePorosityAndPermeability( subRegion ); + updateCompAmount( subRegion ); updatePhaseVolumeFraction( subRegion ); // Now, we initialize and update each constitutive model one by one @@ -883,6 +940,8 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, SolidInternalEnergy const & solidInternalEnergyMaterial = getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); solidInternalEnergyMaterial.saveConvergedState(); + + updateEnergy( subRegion ); } // Step 4.7: if the diffusion and/or dispersion is/are supported, initialize the two models @@ -1251,19 +1310,13 @@ CompositionalMultiphaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM updateSolidInternalEnergyModel( subRegion ); } - // after the update, save the new saturation and phase mobilities + // after the update, save the new saturation arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = subRegion.template getField< fields::flow::phaseVolumeFraction >(); arrayView2d< real64, compflow::USD_PHASE > const phaseVolFrac_n = subRegion.template getField< fields::flow::phaseVolumeFraction_n >(); phaseVolFrac_n.setValues< parallelDevicePolicy<> >( phaseVolFrac ); - arrayView2d< real64 const, compflow::USD_PHASE > const phaseMob = - subRegion.template getField< fields::flow::phaseMobility >(); - arrayView2d< real64, compflow::USD_PHASE > const phaseMob_n = - subRegion.template getField< fields::flow::phaseMobility_n >(); - phaseMob_n.setValues< parallelDevicePolicy<> >( phaseMob ); - } ); } ); } @@ -2381,6 +2434,13 @@ void CompositionalMultiphaseBase::saveConvergedState( ElementSubRegionBase & sub arrayView2d< real64, compflow::USD_COMP > const & compDens_n = subRegion.template getField< fields::flow::globalCompDensity_n >(); compDens_n.setValues< parallelDevicePolicy<> >( compDens ); + + arrayView2d< real64 const, compflow::USD_COMP > const & compAmount = + subRegion.template getField< fields::flow::compAmount >(); + arrayView2d< real64, compflow::USD_COMP > const & compAmount_n = + subRegion.template getField< fields::flow::compAmount_n >(); + compAmount_n.setValues< parallelDevicePolicy<> >( compAmount ); + if( m_isFixedStressPoromechanicsUpdate ) { arrayView2d< real64, compflow::USD_COMP > const & compDens_k = @@ -2438,6 +2498,9 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain ) { GEOS_MARK_FUNCTION; + if( m_keepFlowVariablesConstantDuringInitStep ) + return; + real64 maxDeltaPhaseVolFrac = 0.0; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, @@ -2456,6 +2519,7 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain ) if( m_isThermal ) { updateSolidInternalEnergyModel( subRegion ); + updateEnergy( subRegion ); } } ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 728ab922bfc..594b452f270 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -133,6 +133,18 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ void updateCapPressureModel( ObjectManagerBase & dataGroup ) const; + /** + * @brief Update components mass/moles + * @param subRegion the subregion storing the required fields + */ + void updateCompAmount( ElementSubRegionBase & subRegion ) const; + + /** + * @brief Update energy + * @param subRegion the subregion storing the required fields + */ + void updateEnergy( ElementSubRegionBase & subRegion ) const; + /** * @brief Update all relevant solid internal energy models using current values of temperature * @param dataGroup the group storing the required fields @@ -145,7 +157,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const = 0; - real64 updateFluidState( ObjectManagerBase & dataGroup ) const; + real64 updateFluidState( ElementSubRegionBase & subRegion ) const; virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override final; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp index 847126b111b..2c687f7dfb8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp @@ -119,6 +119,7 @@ DECLARE_FIELD( dPhaseMobility, NO_WRITE, "Derivative of phase volume fraction with respect to pressure, temperature, global component density" ); +// this is needed for time step selector DECLARE_FIELD( phaseVolumeFraction_n, "phaseVolumeFraction_n", array2dLayoutPhase, @@ -127,14 +128,6 @@ DECLARE_FIELD( phaseVolumeFraction_n, WRITE_AND_READ, "Phase volume fraction at the previous converged time step" ); -DECLARE_FIELD( phaseMobility_n, - "phaseMobility_n", - array2dLayoutPhase, - 0, - NOPLOT, - WRITE_AND_READ, - "Phase mobility at the previous converged time step" ); - DECLARE_FIELD( phaseOutflux, "phaseOutflux", array2dLayoutPhase, @@ -175,6 +168,22 @@ DECLARE_FIELD( globalCompDensityScalingFactor, NO_WRITE, "Scaling factors for global component densities" ); +DECLARE_FIELD( compAmount, + "compAmount", + array2dLayoutComp, + 0, + LEVEL_0, + WRITE_AND_READ, + "Component amount" ); + +DECLARE_FIELD( compAmount_n, + "compAmount_n", + array2dLayoutComp, + 0, + LEVEL_0, + WRITE_AND_READ, + "Component amount at the previous converged time step" ); + } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index 3d38dd95302..a116226a893 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -160,6 +160,11 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies ) { subRegion.registerField< fields::flow::temperature_k >( getName() ); // needed for the fixed-stress porosity update } + if( m_isThermal ) + { + subRegion.registerField< fields::flow::energy >( getName() ); + subRegion.registerField< fields::flow::energy_n >( getName() ); + } } ); elemManager.forElementSubRegionsComplete< SurfaceElementSubRegion >( [&]( localIndex const, @@ -218,6 +223,13 @@ void FlowSolverBase::saveConvergedState( ElementSubRegionBase & subRegion ) cons arrayView1d< real64 > const temp_n = subRegion.template getField< fields::flow::temperature_n >(); temp_n.setValues< parallelDevicePolicy<> >( temp ); + if( m_isThermal ) + { + arrayView1d< real64 const > const energy = subRegion.template getField< fields::flow::energy >(); + arrayView1d< real64 > const energy_n = subRegion.template getField< fields::flow::energy_n >(); + energy_n.setValues< parallelDevicePolicy<> >( energy ); + } + if( m_isFixedStressPoromechanicsUpdate ) { arrayView1d< real64 > const pres_k = subRegion.template getField< fields::flow::pressure_k >(); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp index 0b17130255f..64ca6c6a682 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp @@ -96,14 +96,6 @@ DECLARE_FIELD( temperature, WRITE_AND_READ, "Temperature" ); -DECLARE_FIELD( faceTemperature, - "faceTemperature", - array1d< real64 >, - 0, - LEVEL_0, - WRITE_AND_READ, - "Face temperature" ); - DECLARE_FIELD( temperature_n, "temperature_n", array1d< real64 >, @@ -116,7 +108,7 @@ DECLARE_FIELD( temperature_k, "temperature_k", array1d< real64 >, 0, - LEVEL_0, + NOPLOT, NO_WRITE, "Temperature at the previous sequential iteration" ); @@ -128,6 +120,14 @@ DECLARE_FIELD( initialTemperature, WRITE_AND_READ, "Initial temperature" ); +DECLARE_FIELD( faceTemperature, + "faceTemperature", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Face temperature" ); + DECLARE_FIELD( netToGross, "netToGross", array1d< real64 >, @@ -240,6 +240,38 @@ DECLARE_FIELD( temperatureScalingFactor, NO_WRITE, "Scaling factors for temperature" ); +DECLARE_FIELD( mass, + "mass", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Mass" ); + +DECLARE_FIELD( mass_n, + "mass_n", + array1d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Mass at the previous converged time step" ); + +DECLARE_FIELD( energy, + "energy", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Energy" ); + +DECLARE_FIELD( energy_n, + "energy_n", + array1d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Energy at the previous converged time step" ); + } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp index c5863e9bcad..f0a3c194c5b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp @@ -587,21 +587,17 @@ class ElementBasedAssemblyKernel m_dofNumber( subRegion.getReference< array1d< globalIndex > >( dofKey ) ), m_elemGhostRank( subRegion.ghostRank() ), m_volume( subRegion.getElementVolume() ), - m_porosity_n( solid.getPorosity_n() ), m_porosity( solid.getPorosity() ), m_dPoro_dPres( solid.getDporosity_dPressure() ), m_dCompFrac_dCompDens( subRegion.getField< fields::flow::dGlobalCompFraction_dGlobalCompDensity >() ), - m_phaseVolFrac_n( subRegion.getField< fields::flow::phaseVolumeFraction_n >() ), m_phaseVolFrac( subRegion.getField< fields::flow::phaseVolumeFraction >() ), m_dPhaseVolFrac( subRegion.getField< fields::flow::dPhaseVolumeFraction >() ), - m_phaseDens_n( fluid.phaseDensity_n() ), m_phaseDens( fluid.phaseDensity() ), m_dPhaseDens( fluid.dPhaseDensity() ), - m_phaseCompFrac_n( fluid.phaseCompFraction_n() ), m_phaseCompFrac( fluid.phaseCompFraction() ), m_dPhaseCompFrac( fluid.dPhaseCompFraction() ), m_compDens( subRegion.getField< fields::flow::globalCompDensity >() ), - m_compDens_n( subRegion.getField< fields::flow::globalCompDensity_n >() ), + m_compAmount_n( subRegion.getField< fields::flow::compAmount_n >() ), m_localMatrix( localMatrix ), m_localRhs( localRhs ), m_kernelFlags( kernelFlags ) @@ -620,9 +616,6 @@ class ElementBasedAssemblyKernel /// Pore volume at time n+1 real64 poreVolume = 0.0; - /// Pore volume at the previous converged time step - real64 poreVolume_n = 0.0; - /// Derivative of pore volume with respect to pressure real64 dPoreVolume_dPres = 0.0; @@ -663,7 +656,6 @@ class ElementBasedAssemblyKernel { // initialize the pore volume stack.poreVolume = m_volume[ei] * m_porosity[ei][0]; - stack.poreVolume_n = m_volume[ei] * m_porosity_n[ei][0]; stack.dPoreVolume_dPres = m_volume[ei] * m_dPoro_dPres[ei][0]; // set row index and degrees of freedom indices for this element @@ -694,7 +686,7 @@ class ElementBasedAssemblyKernel for( integer ic = 0; ic < numComp; ++ic ) { real64 const compAmount = stack.poreVolume * m_compDens[ei][ic]; - real64 const compAmount_n = stack.poreVolume_n * m_compDens_n[ei][ic]; + real64 const compAmount_n = m_compAmount_n[ei][ic]; stack.localResidual[ic] += compAmount - compAmount_n; @@ -714,15 +706,12 @@ class ElementBasedAssemblyKernel // construct the slices for variables accessed multiple times arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > dCompFrac_dCompDens = m_dCompFrac_dCompDens[ei]; - arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac_n = m_phaseVolFrac_n[ei]; arraySlice1d< real64 const, compflow::USD_PHASE - 1 > phaseVolFrac = m_phaseVolFrac[ei]; arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > dPhaseVolFrac = m_dPhaseVolFrac[ei]; - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseDens_n = m_phaseDens_n[ei][0]; arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseDens = m_phaseDens[ei][0]; arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseDens = m_dPhaseDens[ei][0]; - arraySlice2d< real64 const, multifluid::USD_PHASE_COMP - 2 > phaseCompFrac_n = m_phaseCompFrac_n[ei][0]; arraySlice2d< real64 const, multifluid::USD_PHASE_COMP - 2 > phaseCompFrac = m_phaseCompFrac[ei][0]; arraySlice3d< real64 const, multifluid::USD_PHASE_COMP_DC - 2 > dPhaseCompFrac = m_dPhaseCompFrac[ei][0]; @@ -730,11 +719,16 @@ class ElementBasedAssemblyKernel real64 dPhaseAmount_dC[numComp]{}; real64 dPhaseCompFrac_dC[numComp]{}; + // start with old time step values + for( integer ic = 0; ic < numComp; ++ic ) + { + stack.localResidual[ic] = -m_compAmount_n[ei][ic]; + } + // sum contributions to component accumulation from each phase for( integer ip = 0; ip < m_numPhases; ++ip ) { real64 const phaseAmount = stack.poreVolume * phaseVolFrac[ip] * phaseDens[ip]; - real64 const phaseAmount_n = stack.poreVolume_n * phaseVolFrac_n[ip] * phaseDens_n[ip]; real64 const dPhaseAmount_dP = stack.dPoreVolume_dPres * phaseVolFrac[ip] * phaseDens[ip] + stack.poreVolume * ( dPhaseVolFrac[ip][Deriv::dP] * phaseDens[ip] @@ -754,12 +748,11 @@ class ElementBasedAssemblyKernel for( integer ic = 0; ic < numComp; ++ic ) { real64 const phaseCompAmount = phaseAmount * phaseCompFrac[ip][ic]; - real64 const phaseCompAmount_n = phaseAmount_n * phaseCompFrac_n[ip][ic]; real64 const dPhaseCompAmount_dP = dPhaseAmount_dP * phaseCompFrac[ip][ic] + phaseAmount * dPhaseCompFrac[ip][ic][Deriv::dP]; - stack.localResidual[ic] += phaseCompAmount - phaseCompAmount_n; + stack.localResidual[ic] += phaseCompAmount; stack.localJacobian[ic][0] += dPhaseCompAmount_dP; // jc - index of component w.r.t. whose compositional var the derivative is being taken @@ -778,7 +771,7 @@ class ElementBasedAssemblyKernel // call the lambda in the phase loop to allow the reuse of the phase amounts and their derivatives // possible use: assemble the derivatives wrt temperature, and the accumulation term of the energy equation for this phase - phaseAmountKernelOp( ip, phaseAmount, phaseAmount_n, dPhaseAmount_dP, dPhaseAmount_dC ); + phaseAmountKernelOp( ip, phaseAmount, dPhaseAmount_dP, dPhaseAmount_dC ); } } @@ -911,7 +904,6 @@ class ElementBasedAssemblyKernel arrayView1d< real64 const > const m_volume; /// Views on the porosity - arrayView2d< real64 const > const m_porosity_n; arrayView2d< real64 const > const m_porosity; arrayView2d< real64 const > const m_dPoro_dPres; @@ -919,23 +911,22 @@ class ElementBasedAssemblyKernel arrayView3d< real64 const, compflow::USD_COMP_DC > const m_dCompFrac_dCompDens; /// Views on the phase volume fractions - arrayView2d< real64 const, compflow::USD_PHASE > const m_phaseVolFrac_n; arrayView2d< real64 const, compflow::USD_PHASE > const m_phaseVolFrac; arrayView3d< real64 const, compflow::USD_PHASE_DC > const m_dPhaseVolFrac; /// Views on the phase densities - arrayView3d< real64 const, multifluid::USD_PHASE > const m_phaseDens_n; arrayView3d< real64 const, multifluid::USD_PHASE > const m_phaseDens; arrayView4d< real64 const, multifluid::USD_PHASE_DC > const m_dPhaseDens; /// Views on the phase component fraction - arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const m_phaseCompFrac_n; arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const m_phaseCompFrac; arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > const m_dPhaseCompFrac; - // Views on component densities + // View on component densities arrayView2d< real64 const, compflow::USD_COMP > m_compDens; - arrayView2d< real64 const, compflow::USD_COMP > m_compDens_n; + + // View on component amount (mass/moles) from previous time step + arrayView2d< real64 const, compflow::USD_COMP > m_compAmount_n; /// View on the local CRS matrix CRSMatrixView< real64, globalIndex const > const m_localMatrix; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 1b344dc0130..4c2538c4d67 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -85,6 +85,9 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< mobility >( getName() ); subRegion.registerField< dMobility_dPressure >( getName() ); + subRegion.registerField< fields::flow::mass >( getName() ); + subRegion.registerField< fields::flow::mass_n >( getName() ); + if( m_isThermal ) { subRegion.registerField< dMobility_dTemperature >( getName() ); @@ -244,6 +247,60 @@ void SinglePhaseBase::updateFluidModel( ObjectManagerBase & dataGroup ) const } ); } +void SinglePhaseBase::updateMass( ElementSubRegionBase & subRegion ) const +{ + GEOS_MARK_FUNCTION; + + arrayView1d< real64 > const mass = subRegion.getField< fields::flow::mass >(); + arrayView1d< real64 > const mass_n = subRegion.getField< fields::flow::mass_n >(); + + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); + arrayView2d< real64 const > const porosity_n = porousSolid.getPorosity_n(); + + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView1d< real64 > const deltaVolume = subRegion.getField< fields::flow::deltaVolume >(); + + SingleFluidBase & fluid = + getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) ); + arrayView2d< real64 const > const density = fluid.density(); + arrayView2d< real64 const > const density_n = fluid.density_n(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + mass[ei] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * density[ei][0]; + if( isZero( mass_n[ei] ) ) // this is a hack for hydrofrac cases + mass_n[ei] = porosity_n[ei][0] * volume[ei] * density_n[ei][0]; // initialize newly created element mass + } ); +} + +void SinglePhaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const +{ + GEOS_MARK_FUNCTION; + + arrayView1d< real64 > const energy = subRegion.getField< fields::flow::energy >(); + + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); + arrayView2d< real64 const > const rockInternalEnergy = porousSolid.getInternalEnergy(); + + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView1d< real64 > const deltaVolume = subRegion.getField< fields::flow::deltaVolume >(); + + SingleFluidBase & fluid = + getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) ); + arrayView2d< real64 const > const density = fluid.density(); + arrayView2d< real64 const > const fluidInternalEnergy = fluid.internalEnergy(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + energy[ei] = ( volume[ei] + deltaVolume[ei] ) * + ( porosity[ei][0] * density[ei][0] * fluidInternalEnergy[ei][0] + ( 1.0 - porosity[ei][0] ) * rockInternalEnergy[ei][0] ); + } ); +} + void SinglePhaseBase::updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const { arrayView1d< real64 const > const temp = dataGroup.getField< fields::flow::temperature >(); @@ -271,9 +328,10 @@ void SinglePhaseBase::updateThermalConductivity( ElementSubRegionBase & subRegio conductivityMaterial.update( porosity ); } -void SinglePhaseBase::updateFluidState( ObjectManagerBase & subRegion ) const +void SinglePhaseBase::updateFluidState( ElementSubRegionBase & subRegion ) const { updateFluidModel( subRegion ); + updateMass( subRegion ); updateMobility( subRegion ); } @@ -388,7 +446,6 @@ void SinglePhaseBase::initializePostInitialConditionsPreSubGroups() SolidInternalEnergy const & solidInternalEnergyMaterial = getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); solidInternalEnergyMaterial.saveConvergedState(); - } } ); @@ -409,16 +466,21 @@ void SinglePhaseBase::initializePostInitialConditionsPreSubGroups() } ); } ); - // Save initial pressure field mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, ElementSubRegionBase & subRegion ) { + // Save initial pressure field arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >(); arrayView1d< real64 > const initPres = subRegion.getField< fields::flow::initialPressure >(); arrayView1d< real64 const > const & temp = subRegion.template getField< fields::flow::temperature >(); arrayView1d< real64 > const initTemp = subRegion.template getField< fields::flow::initialTemperature >(); initPres.setValues< parallelDevicePolicy<> >( pres ); initTemp.setValues< parallelDevicePolicy<> >( temp ); + + // finally update mass and energy + updateMass( subRegion ); + if( m_isThermal ) + updateEnergy( subRegion ); } ); } ); @@ -628,6 +690,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_ { updateSolidInternalEnergyModel( subRegion ); updateThermalConductivity( subRegion ); + updateEnergy( subRegion ); } } ); @@ -1193,6 +1256,9 @@ void SinglePhaseBase::updateState( DomainPartition & domain ) { GEOS_MARK_FUNCTION; + if( m_keepFlowVariablesConstantDuringInitStep ) + return; + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) @@ -1206,6 +1272,7 @@ void SinglePhaseBase::updateState( DomainPartition & domain ) if( m_isThermal ) { updateSolidInternalEnergyModel( subRegion ); + updateEnergy( subRegion ); } } ); } ); @@ -1237,6 +1304,7 @@ void SinglePhaseBase::resetStateToBeginningOfStep( DomainPartition & domain ) if( m_isThermal ) { updateSolidInternalEnergyModel( subRegion ); + updateEnergy( subRegion ); } } ); } ); @@ -1324,4 +1392,13 @@ bool SinglePhaseBase::checkSystemSolution( DomainPartition & domain, return (m_allowNegativePressure || numNegativePressures == 0) ? 1 : 0; } +void SinglePhaseBase::saveConvergedState( ElementSubRegionBase & subRegion ) const +{ + FlowSolverBase::saveConvergedState( subRegion ); + + arrayView1d< real64 const > const mass = subRegion.template getField< fields::flow::mass >(); + arrayView1d< real64 > const mass_n = subRegion.template getField< fields::flow::mass_n >(); + mass_n.setValues< parallelDevicePolicy<> >( mass ); +} + } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index 56dbab8979e..4a61817a2f4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -272,10 +272,10 @@ class SinglePhaseBase : public FlowSolverBase /** * @brief Function to update all constitutive state and dependent variables - * @param dataGroup group that contains the fields + * @param subRegion subregion that contains the fields */ void - updateFluidState( ObjectManagerBase & subRegion ) const; + updateFluidState( ElementSubRegionBase & subRegion ) const; /** @@ -285,6 +285,20 @@ class SinglePhaseBase : public FlowSolverBase virtual void updateFluidModel( ObjectManagerBase & dataGroup ) const; + /** + * @brief Function to update fluid mass + * @param subRegion subregion that contains the fields + */ + void + updateMass( ElementSubRegionBase & subRegion ) const; + + /** + * @brief Function to update energy + * @param subRegion subregion that contains the fields + */ + void + updateEnergy( ElementSubRegionBase & subRegion ) const; + /** * @brief Update all relevant solid internal energy models using current values of temperature * @param dataGroup the group storing the required fields @@ -354,6 +368,11 @@ class SinglePhaseBase : public FlowSolverBase virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; + /** + * @brief Utility function to save the converged state + * @param[in] subRegion the element subRegion + */ + virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override; /** * @brief Structure holding views into fluid properties used by the base solver. diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp index 8a101bea087..2d908b6dc10 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBaseKernels.hpp @@ -151,12 +151,11 @@ class ElementBasedAssemblyKernel m_elemGhostRank( subRegion.ghostRank() ), m_volume( subRegion.getElementVolume() ), m_deltaVolume( subRegion.template getField< fields::flow::deltaVolume >() ), - m_porosity_n( solid.getPorosity_n() ), - m_porosityNew( solid.getPorosity() ), + m_porosity( solid.getPorosity() ), m_dPoro_dPres( solid.getDporosity_dPressure() ), - m_density_n( fluid.density_n() ), m_density( fluid.density() ), m_dDensity_dPres( fluid.dDensity_dPressure() ), + m_mass_n( subRegion.template getField< fields::flow::mass_n >() ), m_localMatrix( localMatrix ), m_localRhs( localRhs ) {} @@ -174,9 +173,6 @@ class ElementBasedAssemblyKernel /// Pore volume at time n+1 real64 poreVolume = 0.0; - /// Pore volume at the previous converged time step - real64 poreVolume_n = 0.0; - /// Derivative of pore volume with respect to pressure real64 dPoreVolume_dPres = 0.0; @@ -216,8 +212,7 @@ class ElementBasedAssemblyKernel StackVariables & stack ) const { // initialize the pore volume - stack.poreVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * m_porosityNew[ei][0]; - stack.poreVolume_n = m_volume[ei] * m_porosity_n[ei][0]; + stack.poreVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * m_porosity[ei][0]; stack.dPoreVolume_dPres = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dPres[ei][0]; // set row index and degrees of freedom indices for this element @@ -242,7 +237,7 @@ class ElementBasedAssemblyKernel FUNC && kernelOp = NoOpFunc{} ) const { // Residual contribution is mass conservation in the cell - stack.localResidual[0] = stack.poreVolume * m_density[ei][0] - stack.poreVolume_n * m_density_n[ei][0]; + 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; @@ -314,15 +309,16 @@ class ElementBasedAssemblyKernel arrayView1d< real64 const > const m_deltaVolume; /// Views on the porosity - arrayView2d< real64 const > const m_porosity_n; - arrayView2d< real64 const > const m_porosityNew; + arrayView2d< real64 const > const m_porosity; arrayView2d< real64 const > const m_dPoro_dPres; /// Views on density - arrayView2d< real64 const > const m_density_n; arrayView2d< real64 const > const m_density; arrayView2d< real64 const > const m_dDensity_dPres; + /// View on mass + arrayView1d< real64 const > const m_mass_n; + /// View on the local CRS matrix CRSMatrixView< real64, globalIndex const > const m_localMatrix; /// View on the local RHS @@ -374,7 +370,7 @@ class SurfaceElementBasedAssemblyKernel : public ElementBasedAssemblyKernel< Sur { Base::computeAccumulation( ei, stack, [&] () { - if( Base::m_volume[ei] * Base::m_density_n[ei][0] > 1.1 * m_creationMass[ei] ) + if( Base::m_mass_n[ei] > 1.1 * m_creationMass[ei] ) { stack.localResidual[0] += m_creationMass[ei] * 0.25; } @@ -470,24 +466,20 @@ class ResidualNormKernel : public solverBaseKernels::ResidualNormKernelBase< 1 > arrayView1d< globalIndex const > const & dofNumber, arrayView1d< localIndex const > const & ghostRank, ElementSubRegionBase const & subRegion, - constitutive::SingleFluidBase const & fluid, - constitutive::CoupledSolidBase const & solid, real64 const minNormalizer ) : Base( rankOffset, localResidual, dofNumber, ghostRank, minNormalizer ), - m_volume( subRegion.getElementVolume() ), - m_porosity_n( solid.getPorosity_n() ), - m_density_n( fluid.density_n() ) + m_mass_n( subRegion.template getField< fields::flow::mass_n >() ) {} GEOS_HOST_DEVICE virtual void computeLinf( localIndex const ei, LinfStackVariables & stack ) const override { - real64 const massNormalizer = LvArray::math::max( m_minNormalizer, m_density_n[ei][0] * m_porosity_n[ei][0] * m_volume[ei] ); + real64 const massNormalizer = LvArray::math::max( m_minNormalizer, m_mass_n[ei] ); real64 const valMass = LvArray::math::abs( m_localResidual[stack.localRow] ) / massNormalizer; if( valMass > stack.localValue[0] ) { @@ -499,7 +491,7 @@ class ResidualNormKernel : public solverBaseKernels::ResidualNormKernelBase< 1 > virtual void computeL2( localIndex const ei, L2StackVariables & stack ) const override { - real64 const massNormalizer = LvArray::math::max( m_minNormalizer, m_density_n[ei][0] * m_porosity_n[ei][0] * m_volume[ei] ); + real64 const massNormalizer = LvArray::math::max( m_minNormalizer, m_mass_n[ei] ); stack.localValue[0] += m_localResidual[stack.localRow] * m_localResidual[stack.localRow]; stack.localNormalizer[0] += massNormalizer; } @@ -507,14 +499,8 @@ class ResidualNormKernel : public solverBaseKernels::ResidualNormKernelBase< 1 > protected: - /// View on the volume - arrayView1d< real64 const > const m_volume; - - /// View on porosity at the previous converged time step - arrayView2d< real64 const > const m_porosity_n; - - /// View on total mass/molar density at the previous converged time step - arrayView2d< real64 const > const m_density_n; + /// View on mass at the previous converged time step + arrayView1d< real64 const > const m_mass_n; }; @@ -545,8 +531,6 @@ class ResidualNormKernelFactory string const dofKey, arrayView1d< real64 const > const & localResidual, ElementSubRegionBase const & subRegion, - constitutive::SingleFluidBase const & fluid, - constitutive::CoupledSolidBase const & solid, real64 const minNormalizer, real64 (& residualNorm)[1], real64 (& residualNormalizer)[1] ) @@ -554,7 +538,7 @@ class ResidualNormKernelFactory arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); - ResidualNormKernel kernel( rankOffset, localResidual, dofNumber, ghostRank, subRegion, fluid, solid, minNormalizer ); + ResidualNormKernel kernel( rankOffset, localResidual, dofNumber, ghostRank, subRegion, minNormalizer ); if( normType == solverBaseKernels::NormType::Linf ) { ResidualNormKernel::launchLinf< POLICY >( subRegion.size(), kernel, residualNorm ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index e207d50ad99..565a4621fbe 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -150,13 +150,13 @@ real64 SinglePhaseFVM< BASE >::calculateResidualNorm( real64 const & GEOS_UNUSED string const & fluidName = subRegion.template getReference< string >( BASE::viewKeyStruct::fluidNamesString() ); SingleFluidBase const & fluid = SolverBase::getConstitutiveModel< SingleFluidBase >( subRegion, fluidName ); - string const & solidName = subRegion.template getReference< string >( BASE::viewKeyStruct::solidNamesString() ); - CoupledSolidBase const & solid = SolverBase::getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - // step 1: compute the norm in the subRegion if( m_isThermal ) { + string const & solidName = subRegion.template getReference< string >( BASE::viewKeyStruct::solidNamesString() ); + CoupledSolidBase const & solid = SolverBase::getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); + string const & solidInternalEnergyName = subRegion.template getReference< string >( BASE::viewKeyStruct::solidInternalEnergyNamesString() ); SolidInternalEnergy const & solidInternalEnergy = SolverBase::getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); @@ -185,8 +185,6 @@ real64 SinglePhaseFVM< BASE >::calculateResidualNorm( real64 const & GEOS_UNUSED dofKey, localRhs, subRegion, - fluid, - solid, m_nonlinearSolverParameters.m_minNormalizer, subRegionFlowResidualNorm, subRegionFlowResidualNormalizer ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp index a21e8cc461e..1ce44237606 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp @@ -471,9 +471,6 @@ real64 SinglePhaseHybridFVM::calculateResidualNorm( real64 const & GEOS_UNUSED_P defaultViscosity += fluid.defaultViscosity(); subRegionCounter++; - string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() ); - CoupledSolidBase const & solid = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - // step 1.1: compute the norm in the subRegion singlePhaseBaseKernels:: @@ -483,8 +480,6 @@ real64 SinglePhaseHybridFVM::calculateResidualNorm( real64 const & GEOS_UNUSED_P elemDofKey, localRhs, subRegion, - fluid, - solid, m_nonlinearSolverParameters.m_minNormalizer, subRegionResidualNorm, subRegionResidualNormalizer ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp index 0fec20002a4..8df182ad30e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp @@ -159,17 +159,13 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK using Base::m_dofNumber; using Base::m_elemGhostRank; using Base::m_volume; - using Base::m_porosity_n; using Base::m_porosity; using Base::m_dPoro_dPres; using Base::m_dCompFrac_dCompDens; - using Base::m_phaseVolFrac_n; using Base::m_phaseVolFrac; using Base::m_dPhaseVolFrac; - using Base::m_phaseDens_n; using Base::m_phaseDens; using Base::m_dPhaseDens; - using Base::m_phaseCompFrac_n; using Base::m_phaseCompFrac; using Base::m_dPhaseCompFrac; using Base::m_localMatrix; @@ -197,12 +193,11 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK BitFlags< isothermalCompositionalMultiphaseBaseKernels::ElementBasedAssemblyKernelFlags > const kernelFlags ) : Base( numPhases, rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs, kernelFlags ), m_dPoro_dTemp( solid.getDporosity_dTemperature() ), - m_phaseInternalEnergy_n( fluid.phaseInternalEnergy_n() ), m_phaseInternalEnergy( fluid.phaseInternalEnergy() ), m_dPhaseInternalEnergy( fluid.dPhaseInternalEnergy() ), - m_rockInternalEnergy_n( solid.getInternalEnergy_n() ), m_rockInternalEnergy( solid.getInternalEnergy() ), - m_dRockInternalEnergy_dTemp( solid.getDinternalEnergy_dTemperature() ) + m_dRockInternalEnergy_dTemp( solid.getDinternalEnergy_dTemperature() ), + m_energy_n( subRegion.getField< fields::flow::energy_n >() ) {} struct StackVariables : public Base::StackVariables @@ -214,9 +209,6 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK : Base::StackVariables() {} - using Base::StackVariables::poreVolume; - using Base::StackVariables::poreVolume_n; - using Base::StackVariables::dPoreVolume_dPres; using Base::StackVariables::localRow; using Base::StackVariables::dofIndices; using Base::StackVariables::localResidual; @@ -230,9 +222,6 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK /// Solid energy at time n+1 real64 solidEnergy = 0.0; - /// Solid energy at the previous converged time step - real64 solidEnergy_n = 0.0; - /// Derivative of solid internal energy with respect to pressure real64 dSolidEnergy_dPres = 0.0; @@ -257,13 +246,11 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK // initialize the solid volume real64 const solidVolume = m_volume[ei] * ( 1.0 - m_porosity[ei][0] ); - real64 const solidVolume_n = m_volume[ei] * ( 1.0 - m_porosity_n[ei][0] ); real64 const dSolidVolume_dPres = -m_volume[ei] * m_dPoro_dPres[ei][0]; real64 const dSolidVolume_dTemp = -stack.dPoreVolume_dTemp; // initialize the solid internal energy stack.solidEnergy = solidVolume * m_rockInternalEnergy[ei][0]; - stack.solidEnergy_n = solidVolume_n * m_rockInternalEnergy_n[ei][0]; stack.dSolidEnergy_dPres = dSolidVolume_dPres * m_rockInternalEnergy[ei][0]; stack.dSolidEnergy_dTemp = solidVolume * m_dRockInternalEnergy_dTemp[ei][0] + dSolidVolume_dTemp * m_rockInternalEnergy[ei][0]; @@ -281,9 +268,11 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK { using Deriv = multifluid::DerivativeOffset; + // start with old time step value + stack.localResidual[numEqn-1] = -m_energy_n[ei]; + Base::computeAccumulation( ei, stack, [&] ( integer const ip, real64 const & phaseAmount, - real64 const & phaseAmount_n, real64 const & dPhaseAmount_dP, real64 const (&dPhaseAmount_dC)[numComp] ) { @@ -302,7 +291,6 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseDens = m_dPhaseDens[ei][0]; arraySlice2d< real64 const, multifluid::USD_PHASE_COMP - 2 > phaseCompFrac = m_phaseCompFrac[ei][0]; arraySlice3d< real64 const, multifluid::USD_PHASE_COMP_DC - 2 > dPhaseCompFrac = m_dPhaseCompFrac[ei][0]; - arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseInternalEnergy_n = m_phaseInternalEnergy_n[ei][0]; arraySlice1d< real64 const, multifluid::USD_PHASE - 2 > phaseInternalEnergy = m_phaseInternalEnergy[ei][0]; arraySlice2d< real64 const, multifluid::USD_PHASE_DC - 2 > dPhaseInternalEnergy = m_dPhaseInternalEnergy[ei][0]; @@ -319,14 +307,13 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK // Step 2: assemble the phase-dependent part of the accumulation term of the energy equation real64 const phaseEnergy = phaseAmount * phaseInternalEnergy[ip]; - real64 const phaseEnergy_n = phaseAmount_n * phaseInternalEnergy_n[ip]; real64 const dPhaseEnergy_dP = dPhaseAmount_dP * phaseInternalEnergy[ip] + phaseAmount * dPhaseInternalEnergy[ip][Deriv::dP]; real64 const dPhaseEnergy_dT = dPhaseAmount_dT * phaseInternalEnergy[ip] + phaseAmount * dPhaseInternalEnergy[ip][Deriv::dT]; // local accumulation - stack.localResidual[numEqn-1] += phaseEnergy - phaseEnergy_n; + stack.localResidual[numEqn-1] += phaseEnergy; // derivatives w.r.t. pressure and temperature stack.localJacobian[numEqn-1][0] += dPhaseEnergy_dP; @@ -344,7 +331,7 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK // Step 3: assemble the solid part of the accumulation term // local accumulation and derivatives w.r.t. pressure and temperature - stack.localResidual[numEqn-1] += stack.solidEnergy - stack.solidEnergy_n; + stack.localResidual[numEqn-1] += stack.solidEnergy; stack.localJacobian[numEqn-1][0] += stack.dSolidEnergy_dPres; stack.localJacobian[numEqn-1][numDof-1] += stack.dSolidEnergy_dTemp; @@ -396,15 +383,16 @@ class ElementBasedAssemblyKernel : public isothermalCompositionalMultiphaseBaseK arrayView2d< real64 const > const m_dPoro_dTemp; /// Views on phase internal energy - arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseInternalEnergy_n; arrayView3d< real64 const, multifluid::USD_PHASE > m_phaseInternalEnergy; arrayView4d< real64 const, multifluid::USD_PHASE_DC > m_dPhaseInternalEnergy; /// Views on rock internal energy - arrayView2d< real64 const > m_rockInternalEnergy_n; arrayView2d< real64 const > m_rockInternalEnergy; arrayView2d< real64 const > m_dRockInternalEnergy_dTemp; + /// Views on energy + arrayView1d< real64 const > m_energy_n; + }; /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ThermalSinglePhaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ThermalSinglePhaseBaseKernels.hpp index 12047977050..bc006a230c6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ThermalSinglePhaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ThermalSinglePhaseBaseKernels.hpp @@ -120,10 +120,8 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs using Base::m_elemGhostRank; using Base::m_volume; using Base::m_deltaVolume; - using Base::m_porosity_n; - using Base::m_porosityNew; + using Base::m_porosity; using Base::m_dPoro_dPres; - using Base::m_density_n; using Base::m_density; using Base::m_dDensity_dPres; using Base::m_localMatrix; @@ -149,13 +147,12 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs : Base( rankOffset, dofKey, subRegion, fluid, solid, localMatrix, localRhs ), m_dDensity_dTemp( fluid.dDensity_dTemperature() ), m_dPoro_dTemp( solid.getDporosity_dTemperature() ), - m_internalEnergy_n( fluid.internalEnergy_n() ), m_internalEnergy( fluid.internalEnergy() ), m_dInternalEnergy_dPres( fluid.dInternalEnergy_dPressure() ), m_dInternalEnergy_dTemp( fluid.dInternalEnergy_dTemperature() ), - m_rockInternalEnergy_n( solid.getInternalEnergy_n() ), m_rockInternalEnergy( solid.getInternalEnergy() ), - m_dRockInternalEnergy_dTemp( solid.getDinternalEnergy_dTemperature() ) + m_dRockInternalEnergy_dTemp( solid.getDinternalEnergy_dTemperature() ), + m_energy_n( subRegion.template getField< fields::flow::energy_n >() ) {} /** @@ -172,7 +169,6 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs {} using Base::StackVariables::poreVolume; - using Base::StackVariables::poreVolume_n; using Base::StackVariables::dPoreVolume_dPres; using Base::StackVariables::localRow; using Base::StackVariables::dofIndices; @@ -187,9 +183,6 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs /// Solid energy at time n+1 real64 solidEnergy = 0.0; - /// Solid energy at the previous converged time step - real64 solidEnergy_n = 0.0; - /// Derivative of solid internal energy with respect to pressure real64 dSolidEnergy_dPres = 0.0; @@ -212,14 +205,12 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs stack.dPoreVolume_dTemp = ( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dTemp[ei][0]; // initialize the solid volume - real64 const solidVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * ( 1.0 - m_porosityNew[ei][0] ); - real64 const solidVolume_n = m_volume[ei] * ( 1.0 - m_porosity_n[ei][0] ); + real64 const solidVolume = ( m_volume[ei] + m_deltaVolume[ei] ) * ( 1.0 - m_porosity[ei][0] ); real64 const dSolidVolume_dPres = -( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dPres[ei][0]; real64 const dSolidVolume_dTemp = -( m_volume[ei] + m_deltaVolume[ei] ) * m_dPoro_dTemp[ei][0]; // initialize the solid internal energy stack.solidEnergy = solidVolume * m_rockInternalEnergy[ei][0]; - stack.solidEnergy_n = solidVolume_n * m_rockInternalEnergy_n[ei][0]; stack.dSolidEnergy_dPres = dSolidVolume_dPres * m_rockInternalEnergy[ei][0]; stack.dSolidEnergy_dTemp = solidVolume * m_dRockInternalEnergy_dTemp[ei][0] + dSolidVolume_dTemp * m_rockInternalEnergy[ei][0]; } @@ -235,6 +226,8 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs void computeAccumulation( localIndex const ei, StackVariables & stack ) const { + stack.localResidual[numEqn-1] = -m_energy_n[ei]; + Base::computeAccumulation( ei, stack, [&] () { // Step 1: assemble the derivatives of the mass balance equation w.r.t temperature @@ -242,7 +235,6 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs // 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]; - real64 const fluidEnergy_n = stack.poreVolume_n * m_density_n[ei][0] * m_internalEnergy_n[ei][0]; real64 const dFluidEnergy_dP = stack.dPoreVolume_dPres * m_density[ei][0] * m_internalEnergy[ei][0] + stack.poreVolume * m_dDensity_dPres[ei][0] * m_internalEnergy[ei][0] @@ -253,7 +245,7 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs + stack.dPoreVolume_dTemp * m_density[ei][0] * m_internalEnergy[ei][0]; // local accumulation - stack.localResidual[numEqn-1] = fluidEnergy - fluidEnergy_n; + stack.localResidual[numEqn-1] += fluidEnergy; // derivatives w.r.t. pressure and temperature stack.localJacobian[numEqn-1][0] = dFluidEnergy_dP; @@ -261,7 +253,7 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs } ); // Step 3: assemble the solid part of the accumulation term of the energy equation - stack.localResidual[numEqn-1] += stack.solidEnergy - stack.solidEnergy_n; + stack.localResidual[numEqn-1] += stack.solidEnergy; stack.localJacobian[numEqn-1][0] += stack.dSolidEnergy_dPres; stack.localJacobian[numEqn-1][numDof-1] += stack.dSolidEnergy_dTemp; } @@ -297,15 +289,16 @@ class ElementBasedAssemblyKernel : public singlePhaseBaseKernels::ElementBasedAs arrayView2d< real64 const > const m_dPoro_dTemp; /// Views on fluid internal energy - arrayView2d< real64 const > const m_internalEnergy_n; arrayView2d< real64 const > const m_internalEnergy; arrayView2d< real64 const > const m_dInternalEnergy_dPres; arrayView2d< real64 const > const m_dInternalEnergy_dTemp; /// Views on rock internal energy - arrayView2d< real64 const > m_rockInternalEnergy_n; - arrayView2d< real64 const > m_rockInternalEnergy; - arrayView2d< real64 const > m_dRockInternalEnergy_dTemp; + arrayView2d< real64 const > const m_rockInternalEnergy; + arrayView2d< real64 const > const m_dRockInternalEnergy_dTemp; + + /// View on energy + arrayView1d< real64 const > const m_energy_n; }; @@ -352,7 +345,7 @@ class SurfaceElementBasedAssemblyKernel : public ElementBasedAssemblyKernel< Sur Base::StackVariables & stack ) const { Base::computeAccumulation( ei, stack ); - if( Base::m_volume[ei] * Base::m_density_n[ei][0] > 1.1 * m_creationMass[ei] ) + if( Base::m_mass_n[ei] > 1.1 * m_creationMass[ei] ) { stack.localResidual[0] += m_creationMass[ei] * 0.25; } diff --git a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp index 6dafaf94948..e884fb94cfd 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp @@ -90,7 +90,7 @@ class CompositionalMultiphaseReservoirAndWells : public CoupledReservoirAndWells void keepFlowVariablesConstantDuringInitStep( bool const keepFlowVariablesConstantDuringInitStep ) { flowSolver()->keepFlowVariablesConstantDuringInitStep( keepFlowVariablesConstantDuringInitStep ); } - real64 updateFluidState( ObjectManagerBase & subRegion ) const + real64 updateFluidState( ElementSubRegionBase & subRegion ) const { return flowSolver()->updateFluidState( subRegion ); } void updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const { flowSolver()->updatePorosityAndPermeability( subRegion ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp index 36eaa63652f..d51345298cd 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp @@ -539,10 +539,10 @@ class CoupledSolver : public SolverBase GEOS_UNUSED_VAR( domain, solverType ); } - bool checkSequentialConvergence( int const & iter, - real64 const & time_n, - real64 const & dt, - DomainPartition & domain ) + virtual bool checkSequentialConvergence( int const & iter, + real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) { NonlinearSolverParameters const & params = getNonlinearSolverParameters(); bool isConverged = true; diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index b5c05483d87..e0b93cb71bd 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -186,6 +186,25 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER this->setupCoupling( domain, dofManager ); } + virtual bool checkSequentialConvergence( int const & iter, + real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) override + { + // always force outer loop for initialization + auto & subcycling = this->getNonlinearSolverParameters().m_subcyclingOption; + auto const subcycling_orig = subcycling; + if( m_performStressInitialization ) + subcycling = 1; + + bool isConverged = Base::checkSequentialConvergence( iter, time_n, dt, domain ); + + // restore original + subcycling = subcycling_orig; + + return isConverged; + } + /** * @brief accessor for the pointer to the solid mechanics solver * @return a pointer to the solid mechanics solver diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp index 04891c91dd3..2d1c5057c85 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp @@ -82,7 +82,7 @@ class SinglePhaseReservoirAndWells : public CoupledReservoirAndWellsBase< SINGLE void keepFlowVariablesConstantDuringInitStep( bool const keepFlowVariablesConstantDuringInitStep ) { flowSolver()->keepFlowVariablesConstantDuringInitStep( keepFlowVariablesConstantDuringInitStep ); } - void updateFluidState( ObjectManagerBase & subRegion ) const + void updateFluidState( ElementSubRegionBase & subRegion ) const { flowSolver()->updateFluidState( subRegion ); } void updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const { flowSolver()->updatePorosityAndPermeability( subRegion ); }