From 470921ff28312a4bf0039bc2f9acb50941541457 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 3 Dec 2024 09:36:54 -0600 Subject: [PATCH] fix: detangle and unify flow initialization, fix netToGross bug (#3393) --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 6 +- .../fluidFlow/CompositionalMultiphaseBase.cpp | 247 +++++++----------- .../fluidFlow/CompositionalMultiphaseBase.hpp | 8 +- .../fluidFlow/FlowSolverBase.cpp | 107 +++++++- .../fluidFlow/FlowSolverBase.hpp | 73 +++--- .../fluidFlow/SinglePhaseBase.cpp | 139 ++++------ .../fluidFlow/SinglePhaseBase.hpp | 14 +- .../fluidFlow/SinglePhaseHybridFVM.cpp | 17 -- .../fluidFlow/SinglePhaseHybridFVM.hpp | 2 - 10 files changed, 296 insertions(+), 319 deletions(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 0ebe2652f98..92b83368c51 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3381-9063-399123c + baseline: integratedTests/baseline_integratedTests-pr3393-9089-c187f5d allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 1fbd4cbfa9d..46193312093 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,11 +6,14 @@ 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 #3393 (2024-12-2) +===================== +Fix netToGross bug. + PR #3381 (2024-12-01) ===================== A few baseline diffs for order FaceElementSubRegion::m_toFacesRelation map. Not sure why this was changed by this PR, but the previous order seems incorrect for a couple of cases. - PR #2957 (2024-11-27) ===================== Added ExternalDataRepository. @@ -19,7 +22,6 @@ PR #3448 (2024-11-21) ===================== Switched the FaceElementSubRegion::m_toFacesRelation and FaceElementSubRegion::m_2dElemToElems back to array2d instead of ArrayOfArray. This results in a reordering m_toFacesRelation back to the "correct" assumed order of "original face first". This fixes a bug that failed to remove the CellStencil entry when a FaceElement splits two cells. - PR #2637 (2024-11-21) ===================== Added numberOfTargetProcesses. diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 08e7176fdbe..5e59ba4bc2c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -831,30 +831,25 @@ real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & sub } void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, - DomainPartition & domain, arrayView1d< string const > const & regionNames ) { GEOS_MARK_FUNCTION; integer const numComp = m_numComponents; - // 1. Compute hydrostatic equilibrium in the regions for which corresponding field specification tag has been specified - computeHydrostaticEquilibrium(); - mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, ElementSubRegionBase & subRegion ) { - // 2. Assume global component fractions have been prescribed. + // Assume global component fractions have been prescribed. // Initialize constitutive state to get fluid density. updateFluidModel( subRegion ); - // 3. Back-calculate global component densities from fractions and total fluid density + // Back-calculate global component densities from fractions and total fluid density // in order to initialize the primary solution variables - string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); + string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); arrayView2d< real64 const, multifluid::USD_FLUID > const totalDens = fluid.totalDensity(); - arrayView2d< real64 const, compflow::USD_COMP > const compFrac = subRegion.getField< fields::flow::globalCompFraction >(); arrayView2d< real64, compflow::USD_COMP > const compDens = @@ -867,10 +862,13 @@ 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 - chopNegativeDensities( domain ); + // with initial component densities defined - check if they need to be corrected to avoid zero diags etc + if( m_allowCompDensChopping ) + { + chopNegativeDensities( subRegion ); + } + } ); // for some reason CUDA does not want the host_device lambda to be defined inside the generic lambda // I need the exact type of the subRegion for updateSolidflowProperties to work well. @@ -878,115 +876,65 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, auto & subRegion ) { - // 4. Initialize/update dependent state quantities + // Initialize/update dependent state quantities - // 4.1 Update the constitutive models that only depend on - // - the primary variables - // - the fluid constitutive quantities (as they have already been updated) - // We postpone the other constitutive models for now - // In addition, to avoid multiplying permeability/porosity bay netToGross in the assembly kernel, we do it once and for all here - arrayView1d< real64 const > const netToGross = subRegion.template getField< fields::flow::netToGross >(); - CoupledSolidBase const & porousSolid = - getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); - PermeabilityBase const & permeabilityModel = - getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) ); - permeabilityModel.scaleHorizontalPermeability( netToGross ); - 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 + // Update the constitutive models that only depend on + // - the primary variables + // - the fluid constitutive quantities (as they have already been updated) + // We postpone the other constitutive models for now - // 4.2 Save the computed porosity into the old porosity - // - // Note: - // - This must be called after updatePorosityAndPermeability - // - This step depends on porosity - string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ); - CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); - porousMaterial.initializeState(); - - // 4.3 Initialize/update the relative permeability model using the initial phase volume fraction - // This is needed to handle relative permeability hysteresis - // Also, initialize the fluid model - // - // Note: - // - This must be called after updatePhaseVolumeFraction - // - This step depends on phaseVolFraction + // Now, we initialize and update each constitutive model one by one // initialized phase volume fraction arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = subRegion.template getField< fields::flow::phaseVolumeFraction >(); - string const & relpermName = subRegion.template getReference< string >( viewKeyStruct::relPermNamesString() ); - RelativePermeabilityBase & relPermMaterial = - getConstitutiveModel< RelativePermeabilityBase >( subRegion, relpermName ); - relPermMaterial.saveConvergedPhaseVolFractionState( phaseVolFrac ); // this needs to happen before calling updateRelPermModel + // Initialize/update the relative permeability model using the initial phase volume fraction + // Note: + // - This must be called after updatePhaseVolumeFraction + // - This step depends on phaseVolFraction + RelativePermeabilityBase & relPerm = + getConstitutiveModel< RelativePermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::relPermNamesString() ) ); + relPerm.saveConvergedPhaseVolFractionState( phaseVolFrac ); // this needs to happen before calling updateRelPermModel updateRelPermModel( subRegion ); - relPermMaterial.saveConvergedState(); // this needs to happen after calling updateRelPermModel + relPerm.saveConvergedState(); // this needs to happen after calling updateRelPermModel string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); - MultiFluidBase & fluidMaterial = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); - fluidMaterial.initializeState(); + MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); + fluid.initializeState(); + + // Update the phase mobility + // Note: + // - This must be called after updateRelPermModel + // - This step depends phaseRelPerm + updatePhaseMobility( subRegion ); - // 4.4 Then, we initialize/update the capillary pressure model - // + // Initialize/update the capillary pressure model // Note: - // - This must be called after updatePorosityAndPermeability - // - This step depends on porosity and permeability + // - This must be called after updatePorosityAndPermeability and updatePhaseVolumeFraction + // - This step depends on porosity, permeability, and phaseVolFraction if( m_hasCapPressure ) { // initialized porosity - arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); - string const & permName = subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ); - PermeabilityBase const & permeabilityMaterial = - getConstitutiveModel< PermeabilityBase >( subRegion, permName ); // initialized permeability + PermeabilityBase const & permeabilityMaterial = + getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) ); arrayView3d< real64 const > const permeability = permeabilityMaterial.permeability(); - string const & capPressureName = subRegion.template getReference< string >( viewKeyStruct::capPressureNamesString() ); - CapillaryPressureBase const & capPressureMaterial = - getConstitutiveModel< CapillaryPressureBase >( subRegion, capPressureName ); - capPressureMaterial.initializeRockState( porosity, permeability ); // this needs to happen before calling updateCapPressureModel + CapillaryPressureBase const & capPressure = + getConstitutiveModel< CapillaryPressureBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::capPressureNamesString() ) ); + capPressure.initializeRockState( porosity, permeability ); // this needs to happen before calling updateCapPressureModel updateCapPressureModel( subRegion ); } - // 4.5 Update the phase mobility - // - // Note: - // - This must be called after updateRelPermModel - // - This step depends phaseRelPerm - updatePhaseMobility( subRegion ); - - // 4.6 We initialize the rock thermal quantities: conductivity and solid internal energy - // - // Note: - // - This must be called after updatePorosityAndPermeability and updatePhaseVolumeFraction - // - This step depends on porosity and phaseVolFraction - if( m_isThermal ) - { - // initialized porosity - arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); - - string const & thermalConductivityName = subRegion.template getReference< string >( viewKeyStruct::thermalConductivityNamesString() ); - MultiPhaseThermalConductivityBase const & conductivityMaterial = - getConstitutiveModel< MultiPhaseThermalConductivityBase >( subRegion, thermalConductivityName ); - conductivityMaterial.initializeRockFluidState( porosity, phaseVolFrac ); - // note that there is nothing to update here because thermal conductivity is explicit for now - - updateSolidInternalEnergyModel( subRegion ); - string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString() ); - 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 + // If the diffusion and/or dispersion is/are supported, initialize the two models if( m_hasDiffusion ) { string const & diffusionName = subRegion.template getReference< string >( viewKeyStruct::diffusionNamesString() ); @@ -1004,24 +952,42 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, } } ); +} - // 5. Save initial pressure - mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) +void CompositionalMultiphaseBase::initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, + SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) { - 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.getField< fields::flow::temperature >(); - arrayView1d< real64 > const initTemp = subRegion.template getField< fields::flow::initialTemperature >(); - initPres.setValues< parallelDevicePolicy<> >( pres ); - initTemp.setValues< parallelDevicePolicy<> >( temp ); + // initialized porosity + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); + + // initialized phase volume fraction + arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = + subRegion.template getField< fields::flow::phaseVolumeFraction >(); + + string const & thermalConductivityName = subRegion.template getReference< string >( viewKeyStruct::thermalConductivityNamesString()); + MultiPhaseThermalConductivityBase const & conductivityMaterial = + getConstitutiveModel< MultiPhaseThermalConductivityBase >( subRegion, thermalConductivityName ); + conductivityMaterial.initializeRockFluidState( porosity, phaseVolFrac ); + // note that there is nothing to update here because thermal conductivity is explicit for now + + updateSolidInternalEnergyModel( subRegion ); + string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString()); + SolidInternalEnergy const & solidInternalEnergyMaterial = + getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); + solidInternalEnergyMaterial.saveConvergedState(); + + updateEnergy( subRegion ); } ); } -void CompositionalMultiphaseBase::computeHydrostaticEquilibrium() +void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) { FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); - DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); integer const numComps = m_numComponents; integer const numPhases = m_numPhases; @@ -1281,8 +1247,7 @@ void CompositionalMultiphaseBase::initializePostInitialConditionsPreSubGroups() arrayView1d< string const > const & regionNames ) { FieldIdentifiers fieldsToBeSync; - fieldsToBeSync.addElementFields( { fields::flow::pressure::key(), - fields::flow::globalCompDensity::key() }, + fieldsToBeSync.addElementFields( { fields::flow::globalCompDensity::key() }, regionNames ); CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); @@ -1295,35 +1260,10 @@ void CompositionalMultiphaseBase::initializePostInitialConditionsPreSubGroups() string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); fluid.setMassFlag( m_useMass ); - - saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes - updatePorosityAndPermeability( subRegion ); - - CoupledSolidBase const & porousSolid = - getConstitutiveModel< CoupledSolidBase >( subRegion, - subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); - porousSolid.initializeState(); } ); - - // Initialize primary variables from applied initial conditions - initializeFluidState( mesh, domain, regionNames ); - - mesh.getElemManager().forElementRegions< SurfaceElementRegion >( regionNames, - [&]( localIndex const, - SurfaceElementRegion & region ) - { - region.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) - { - subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key() ). - setApplyDefaultValue( region.getDefaultAperture() ); - } ); - } ); - } ); - // report to the user if some pore volumes are very small - // note: this function is here because: 1) porosity has been initialized and 2) NTG has been applied - validatePoreVolumes( domain ); + initialize( domain ); } void @@ -2071,9 +2011,6 @@ void CompositionalMultiphaseBase::chopNegativeDensities( DomainPartition & domai using namespace isothermalCompositionalMultiphaseBaseKernels; - integer const numComp = m_numComponents; - real64 const minCompDens = m_minCompDens; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) @@ -2082,25 +2019,33 @@ void CompositionalMultiphaseBase::chopNegativeDensities( DomainPartition & domai [&]( localIndex const, ElementSubRegionBase & subRegion ) { - arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + chopNegativeDensities( subRegion ); + } ); + } ); +} - arrayView2d< real64, compflow::USD_COMP > const compDens = - subRegion.getField< fields::flow::globalCompDensity >(); +void CompositionalMultiphaseBase::chopNegativeDensities( ElementSubRegionBase & subRegion ) +{ + integer const numComp = m_numComponents; + real64 const minCompDens = m_minCompDens; - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + arrayView2d< real64, compflow::USD_COMP > const compDens = + subRegion.getField< fields::flow::globalCompDensity >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + if( ghostRank[ei] < 0 ) + { + for( integer ic = 0; ic < numComp; ++ic ) { - if( ghostRank[ei] < 0 ) + if( compDens[ei][ic] < minCompDens ) { - for( integer ic = 0; ic < numComp; ++ic ) - { - if( compDens[ei][ic] < minCompDens ) - { - compDens[ei][ic] = minCompDens; - } - } + compDens[ei][ic] = minCompDens; } - } ); - } ); + } + } } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index cfd06428707..278d8d6b453 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -282,12 +282,14 @@ class CompositionalMultiphaseBase : public FlowSolverBase * from prescribed intermediate values (i.e. global densities from global fractions) * and any applicable hydrostatic equilibration of the domain */ - void initializeFluidState( MeshLevel & mesh, DomainPartition & domain, arrayView1d< string const > const & regionNames ); + virtual void initializeFluidState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + + virtual void initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ - void computeHydrostaticEquilibrium(); + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; /** * @brief Function to perform the Application of Dirichlet type BC's @@ -362,6 +364,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ void chopNegativeDensities( DomainPartition & domain ); + void chopNegativeDensities( ElementSubRegionBase & subRegion ); + virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, DomainPartition & domain ) override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index f803e7c7e1e..bfbbbf368cd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -286,16 +286,6 @@ void FlowSolverBase::saveSequentialIterationState( DomainPartition & domain ) m_sequentialTempChange = m_isThermal ? MpiWrapper::max( maxTempChange ) : 0.0; } -void FlowSolverBase::enableFixedStressPoromechanicsUpdate() -{ - m_isFixedStressPoromechanicsUpdate = true; -} - -void FlowSolverBase::enableJumpStabilization() -{ - m_isJumpStabilized = true; -} - void FlowSolverBase::setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const { PhysicsSolverBase::setConstitutiveNamesCallSuper( subRegion ); @@ -455,6 +445,12 @@ void FlowSolverBase::initializePostInitialConditionsPreSubGroups() arrayView1d< string const > const & regionNames ) { precomputeData( mesh, regionNames ); + + FieldIdentifiers fieldsToBeSync; + fieldsToBeSync.addElementFields( { fields::flow::pressure::key(), fields::flow::temperature::key() }, + regionNames ); + + CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); } ); } @@ -491,6 +487,97 @@ void FlowSolverBase::precomputeData( MeshLevel & mesh, } } +void FlowSolverBase::initialize( DomainPartition & domain ) +{ + // Compute hydrostatic equilibrium in the regions for which corresponding field specification tag has been specified + computeHydrostaticEquilibrium( domain ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + initializePorosityAndPermeability( mesh, regionNames ); + initializeHydraulicAperture( mesh, regionNames ); + + // Initialize primary variables from applied initial conditions + initializeFluidState( mesh, regionNames ); + + // Initialize the rock thermal quantities: conductivity and solid internal energy + // Note: + // - This must be called after updatePorosityAndPermeability and updatePhaseVolumeFraction + // - This step depends on porosity and phaseVolFraction + if( m_isThermal ) + { + initializeThermalState( mesh, regionNames ); + } + + // Save initial pressure and temperature fields + saveInitialPressureAndTemperature( mesh, regionNames ); + } ); + + // report to the user if some pore volumes are very small + // note: this function is here because: 1) porosity has been initialized and 2) NTG has been applied + validatePoreVolumes( domain ); +} + +void FlowSolverBase::initializePorosityAndPermeability( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + // Update porosity and permeability + // In addition, to avoid multiplying permeability/porosity bay netToGross in the assembly kernel, we do it once and for all here + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // Apply netToGross to reference porosity and horizontal permeability + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + PermeabilityBase const & permeability = + getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) ); + arrayView1d< real64 const > const netToGross = subRegion.template getField< fields::flow::netToGross >(); + porousSolid.scaleReferencePorosity( netToGross ); + permeability.scaleHorizontalPermeability( netToGross ); + + // in some initializeState versions it uses newPorosity, so let's run updatePorosityAndPermeability to compute something + saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes + updatePorosityAndPermeability( subRegion ); + porousSolid.initializeState(); + + // run final update + saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes + updatePorosityAndPermeability( subRegion ); + + // Save the computed porosity into the old porosity + // Note: + // - This must be called after updatePorosityAndPermeability + // - This step depends on porosity + porousSolid.saveConvergedState(); + } ); +} + +void FlowSolverBase::initializeHydraulicAperture( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) +{ + mesh.getElemManager().forElementRegions< SurfaceElementRegion >( regionNames, + [&]( localIndex const, + SurfaceElementRegion & region ) + { + region.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) + { subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key()).setApplyDefaultValue( region.getDefaultAperture()); } ); + } ); +} + +void FlowSolverBase::saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) +{ + mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + 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.getField< fields::flow::temperature >(); + arrayView1d< real64 > const initTemp = subRegion.template getField< fields::flow::initialTemperature >(); + initPres.setValues< parallelDevicePolicy<> >( pres ); + initTemp.setValues< parallelDevicePolicy<> >( temp ); + } ); +} + void FlowSolverBase::updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const { GEOS_MARK_FUNCTION; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 6b65e91e979..1d46d170e37 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -100,9 +100,9 @@ class FlowSolverBase : public PhysicsSolverBase */ void updateStencilWeights( DomainPartition & domain ) const; - void enableFixedStressPoromechanicsUpdate(); + void enableFixedStressPoromechanicsUpdate() { m_isFixedStressPoromechanicsUpdate = true; } - void enableJumpStabilization(); + void enableJumpStabilization() { m_isJumpStabilized = true; } void updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const; @@ -114,39 +114,12 @@ class FlowSolverBase : public PhysicsSolverBase */ virtual void saveSequentialIterationState( DomainPartition & domain ) override; - /** - * @brief For each equilibrium initial condition, loop over all the target cells and compute the min/max elevation - * @param[in] domain the domain partition - * @param[in] equilNameToEquilId the map from the name of the initial condition to the initial condition index (used in min/maxElevation) - * @param[out] maxElevation the max elevation for each initial condition - * @param[out] minElevation the min elevation for each initial condition - */ - void findMinMaxElevationInEquilibriumTarget( DomainPartition & domain, // cannot be const... - std::map< string, localIndex > const & equilNameToEquilId, - arrayView1d< real64 > const & maxElevation, - arrayView1d< real64 > const & minElevation ) const; - - /** - * @brief For each source flux boundary condition, loop over all the target cells and sum the owned cells - * @param[in] time the time at the beginning of the time step - * @param[in] dt the time step size - * @param[in] domain the domain partition - * @param[in] bcNameToBcId the map from the name of the boundary condition to the boundary condition index - * @param[out] bcAllSetsSize the total number of owned cells for each source flux boundary condition - */ - void computeSourceFluxSizeScalingFactor( real64 const & time, - real64 const & dt, - DomainPartition & domain, // cannot be const... - std::map< string, localIndex > const & bcNameToBcId, - arrayView1d< globalIndex > const & bcAllSetsSize ) const; - integer & isThermal() { return m_isThermal; } /** * @return The unit in which we evaluate the amount of fluid per element (Mass or Mole). */ - virtual units::Unit getMassUnit() const - { return units::Unit::Mass; } + virtual units::Unit getMassUnit() const { return units::Unit::Mass; } /** * @brief Function to activate the flag allowing negative pressure @@ -173,6 +146,36 @@ class FlowSolverBase : public PhysicsSolverBase real64 const & timeAtBeginningOfStep, real64 const & dt ); + virtual void initializeFluidState( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } + + virtual void initializeThermalState( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } + + /** + * @brief For each equilibrium initial condition, loop over all the target cells and compute the min/max elevation + * @param[in] domain the domain partition + * @param[in] equilNameToEquilId the map from the name of the initial condition to the initial condition index (used in min/maxElevation) + * @param[out] maxElevation the max elevation for each initial condition + * @param[out] minElevation the min elevation for each initial condition + */ + void findMinMaxElevationInEquilibriumTarget( DomainPartition & domain, // cannot be const... + std::map< string, localIndex > const & equilNameToEquilId, + arrayView1d< real64 > const & maxElevation, + arrayView1d< real64 > const & minElevation ) const; + + /** + * @brief For each source flux boundary condition, loop over all the target cells and sum the owned cells + * @param[in] time the time at the beginning of the time step + * @param[in] dt the time step size + * @param[in] domain the domain partition + * @param[in] bcNameToBcId the map from the name of the boundary condition to the boundary condition index + * @param[out] bcAllSetsSize the total number of owned cells for each source flux boundary condition + */ + void computeSourceFluxSizeScalingFactor( real64 const & time, + real64 const & dt, + DomainPartition & domain, // cannot be const... + std::map< string, localIndex > const & bcNameToBcId, + arrayView1d< globalIndex > const & bcAllSetsSize ) const; + protected: /** @@ -207,6 +210,16 @@ class FlowSolverBase : public PhysicsSolverBase virtual void initializePostInitialConditionsPreSubGroups() override; + void initialize( DomainPartition & domain ); + + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); } + + void initializePorosityAndPermeability( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + + void initializeHydraulicAperture( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); + + void saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); + virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; /// the number of Degrees of Freedom per cell diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index f1a77c2001d..aca05df65e3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -404,105 +404,12 @@ void SinglePhaseBase::initializePostInitialConditionsPreSubGroups() DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - FieldIdentifiers fieldsToBeSync; - fieldsToBeSync.addElementFields( { fields::flow::pressure::key() }, - regionNames ); - - CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); - - // Moved the following part from ImplicitStepSetup to here since it only needs to be initialized once - // They will be updated in applySystemSolution and ImplicitStepComplete, respectively - mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, - auto & subRegion ) - { - // Compute hydrostatic equilibrium in the regions for which corresponding field specification tag has been specified - computeHydrostaticEquilibrium(); - - // 1. update porosity, permeability, and density/viscosity - // In addition, to avoid multiplying permeability/porosity bay netToGross in the assembly kernel, we do it once and for all here - arrayView1d< real64 const > const netToGross = subRegion.template getField< fields::flow::netToGross >(); - CoupledSolidBase const & porousSolid = - getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); - PermeabilityBase const & permeabilityModel = - getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) ); - permeabilityModel.scaleHorizontalPermeability( netToGross ); - porousSolid.scaleReferencePorosity( netToGross ); - saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes - updatePorosityAndPermeability( subRegion ); - - SingleFluidBase const & fluid = - getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) ); - updateFluidState( subRegion ); - - // 2. save the initial density (for use in the single-phase poromechanics solver to compute the deltaBodyForce) - fluid.initializeState(); - - // 3. save the initial/old porosity - porousSolid.initializeState(); - - // 4. initialize the rock thermal quantities: conductivity and solid internal energy - if( m_isThermal ) - { - // initialized porosity - arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); - - string const & thermalConductivityName = subRegion.template getReference< string >( viewKeyStruct::thermalConductivityNamesString() ); - SinglePhaseThermalConductivityBase const & conductivityMaterial = - getConstitutiveModel< SinglePhaseThermalConductivityBase >( subRegion, thermalConductivityName ); - conductivityMaterial.initializeRockFluidState( porosity ); - // note that there is nothing to update here because thermal conductivity is explicit for now - - updateSolidInternalEnergyModel( subRegion ); - string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString() ); - SolidInternalEnergy const & solidInternalEnergyMaterial = - getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); - solidInternalEnergyMaterial.saveConvergedState(); - } - } ); - - mesh.getElemManager().forElementRegions< SurfaceElementRegion >( regionNames, - [&]( localIndex const, - SurfaceElementRegion & region ) - { - region.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) - { - subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key() ). - setApplyDefaultValue( region.getDefaultAperture() ); - } ); - } ); - - 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 ); - } ); - } ); - - // report to the user if some pore volumes are very small - // note: this function is here because: 1) porosity has been initialized and 2) NTG has been applied - validatePoreVolumes( domain ); + FlowSolverBase::initialize( domain ); } -void SinglePhaseBase::computeHydrostaticEquilibrium() +void SinglePhaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) { FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); - DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); real64 const gravVector[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( gravityVector() ); @@ -676,6 +583,48 @@ void SinglePhaseBase::computeHydrostaticEquilibrium() } ); } +void SinglePhaseBase::initializeFluidState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + SingleFluidBase const & fluid = + getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString())); + updateFluidState( subRegion ); + + // 2. save the initial density (for use in the single-phase poromechanics solver to compute the deltaBodyForce) + fluid.initializeState(); + + updateMass( subRegion ); + } ); +} + +void SinglePhaseBase::initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +{ + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // initialized porosity + CoupledSolidBase const & porousSolid = + getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); + arrayView2d< real64 const > const porosity = porousSolid.getPorosity(); + + string const & thermalConductivityName = subRegion.template getReference< string >( viewKeyStruct::thermalConductivityNamesString()); + SinglePhaseThermalConductivityBase const & conductivityMaterial = + getConstitutiveModel< SinglePhaseThermalConductivityBase >( subRegion, thermalConductivityName ); + conductivityMaterial.initializeRockFluidState( porosity ); + // note that there is nothing to update here because thermal conductivity is explicit for now + + updateSolidInternalEnergyModel( subRegion ); + string const & solidInternalEnergyName = subRegion.template getReference< string >( viewKeyStruct::solidInternalEnergyNamesString()); + SolidInternalEnergy const & solidInternalEnergyMaterial = + getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName ); + solidInternalEnergyMaterial.saveConvergedState(); + + updateEnergy( subRegion ); + } ); +} + void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_n ), real64 const & GEOS_UNUSED_PARAM( dt ), DomainPartition & domain ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index bf1cd19eadd..ad175412dc1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -335,18 +335,14 @@ class SinglePhaseBase : public FlowSolverBase virtual void initializePostInitialConditionsPreSubGroups() override; - /** - * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables - */ - void computeHydrostaticEquilibrium(); + virtual void initializeFluidState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + + virtual void initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** - * @brief Update the cell-wise pressure gradient + * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ - virtual void updatePressureGradient( DomainPartition & domain ) - { - GEOS_UNUSED_VAR( domain ); - } + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; /** * @brief Function to fix the initial state during the initialization step in coupled problems diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp index 233eb1d1363..f2022bc9e40 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp @@ -656,22 +656,5 @@ void SinglePhaseHybridFVM::resetStateToBeginningOfStep( DomainPartition & domain } ); } -void SinglePhaseHybridFVM::updatePressureGradient( DomainPartition & domain ) -{ - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - FaceManager & faceManager = mesh.getFaceManager(); - - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - auto & subRegion ) - { - singlePhaseHybridFVMKernels::AveragePressureGradientKernelFactory::createAndLaunch< parallelHostPolicy >( subRegion, - faceManager ); - } ); - } ); -} - REGISTER_CATALOG_ENTRY( PhysicsSolverBase, SinglePhaseHybridFVM, string const &, Group * const ) } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp index a1ae07fb5b5..f8427490c60 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp @@ -164,8 +164,6 @@ class SinglePhaseHybridFVM : public SinglePhaseBase real64 const & dt, DomainPartition & domain ) override; - virtual void updatePressureGradient( DomainPartition & domain ) override final; - /** * @brief Function to perform the application of Dirichlet BCs on faces * @param[in] time_n current time