Add the function to update the pressure gradient in the hybrid solver (…

…#2110)
GEOS-DEV · May 14, 2024 · a1fec52 · a1fec52
1 parent 44684e1
commit a1fec52
Show file tree

Hide file tree

Showing 13 changed files with 312 additions and 2 deletions.
diff --git a/.gitmodules b/.gitmodules
@@ -12,4 +12,4 @@
 	url = ../../GEOS-DEV/hdf5_interface.git
 [submodule "scripts/uberenv"]
 	path = scripts/uberenv
-	url = ../../LLNL/uberenv.git
+	url = git@github.com:LLNL/uberenv.git
diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
@@ -1,7 +1,7 @@
 ---
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3060-4956-54ac336
+  baseline: integratedTests/baseline_integratedTests-pr2110-4992-9512556
 
 allow_fail:
   all: ''

diff --git 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 #2110 (2024-05-13)
+=====================
+new field to store pressure gradient cell-wise.
+
 PR #3060 (2024-05-13)
 ======================
 

diff --git a/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackFunctions.h b/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackFunctions.h
@@ -162,6 +162,19 @@ void GEOS_dgetrs( char const * TRANS,
                   int const * LDB,
                   int * INFO );
 
+#define GEOS_dgels FORTRAN_MANGLE( dgels )
+void GEOS_dgels( char const * TRANS,
+                 int const * M,
+                 int const * N,
+                 int const * NRHS,
+                 double * A,
+                 int const * LDA,
+                 double * B,
+                 int const * LDB,
+                 double * WORK,
+                 int const * LWORK,
+                 int * INFO );
+
 }
 #endif
 

diff --git a/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.cpp b/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.cpp
@@ -349,6 +349,36 @@ void matrixInverse( arraySlice2d< real64 const, USD > const & A,
   matrixInverse( A, Ainv, detA );
 }
 
+template< int USD >
+void matrixLeastSquaresSolutionSolve( arraySlice2d< real64, USD > const & A,
+                                      arraySlice1d< real64 > const & B,
+                                      arraySlice1d< real64 > const & X )
+{
+  GEOS_ASSERT_MSG( A.size( 1 ) == X.size() && A.size( 0 ) == B.size(),
+                   "Matrix, unknown vector and rhs vector not compatible" );
+
+  GEOS_ASSERT_MSG( X.size() <= B.size(),
+                   "Matrix, unknown vector and rhs vector not compatible" );
+
+  int const M = LvArray::integerConversion< int >( A.size( 0 ) );
+  int const N = LvArray::integerConversion< int >( A.size( 1 ) );
+  int const NRHS = 1;
+
+  int const LWORK = N + N;
+  array1d< double > WORK( LWORK );
+
+  int INFO = 0;
+
+  GEOS_dgels( "N", &M, &N, &NRHS, A.dataIfContiguous(), &M, B.dataIfContiguous(), &M, WORK.data(), &LWORK, &INFO );
+
+  for( int i = 0; i < N; ++i )
+  {
+    X[i] = B[i];
+  }
+
+  GEOS_ERROR_IF( INFO != 0, "The algorithm computing matrix linear system failed to converge." );
+}
+
 } // namespace detail
 
 real64 BlasLapackLA::determinant( arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & A )
@@ -911,6 +941,37 @@ void BlasLapackLA::solveLinearSystem( MatRowMajor< real64 const > const & A,
   solveLinearSystem( AT.toSliceConst(), rhs, solution );
 }
 
+void BlasLapackLA::matrixLeastSquaresSolutionSolve( arraySlice2d< real64 const, MatrixLayout::ROW_MAJOR > const & A,
+                                                    arraySlice1d< real64 const > const & B,
+                                                    arraySlice1d< real64 > const & X )
+{
+  array2d< real64, MatrixLayout::COL_MAJOR_PERM > AT( A.size( 0 ), A.size( 1 ) );
+
+  // convert A to a row major format
+  for( int i = 0; i < A.size( 0 ); ++i )
+  {
+    for( int j = 0; j < A.size( 1 ); ++j )
+    {
+      AT( i, j ) = A( i, j );
+    }
+  }
+
+  matrixLeastSquaresSolutionSolve( AT.toSliceConst(), B, X );
+}
+
+void BlasLapackLA::matrixLeastSquaresSolutionSolve( arraySlice2d< real64 const, MatrixLayout::COL_MAJOR > const & A,
+                                                    arraySlice1d< real64 const > const & B,
+                                                    arraySlice1d< real64 > const & X )
+{
+  // make a copy of A, since dgels modifies the components in A
+  array2d< real64, MatrixLayout::COL_MAJOR_PERM > ACOPY( A.size( 0 ), A.size( 1 ) );
+  BlasLapackLA::matrixCopy( A, ACOPY );
+
+  // make a copy of B, since dgels modifies the components in B
+  array1d< real64 > BCOPY( B.size() );
+  BlasLapackLA::vectorCopy( B, BCOPY );
 
+  detail::matrixLeastSquaresSolutionSolve( ACOPY.toSlice(), BCOPY.toSlice(), X );
+}
 
 } // end geosx namespace
