From f181983859162d92a53b2d1b6782be11034a8762 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 10 Oct 2024 14:23:00 -0500 Subject: [PATCH 01/27] un..ck flow initialization a bit --- .../fluidFlow/CompositionalMultiphaseBase.cpp | 257 +++++++----------- .../fluidFlow/CompositionalMultiphaseBase.hpp | 8 +- .../fluidFlow/FlowSolverBase.cpp | 89 ++++++ .../fluidFlow/FlowSolverBase.hpp | 14 + .../fluidFlow/SinglePhaseBase.hpp | 6 +- 5 files changed, 211 insertions(+), 163 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index a8d6f99ecec..44b564052f5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -819,31 +819,29 @@ real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & sub return maxDeltaPhaseVolFrac; } -void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, - DomainPartition & domain, - arrayView1d< string const > const & regionNames ) +void CompositionalMultiphaseBase::initializeFluid( MeshLevel & mesh, + 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. + // set mass fraction flag on fluid models + string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); + MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); + fluid.setMassFlag( m_useMass ); + + // 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() ); - 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 = @@ -856,10 +854,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. @@ -867,115 +868,64 @@ 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) + // 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 - // 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 - // initialized phase volume fraction arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = subRegion.template getField< fields::flow::phaseVolumeFraction >(); + // Initialize/update the relative permeability model using the initial phase volume fraction + // Note: + // - This must be called after updatePhaseVolumeFraction + // - This step depends on phaseVolFraction string const & relpermName = subRegion.template getReference< string >( viewKeyStruct::relPermNamesString() ); - RelativePermeabilityBase & relPermMaterial = - getConstitutiveModel< RelativePermeabilityBase >( subRegion, relpermName ); + RelativePermeabilityBase & relPermMaterial = getConstitutiveModel< RelativePermeabilityBase >( subRegion, relpermName ); relPermMaterial.saveConvergedPhaseVolFractionState( phaseVolFrac ); // this needs to happen before calling updateRelPermModel updateRelPermModel( subRegion ); relPermMaterial.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() ); @@ -993,24 +943,42 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, } } ); +} - // 5. Save initial pressure - mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) +void CompositionalMultiphaseBase::initializeThermal( 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; @@ -1270,49 +1238,13 @@ 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 ); - - mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, - [&]( localIndex const, - auto & subRegion ) - { - // set mass fraction flag on fluid models - 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 @@ -2043,9 +1975,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 ) @@ -2054,25 +1983,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 55dcb3460b6..1e64811dd6f 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 ); + void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + + void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); /** * @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 5c54deae9b5..a9533947b4a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -455,6 +455,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 +497,89 @@ 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 + initializeFluid( 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 ) + { + initializeThermal( 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 ) + { + 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 ); + + saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes + + updatePorosityAndPermeability( subRegion ); + + // save the initial/old porosity + porousSolid.initializeState(); + } ); +} + +void FlowSolverBase::initializeHydraulicAperture( MeshLevel & mesh, const arrayView1d< const string > & 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()); } ); + } ); +} + +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 883661f25a0..89ae8b618c9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -190,6 +190,20 @@ class FlowSolverBase : public SolverBase virtual void precomputeData( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + 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 ); + + virtual void initializeFluid( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) {GEOS_UNUSED_VAR( mesh, regionNames );} + + virtual void initializeThermal( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) {GEOS_UNUSED_VAR( mesh, regionNames );} + + void saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); + virtual void initializePreSubGroups() override; virtual void initializePostInitialConditionsPreSubGroups() override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index 59ad25124b4..388d9830020 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -338,7 +338,11 @@ class SinglePhaseBase : public FlowSolverBase /** * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ - void computeHydrostaticEquilibrium(); + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; + + void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + + void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); /** * @brief Update the cell-wise pressure gradient From 053377b24c905466e0c4a483ab1b028ee553b1d0 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 10 Oct 2024 15:40:59 -0500 Subject: [PATCH 02/27] missing file --- .../fluidFlow/SinglePhaseBase.cpp | 152 ++++++------------ 1 file changed, 51 insertions(+), 101 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 05dbc006843..f414cf54df9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -399,105 +399,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< FaceElementSubRegion >( [&]( FaceElementSubRegion & 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() ); @@ -671,6 +578,48 @@ void SinglePhaseBase::computeHydrostaticEquilibrium() } ); } +void SinglePhaseBase::initializeFluid( 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::initializeThermal( 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 ) @@ -711,7 +660,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_ CoupledSolidBase const & porousSolid = getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::solidNamesString() ) ); porousSolid.saveConvergedState(); - saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes + saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes updatePorosityAndPermeability( subRegion ); updateFluidState( subRegion ); @@ -765,11 +714,11 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); if( m_keepFlowVariablesConstantDuringInitStep ) { - porousSolid.ignoreConvergedState(); // newPorosity <- porosity_n + porousSolid.ignoreConvergedState(); // newPorosity <- porosity_n } else { - porousSolid.saveConvergedState(); // porosity_n <- porosity + porousSolid.saveConvergedState(); // porosity_n <- porosity } if( m_isThermal ) @@ -1125,7 +1074,8 @@ void SinglePhaseBase::applySourceFluxBC( real64 const time_n, // add the value to the mass balance equation globalIndex const massRowIndex = dofNumber[ei] - rankOffset; globalIndex const energyRowIndex = massRowIndex + 1; - real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor here! + real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor + // here! localRhs[massRowIndex] += rhsValue; massProd += rhsValue; //add the value to the energy balance equation if the flux is positive (i.e., it's a producer) @@ -1229,7 +1179,7 @@ void SinglePhaseBase::keepFlowVariablesConstantDuringInitStep( real64 const time rankOffset, localMatrix, rhsValue, - pres[ei], // freeze the current pressure value + pres[ei], // freeze the current pressure value pres[ei] ); localRhs[localRow] = rhsValue; @@ -1240,7 +1190,7 @@ void SinglePhaseBase::keepFlowVariablesConstantDuringInitStep( real64 const time rankOffset, localMatrix, rhsValue, - temp[ei], // freeze the current temperature value + temp[ei], // freeze the current temperature value temp[ei] ); localRhs[localRow + 1] = rhsValue; } From ec91135a243735f6929fb4b0e50556ea5b9cdd7d Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 10 Oct 2024 16:04:40 -0500 Subject: [PATCH 03/27] build fix --- .../physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp | 4 ++-- .../physicsSolvers/fluidFlow/SinglePhaseBase.hpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 1e64811dd6f..d58a4b7e364 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -282,9 +282,9 @@ class CompositionalMultiphaseBase : public FlowSolverBase * from prescribed intermediate values (i.e. global densities from global fractions) * and any applicable hydrostatic equilibration of the domain */ - void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index 388d9830020..98b38548e0e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -340,9 +340,9 @@ class SinglePhaseBase : public FlowSolverBase */ virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; - void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** * @brief Update the cell-wise pressure gradient From 6018d8c7542ddd8d2f13239cbfbcdb3d04dd0d77 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 11 Oct 2024 14:03:29 -0500 Subject: [PATCH 04/27] bug fix and functions cleanup --- .../fluidFlow/CompositionalMultiphaseBase.cpp | 27 ++- .../fluidFlow/CompositionalMultiphaseBase.hpp | 92 +++++----- .../fluidFlow/CompositionalMultiphaseFVM.hpp | 3 +- .../CompositionalMultiphaseHybridFVM.hpp | 9 +- .../fluidFlow/FlowSolverBase.cpp | 15 +- .../fluidFlow/FlowSolverBase.hpp | 61 ++++--- .../fluidFlow/SinglePhaseBase.hpp | 167 +++++++++--------- .../fluidFlow/SinglePhaseFVM.hpp | 2 + .../fluidFlow/SinglePhaseHybridFVM.cpp | 17 -- .../fluidFlow/SinglePhaseHybridFVM.hpp | 4 +- .../fluidFlow/SinglePhaseProppantBase.hpp | 4 +- 11 files changed, 189 insertions(+), 212 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 44b564052f5..40c8af293ed 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -869,13 +869,14 @@ void CompositionalMultiphaseBase::initializeFluid( MeshLevel & mesh, auto & subRegion ) { // Initialize/update dependent state quantities - + + updateCompAmount( subRegion ); + updatePhaseVolumeFraction( subRegion ); + // 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 - updateCompAmount( subRegion ); - updatePhaseVolumeFraction( subRegion ); // Now, we initialize and update each constitutive model one by one @@ -887,11 +888,11 @@ void CompositionalMultiphaseBase::initializeFluid( MeshLevel & mesh, // Note: // - This must be called after updatePhaseVolumeFraction // - This step depends on phaseVolFraction - 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 + 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 & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); @@ -978,6 +979,8 @@ void CompositionalMultiphaseBase::initializeThermal( MeshLevel & mesh, arrayView void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) { + std::cout << "CompositionalMultiphaseBase::computeHydrostaticEquilibrium" << std::endl; + FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); integer const numComps = m_numComponents; @@ -1242,6 +1245,16 @@ void CompositionalMultiphaseBase::initializePostInitialConditionsPreSubGroups() regionNames ); CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false ); + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + auto & subRegion ) + { + // set mass fraction flag on fluid models + string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); + MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); + fluid.setMassFlag( m_useMass ); + } ); } ); initialize( domain ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index d58a4b7e364..ffdcd5a541f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -274,6 +274,45 @@ class CompositionalMultiphaseBase : public FlowSolverBase }; + virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, + DomainPartition & domain ) override; + + /** + * @brief function to set the next time step size + * @param[in] currentDt the current time step size + * @param[in] domain the domain object + * @return the prescribed time step size + */ + real64 setNextDt( real64 const & currentDt, + DomainPartition & domain ) override; + + virtual real64 setNextDtBasedOnCFL( real64 const & currentDt, + DomainPartition & domain ) override; + + void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL ); + + integer useSimpleAccumulation() const { return m_useSimpleAccumulation; } + + integer useTotalMassEquation() const { return m_useTotalMassEquation; } + + virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override; + +protected: + + virtual void postInputInitialization() override; + + virtual void initializePreSubGroups() override; + + virtual void initializePostInitialConditionsPreSubGroups() override; + + /** + * @brief Utility function that checks the consistency of the constitutive models + * @param[in] domain the domain partition + * This function will produce an error if one of the constitutive models + * (fluid, relperm) is incompatible with the reference fluid model. + */ + void validateConstitutiveModels( DomainPartition const & domain ) const; + /** * @brief Initialize all variables from initial conditions * @param domain the domain containing the mesh and fields @@ -291,6 +330,12 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; + /** + * @brief Initialize the aquifer boundary condition (gravity vector, water phase index) + * @param[in] cm reference to the global constitutive model manager + */ + void initializeAquiferBC( constitutive::ConstitutiveManager const & cm ) const; + /** * @brief Function to perform the Application of Dirichlet type BC's * @param time current time @@ -357,7 +402,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const; - /** * @brief Sets all the negative component densities (if any) to zero. * @param domain the physical domain object @@ -366,52 +410,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase void chopNegativeDensities( ElementSubRegionBase & subRegion ); - virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, - DomainPartition & domain ) override; - - void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL ); - - /** - * @brief function to set the next time step size - * @param[in] currentDt the current time step size - * @param[in] domain the domain object - * @return the prescribed time step size - */ - real64 setNextDt( real64 const & currentDt, - DomainPartition & domain ) override; - - virtual real64 setNextDtBasedOnCFL( real64 const & currentDt, - DomainPartition & domain ) override; - - virtual void initializePostInitialConditionsPreSubGroups() override; - - integer useSimpleAccumulation() const { return m_useSimpleAccumulation; } - - integer useTotalMassEquation() const { return m_useTotalMassEquation; } - - virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override; - -protected: - - virtual void postInputInitialization() override; - - virtual void initializePreSubGroups() override; - - - /** - * @brief Utility function that checks the consistency of the constitutive models - * @param[in] domain the domain partition - * This function will produce an error if one of the constitutive models - * (fluid, relperm) is incompatible with the reference fluid model. - */ - void validateConstitutiveModels( DomainPartition const & domain ) const; - - /** - * @brief Initialize the aquifer boundary condition (gravity vector, water phase index) - * @param[in] cm reference to the global constitutive model manager - */ - void initializeAquiferBC( constitutive::ConstitutiveManager const & cm ) const; - /** * @brief Utility function that encapsulates the call to FieldSpecificationBase::applyFieldValue in BC application * @param[in] time_n the time at the beginning of the step diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp index 93b54d898a7..1e5fc325c7b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp @@ -178,8 +178,7 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase virtual void postInputInitialization() override; - virtual void - initializePreSubGroups() override; + virtual void initializePreSubGroups() override; struct DBCParameters { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp index 5562f6202c9..9fe8e8b702a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp @@ -141,9 +141,6 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const override; - virtual void - updatePhaseMobility( ObjectManagerBase & dataGroup ) const override; - virtual void applyAquiferBC( real64 const time, real64 const dt, @@ -164,15 +161,17 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase static constexpr char const * faceDofFieldString() { return "faceCenteredVariables"; } }; - virtual void initializePostInitialConditionsPreSubGroups() override; +protected: virtual void initializePreSubGroups() override; -protected: + virtual void initializePostInitialConditionsPreSubGroups() override; /// precompute the minGravityCoefficient for the buoyancy term void precomputeData( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const override; + private: /// tolerance used in the computation of the transmissibility matrix diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index a9533947b4a..a6c7e930387 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 { SolverBase::setConstitutiveNamesCallSuper( subRegion ); @@ -550,7 +540,10 @@ void FlowSolverBase::initializePorosityAndPermeability( MeshLevel & mesh, arrayV updatePorosityAndPermeability( subRegion ); - // save the initial/old porosity + // Save the computed porosity into the old porosity + // Note: + // - This must be called after updatePorosityAndPermeability + // - This step depends on porosity porousSolid.initializeState(); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 89ae8b618c9..d13e7940064 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -95,9 +95,9 @@ class FlowSolverBase : public SolverBase */ void updateStencilWeights( DomainPartition & domain ) const; - void enableFixedStressPoromechanicsUpdate(); + void enableFixedStressPoromechanicsUpdate() { m_isFixedStressPoromechanicsUpdate = true; } - void enableJumpStabilization(); + void enableJumpStabilization() { m_isJumpStabilized = true; } void updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const; @@ -109,39 +109,12 @@ class FlowSolverBase : public SolverBase */ 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 @@ -190,9 +163,21 @@ class FlowSolverBase : public SolverBase virtual void precomputeData( MeshLevel & mesh, arrayView1d< string const > const & 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; + void initialize( DomainPartition & domain ); - virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) {GEOS_UNUSED_VAR( domain );} + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) {GEOS_UNUSED_VAR( domain ); std::cout << "computeHydrostaticEquilibrium" << std::endl;} void initializePorosityAndPermeability( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); @@ -210,6 +195,20 @@ class FlowSolverBase : public SolverBase virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; + /** + * @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; + /// the number of Degrees of Freedom per cell integer m_numDofPerCell; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index 98b38548e0e..75bc31dde77 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -229,67 +229,49 @@ class SinglePhaseBase : public FlowSolverBase static constexpr char const * elemDofFieldString() { return "singlePhaseVariables"; } }; + virtual void + updateState ( DomainPartition & domain ) override final; + /** - * @brief Function to perform the Application of Dirichlet type BC's - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector + * @brief Function to update all constitutive state and dependent variables + * @param subRegion subregion that contains the fields */ - void - applyDirichletBC( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; + real64 + updateFluidState( ElementSubRegionBase & subRegion ) const; /** - * @brief Apply source flux boundary conditions to the system - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector + * @brief Update all relevant solid internal energy models using current values of temperature + * @param dataGroup the group storing the required fields */ - void - applySourceFluxBC( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; + void updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const; + +protected: + + virtual void initializePreSubGroups() override; + + virtual void initializePostInitialConditionsPreSubGroups() override; /** - * @brief Apply aquifer boundary conditions to the system - * @param time current time - * @param dt time step - * @param domain the domain - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector + * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ - virtual void - applyAquiferBC( real64 const time, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const = 0; + virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; - virtual void - updateState ( DomainPartition & domain ) override final; + virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + + virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** - * @brief Function to update all constitutive state and dependent variables - * @param subRegion subregion that contains the fields + * @brief Checks constitutive models for consistency + * @param[in] domain the domain partition */ - real64 - updateFluidState( ElementSubRegionBase & subRegion ) const; + virtual void validateConstitutiveModels( DomainPartition & domain ) const; + /** + * @brief Initialize the aquifer boundary condition (gravity vector, water phase index) + */ + void initializeAquiferBC() const; + + virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; /** * @brief Function to update all constitutive models @@ -312,12 +294,6 @@ class SinglePhaseBase : public FlowSolverBase 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 - */ - void updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const; - /** * @brief Update thermal conductivity * @param subRegion the group storing the required fields @@ -331,26 +307,62 @@ class SinglePhaseBase : public FlowSolverBase void updateMobility( ObjectManagerBase & dataGroup ) const; - virtual void initializePreSubGroups() override; - - virtual void initializePostInitialConditionsPreSubGroups() override; - /** - * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables + * @brief Utility function to save the converged state + * @param[in] subRegion the element subRegion */ - virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; + virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override; - virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /** + * @brief Function to perform the Application of Dirichlet type BC's + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void + applyDirichletBC( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; - virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /** + * @brief Apply source flux boundary conditions to the system + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void + applySourceFluxBC( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; /** - * @brief Update the cell-wise pressure gradient + * @brief Apply aquifer boundary conditions to the system + * @param time current time + * @param dt time step + * @param domain the domain + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector */ - virtual void updatePressureGradient( DomainPartition & domain ) - { - GEOS_UNUSED_VAR( domain ); - } + virtual void + applyAquiferBC( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const = 0; /** * @brief Function to fix the initial state during the initialization step in coupled problems @@ -370,27 +382,6 @@ class SinglePhaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const; -protected: - - /** - * @brief Checks constitutive models for consistency - * @param[in] domain the domain partition - */ - virtual void validateConstitutiveModels( DomainPartition & domain ) const; - - /** - * @brief Initialize the aquifer boundary condition (gravity vector, water phase index) - */ - void initializeAquiferBC() const; - - 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/SinglePhaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp index 882ae4c7e8b..25e0ea546ab 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp @@ -191,6 +191,8 @@ class SinglePhaseFVM : public BASE /**@}*/ +protected: + virtual void applyAquiferBC( real64 const time, real64 const dt, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp index 6cbabf723f7..61873d5d9df 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp @@ -652,22 +652,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( SolverBase, SinglePhaseHybridFVM, string const &, Group * const ) } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp index 21c1bcaf09d..bec559ad6bf 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp @@ -151,6 +151,8 @@ class SinglePhaseHybridFVM : public SinglePhaseBase arrayView1d< real64 > const & localRhs, CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; +protected: + virtual void applyAquiferBC( real64 const time, real64 const dt, @@ -164,8 +166,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 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp index ab07e60f9aa..8b0ae10a990 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp @@ -55,14 +55,14 @@ class SinglePhaseProppantBase : public SinglePhaseBase */ virtual ~SinglePhaseProppantBase(); - virtual void updateFluidModel( ObjectManagerBase & dataGroup ) const override; - virtual void updatePorosityAndPermeability( SurfaceElementSubRegion & subRegion ) const override; protected: virtual void validateConstitutiveModels( DomainPartition & domain ) const override; + virtual void updateFluidModel( ObjectManagerBase & dataGroup ) const override; + virtual FluidPropViews getFluidProperties( constitutive::ConstitutiveBase const & fluid ) const override; private: From c1be658acf33d3977ade6d5f7b2d4766502e8179 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 11 Oct 2024 14:09:29 -0500 Subject: [PATCH 05/27] code style --- .../physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 40c8af293ed..55541fb65d3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -869,10 +869,10 @@ void CompositionalMultiphaseBase::initializeFluid( MeshLevel & mesh, auto & subRegion ) { // Initialize/update dependent state quantities - + updateCompAmount( subRegion ); updatePhaseVolumeFraction( subRegion ); - + // Update the constitutive models that only depend on // - the primary variables // - the fluid constitutive quantities (as they have already been updated) From 1313a67493f0d45be8a916425d01fcbb89a5eada Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 11 Oct 2024 20:40:16 -0500 Subject: [PATCH 06/27] try this --- .../fluidFlow/CompositionalMultiphaseBase.cpp | 9 ++------- .../physicsSolvers/fluidFlow/FlowSolverBase.cpp | 9 +++++++-- .../physicsSolvers/fluidFlow/FlowSolverBase.hpp | 6 +++--- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 55541fb65d3..3e3aeb4eacc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -830,17 +830,14 @@ void CompositionalMultiphaseBase::initializeFluid( MeshLevel & mesh, [&]( localIndex const, ElementSubRegionBase & subRegion ) { - // set mass fraction flag on fluid models - string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); - MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); - fluid.setMassFlag( m_useMass ); - // Assume global component fractions have been prescribed. // Initialize constitutive state to get fluid density. updateFluidModel( subRegion ); // Back-calculate global component densities from fractions and total fluid density // in order to initialize the primary solution variables + string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); + MultiFluidBase & 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 >(); @@ -979,8 +976,6 @@ void CompositionalMultiphaseBase::initializeThermal( MeshLevel & mesh, arrayView void CompositionalMultiphaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) { - std::cout << "CompositionalMultiphaseBase::computeHydrostaticEquilibrium" << std::endl; - FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); integer const numComps = m_numComponents; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index a6c7e930387..7ff0bd4232b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -527,24 +527,29 @@ void FlowSolverBase::initializePorosityAndPermeability( MeshLevel & mesh, arrayV 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.initializeState(); + porousSolid.saveConvergedState(); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index d13e7940064..ae5611d3717 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -177,15 +177,15 @@ class FlowSolverBase : public SolverBase void initialize( DomainPartition & domain ); - virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) {GEOS_UNUSED_VAR( domain ); std::cout << "computeHydrostaticEquilibrium" << std::endl;} + 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 ); - virtual void initializeFluid( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) {GEOS_UNUSED_VAR( mesh, regionNames );} + virtual void initializeFluid( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } - virtual void initializeThermal( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) {GEOS_UNUSED_VAR( mesh, regionNames );} + virtual void initializeThermal( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } void saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); From 81b57c217d101b172d7684b40bca13c8be0095ec Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Sat, 12 Oct 2024 11:27:13 -0500 Subject: [PATCH 07/27] unify poromech init again --- .../PoroElastic_staircase_co2_3d_fim_new.xml | 152 ++++++++++++++++++ .../fluidFlow/FlowSolverBase.hpp | 4 +- .../CoupledReservoirAndWellsBase.hpp | 2 + .../PoromechanicsInitialization.cpp | 5 +- .../multiphysics/PoromechanicsSolver.hpp | 51 +++--- .../SolidMechanicsLagrangianFEM.cpp | 5 +- .../SolidMechanicsLagrangianFEM.hpp | 13 ++ 7 files changed, 202 insertions(+), 30 deletions(-) create mode 100644 inputFiles/poromechanics/PoroElastic_staircase_co2_3d_fim_new.xml diff --git a/inputFiles/poromechanics/PoroElastic_staircase_co2_3d_fim_new.xml b/inputFiles/poromechanics/PoroElastic_staircase_co2_3d_fim_new.xml new file mode 100644 index 00000000000..ab47b718968 --- /dev/null +++ b/inputFiles/poromechanics/PoroElastic_staircase_co2_3d_fim_new.xml @@ -0,0 +1,152 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index ae5611d3717..b990c8abce2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -133,6 +133,8 @@ class FlowSolverBase : public SolverBase void enableLaggingFractureStencilWeightsUpdate(){ m_isLaggingFractureStencilWeightsUpdate = 1; }; + void initialize( DomainPartition & domain ); + protected: /** @@ -175,8 +177,6 @@ class FlowSolverBase : public SolverBase arrayView1d< real64 > const & maxElevation, arrayView1d< real64 > const & minElevation ) const; - void initialize( DomainPartition & domain ); - virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); } void initializePorosityAndPermeability( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp index af10776cdf4..3359cebbf6f 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp @@ -223,6 +223,8 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL } } + void initialize( DomainPartition & domain ) const { return reservoirSolver()->initialize( domain ); } + void assembleFluxTerms( real64 const dt, DomainPartition const & domain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index e86c9553544..9d030c30063 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -84,8 +84,9 @@ postInputInitialization() TasksManager & tasksManager = problemManager.getGroup< TasksManager >( "Tasks" ); GEOS_THROW_IF( !tasksManager.hasGroup( m_solidMechanicsStatisticsName ), - GEOS_FMT( "{}: statistics task named {} not found", + GEOS_FMT( "{}: {} task named {} not found", getWrapperDataContext( viewKeyStruct::solidMechanicsStatisticsNameString() ), + SolidMechanicsStatistics::catalogName(), m_solidMechanicsStatisticsName ), InputError ); @@ -132,6 +133,7 @@ execute( real64 const time_n, namespace { +typedef PoromechanicsInitialization< SolidMechanicsLagrangianFEM > MechanicsInitialization; typedef PoromechanicsInitialization< MultiphasePoromechanics<> > MultiphasePoromechanicsInitialization; typedef PoromechanicsInitialization< MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> > > MultiphaseReservoirPoromechanicsInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanics<> > SinglePhasePoromechanicsInitialization; @@ -139,6 +141,7 @@ typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFractures typedef PoromechanicsInitialization< SinglePhasePoromechanicsEmbeddedFractures > SinglePhasePoromechanicsEmbeddedFracturesInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanics< SinglePhaseReservoirAndWells<> > > SinglePhaseReservoirPoromechanicsInitialization; typedef PoromechanicsInitialization< HydrofractureSolver< SinglePhasePoromechanics<> > > HydrofractureInitialization; +REGISTER_CATALOG_ENTRY( TaskBase, MechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, MultiphaseReservoirPoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsInitialization, string const &, Group * const ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index b2fd6d4a603..05a922d6cfb 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -80,18 +80,14 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER PoromechanicsSolver( const string & name, dataRepository::Group * const parent ) : Base( name, parent ), - m_isThermal( 0 ) + m_isThermal( 0 ), + m_performStressInitialization( false ) { this->registerWrapper( viewKeyStruct::isThermalString(), &m_isThermal ). setApplyDefaultValue( 0 ). setInputFlag( dataRepository::InputFlags::OPTIONAL ). setDescription( "Flag indicating whether the problem is thermal or not. Set isThermal=\"1\" to enable the thermal coupling" ); - this->registerWrapper( viewKeyStruct::performStressInitializationString(), &m_performStressInitialization ). - setApplyDefaultValue( false ). - setInputFlag( dataRepository::InputFlags::FALSE ). - setDescription( "Flag to indicate that the solver is going to perform stress initialization" ); - this->registerWrapper( viewKeyStruct::stabilizationTypeString(), &m_stabilizationType ). setInputFlag( dataRepository::InputFlags::OPTIONAL ). setDescription( "StabilizationType. Options are:\n" + @@ -118,6 +114,10 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER GEOS_FMT( "{} {}: The attribute `{}` of the flow solver `{}` must be set to 1 since the poromechanics solver is thermal", this->getCatalogName(), this->getName(), FlowSolverBase::viewKeyStruct::isThermalString(), this->flowSolver()->getName() ), InputError ); + + DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" ); + flowSolver()->initialize( domain ); + updateBulkDensity( domain ); } virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override final @@ -219,7 +219,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER setRestartFlags( dataRepository::RestartFlags::NO_WRITE ). setSizedFromParent( 0 ); - if( this->getNonlinearSolverParameters().m_couplingType == NonlinearSolverParameters::CouplingType::Sequential ) + //if( this->getNonlinearSolverParameters().m_couplingType == NonlinearSolverParameters::CouplingType::Sequential ) { // register the bulk density for use in the solid mechanics solver // ideally we would resize it here as well, but the solid model name is not available yet (see below) @@ -314,9 +314,6 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER /// Flag to indicate that the simulation is thermal constexpr static char const * isThermalString() { return "isThermal"; } - /// Flag to indicate that the solver is going to perform stress initialization - constexpr static char const * performStressInitializationString() { return "performStressInitialization"; } - /// Type of pressure stabilization constexpr static char const * stabilizationTypeString() {return "stabilizationType"; } @@ -521,20 +518,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER /// After the flow solver if( solverType == static_cast< integer >( SolverType::Flow ) ) { - this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - auto & subRegion ) - { - // update the bulk density - // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the - // pressure/temperature terms - updateBulkDensity( subRegion ); - } ); - } ); + updateBulkDensity( domain ); } /// After the solid mechanics solver @@ -614,6 +598,23 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER } ); } + void updateBulkDensity( DomainPartition & domain ) + { + this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // update the bulk density + // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the + // pressure/temperature terms + updateBulkDensity( subRegion ); + } ); + } ); + } + virtual void updateBulkDensity( ElementSubRegionBase & subRegion ) = 0; virtual void validateNonlinearAcceleration() override @@ -628,7 +629,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER integer m_isThermal; /// Flag to indicate that the solver is going to perform stress initialization - integer m_performStressInitialization; + bool m_performStressInitialization; /// Type of stabilization used stabilization::StabilizationType m_stabilizationType; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 0ad7c46f776..e64e4b8428c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -60,7 +60,8 @@ SolidMechanicsLagrangianFEM::SolidMechanicsLagrangianFEM( const string & name, m_maxNumResolves( 10 ), m_strainTheory( 0 ), m_iComm( CommunicationTools::getInstance().getCommID() ), - m_isFixedStressPoromechanicsUpdate( false ) + m_isFixedStressPoromechanicsUpdate( false ), + m_performStressInitialization( false ) { registerWrapper( viewKeyStruct::newmarkGammaString(), &m_newmarkGamma ). @@ -1055,7 +1056,7 @@ void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - if( m_isFixedStressPoromechanicsUpdate ) + if( m_isFixedStressPoromechanicsUpdate || m_performStressInitialization ) { set< string > poromechanicsRegions; set< string > mechanicsRegions; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index c59227142b1..5af604d4076 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -279,6 +279,15 @@ class SolidMechanicsLagrangianFEM : public SolverBase return m_rigidBodyModes; } + /* + * @brief Utility function to set the stress initialization flag + * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise + */ + void setStressInitialization( integer const performStressInitialization ) + { + m_performStressInitialization = performStressInitialization; + } + protected: virtual void postInputInitialization() override; @@ -295,7 +304,11 @@ class SolidMechanicsLagrangianFEM : public SolverBase integer m_maxNumResolves; integer m_strainTheory; MPI_iCommData m_iComm; + + /// Flag to indicate that the solver is running in fixed stress (sequential) mode bool m_isFixedStressPoromechanicsUpdate; + /// Flag to indicate that the solver is going to perform stress initialization + bool m_performStressInitialization; /// Rigid body modes array1d< ParallelVector > m_rigidBodyModes; From 035a98637860834102300e7668035f82c8a1a77e Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 14 Nov 2024 13:24:15 -0600 Subject: [PATCH 08/27] Update CompositionalMultiphaseBase.cpp --- .../physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 513a7337618..1c4fb095738 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -848,7 +848,7 @@ void CompositionalMultiphaseBase::initializeFluid( MeshLevel & mesh, // Back-calculate global component densities from fractions and total fluid density // in order to initialize the primary solution variables string const & fluidName = subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ); - MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); + 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 >(); From d9ed44f06fddb12611efbf0e7de98adbf944718b Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 15 Nov 2024 08:29:29 -0600 Subject: [PATCH 09/27] cuda build --- .../fluidFlow/CompositionalMultiphaseBase.hpp | 143 +++++++++--------- .../CompositionalMultiphaseHybridFVM.hpp | 4 +- .../fluidFlow/FlowSolverBase.hpp | 8 +- .../fluidFlow/SinglePhaseBase.hpp | 4 +- 4 files changed, 80 insertions(+), 79 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index d5c91e8376e..6dd01a081f3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -68,6 +68,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual void registerDataOnMesh( Group & meshBodies ) override; + virtual void initializePostInitialConditionsPreSubGroups() override; + /** * @defgroup Solver Interface Functions * @@ -166,6 +168,12 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual void updateState( DomainPartition & domain ) override final; + /** + * @brief Sets all the negative component densities (if any) to zero. + * @param domain the physical domain object + */ + void chopNegativeDensities( DomainPartition & domain ); + /** * @brief Getter for the number of fluid components (species) * @return the number of components @@ -234,6 +242,39 @@ class CompositionalMultiphaseBase : public FlowSolverBase DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const = 0; + + /** + * @brief Apply source flux boundary conditions to the system + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void applySourceFluxBC( real64 const time, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + + /** + * @brief Function to perform the Application of Dirichlet type BC's + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void applyDirichletBC( real64 const time, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + /**@}*/ struct viewKeyStruct : FlowSolverBase::viewKeyStruct @@ -297,22 +338,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override; -protected: - - virtual void postInputInitialization() override; - - virtual void initializePreSubGroups() override; - - virtual void initializePostInitialConditionsPreSubGroups() override; - - /** - * @brief Utility function that checks the consistency of the constitutive models - * @param[in] domain the domain partition - * This function will produce an error if one of the constitutive models - * (fluid, relperm) is incompatible with the reference fluid model. - */ - void validateConstitutiveModels( DomainPartition const & domain ) const; - /** * @brief Initialize all variables from initial conditions * @param domain the domain containing the mesh and fields @@ -323,50 +348,50 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - /** * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; /** - * @brief Initialize the aquifer boundary condition (gravity vector, water phase index) - * @param[in] cm reference to the global constitutive model manager + * @brief Function to fix the initial state during the initialization step in coupled problems + * @param[in] time current time + * @param[in] dt time step + * @param[in] dofManager degree-of-freedom manager associated with the linear system + * @param[in] domain the domain + * @param[in] localMatrix local system matrix + * @param[in] localRhs local system right-hand side vector + * @detail This function is meant to be called when the flag m_keepVariablesConstantDuringInitStep is on + * The main use case is the initialization step in coupled problems during which we solve an elastic problem for a fixed pressure */ - void initializeAquiferBC( constitutive::ConstitutiveManager const & cm ) const; + void keepVariablesConstantDuringInitStep( real64 const time, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + +protected: + + virtual void postInputInitialization() override; + + virtual void initializePreSubGroups() override; /** - * @brief Function to perform the Application of Dirichlet type BC's - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector + * @brief Utility function that checks the consistency of the constitutive models + * @param[in] domain the domain partition + * This function will produce an error if one of the constitutive models + * (fluid, relperm) is incompatible with the reference fluid model. */ - void applyDirichletBC( real64 const time, - real64 const dt, - DofManager const & dofManager, - DomainPartition & domain, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; + void validateConstitutiveModels( DomainPartition const & domain ) const; + + virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** - * @brief Apply source flux boundary conditions to the system - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector + * @brief Initialize the aquifer boundary condition (gravity vector, water phase index) + * @param[in] cm reference to the global constitutive model manager */ - void applySourceFluxBC( real64 const time, - real64 const dt, - DofManager const & dofManager, - DomainPartition & domain, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; + void initializeAquiferBC( constitutive::ConstitutiveManager const & cm ) const; /** * @brief Apply aquifer boundary conditions to the system @@ -384,30 +409,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const = 0; - /** - * @brief Function to fix the initial state during the initialization step in coupled problems - * @param[in] time current time - * @param[in] dt time step - * @param[in] dofManager degree-of-freedom manager associated with the linear system - * @param[in] domain the domain - * @param[in] localMatrix local system matrix - * @param[in] localRhs local system right-hand side vector - * @detail This function is meant to be called when the flag m_keepVariablesConstantDuringInitStep is on - * The main use case is the initialization step in coupled problems during which we solve an elastic problem for a fixed pressure - */ - void keepVariablesConstantDuringInitStep( real64 const time, - real64 const dt, - DofManager const & dofManager, - DomainPartition & domain, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; - - /** - * @brief Sets all the negative component densities (if any) to zero. - * @param domain the physical domain object - */ - void chopNegativeDensities( DomainPartition & domain ); - void chopNegativeDensities( ElementSubRegionBase & subRegion ); /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp index 4a54209df59..e6f0c8f5bf5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp @@ -75,6 +75,8 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase virtual void registerDataOnMesh( Group & MeshBodies ) override; + virtual void initializePostInitialConditionsPreSubGroups() override; + /** * @defgroup Solver Interface Functions * @@ -165,8 +167,6 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase virtual void initializePreSubGroups() override; - virtual void initializePostInitialConditionsPreSubGroups() override; - /// precompute the minGravityCoefficient for the buoyancy term void precomputeData( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index fb205f475e3..9fee04b2672 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -70,6 +70,8 @@ class FlowSolverBase : public PhysicsSolverBase virtual void registerDataOnMesh( Group & MeshBodies ) override; + virtual void initializePostInitialConditionsPreSubGroups() override; + struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct { // misc inputs @@ -146,6 +148,8 @@ class FlowSolverBase : public PhysicsSolverBase real64 const & timeAtBeginningOfStep, real64 const & dt ); + virtual void initializeFluid( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } + protected: /** @@ -196,16 +200,12 @@ class FlowSolverBase : public PhysicsSolverBase void initializeHydraulicAperture( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); - virtual void initializeFluid( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } - virtual void initializeThermal( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } void saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); virtual void initializePreSubGroups() override; - virtual void initializePostInitialConditionsPreSubGroups() override; - virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; /** diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index e4af3f705cc..4cd41cb7626 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -245,6 +245,8 @@ class SinglePhaseBase : public FlowSolverBase */ void updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const; + virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + protected: virtual void initializePreSubGroups() override; @@ -256,8 +258,6 @@ class SinglePhaseBase : public FlowSolverBase */ virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; - virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** From 62b27fcae956c8ee443715610adde7f06634a2b0 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 15 Nov 2024 08:34:39 -0600 Subject: [PATCH 10/27] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 00502a919e4..4d81ecd0048 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3339-8707-7c55c70 + baseline: integratedTests/baseline_integratedTests-pr3393-8760-b704634 allow_fail: all: '' streak: '' From 999f0b498e13e5864f3251636557b6b2002390f3 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 15 Nov 2024 08:35:17 -0600 Subject: [PATCH 11/27] Update BASELINE_NOTES.md --- BASELINE_NOTES.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index f4dcac6a76e..b79e6920179 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,9 @@ 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-11-15) +===================== +Fix netToGross bug. PR #3339 (2024-11-14) ===================== From e7ca59fb0c544079509a8ad817912c256a144650 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 15 Nov 2024 08:53:56 -0600 Subject: [PATCH 12/27] cuda --- .../fluidFlow/CompositionalMultiphaseBase.hpp | 4 +- .../fluidFlow/FlowSolverBase.hpp | 52 ++++----- .../fluidFlow/SinglePhaseBase.hpp | 106 +++++++++--------- 3 files changed, 81 insertions(+), 81 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 6dd01a081f3..02a73e178b9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -174,6 +174,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ void chopNegativeDensities( DomainPartition & domain ); + void chopNegativeDensities( ElementSubRegionBase & subRegion ); + /** * @brief Getter for the number of fluid components (species) * @return the number of components @@ -409,8 +411,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const = 0; - void chopNegativeDensities( ElementSubRegionBase & subRegion ); - /** * @brief Utility function that encapsulates the call to FieldSpecificationBase::applyFieldValue in BC application * @param[in] time_n the time at the beginning of the step diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 9fee04b2672..0888252bf9c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -150,6 +150,32 @@ class FlowSolverBase : public PhysicsSolverBase virtual void initializeFluid( 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: /** @@ -180,18 +206,6 @@ class FlowSolverBase : public PhysicsSolverBase virtual void precomputeData( MeshLevel & mesh, arrayView1d< string const > const & 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; - void initialize( DomainPartition & domain ); virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); } @@ -208,20 +222,6 @@ class FlowSolverBase : public PhysicsSolverBase virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; - /** - * @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; - /// the number of Degrees of Freedom per cell integer m_numDofPerCell; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index 4cd41cb7626..b68b35e9af2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -224,6 +224,23 @@ class SinglePhaseBase : public FlowSolverBase arrayView1d< real64 > const & localRhs, CRSMatrixView< real64, localIndex const > const & dR_dAper ) = 0; + /** + * @brief Apply source flux boundary conditions to the system + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void + applySourceFluxBC( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + struct viewKeyStruct : FlowSolverBase::viewKeyStruct { static constexpr char const * elemDofFieldString() { return "singlePhaseVariables"; } @@ -247,17 +264,49 @@ class SinglePhaseBase : public FlowSolverBase virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; -protected: - - virtual void initializePreSubGroups() override; + /** + * @brief Function to update fluid mass + * @param subRegion subregion that contains the fields + */ + void + updateMass( ElementSubRegionBase & subRegion ) const; - virtual void initializePostInitialConditionsPreSubGroups() override; + /** + * @brief Function to update energy + * @param subRegion subregion that contains the fields + */ + void + updateEnergy( ElementSubRegionBase & subRegion ) const; /** * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; + /** + * @brief Function to fix the initial state during the initialization step in coupled problems + * @param[in] time current time + * @param[in] dt time step + * @param[in] dofManager degree-of-freedom manager associated with the linear system + * @param[in] domain the domain + * @param[in] localMatrix local system matrix + * @param[in] localRhs local system right-hand side vector + * @detail This function is meant to be called when the flag m_keepVariablesConstantDuringInitStep is on + * The main use case is the initialization step in coupled problems during which we solve an elastic problem for a fixed pressure + */ + void keepVariablesConstantDuringInitStep( real64 const time, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + +protected: + + virtual void initializePreSubGroups() override; + + virtual void initializePostInitialConditionsPreSubGroups() override; + virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; /** @@ -280,20 +329,6 @@ 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 thermal conductivity * @param subRegion the group storing the required fields @@ -330,23 +365,6 @@ class SinglePhaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const; - /** - * @brief Apply source flux boundary conditions to the system - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector - */ - void - applySourceFluxBC( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; - /** * @brief Apply aquifer boundary conditions to the system * @param time current time @@ -364,24 +382,6 @@ class SinglePhaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const = 0; - /** - * @brief Function to fix the initial state during the initialization step in coupled problems - * @param[in] time current time - * @param[in] dt time step - * @param[in] dofManager degree-of-freedom manager associated with the linear system - * @param[in] domain the domain - * @param[in] localMatrix local system matrix - * @param[in] localRhs local system right-hand side vector - * @detail This function is meant to be called when the flag m_keepVariablesConstantDuringInitStep is on - * The main use case is the initialization step in coupled problems during which we solve an elastic problem for a fixed pressure - */ - void keepVariablesConstantDuringInitStep( real64 const time, - real64 const dt, - DofManager const & dofManager, - DomainPartition & domain, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; - /** * @brief Structure holding views into fluid properties used by the base solver. */ From f1bbcfe782d5d3dec0706027cceb105de73b2728 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 21 Nov 2024 08:31:38 -0600 Subject: [PATCH 13/27] Update FlowSolverBase.cpp --- src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index bdb15b4e0a7..b6375094eb5 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -559,7 +559,7 @@ void FlowSolverBase::initializeHydraulicAperture( MeshLevel & mesh, const arrayV [&]( localIndex const, SurfaceElementRegion & region ) { - region.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) + region.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) { subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key()).setApplyDefaultValue( region.getDefaultAperture()); } ); } ); } From 7440d13fc47fe4747e7c38f53ef7309db64e75f3 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 2 Dec 2024 15:06:34 -0600 Subject: [PATCH 14/27] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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: '' From 74b5a7b7510c0accc690b3b7868fb420ee761aa7 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 2 Dec 2024 15:06:56 -0600 Subject: [PATCH 15/27] Update BASELINE_NOTES.md --- BASELINE_NOTES.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index a4c76d4f455..46193312093 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,7 +6,7 @@ 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 (XXXX-XX-XX) +PR #3393 (2024-12-2) ===================== Fix netToGross bug. From b95d65a5af869e5932c78bd13a9a1a18b3a3754a Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 2 Dec 2024 16:53:29 -0600 Subject: [PATCH 16/27] revert to minimize changes --- .../fluidFlow/CompositionalMultiphaseBase.hpp | 127 ++++++++-------- .../fluidFlow/CompositionalMultiphaseFVM.hpp | 3 +- .../CompositionalMultiphaseHybridFVM.hpp | 11 +- .../fluidFlow/FlowSolverBase.hpp | 8 +- .../fluidFlow/SinglePhaseBase.hpp | 139 +++++++++--------- .../fluidFlow/SinglePhaseFVM.hpp | 2 - .../fluidFlow/SinglePhaseHybridFVM.hpp | 2 - .../fluidFlow/SinglePhaseProppantBase.hpp | 4 +- 8 files changed, 147 insertions(+), 149 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 02a73e178b9..ae5f259c14c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -168,14 +168,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual void updateState( DomainPartition & domain ) override final; - /** - * @brief Sets all the negative component densities (if any) to zero. - * @param domain the physical domain object - */ - void chopNegativeDensities( DomainPartition & domain ); - - void chopNegativeDensities( ElementSubRegionBase & subRegion ); - /** * @brief Getter for the number of fluid components (species) * @return the number of components @@ -244,39 +236,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const = 0; - - /** - * @brief Apply source flux boundary conditions to the system - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector - */ - void applySourceFluxBC( real64 const time, - real64 const dt, - DofManager const & dofManager, - DomainPartition & domain, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; - - /** - * @brief Function to perform the Application of Dirichlet type BC's - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector - */ - void applyDirichletBC( real64 const time, - real64 const dt, - DofManager const & dofManager, - DomainPartition & domain, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; - /**@}*/ struct viewKeyStruct : FlowSolverBase::viewKeyStruct @@ -317,29 +276,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase }; - virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, - DomainPartition & domain ) override; - - /** - * @brief function to set the next time step size - * @param[in] currentDt the current time step size - * @param[in] domain the domain object - * @return the prescribed time step size - */ - real64 setNextDt( real64 const & currentDt, - DomainPartition & domain ) override; - - virtual real64 setNextDtBasedOnCFL( real64 const & currentDt, - DomainPartition & domain ) override; - - void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL ); - - integer useSimpleAccumulation() const { return m_useSimpleAccumulation; } - - integer useTotalMassEquation() const { return m_useTotalMassEquation; } - - virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override; - /** * @brief Initialize all variables from initial conditions * @param domain the domain containing the mesh and fields @@ -355,6 +291,38 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; + /** + * @brief Apply source flux boundary conditions to the system + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void applySourceFluxBC( real64 const time, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + + /** + * @brief Function to perform the Application of Dirichlet type BC's + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void applyDirichletBC( real64 const time, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + /** * @brief Function to fix the initial state during the initialization step in coupled problems * @param[in] time current time @@ -373,6 +341,37 @@ class CompositionalMultiphaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const; + /** + * @brief Sets all the negative component densities (if any) to zero. + * @param domain the physical domain object + */ + void chopNegativeDensities( DomainPartition & domain ); + + void chopNegativeDensities( ElementSubRegionBase & subRegion ); + + virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, + DomainPartition & domain ) override; + + /** + * @brief function to set the next time step size + * @param[in] currentDt the current time step size + * @param[in] domain the domain object + * @return the prescribed time step size + */ + real64 setNextDt( real64 const & currentDt, + DomainPartition & domain ) override; + + virtual real64 setNextDtBasedOnCFL( real64 const & currentDt, + DomainPartition & domain ) override; + + void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL ); + + integer useSimpleAccumulation() const { return m_useSimpleAccumulation; } + + integer useTotalMassEquation() const { return m_useTotalMassEquation; } + + virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override; + protected: virtual void postInputInitialization() override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp index 61e2a697f53..48efd363a96 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp @@ -191,7 +191,8 @@ class CompositionalMultiphaseFVM : public CompositionalMultiphaseBase virtual void postInputInitialization() override; - virtual void initializePreSubGroups() override; + virtual void + initializePreSubGroups() override; struct DBCParameters { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp index e6f0c8f5bf5..45bea1e810e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp @@ -75,8 +75,6 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase virtual void registerDataOnMesh( Group & MeshBodies ) override; - virtual void initializePostInitialConditionsPreSubGroups() override; - /** * @defgroup Solver Interface Functions * @@ -143,6 +141,9 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const override; + virtual void + updatePhaseMobility( ObjectManagerBase & dataGroup ) const override; + virtual void applyAquiferBC( real64 const time, real64 const dt, @@ -163,15 +164,15 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase static constexpr char const * faceDofFieldString() { return "faceCenteredVariables"; } }; -protected: + virtual void initializePostInitialConditionsPreSubGroups() override; virtual void initializePreSubGroups() override; +protected: + /// precompute the minGravityCoefficient for the buoyancy term void precomputeData( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const override; - private: /// tolerance used in the computation of the transmissibility matrix diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 0888252bf9c..225d2368f62 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -70,8 +70,6 @@ class FlowSolverBase : public PhysicsSolverBase virtual void registerDataOnMesh( Group & MeshBodies ) override; - virtual void initializePostInitialConditionsPreSubGroups() override; - struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct { // misc inputs @@ -206,6 +204,10 @@ class FlowSolverBase : public PhysicsSolverBase virtual void precomputeData( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); + virtual void initializePreSubGroups() override; + + virtual void initializePostInitialConditionsPreSubGroups() override; + void initialize( DomainPartition & domain ); virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); } @@ -218,8 +220,6 @@ class FlowSolverBase : public PhysicsSolverBase void saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); - virtual void initializePreSubGroups() override; - virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; /// the number of Degrees of Freedom per cell diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index b68b35e9af2..1ff207bfbbd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -224,6 +224,28 @@ class SinglePhaseBase : public FlowSolverBase arrayView1d< real64 > const & localRhs, CRSMatrixView< real64, localIndex const > const & dR_dAper ) = 0; + struct viewKeyStruct : FlowSolverBase::viewKeyStruct + { + static constexpr char const * elemDofFieldString() { return "singlePhaseVariables"; } + }; + + /** + * @brief Function to perform the Application of Dirichlet type BC's + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void + applyDirichletBC( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + /** * @brief Apply source flux boundary conditions to the system * @param time current time @@ -241,10 +263,22 @@ class SinglePhaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const; - struct viewKeyStruct : FlowSolverBase::viewKeyStruct - { - static constexpr char const * elemDofFieldString() { return "singlePhaseVariables"; } - }; + /** + * @brief Apply aquifer boundary conditions to the system + * @param time current time + * @param dt time step + * @param domain the domain + * @param dofManager degree-of-freedom manager associated with the linear system + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + virtual void + applyAquiferBC( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const = 0; virtual void updateState ( DomainPartition & domain ) override final; @@ -256,13 +290,13 @@ class SinglePhaseBase : public FlowSolverBase real64 updateFluidState( ElementSubRegionBase & subRegion ) const; + /** - * @brief Update all relevant solid internal energy models using current values of temperature - * @param dataGroup the group storing the required fields + * @brief Function to update all constitutive models + * @param dataGroup group that contains the fields */ - void updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const; - - virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + virtual void + updateFluidModel( ObjectManagerBase & dataGroup ) const; /** * @brief Function to update fluid mass @@ -278,6 +312,33 @@ class SinglePhaseBase : public FlowSolverBase 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 + */ + void updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const; + + /** + * @brief Update thermal conductivity + * @param subRegion the group storing the required fields + */ + void updateThermalConductivity( ElementSubRegionBase & subRegion ) const; + + /** + * @brief Function to update fluid mobility + * @param dataGroup group that contains the fields + */ + void + updateMobility( ObjectManagerBase & dataGroup ) const; + + virtual void initializePreSubGroups() override; + + virtual void initializePostInitialConditionsPreSubGroups() override; + + virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + + virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + /** * @brief Compute the hydrostatic equilibrium using the compositions and temperature input tables */ @@ -303,12 +364,6 @@ class SinglePhaseBase : public FlowSolverBase protected: - virtual void initializePreSubGroups() override; - - virtual void initializePostInitialConditionsPreSubGroups() override; - - virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - /** * @brief Checks constitutive models for consistency * @param[in] domain the domain partition @@ -322,66 +377,12 @@ class SinglePhaseBase : public FlowSolverBase virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; - /** - * @brief Function to update all constitutive models - * @param dataGroup group that contains the fields - */ - virtual void - updateFluidModel( ObjectManagerBase & dataGroup ) const; - - /** - * @brief Update thermal conductivity - * @param subRegion the group storing the required fields - */ - void updateThermalConductivity( ElementSubRegionBase & subRegion ) const; - - /** - * @brief Function to update fluid mobility - * @param dataGroup group that contains the fields - */ - void - updateMobility( ObjectManagerBase & dataGroup ) const; - /** * @brief Utility function to save the converged state * @param[in] subRegion the element subRegion */ virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override; - /** - * @brief Function to perform the Application of Dirichlet type BC's - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector - */ - void - applyDirichletBC( real64 const time_n, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; - - /** - * @brief Apply aquifer boundary conditions to the system - * @param time current time - * @param dt time step - * @param domain the domain - * @param dofManager degree-of-freedom manager associated with the linear system - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector - */ - virtual void - applyAquiferBC( real64 const time, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const = 0; - /** * @brief Structure holding views into fluid properties used by the base solver. */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp index 67d3149f3b1..b93ce7e6b33 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.hpp @@ -191,8 +191,6 @@ class SinglePhaseFVM : public BASE /**@}*/ -protected: - virtual void applyAquiferBC( real64 const time, real64 const dt, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp index 6e8bbcf23c9..f8427490c60 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.hpp @@ -151,8 +151,6 @@ class SinglePhaseHybridFVM : public SinglePhaseBase arrayView1d< real64 > const & localRhs, CRSMatrixView< real64, localIndex const > const & dR_dAper ) override final; -protected: - virtual void applyAquiferBC( real64 const time, real64 const dt, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp index 8b0ae10a990..ab07e60f9aa 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseProppantBase.hpp @@ -55,14 +55,14 @@ class SinglePhaseProppantBase : public SinglePhaseBase */ virtual ~SinglePhaseProppantBase(); + virtual void updateFluidModel( ObjectManagerBase & dataGroup ) const override; + virtual void updatePorosityAndPermeability( SurfaceElementSubRegion & subRegion ) const override; protected: virtual void validateConstitutiveModels( DomainPartition & domain ) const override; - virtual void updateFluidModel( ObjectManagerBase & dataGroup ) const override; - virtual FluidPropViews getFluidProperties( constitutive::ConstitutiveBase const & fluid ) const override; private: From aded90fabf2850cbae6a2cfca9dcebede6e8b7fe Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 2 Dec 2024 17:06:49 -0600 Subject: [PATCH 17/27] revert more --- .../fluidFlow/CompositionalMultiphaseBase.cpp | 6 +- .../fluidFlow/CompositionalMultiphaseBase.hpp | 60 ++++++++++--------- .../fluidFlow/FlowSolverBase.cpp | 4 +- .../fluidFlow/FlowSolverBase.hpp | 6 +- .../fluidFlow/SinglePhaseBase.cpp | 17 +++--- .../fluidFlow/SinglePhaseBase.hpp | 4 +- 6 files changed, 49 insertions(+), 48 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 1c4fb095738..5e59ba4bc2c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -830,8 +830,8 @@ real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & sub return maxDeltaPhaseVolFrac; } -void CompositionalMultiphaseBase::initializeFluid( MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) +void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) { GEOS_MARK_FUNCTION; @@ -954,7 +954,7 @@ void CompositionalMultiphaseBase::initializeFluid( MeshLevel & mesh, } ); } -void CompositionalMultiphaseBase::initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +void CompositionalMultiphaseBase::initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index ae5f259c14c..278d8d6b453 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -68,8 +68,6 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual void registerDataOnMesh( Group & meshBodies ) override; - virtual void initializePostInitialConditionsPreSubGroups() override; - /** * @defgroup Solver Interface Functions * @@ -284,13 +282,31 @@ class CompositionalMultiphaseBase : public FlowSolverBase * from prescribed intermediate values (i.e. global densities from global fractions) * and any applicable hydrostatic equilibration of the domain */ - virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + 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 */ virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) override; + /** + * @brief Function to perform the Application of Dirichlet type BC's + * @param time current time + * @param dt time step + * @param dofManager degree-of-freedom manager associated with the linear system + * @param domain the domain + * @param localMatrix local system matrix + * @param localRhs local system right-hand side vector + */ + void applyDirichletBC( real64 const time, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const; + /** * @brief Apply source flux boundary conditions to the system * @param time current time @@ -308,7 +324,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase arrayView1d< real64 > const & localRhs ) const; /** - * @brief Function to perform the Application of Dirichlet type BC's + * @brief Apply aquifer boundary conditions to the system * @param time current time * @param dt time step * @param dofManager degree-of-freedom manager associated with the linear system @@ -316,12 +332,12 @@ class CompositionalMultiphaseBase : public FlowSolverBase * @param localMatrix local system matrix * @param localRhs local system right-hand side vector */ - void applyDirichletBC( real64 const time, - real64 const dt, - DofManager const & dofManager, - DomainPartition & domain, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const; + virtual void applyAquiferBC( real64 const time, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) const = 0; /** * @brief Function to fix the initial state during the initialization step in coupled problems @@ -341,6 +357,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) const; + /** * @brief Sets all the negative component densities (if any) to zero. * @param domain the physical domain object @@ -352,6 +369,8 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual real64 setNextDtBasedOnStateChange( real64 const & currentDt, DomainPartition & domain ) override; + void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL ); + /** * @brief function to set the next time step size * @param[in] currentDt the current time step size @@ -364,7 +383,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual real64 setNextDtBasedOnCFL( real64 const & currentDt, DomainPartition & domain ) override; - void computeCFLNumbers( DomainPartition & domain, real64 const & dt, real64 & maxPhaseCFL, real64 & maxCompCFL ); + virtual void initializePostInitialConditionsPreSubGroups() override; integer useSimpleAccumulation() const { return m_useSimpleAccumulation; } @@ -378,6 +397,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase virtual void initializePreSubGroups() override; + /** * @brief Utility function that checks the consistency of the constitutive models * @param[in] domain the domain partition @@ -386,30 +406,12 @@ class CompositionalMultiphaseBase : public FlowSolverBase */ void validateConstitutiveModels( DomainPartition const & domain ) const; - virtual void initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - /** * @brief Initialize the aquifer boundary condition (gravity vector, water phase index) * @param[in] cm reference to the global constitutive model manager */ void initializeAquiferBC( constitutive::ConstitutiveManager const & cm ) const; - /** - * @brief Apply aquifer boundary conditions to the system - * @param time current time - * @param dt time step - * @param dofManager degree-of-freedom manager associated with the linear system - * @param domain the domain - * @param localMatrix local system matrix - * @param localRhs local system right-hand side vector - */ - virtual void applyAquiferBC( real64 const time, - real64 const dt, - DofManager const & dofManager, - DomainPartition & domain, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) const = 0; - /** * @brief Utility function that encapsulates the call to FieldSpecificationBase::applyFieldValue in BC application * @param[in] time_n the time at the beginning of the step diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index b6375094eb5..bfbbbf368cd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -500,7 +500,7 @@ void FlowSolverBase::initialize( DomainPartition & domain ) initializeHydraulicAperture( mesh, regionNames ); // Initialize primary variables from applied initial conditions - initializeFluid( mesh, regionNames ); + initializeFluidState( mesh, regionNames ); // Initialize the rock thermal quantities: conductivity and solid internal energy // Note: @@ -508,7 +508,7 @@ void FlowSolverBase::initialize( DomainPartition & domain ) // - This step depends on porosity and phaseVolFraction if( m_isThermal ) { - initializeThermal( mesh, regionNames ); + initializeThermalState( mesh, regionNames ); } // Save initial pressure and temperature fields diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 225d2368f62..1d46d170e37 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -146,7 +146,9 @@ class FlowSolverBase : public PhysicsSolverBase real64 const & timeAtBeginningOfStep, real64 const & dt ); - virtual void initializeFluid( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } + 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 @@ -216,8 +218,6 @@ class FlowSolverBase : public PhysicsSolverBase void initializeHydraulicAperture( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); - virtual void initializeThermal( MeshLevel & mesh, const arrayView1d< const string > & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } - void saveInitialPressureAndTemperature( MeshLevel & mesh, const arrayView1d< const string > & regionNames ); virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 142f338b5bb..aca05df65e3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -583,7 +583,7 @@ void SinglePhaseBase::computeHydrostaticEquilibrium( DomainPartition & domain ) } ); } -void SinglePhaseBase::initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +void SinglePhaseBase::initializeFluidState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, auto & subRegion ) @@ -599,7 +599,7 @@ void SinglePhaseBase::initializeFluid( MeshLevel & mesh, arrayView1d< string con } ); } -void SinglePhaseBase::initializeThermal( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) +void SinglePhaseBase::initializeThermalState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { mesh.getElemManager().forElementSubRegions< CellElementSubRegion, SurfaceElementSubRegion >( regionNames, [&]( localIndex const, auto & subRegion ) @@ -665,7 +665,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_ CoupledSolidBase const & porousSolid = getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::solidNamesString() ) ); porousSolid.saveConvergedState(); - saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes + saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes updatePorosityAndPermeability( subRegion ); updateFluidState( subRegion ); @@ -720,11 +720,11 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) ); if( m_keepVariablesConstantDuringInitStep ) { - porousSolid.ignoreConvergedState(); // newPorosity <- porosity_n + porousSolid.ignoreConvergedState(); // newPorosity <- porosity_n } else { - porousSolid.saveConvergedState(); // porosity_n <- porosity + porousSolid.saveConvergedState(); // porosity_n <- porosity } if( m_isThermal ) @@ -1079,8 +1079,7 @@ void SinglePhaseBase::applySourceFluxBC( real64 const time_n, // add the value to the mass balance equation globalIndex const massRowIndex = dofNumber[ei] - rankOffset; globalIndex const energyRowIndex = massRowIndex + 1; - real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor - // here! + real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor here! localRhs[massRowIndex] += rhsValue; massProd += rhsValue; //add the value to the energy balance equation if the flux is positive (i.e., it's a producer) @@ -1184,7 +1183,7 @@ void SinglePhaseBase::keepVariablesConstantDuringInitStep( real64 const time, rankOffset, localMatrix, rhsValue, - pres[ei], // freeze the current pressure value + pres[ei], // freeze the current pressure value pres[ei] ); localRhs[localRow] = rhsValue; @@ -1195,7 +1194,7 @@ void SinglePhaseBase::keepVariablesConstantDuringInitStep( real64 const time, rankOffset, localMatrix, rhsValue, - temp[ei], // freeze the current temperature value + temp[ei], // freeze the current temperature value temp[ei] ); localRhs[localRow + 1] = rhsValue; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index 1ff207bfbbd..ad175412dc1 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -335,9 +335,9 @@ class SinglePhaseBase : public FlowSolverBase virtual void initializePostInitialConditionsPreSubGroups() override; - virtual void initializeFluid( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; + virtual void initializeFluidState( MeshLevel & mesh, arrayView1d< string const > const & regionNames ) override; - virtual void initializeThermal( 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 From a75dcfa9be2490c33bf55937617291a893e748f8 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 3 Dec 2024 14:20:50 -0600 Subject: [PATCH 18/27] hope it works --- .../fluidFlow/FlowSolverBase.cpp | 2 +- .../fluidFlow/FlowSolverBase.hpp | 4 +- .../CoupledReservoirAndWellsBase.hpp | 2 - .../PoromechanicsInitialization.cpp | 83 ++++++++++++++----- .../PoromechanicsInitialization.hpp | 29 ++++--- .../multiphysics/PoromechanicsSolver.hpp | 7 +- .../SolidMechanicsLagrangianFEM.hpp | 2 +- .../SolidMechanicsStateReset.hpp | 2 +- 8 files changed, 87 insertions(+), 44 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index bfbbbf368cd..db576cbd0fa 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -92,7 +92,7 @@ FlowSolverBase::FlowSolverBase( string const & name, m_numDofPerCell( 0 ), m_isThermal( 0 ), m_keepVariablesConstantDuringInitStep( 0 ), - m_isFixedStressPoromechanicsUpdate( false ), + m_isFixedStressPoromechanicsUpdate( 0 ), m_isJumpStabilized( false ), m_isLaggingFractureStencilWeightsUpdate( 0 ) { diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 93daa5f6a3d..1d46d170e37 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -138,8 +138,6 @@ class FlowSolverBase : public PhysicsSolverBase void enableLaggingFractureStencilWeightsUpdate(){ m_isLaggingFractureStencilWeightsUpdate = 1; }; - void initialize( DomainPartition & domain ); - real64 sumAquiferFluxes( BoundaryStencil const & stencil, AquiferBoundaryCondition::KernelWrapper const & aquiferBCWrapper, ElementViewConst< arrayView1d< real64 const > > const & pres, @@ -212,6 +210,8 @@ 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 ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp index 3fc2baa27df..264525673e3 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp @@ -223,8 +223,6 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL } } - void initialize( DomainPartition & domain ) const { return reservoirSolver()->initialize( domain ); } - void assembleFluxTerms( real64 const dt, DomainPartition const & domain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index 5a659d4bd04..b43682e88c4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -22,7 +22,13 @@ #include "events/tasks/TasksManager.hpp" #include "physicsSolvers/PhysicsSolverManager.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsStatistics.hpp" -#include "physicsSolvers/multiphysics/PoromechanicsSolver.hpp" +#include "physicsSolvers/multiphysics/MultiphasePoromechanics.hpp" +#include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" +#include "physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp" +#include "physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp" +#include "physicsSolvers/multiphysics/HydrofractureSolver.hpp" +#include "physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp" +#include "physicsSolvers/multiphysics/CompositionalMultiphaseReservoirAndWells.hpp" #include "physicsSolvers/multiphysics/LogLevelsInfo.hpp" namespace geos @@ -30,15 +36,18 @@ namespace geos using namespace dataRepository; -PoromechanicsInitialization::PoromechanicsInitialization( const string & name, Group * const parent ) : +template< typename POROMECHANICS_SOLVER > +PoromechanicsInitialization< POROMECHANICS_SOLVER >:: +PoromechanicsInitialization( const string & name, + Group * const parent ): TaskBase( name, parent ), - m_solidMechanicsSolverName(), + m_poromechanicsSolverName(), m_solidMechanicsStatistics(), m_solidMechanicsStateResetTask( name, parent ) { enableLogLevelInput(); - registerWrapper( viewKeyStruct::solidMechanicsSolverNameString(), &m_solidMechanicsSolverName ). + registerWrapper( viewKeyStruct::poromechanicsSolverNameString(), &m_poromechanicsSolverName ). setRTTypeName( rtTypes::CustomTypes::groupNameRef ). setInputFlag( InputFlags::REQUIRED ). setDescription( "Name of the poromechanics solver" ); @@ -52,20 +61,25 @@ PoromechanicsInitialization::PoromechanicsInitialization( const string & name, G addLogLevel< logInfo::Initialization >(); } -PoromechanicsInitialization::~PoromechanicsInitialization() {} +template< typename POROMECHANICS_SOLVER > +PoromechanicsInitialization< POROMECHANICS_SOLVER >::~PoromechanicsInitialization() {} -void PoromechanicsInitialization::postInputInitialization() +template< typename POROMECHANICS_SOLVER > +void +PoromechanicsInitialization< POROMECHANICS_SOLVER >:: +postInputInitialization() { Group & problemManager = this->getGroupByPath( "/Problem" ); Group & physicsSolverManager = problemManager.getGroup( "Solvers" ); - GEOS_THROW_IF( !physicsSolverManager.hasGroup( m_solidMechanicsSolverName ), - GEOS_FMT( "{}: solid mechanics solver named {} not found", - getWrapperDataContext( viewKeyStruct::solidMechanicsSolverNameString() ), - m_solidMechanicsSolverName ), + GEOS_THROW_IF( !physicsSolverManager.hasGroup( m_poromechanicsSolverName ), + GEOS_FMT( "{}: {} solver named {} not found", + getWrapperDataContext( viewKeyStruct::poromechanicsSolverNameString() ), + POROMECHANICS_SOLVER::catalogName(), + m_poromechanicsSolverName ), InputError ); - m_solidMechanicsSolver = &physicsSolverManager.getGroup< SolidMechanicsLagrangianFEM >( m_solidMechanicsSolverName ); + m_poromechanicsSolver = &physicsSolverManager.getGroup< POROMECHANICS_SOLVER >( m_poromechanicsSolverName ); if( !m_solidMechanicsStatisticsName.empty()) { @@ -82,28 +96,40 @@ void PoromechanicsInitialization::postInputInitialization() } m_solidMechanicsStateResetTask.setLogLevel( getLogLevel()); - m_solidMechanicsStateResetTask.m_solidSolverName = m_solidMechanicsSolver->getName(); + m_solidMechanicsStateResetTask.m_solidSolverName = m_poromechanicsSolver->solidMechanicsSolver()->getName(); m_solidMechanicsStateResetTask.postInputInitialization(); } -bool PoromechanicsInitialization::execute( real64 const time_n, +template< typename POROMECHANICS_SOLVER > +bool +PoromechanicsInitialization< POROMECHANICS_SOLVER >:: +execute( real64 const time_n, real64 const dt, integer const cycleNumber, integer const eventCounter, real64 const eventProgress, DomainPartition & domain ) { - GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` is set to perform stress initialization during the next time step(s)", - getName(), time_n, m_solidMechanicsSolverName ) ); - m_solidMechanicsSolver->setStressInitialization( true ); + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, + GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` is set to perform stress initialization during the next time step(s)", + getName(), time_n, m_poromechanicsSolverName ) ); + m_poromechanicsSolver->setStressInitialization( true ); m_solidMechanicsStateResetTask.execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); - m_solidMechanicsSolver->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); + if constexpr ( std::is_same_v< POROMECHANICS_SOLVER, HydrofractureSolver<> > ) // special case + { + m_poromechanicsSolver->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); + } + else // default + { + m_poromechanicsSolver->solidMechanicsSolver()->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); + } - GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` has completed stress initialization", - getName(), time_n + dt, m_solidMechanicsSolverName ) ); - m_solidMechanicsSolver->setStressInitialization( false ); + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, + GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` has completed stress initialization", + getName(), time_n + dt, m_poromechanicsSolverName ) ); + m_poromechanicsSolver->setStressInitialization( false ); if( m_solidMechanicsStatistics != nullptr ) { @@ -118,7 +144,22 @@ bool PoromechanicsInitialization::execute( real64 const time_n, namespace { -REGISTER_CATALOG_ENTRY( TaskBase, PoromechanicsInitialization, string const &, Group * const ) +typedef PoromechanicsInitialization< MultiphasePoromechanics<> > MultiphasePoromechanicsInitialization; +typedef PoromechanicsInitialization< MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> > > MultiphaseReservoirPoromechanicsInitialization; +typedef PoromechanicsInitialization< SinglePhasePoromechanics<> > SinglePhasePoromechanicsInitialization; +typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFractures<> > SinglePhasePoromechanicsConformingFracturesInitialization; +typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFractures< SinglePhaseReservoirAndWells<> > > SinglePhaseReservoirPoromechanicsConformingFracturesInitialization; +typedef PoromechanicsInitialization< SinglePhasePoromechanicsEmbeddedFractures > SinglePhasePoromechanicsEmbeddedFracturesInitialization; +typedef PoromechanicsInitialization< SinglePhasePoromechanics< SinglePhaseReservoirAndWells<> > > SinglePhaseReservoirPoromechanicsInitialization; +typedef PoromechanicsInitialization< HydrofractureSolver< SinglePhasePoromechanics<> > > HydrofractureInitialization; +REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, MultiphaseReservoirPoromechanicsInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsConformingFracturesInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, SinglePhaseReservoirPoromechanicsConformingFracturesInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsEmbeddedFracturesInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, SinglePhaseReservoirPoromechanicsInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, HydrofractureInitialization, string const &, Group * const ) } } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp index dbfb51d96f8..58e41ecb355 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp @@ -21,11 +21,11 @@ #define SRC_CORECOMPONENTS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSINITIALIZATION_HPP_ #include "events/tasks/TaskBase.hpp" +#include "physicsSolvers/multiphysics/HydrofractureSolver.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsStateReset.hpp" namespace geos { -class SolidMechanicsLagrangianFEM; class SolidMechanicsStatistics; /** @@ -39,6 +39,7 @@ class SolidMechanicsStatistics; * (4) run a normal simulation * */ +template< typename POROMECHANICS_SOLVER > class PoromechanicsInitialization : public TaskBase { public: @@ -48,7 +49,8 @@ class PoromechanicsInitialization : public TaskBase * @param[in] name the name of the task coming from the xml * @param[in] parent the parent group of the task */ - PoromechanicsInitialization( const string & name, Group * const parent ); + PoromechanicsInitialization( const string & name, + Group * const parent ); /// Destructor for the class ~PoromechanicsInitialization() override; @@ -56,7 +58,14 @@ class PoromechanicsInitialization : public TaskBase /// Accessor for the catalog name static string catalogName() { - return "PoromechanicsInitialization"; + if constexpr ( std::is_same_v< POROMECHANICS_SOLVER, HydrofractureSolver<> > ) // special case + { + return "HydrofractureInitialization"; + } + else // default + { + return "PoromechanicsInitialization"; + } } /** @@ -82,24 +91,22 @@ class PoromechanicsInitialization : public TaskBase */ struct viewKeyStruct { - /// String for the solid mechanics solver name - constexpr static char const * solidMechanicsSolverNameString() { return "solidMechanicsSolverName"; } + /// String for the poromechanics solver name + constexpr static char const * poromechanicsSolverNameString() { return "poromechanicsSolverName"; } /// String for the solid mechanics statistics name constexpr static char const * solidMechanicsStatisticsNameString() { return "solidMechanicsStatisticsName"; } }; void postInputInitialization() override; -// void registerDataOnMesh( Group & meshBodies ) override; - - /// Name of the solid mechanics solver - string m_solidMechanicsSolverName; + /// Name of the poromechanics solver + string m_poromechanicsSolverName; /// Name of the solid mechanics statistics string m_solidMechanicsStatisticsName; - /// Pointer to the solid mechanics solver - SolidMechanicsLagrangianFEM * m_solidMechanicsSolver; + /// Pointer to the poromechanics solver + POROMECHANICS_SOLVER * m_poromechanicsSolver; /// Pointer to the solid mechanics statistics SolidMechanicsStatistics * m_solidMechanicsStatistics; diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 72084abf0bc..a28e5aca800 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -114,10 +114,6 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER GEOS_FMT( "{} {}: The attribute `{}` of the flow solver `{}` must be set to 1 since the poromechanics solver is thermal", this->getCatalogName(), this->getName(), FlowSolverBase::viewKeyStruct::isThermalString(), this->flowSolver()->getName() ), InputError ); - - DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" ); - flowSolver()->initialize( domain ); - updateBulkDensity( domain ); } virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override final @@ -298,9 +294,10 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER * @brief Utility function to set the stress initialization flag * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise */ - void setStressInitialization( integer const performStressInitialization ) + void setStressInitialization( bool const performStressInitialization ) { m_performStressInitialization = performStressInitialization; + solidMechanicsSolver()->setStressInitialization( performStressInitialization ); } struct viewKeyStruct : Base::viewKeyStruct diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index e2c3e133e0d..7bae91b1b98 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -283,7 +283,7 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase * @brief Utility function to set the stress initialization flag * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise */ - void setStressInitialization( integer const performStressInitialization ) + void setStressInitialization( bool const performStressInitialization ) { m_performStressInitialization = performStressInitialization; } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.hpp index ef3b4b03110..73d3772be1a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.hpp @@ -42,7 +42,7 @@ class SolidMechanicsLagrangianFEM; */ class SolidMechanicsStateReset : public TaskBase { - friend class PoromechanicsInitialization; + template< typename > friend class PoromechanicsInitialization; public: From a9e15a1b3d9ba9da7c66ff4f939759babb46e41f Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 3 Dec 2024 14:26:46 -0600 Subject: [PATCH 19/27] test --- .../physicsSolvers/fluidFlow/FlowSolverBase.cpp | 2 +- .../multiphysics/PoromechanicsInitialization.cpp | 10 ++++------ .../multiphysics/PoromechanicsInitialization.hpp | 10 +--------- 3 files changed, 6 insertions(+), 16 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index db576cbd0fa..bfbbbf368cd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -92,7 +92,7 @@ FlowSolverBase::FlowSolverBase( string const & name, m_numDofPerCell( 0 ), m_isThermal( 0 ), m_keepVariablesConstantDuringInitStep( 0 ), - m_isFixedStressPoromechanicsUpdate( 0 ), + m_isFixedStressPoromechanicsUpdate( false ), m_isJumpStabilized( false ), m_isLaggingFractureStencilWeightsUpdate( 0 ) { diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index b43682e88c4..da276dd2fb3 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -110,9 +110,8 @@ execute( real64 const time_n, real64 const eventProgress, DomainPartition & domain ) { - GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, - GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` is set to perform stress initialization during the next time step(s)", - getName(), time_n, m_poromechanicsSolverName ) ); + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` is set to perform stress initialization during the next time step(s)", + getName(), time_n, m_poromechanicsSolverName ) ); m_poromechanicsSolver->setStressInitialization( true ); m_solidMechanicsStateResetTask.execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); @@ -126,9 +125,8 @@ execute( real64 const time_n, m_poromechanicsSolver->solidMechanicsSolver()->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); } - GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, - GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` has completed stress initialization", - getName(), time_n + dt, m_poromechanicsSolverName ) ); + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` has completed stress initialization", + getName(), time_n + dt, m_poromechanicsSolverName ) ); m_poromechanicsSolver->setStressInitialization( false ); if( m_solidMechanicsStatistics != nullptr ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp index 58e41ecb355..c1c5acaf3dd 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp @@ -21,7 +21,6 @@ #define SRC_CORECOMPONENTS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSINITIALIZATION_HPP_ #include "events/tasks/TaskBase.hpp" -#include "physicsSolvers/multiphysics/HydrofractureSolver.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsStateReset.hpp" namespace geos @@ -58,14 +57,7 @@ class PoromechanicsInitialization : public TaskBase /// Accessor for the catalog name static string catalogName() { - if constexpr ( std::is_same_v< POROMECHANICS_SOLVER, HydrofractureSolver<> > ) // special case - { - return "HydrofractureInitialization"; - } - else // default - { - return "PoromechanicsInitialization"; - } + POROMECHANICS_SOLVER::catalogName() + "Initialization"; } /** From ea3467e6acf9e19ab183886ff5a686037bff5d5a Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 3 Dec 2024 14:27:20 -0600 Subject: [PATCH 20/27] merge fix --- .../physicsSolvers/multiphysics/PoromechanicsInitialization.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp index c1c5acaf3dd..e3d2a07c96d 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp @@ -57,7 +57,7 @@ class PoromechanicsInitialization : public TaskBase /// Accessor for the catalog name static string catalogName() { - POROMECHANICS_SOLVER::catalogName() + "Initialization"; + return POROMECHANICS_SOLVER::catalogName() + "Initialization"; } /** From e635c5842ef8620c9403afaceec59066fa230e39 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 3 Dec 2024 17:01:32 -0600 Subject: [PATCH 21/27] add back init --- .../physicsSolvers/fluidFlow/FlowSolverBase.hpp | 4 ++-- .../multiphysics/CoupledReservoirAndWellsBase.hpp | 2 ++ .../physicsSolvers/multiphysics/PoromechanicsSolver.hpp | 4 ++++ 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 23897ccb899..38765f57471 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -146,6 +146,8 @@ class FlowSolverBase : public PhysicsSolverBase real64 const & timeAtBeginningOfStep, real64 const & dt ); + void initializeState( DomainPartition & domain ); + 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 ); } @@ -210,8 +212,6 @@ class FlowSolverBase : public PhysicsSolverBase virtual void initializePostInitialConditionsPreSubGroups() override; - void initializeState( DomainPartition & domain ); - virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); } void initializePorosityAndPermeability( MeshLevel & mesh, arrayView1d< string const > const & regionNames ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp index 264525673e3..5e41d9c23b7 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp @@ -223,6 +223,8 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL } } + void initializeState( DomainPartition & domain ) const { return reservoirSolver()->initializeState( domain ); } + void assembleFluxTerms( real64 const dt, DomainPartition const & domain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index a28e5aca800..80c0f469ea3 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -114,6 +114,10 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER GEOS_FMT( "{} {}: The attribute `{}` of the flow solver `{}` must be set to 1 since the poromechanics solver is thermal", this->getCatalogName(), this->getName(), FlowSolverBase::viewKeyStruct::isThermalString(), this->flowSolver()->getName() ), InputError ); + + DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" ); + flowSolver()->initializeState( domain ); + updateBulkDensity( domain ); } virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override final From ee92b4b9ef48e81a999999b524879aa461a95127 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 3 Dec 2024 18:05:46 -0600 Subject: [PATCH 22/27] Update PoromechanicsSolver.hpp --- .../physicsSolvers/multiphysics/PoromechanicsSolver.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 80c0f469ea3..5d17eba1959 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -115,6 +115,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER this->getCatalogName(), this->getName(), FlowSolverBase::viewKeyStruct::isThermalString(), this->flowSolver()->getName() ), InputError ); + // the following is needed for proper geomechanics initialization DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" ); flowSolver()->initializeState( domain ); updateBulkDensity( domain ); From 55ec38a30e87cf4f0783c41a0c37215f65bae709 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Wed, 4 Dec 2024 12:28:09 -0600 Subject: [PATCH 23/27] better version --- .../multiphysics/MultiphasePoromechanics.hpp | 2 +- .../PoromechanicsInitialization.cpp | 2 + .../multiphysics/PoromechanicsSolver.hpp | 38 +++++++++---------- .../multiphysics/SinglePhasePoromechanics.hpp | 1 + 4 files changed, 21 insertions(+), 22 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp index 9251d86dac1..4e7fe989439 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp @@ -41,7 +41,7 @@ class MultiphasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHANI using Base::m_stabilizationType; using Base::m_stabilizationRegionNames; using Base::m_stabilizationMultiplier; - using Base::getLogLevel; + using Base::updateBulkDensity; /** * @brief main constructor for MultiphasePoromechanics Objects diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index da276dd2fb3..d3dfcaf14f0 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -122,6 +122,8 @@ execute( real64 const time_n, } else // default { + m_poromechanicsSolver->flowSolver()->initializeState( domain ); + m_poromechanicsSolver->updateBulkDensity( domain ); m_poromechanicsSolver->solidMechanicsSolver()->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 80c0f469ea3..08143425a53 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -114,10 +114,6 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER GEOS_FMT( "{} {}: The attribute `{}` of the flow solver `{}` must be set to 1 since the poromechanics solver is thermal", this->getCatalogName(), this->getName(), FlowSolverBase::viewKeyStruct::isThermalString(), this->flowSolver()->getName() ), InputError ); - - DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" ); - flowSolver()->initializeState( domain ); - updateBulkDensity( domain ); } virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override final @@ -388,6 +384,23 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER } + void updateBulkDensity( DomainPartition & domain ) + { + this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // update the bulk density + // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the + // pressure/temperature terms + updateBulkDensity( subRegion ); + } ); + } ); + } + protected: /* Implementation of Nonlinear Acceleration (Aitken) of averageMeanTotalStressIncrement */ @@ -596,23 +609,6 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER } ); } - void updateBulkDensity( DomainPartition & domain ) - { - this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - auto & subRegion ) - { - // update the bulk density - // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the - // pressure/temperature terms - updateBulkDensity( subRegion ); - } ); - } ); - } - virtual void updateBulkDensity( ElementSubRegionBase & subRegion ) = 0; virtual void validateNonlinearAcceleration() override diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp index f311b786c38..ecf7960827b 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp @@ -41,6 +41,7 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN using Base::m_stabilizationType; using Base::m_stabilizationRegionNames; using Base::m_stabilizationMultiplier; + using Base::updateBulkDensity; /** * @brief main constructor for SinglePhasePoromechanics objects From 104cbb72f1f00772ef4ef7f35f2c2534ab72a235 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 17 Feb 2025 10:24:40 -0600 Subject: [PATCH 24/27] Update PoromechanicsSolver.hpp --- .../physicsSolvers/multiphysics/PoromechanicsSolver.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 7aff2ac05ce..e4ba06121ad 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -388,7 +388,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER { this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) + string_array const & regionNames ) { mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, auto & subRegion ) From e2a7a5ae804be6ea4bbb9828864bb776a40c259c Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 6 Mar 2025 12:42:58 -0600 Subject: [PATCH 25/27] try unify hydrofrac --- .../PoromechanicsInitialization.cpp | 13 +++--------- .../SolidMechanicsLagrangianFEM.cpp | 21 ++----------------- 2 files changed, 5 insertions(+), 29 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index 20ba67e049e..d0145f47eb8 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -117,16 +117,9 @@ execute( real64 const time_n, m_solidMechanicsStateResetTask.execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); - if constexpr ( std::is_same_v< POROMECHANICS_SOLVER, HydrofractureSolver<> > ) // special case - { - m_poromechanicsSolver->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); - } - else // default - { - m_poromechanicsSolver->flowSolver()->initializeState( domain ); - m_poromechanicsSolver->updateBulkDensity( domain ); - m_poromechanicsSolver->solidMechanicsSolver()->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); - } + m_poromechanicsSolver->flowSolver()->initializeState( domain ); + m_poromechanicsSolver->updateBulkDensity( domain ); + m_poromechanicsSolver->solidMechanicsSolver()->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Initialization, GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` has completed stress initialization", getName(), time_n + dt, m_poromechanicsSolverName ) ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index f6b2bd57835..d5b672efeaa 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -1007,6 +1007,8 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain, ParallelVector & solution, bool const setSparsity ) { + GEOS_LOG( "SolidMechanicsLagrangianFEM::setupSystem" ); + GEOS_MARK_FUNCTION; PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); @@ -1022,25 +1024,6 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain, arrayView1d< globalIndex const > const dofNumber = nodeManager.getReference< globalIndex_array >( dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); - if( m_contactRelationName != viewKeyStruct::noContactRelationNameString() ) - { - ElementRegionManager const & elemManager = mesh.getElemManager(); - string_array allFaceElementRegions; - elemManager.forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion const & elemRegion ) - { - allFaceElementRegions.emplace_back( elemRegion.getName() ); - } ); - - finiteElement:: - fillSparsity< FaceElementSubRegion, - solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, - allFaceElementRegions, - this->getDiscretizationName(), - dofNumber, - dofManager.rankOffset(), - sparsityPattern ); - - } finiteElement:: fillSparsity< CellElementSubRegion, solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, From acf4c459c08739e88248f949c9187a32ac030d32 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 21 Mar 2025 17:54:21 -0500 Subject: [PATCH 26/27] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index c5de93442f2..f303a54d8d7 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr2125-10859-1db8623 + baseline: integratedTests/baseline_integratedTests-pr3396-10875-13b594f allow_fail: all: '' From 2c70d5ee33bc38cd92f0cf32570730e668b0dc13 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 21 Mar 2025 17:55:03 -0500 Subject: [PATCH 27/27] Update BASELINE_NOTES.md --- BASELINE_NOTES.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 26c0bdfc084..6a66008eb18 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3396 (2024-03-21) +===================== +Use solid mechanics solver directly to perform poromechanics initialization. + PR #2125 (2024-03-20) ===================== Phase-field nucleation model.