wip: improving conditioning I

src/coreComponents/integrationTests/fluidFlowTests/CMakeLists.txt

@@ -4,7 +4,7 @@ set( gtest_geosx_tests
      testThermalSinglePhaseFlow.cpp
      testTransmissibility.cpp
      testImmiscibleMultiphaseFlow.cpp
-     testSinglePhaseMFDDiscretization.cpp )
+     testSinglePhaseMFDPolyhedral.cpp )
 
 if( ENABLE_PVTPackage )
     list( APPEND gtest_geosx_tests
@@ -31,6 +31,32 @@ foreach(test ${gtest_geosx_tests})
 
   geos_add_test( NAME ${test_name}
                  COMMAND ${test_name} )
+
+ set_tests_properties(${test_name} PROPERTIES ENVIRONMENT "TEST_BINARY_DIR=$<TARGET_FILE_DIR:${test_name}>")
+
+    # Add TEST_BINARY_DIR as a compile definition
+    target_compile_definitions(${test_name} PRIVATE TEST_BINARY_DIR=\"$<TARGET_FILE_DIR:${test_name}>\")
+
+    # --- Copy mesh file to test binary dir ---
+    add_custom_command(
+        TARGET ${test_name} POST_BUILD
+        COMMAND ${CMAKE_COMMAND} -E copy_if_different
+                ${CMAKE_CURRENT_SOURCE_DIR}/polyhedral_voronoi_regular.vtk
+                $<TARGET_FILE_DIR:${test_name}>
+    )
+    add_custom_command(
+        TARGET ${test_name} POST_BUILD
+        COMMAND ${CMAKE_COMMAND} -E copy_if_different
+                ${CMAKE_CURRENT_SOURCE_DIR}/polyhedral_voronoi_lattice.vtk
+                $<TARGET_FILE_DIR:${test_name}>
+    )
+    add_custom_command(
+        TARGET ${test_name} POST_BUILD
+        COMMAND ${CMAKE_COMMAND} -E copy_if_different
+                ${CMAKE_CURRENT_SOURCE_DIR}/polyhedral_voronoi_complex.vtk
+                $<TARGET_FILE_DIR:${test_name}>
+    )
+
 endforeach()
 
 # For some reason, BLT is not setting CUDA language for these source files

src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseMFDDiscretization.cpp → src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseMFDPolyhedral.cpp

@@ -12,7 +12,7 @@
  */
 
 #include <gtest/gtest.h>
-#include "unitTests/fluidFlowTests/testCompFlowUtils.hpp"
+#include "integrationTests/fluidFlowTests/testCompFlowUtils.hpp"
 #include "mainInterface/initialization.hpp"
 #include "mainInterface/ProblemManager.hpp"
 #include "mainInterface/GeosxState.hpp"
@@ -108,7 +108,7 @@ std::string generateXmlInputTPFA( std::string const & meshFile )
 
     <FieldSpecifications>
       <FieldSpecification name="initialPressure" initialCondition="1"
-        setNames="{ all }" objectPath="ElementRegions/Domain" fieldName="pressure" scale="1.0e7"/>    
+        setNames="{ all }" objectPath="ElementRegions/Domain" fieldName="pressure" scale="1.0e7"/>  
       <FieldSpecification name="west_pressure"
         setNames="{ westBC }" objectPath="faceManager" fieldName="pressure" scale="2.0e7"/>
       <FieldSpecification name="east_pressure"
@@ -157,68 +157,68 @@ class TPFAIntegrationTest : public ::testing::TestWithParam< const char * >
   std::string testBinaryDir;
 };
 
-INSTANTIATE_TEST_SUITE_P(
-  MeshFiles,
-  TPFAIntegrationTest,
-  ::testing::Values(
-    "polyhedral_voronoi_complex.vtk",
-    "polyhedral_voronoi_lattice.vtk",
-    "polyhedral_voronoi_regular.vtk"
-    )
-  );
-
-TEST_P( TPFAIntegrationTest, PressureFieldL2Error )
-{
-  ProblemManager & problemManager = state.getProblemManager();
-  DomainPartition & domain = problemManager.getDomainPartition();
-
-  // Retrieve the solver using the PhysicsSolverManager
-  SinglePhaseFVM< SinglePhaseBase > & solver =
-    dynamic_cast< SinglePhaseFVM< SinglePhaseBase > & >( problemManager.getPhysicsSolverManager().getGroup< SinglePhaseFVM< SinglePhaseBase > >( "SinglePhaseFlow" ) );
-
-  // Run the simulation to compute the numerical pressure
-  solver.setupSystem( domain, solver.getDofManager(), solver.getLocalMatrix(), solver.getSystemRhs(), solver.getSystemSolution() );
-  solver.implicitStepSetup( 0.0, MAX_TIME_STEP, domain );
-  solver.solverStep( 0.0, MAX_TIME_STEP, 0, domain );
-  solver.implicitStepComplete( 0.0, MAX_TIME_STEP, domain );
-
-  // Access the mesh and subregion
-  MeshLevel & mesh = domain.getMeshBody( 0 ).getBaseDiscretization();
-  CellElementSubRegion & subRegion = mesh.getElemManager().getRegion( 0 ).getSubRegion< CellElementSubRegion >( 0 );
-
-  // Retrieve pressure field and cell centers
-  arrayView2d< real64 const > centers = subRegion.getElementCenter();
-  arrayView1d< real64 const > volumes = subRegion.getElementVolume();
-  arrayView1d< real64 const > const p_h = subRegion.getField< fields::flow::pressure >();
-
-  // Compute exact pressure and L2 error
-  real64 l2Error = 0.0;
-  real64 totalVolume = 0.0;
-  for( localIndex i = 0; i < subRegion.size(); ++i )
-  {
-    real64 x = centers[i][0];
-    real64 volume = volumes[i];
-    real64 pNumeric = p_h[i] * to_MPA; // Convert pressure to MPa
-    real64 pExact = 20.0 * (1.0 - x) + 10.0 * x;
-    l2Error += (pNumeric - pExact) * (pNumeric - pExact) * volume;
-    totalVolume += volume;
-  }
-
-  l2Error = std::sqrt( l2Error / totalVolume );
-
-  std::string meshFile = GetParam();
-  if( meshFile == "polyhedral_voronoi_regular.vtk" )
-  {
-    // Assert that the L2 error is within machine precision
-    EXPECT_NEAR( l2Error, 0.0, PRESSURE_L2_TOLERANCE );
-  }
-  else
-  {
-    // Assert that the L2 error is not exact
-    EXPECT_GT( l2Error, PRESSURE_L2_TOLERANCE );
-  }
-
-}
+//INSTANTIATE_TEST_SUITE_P(
+//  MeshFiles,
+//  TPFAIntegrationTest,
+//  ::testing::Values(
+//    "polyhedral_voronoi_complex.vtk",
+//    "polyhedral_voronoi_lattice.vtk",
+//    "polyhedral_voronoi_regular.vtk"
+//    )
+//  );
+//
+//TEST_P( TPFAIntegrationTest, PressureFieldL2Error )
+//{
+//  ProblemManager & problemManager = state.getProblemManager();
+//  DomainPartition & domain = problemManager.getDomainPartition();
+//
+//  // Retrieve the solver using the PhysicsSolverManager
+//  SinglePhaseFVM< SinglePhaseBase > & solver =
+//    dynamic_cast< SinglePhaseFVM< SinglePhaseBase > & >( problemManager.getPhysicsSolverManager().getGroup< SinglePhaseFVM< SinglePhaseBase > >( "SinglePhaseFlow" ) );
+//
+//  // Run the simulation to compute the numerical pressure
+//  solver.setupSystem( domain, solver.getDofManager(), solver.getLocalMatrix(), solver.getSystemRhs(), solver.getSystemSolution() );
+//  solver.implicitStepSetup( 0.0, MAX_TIME_STEP, domain );
+//  solver.solverStep( 0.0, MAX_TIME_STEP, 0, domain );
+//  solver.implicitStepComplete( 0.0, MAX_TIME_STEP, domain );
+//
+//  // Access the mesh and subregion
+//  MeshLevel & mesh = domain.getMeshBody( 0 ).getBaseDiscretization();
+//  CellElementSubRegion & subRegion = mesh.getElemManager().getRegion( 0 ).getSubRegion< CellElementSubRegion >( 0 );
+//
+//  // Retrieve pressure field and cell centers
+//  arrayView2d< real64 const > centers = subRegion.getElementCenter();
+//  arrayView1d< real64 const > volumes = subRegion.getElementVolume();
+//  arrayView1d< real64 const > const p_h = subRegion.getField< fields::flow::pressure >();
+//
+//  // Compute exact pressure and L2 error
+//  real64 l2Error = 0.0;
+//  real64 totalVolume = 0.0;
+//  for( localIndex i = 0; i < subRegion.size(); ++i )
+//  {
+//    real64 x = centers[i][0];
+//    real64 volume = volumes[i];
+//    real64 pNumeric = p_h[i] * to_MPA; // Convert pressure to MPa
+//    real64 pExact = 20.0 * (1.0 - x) + 10.0 * x;
+//    l2Error += (pNumeric - pExact) * (pNumeric - pExact) * volume;
+//    totalVolume += volume;
+//  }
+//
+//  l2Error = std::sqrt( l2Error / totalVolume );
+//
+//  std::string meshFile = GetParam();
+//  if( meshFile == "polyhedral_voronoi_regular.vtk" )
+//  {
+//    // Assert that the L2 error is within machine precision
+//    EXPECT_NEAR( l2Error, 0.0, PRESSURE_L2_TOLERANCE );
+//  }
+//  else
+//  {
+//    // Assert that the L2 error is not exact
+//    EXPECT_GT( l2Error, PRESSURE_L2_TOLERANCE );
+//  }
+//
+//}
 
 std::string generateXmlInputMFD( std::string const & innerProductType,
                                  std::string const & meshFile )
@@ -267,7 +267,9 @@ std::string generateXmlInputMFD( std::string const & innerProductType,
 
   <FieldSpecifications>
     <FieldSpecification name="initialPressure" initialCondition="1"
-      setNames="{ all }" objectPath="ElementRegions/Domain" fieldName="pressure" scale="1.0e7"/>    
+      setNames="{ all }" objectPath="ElementRegions/Domain" fieldName="pressure" scale="1.0e7"/>
+    <FieldSpecification name="initialFacePressure" initialCondition="1"
+      setNames="{ all }" objectPath="faceManager" fieldName="facePressure" scale="1.0e7"/>          
     <FieldSpecification name="west_pressure"
       setNames="{ westBC }" objectPath="faceManager" fieldName="bcPressure" scale="2.0e7"/>
     <FieldSpecification name="east_pressure"
@@ -327,7 +329,7 @@ INSTANTIATE_TEST_SUITE_P(
   InnerProductAndMeshes,
   MFDIntegrationTest,
   ::testing::Combine(
-    ::testing::Values( TPFA, QuasiTPFA, QuasiRT, Simple, BdVLM ),
+    ::testing::Values( QuasiTPFA, QuasiRT, Simple, BdVLM ),
     ::testing::Values(
       "polyhedral_voronoi_complex.vtk",
       "polyhedral_voronoi_lattice.vtk",

src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp

@@ -61,7 +61,6 @@ SinglePhaseHybridFVM::SinglePhaseHybridFVM( const string & name,
 void SinglePhaseHybridFVM::registerDataOnMesh( Group & meshBodies )
 {
   SinglePhaseBase::registerDataOnMesh( meshBodies );
-
   forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
                                                     MeshLevel & mesh,
                                                     string_array const & regionNames )
@@ -83,9 +82,8 @@ void SinglePhaseHybridFVM::registerDataOnMesh( Group & meshBodies )
       // primary variables: face pressures at the previous converged time step
       faceManager.registerField< flow::facePressure_n >( getName() );
     }
-    // 3) Register the bc face data
+    // Register the bc face data
     {
-      // primary variables: face pressures at the previous converged time step
       faceManager.registerField< flow::bcPressure >( getName() );
     }
   } );
@@ -115,9 +113,7 @@ void SinglePhaseHybridFVM::initializePostInitialConditionsPreSubGroups()
   GEOS_MARK_FUNCTION;
 
   SinglePhaseBase::initializePostInitialConditionsPreSubGroups();
-
   DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
-
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
                                                                 MeshLevel & mesh,
                                                                 string_array const & regionNames )
@@ -396,14 +392,14 @@ void SinglePhaseHybridFVM::applyFaceDirichletBC( real64 const time_n,
 
       // next, we use the field specification functions to apply the boundary conditions to the system
 
-      // 1. first, populate the face pressure vector at the boundaries of the domain
+      // Populate the face pressure vector at the boundaries of the domain
       fs.applyFieldValue< FieldSpecificationEqual,
                           parallelDevicePolicy<> >( targetSet,
                                                     time_n + dt,
                                                     targetGroup,
                                                     flow::bcPressure::key() );
 
-      // 2. second, modify the residual/jacobian matrix as needed to impose the boundary conditions
+      // Second, modify the residual/jacobian matrix as needed to impose the boundary conditions
       forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
       {
 
@@ -413,12 +409,12 @@ void SinglePhaseHybridFVM::applyFaceDirichletBC( real64 const time_n,
           return;
         }
 
-        // 2.1 get the dof number of this face
+        // get the dof number of this face
         globalIndex const dofIndex = faceDofNumber[kf];
         localIndex const localRow = dofIndex - rankOffset;
         real64 rhsValue;
 
-        // 2.2 apply field value to the matrix/rhs
+        // apply field value to the matrix/rhs
         FieldSpecificationEqual::SpecifyFieldValue( dofIndex,
                                                     rankOffset,
                                                     localMatrix,
@@ -442,7 +438,6 @@ void SinglePhaseHybridFVM::applyAquiferBC( real64 const time,
                                            arrayView1d< real64 > const & localRhs ) const
 {
   GEOS_MARK_FUNCTION;
-
   GEOS_UNUSED_VAR( time, dt, dofManager, domain, localMatrix, localRhs );
 }
 
@@ -451,7 +446,6 @@ void SinglePhaseHybridFVM::saveAquiferConvergedState( real64 const & time,
                                                       DomainPartition & domain )
 {
   GEOS_MARK_FUNCTION;
-
   GEOS_UNUSED_VAR( time, dt, domain );
 }
 

src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernels.hpp

@@ -369,16 +369,16 @@ class ElementBasedAssemblyKernel
                         StackVariables & stack ) const
   {
     // local (cell-centered) mobility and its derivative w.r.t. pressure
-    real64 const localMobility = m_mob[m_er][m_esr][ei];
-    real64 const dLocalMobility_dPres = m_dMob[m_er][m_esr][ei][DerivOffset::dP];
+    real64 const massMobility = m_mob[m_er][m_esr][ei];
+    real64 const dmassMobility_dPres = m_dMob[m_er][m_esr][ei][DerivOffset::dP];
 
     for( integer iFaceLoc = 0; iFaceLoc < NUM_FACE; ++iFaceLoc )
     {
       // now in the following nested loop,
       // we compute the contribution of face jFaceLoc to the one sided total mass flux at face iFaceLoc
       for( integer jFaceLoc = 0; jFaceLoc < NUM_FACE; ++jFaceLoc )
       {
-        // 1) compute the potential diff between the cell center and the face center
+        // Compute the potential diff between the cell center and the face center
         real64 const ccPres = m_elemPres[ei];
         real64 const fPres = m_facePres[m_elemToFaces[ei][jFaceLoc]];
 
@@ -405,18 +405,14 @@ class ElementBasedAssemblyKernel
         real64 const dPotDif_dFacePres = dPresDif_dFacePres;
         real64 const T_ij = stack.transMatrix[iFaceLoc][jFaceLoc];
 
-        // massic factor: rho * lambda
-        real64 const massMobility = ccDens * localMobility;
-        real64 const dmassMobility_dPres = dCcDens_dPres * localMobility + ccDens * dLocalMobility_dPres;
-
-        // T * rho * lambda * (\nabla p - rho * g * \nabla d)
-        stack.massFlux[iFaceLoc] += T_ij * massMobility * potDif;
+        // T * lambda * (\nabla p - rho * g * \nabla d)
+        stack.massFlux[iFaceLoc] += m_dt * T_ij * massMobility * potDif;
 
         // derivatives w.r.t. element-centered pressure
-        stack.dmassFlux_dPres[iFaceLoc] += T_ij * ( dmassMobility_dPres * potDif + massMobility * dPotDif_dPres );
+        stack.dmassFlux_dPres[iFaceLoc] += m_dt * T_ij * ( dmassMobility_dPres * potDif + massMobility * dPotDif_dPres );
 
         // derivatives w.r.t. face-centered pressures
-        stack.dmassFlux_dFacePres[iFaceLoc][jFaceLoc] += T_ij * massMobility * dPotDif_dFacePres;
+        stack.dmassFlux_dFacePres[iFaceLoc][jFaceLoc] += m_dt * T_ij * massMobility * dPotDif_dFacePres;
       }
     }
   }
@@ -439,12 +435,12 @@ class ElementBasedAssemblyKernel
     for( integer iFaceLoc = 0; iFaceLoc < NUM_FACE; ++iFaceLoc )
     {
       // accumulate the mass flux divergence and its derivatives using the actual mass flux
-      stack.divMassFluxes += m_dt * stack.massFlux[iFaceLoc];
-      stack.dDivMassFluxes_dElemVars[0] += m_dt * stack.dmassFlux_dPres[iFaceLoc];
+      stack.divMassFluxes += stack.massFlux[iFaceLoc];
+      stack.dDivMassFluxes_dElemVars[0] += stack.dmassFlux_dPres[iFaceLoc];
 
       for( integer jFaceLoc = 0; jFaceLoc < NUM_FACE; ++jFaceLoc )
       {
-        stack.dDivMassFluxes_dFaceVars[jFaceLoc] += m_dt * stack.dmassFlux_dFacePres[iFaceLoc][jFaceLoc];
+        stack.dDivMassFluxes_dFaceVars[jFaceLoc] += stack.dmassFlux_dFacePres[iFaceLoc][jFaceLoc];
       }
 
       // collect the relevant dof numbers (always local)
@@ -470,7 +466,6 @@ class ElementBasedAssemblyKernel
     real64 const perm[ 3 ] = { m_elemPerm[ei][0][0], m_elemPerm[ei][0][1], m_elemPerm[ei][0][2] };
 
     // recompute the local transmissibility matrix at each iteration
-    // we can decide later to precompute transMatrix if needed
     IP::template compute< NUM_FACE >( m_nodePosition,
                                       m_transMultiplier,
                                       m_faceToNodes,
@@ -481,14 +476,11 @@ class ElementBasedAssemblyKernel
                                       m_lengthTolerance,
                                       stack.transMatrix );
 
-    // compute the one-sided mass fluxes and their derivatives
-    computeMassFlux( ei, stack );
-
-    // at this point, we know the local flow direction in the element
-    // so we can upwind the transport coefficients (mobilities) at the one sided faces
     if( m_elemGhostRank[ei] < 0 )
     {
-
+      // compute the one-sided mass fluxes and their derivatives
+      computeMassFlux( ei, stack );
+
       /*
        * perform assembly in this element in two steps:
        * 1) mass conservation equations