fix: detangle and unify flow initialization, fix netToGross bug

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
@@ -282,12 +321,20 @@ class CompositionalMultiphaseBase : public FlowSolverBase
    * from prescribed intermediate values (i.e. global densities from global fractions)
    * and any applicable hydrostatic equilibration of the domain
    */
-  void initializeFluidState( MeshLevel & mesh, DomainPartition & domain, arrayView1d< string const > const & regionNames );
+  virtual void 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
    */
-  void computeHydrostaticEquilibrium();
+  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
@@ -355,58 +402,13 @@ 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 );
 
-  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;
+  void chopNegativeDensities( ElementSubRegionBase & subRegion );
 
   /**
    * @brief Utility function that encapsulates the call to FieldSpecificationBase::applyFieldValue in BC application

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
   {

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

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp

@@ -286,16 +286,6 @@ void FlowSolverBase::saveSequentialIterationState( DomainPartition & domain )
   m_sequentialTempChange = m_isThermal ? MpiWrapper::max( maxTempChange ) : 0.0;
 }
 
-void FlowSolverBase::enableFixedStressPoromechanicsUpdate()
-{
-  m_isFixedStressPoromechanicsUpdate = true;
-}
-
-void FlowSolverBase::enableJumpStabilization()
-{
-  m_isJumpStabilized = true;
-}
-
 void FlowSolverBase::setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const
 {
   PhysicsSolverBase::setConstitutiveNamesCallSuper( subRegion );
@@ -455,6 +445,12 @@ void FlowSolverBase::initializePostInitialConditionsPreSubGroups()
                                                                 arrayView1d< string const > const & regionNames )
   {
     precomputeData( mesh, regionNames );
+
+    FieldIdentifiers fieldsToBeSync;
+    fieldsToBeSync.addElementFields( { fields::flow::pressure::key(), fields::flow::temperature::key() },
+                                     regionNames );
+
+    CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, mesh, domain.getNeighbors(), false );
   } );
 }
 
@@ -491,6 +487,97 @@ void FlowSolverBase::precomputeData( MeshLevel & mesh,
   }
 }
 
+void FlowSolverBase::initialize( DomainPartition & domain )
+{
+  // Compute hydrostatic equilibrium in the regions for which corresponding field specification tag has been specified
+  computeHydrostaticEquilibrium( domain );
+
+  forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+                                                               MeshLevel & mesh,
+                                                               arrayView1d< string const > const & regionNames )
+  {
+    initializePorosityAndPermeability( mesh, regionNames );
+    initializeHydraulicAperture( mesh, regionNames );
+
+    // Initialize primary variables from applied initial conditions
+    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 )
+  {
+    // Apply netToGross to reference porosity and horizontal permeability
+    CoupledSolidBase const & porousSolid =
+      getConstitutiveModel< CoupledSolidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ) );
+    PermeabilityBase const & permeability =
+      getConstitutiveModel< PermeabilityBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::permeabilityNamesString() ) );
+    arrayView1d< real64 const > const netToGross = subRegion.template getField< fields::flow::netToGross >();
+    porousSolid.scaleReferencePorosity( netToGross );
+    permeability.scaleHorizontalPermeability( netToGross );
+
+    // in some initializeState versions it uses newPorosity, so let's run updatePorosityAndPermeability to compute something
+    saveConvergedState( subRegion );   // necessary for a meaningful porosity update in sequential schemes
+    updatePorosityAndPermeability( subRegion );
+    porousSolid.initializeState();
+
+    // run final update
+    saveConvergedState( subRegion );   // necessary for a meaningful porosity update in sequential schemes
+    updatePorosityAndPermeability( subRegion );
+
+    // Save the computed porosity into the old porosity
+    // Note:
+    // - This must be called after updatePorosityAndPermeability
+    // - This step depends on porosity
+    porousSolid.saveConvergedState();
+  } );
+}
+
+void FlowSolverBase::initializeHydraulicAperture( MeshLevel & mesh, const arrayView1d< const string > & regionNames )
+{
+  mesh.getElemManager().forElementRegions< SurfaceElementRegion >( regionNames,
+                                                                   [&]( localIndex const,
+                                                                        SurfaceElementRegion & region )
+  {
+    region.forElementSubRegions< 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;