diff --git a/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp b/src/coreComponents/denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp
@@ -596,6 +596,23 @@ struct BlasLapackLA
   static void matrixEigenvalues( MatColMajor< real64 const > const & A,
                                  Vec< std::complex< real64 > > const & lambda );
 
+  /**
+   * @brief Computes the least squares solution of B - AX
+   *
+   * @param [in]    A GEOSX array2d.
+   * @param [in]    B GEOSX array1d.
+   * @param [out]   X GEOSX array1d.
+   */
+  static void matrixLeastSquaresSolutionSolve( MatRowMajor< real64 const > const & A,
+                                               Vec< real64 const > const & B,
+                                               Vec< real64 > const & X );
+
+  /**
+   * @copydoc matrixLeastSquaresSolutionSolve
+   */
+  static void matrixLeastSquaresSolutionSolve( MatColMajor< real64 const > const & A,
+                                               Vec< real64 const > const & B,
+                                               Vec< real64 > const & X );
 };
 
 }

diff --git a/src/coreComponents/denseLinearAlgebra/unitTests/testBlasLapack.cpp b/src/coreComponents/denseLinearAlgebra/unitTests/testBlasLapack.cpp
@@ -1139,6 +1139,47 @@ void matrix_svd_test()
   }
 }
 
+template< typename LAI >
+void matrix_linear_system_least_square_solve_test()
+{
+  INDEX_TYPE M = 3;
+  INDEX_TYPE N = 2;
+
+  array2d< real64 > A( M, N );
+  array1d< real64 > B( M );
+  array1d< real64 > X( N );
+
+  array1d< real64 > vecResult( N );
+
+  // Assign component values to A
+  A( 0, 0 ) = 0.0;
+  A( 0, 1 ) = 1.0;
+  A( 1, 0 ) = 1.0;
+  A( 1, 1 ) = 1.0;
+  A( 2, 0 ) = 2.0;
+  A( 2, 1 ) = 1.0;
+
+  // Assign component values to B
+  B( 0 ) = 6.0;
+  B( 1 ) = 0.0;
+  B( 2 ) = 0.0;
+
+  // Compute X
+  LAI::matrixLeastSquaresSolutionSolve( A, B, X );
+
+  // Assign component values to the reference solution
+  vecResult( 0 ) = -3.0;
+  vecResult( 1 ) = 5.0;
+
+  // Check
+  for( INDEX_TYPE i = 0; i < N; ++i )
+  {
+    EXPECT_NEAR( vecResult( i ),
+                 X( i ),
+                 machinePrecision );
+  }
+}
+
 TEST( Array1D, vectorNorm1 )
 {
   vector_norm1_test< BlasLapackLA >();
@@ -1264,6 +1305,12 @@ TEST( DenseLAInterface, matrixSVD )
   matrix_svd_test< BlasLapackLA >();
 }
 
+TEST( Array2D, matrixLinearSystemSolve )
+{
+  matrix_linear_system_least_square_solve_test< BlasLapackLA >();
+}
+
+
 int main( int argc, char * * argv )
 {
   ::testing::InitGoogleTest( &argc, argv );

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp
@@ -88,6 +88,14 @@ DECLARE_FIELD( facePressure_n,
                NO_WRITE,
                "Face pressure at the previous converged time step" );
 
+DECLARE_FIELD( pressureGradient,
+               "pressureGradient",
+               array2d< real64 >,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Pressure gradient" );
+
 DECLARE_FIELD( temperature,
                "temperature",
                array1d< real64 >,

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -92,6 +92,9 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies )
         subRegion.registerField< dMobility_dTemperature >( getName() );
       }
 
+      subRegion.registerField< pressureGradient >( getName() ).
+        reference().resizeDimension< 1 >( 3 );
+
     } );
 
     FaceManager & faceManager = mesh.getFaceManager();

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
@@ -323,6 +323,14 @@ class SinglePhaseBase : public FlowSolverBase
    */
   void computeHydrostaticEquilibrium();
 
+  /**
+   * @brief Update the cell-wise pressure gradient
+   */
+  virtual void updatePressureGradient( DomainPartition & domain )
+  {
+    GEOS_UNUSED_VAR( domain );
+  }
+
   /**
    * @brief Function to fix the initial state during the initialization step in coupled problems
    * @param[in] time current time

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVM.cpp
@@ -622,5 +622,22 @@ 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
@@ -156,6 +156,8 @@ 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