From a1fec52ffe4dfb1c24dab2dc3b93272f211d66f8 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Mon, 13 May 2024 19:56:08 -0700 Subject: [PATCH] Add the function to update the pressure gradient in the hybrid solver (#2110) --- .gitmodules | 2 +- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 + .../blaslapack/BlasLapackFunctions.h | 13 ++ .../interfaces/blaslapack/BlasLapackLA.cpp | 61 ++++++++ .../interfaces/blaslapack/BlasLapackLA.hpp | 17 +++ .../unitTests/testBlasLapack.cpp | 47 +++++++ .../fluidFlow/FlowSolverBaseFields.hpp | 8 ++ .../fluidFlow/SinglePhaseBase.cpp | 3 + .../fluidFlow/SinglePhaseBase.hpp | 8 ++ .../fluidFlow/SinglePhaseHybridFVM.cpp | 17 +++ .../fluidFlow/SinglePhaseHybridFVM.hpp | 2 + .../fluidFlow/SinglePhaseHybridFVMKernels.hpp | 130 ++++++++++++++++++ 13 files changed, 312 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index 029ba71a9ab..b3a951ebb50 100644 --- 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 index c2ca9d133b8..7756c489c70 100644 --- 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 index 9e9f73767e9..311616f1976 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 #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 index 5f61325f69c..fe66dcfc2ce 100644 --- 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 index 31fb564cf89..d05bbf8aef2 100644 --- 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 index dbf7eec9f6e..adc72007e49 100644 --- 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 index 8198473c2fc..d819072172a 100644 --- 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 index 64ca6c6a682..53f7ea8704f 100644 --- 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 index 2e0d864eb3f..2a34feff229 100644 --- 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 index 88fc2523de3..3ee103fae33 100644 --- 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 index c611af92cff..6215d942f77 100644 --- 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 index 7ba6c5d6dd6..65ceae5786e 100644 --- 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 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp index 4c6bc274f4c..f6c59f63df9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.hpp @@ -30,6 +30,8 @@ #include "finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp" #include "finiteVolume/mimeticInnerProducts/QuasiRTInnerProduct.hpp" #include "finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp" +#include "linearAlgebra/interfaces/InterfaceTypes.hpp" +#include "denseLinearAlgebra/interfaces/blaslapack/BlasLapackLA.hpp" #include "mesh/MeshLevel.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/HybridFVMHelperKernels.hpp" @@ -43,6 +45,7 @@ namespace geos namespace singlePhaseHybridFVMKernels { +/******************************** AssemblerKernelHelper ********************************/ /******************************** Kernel switches ********************************/ namespace internal @@ -81,6 +84,133 @@ void kernelLaunchSelectorFaceSwitch( T value, LAMBDA && lambda ) } // namespace internal +/******************************** AveragePressureGradientKernel ********************************/ + +template< integer NUM_FACES > +class AveragePressureGradientKernel +{ +public: + AveragePressureGradientKernel( CellElementSubRegion & subRegion, + FaceManager const & faceManager ) + : + // get the face-centered pressures + m_facePressure( faceManager.getField< fields::flow::facePressure >() ), + m_faceCenter( faceManager.faceCenter() ), + m_pres( subRegion.template getField< fields::flow::pressure >() ), + m_elemCenter( subRegion.getElementCenter() ), + m_elemsToFaces( subRegion.faceList() ), + m_presGradient( subRegion.template getField< fields::flow::pressureGradient >() ) + {} + + /** + * @struct StackVariables + * @brief Kernel variables (dof numbers, jacobian and residual) located on the stack + */ + struct StackVariables + { + + GEOS_HOST_DEVICE + StackVariables(): + coordinates( NUM_FACES+1, 4 ), + pressures( NUM_FACES+1 ), + presGradientLocal( 4 ) + {} + + stackArray2d< real64, (NUM_FACES + 1) * 4 > coordinates; + stackArray1d< real64, NUM_FACES + 1 > pressures; + stackArray1d< real64, 4 > presGradientLocal; + }; + + + inline + void compute( localIndex const elemIndex, + StackVariables stack ) const + { + + for( integer dim=0; dim<3; ++dim ) + { + stack.coordinates( 0, dim ) = m_elemCenter( elemIndex, dim ); + } + stack.coordinates( 0, 3 ) = 1.0; + stack.pressures[0] = m_pres[elemIndex]; + + for( integer fi=0; fi + static void + launch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + GEOS_MARK_FUNCTION; + + forAll< POLICY >( numElems, [=] ( localIndex const ei ) + { + typename KERNEL_TYPE::StackVariables stack; + kernelComponent.compute( ei, stack ); + } ); + } + +private: + /// face pressure + arrayView1d< real64 const > const m_facePressure; + + /// the face center coordinates + arrayView2d< real64 const > const m_faceCenter; + + /// the cell-centered pressures + arrayView1d< real64 const > const m_pres; + + /// the cell center coordinates + arrayView2d< real64 const > const m_elemCenter; + + /// the elements to faces map + arrayView2d< localIndex const > const m_elemsToFaces; + + /// pressure gradient in the cell + arrayView2d< real64 > const m_presGradient; + +}; + +class AveragePressureGradientKernelFactory +{ +public: + + template< typename POLICY > + static void createAndLaunch( CellElementSubRegion & subRegion, + FaceManager const & faceManager ) + { + internal::kernelLaunchSelectorFaceSwitch( subRegion.numFacesPerElement(), [&] ( auto NUM_FACES ) + { + AveragePressureGradientKernel< NUM_FACES > kernel( subRegion, faceManager ); + AveragePressureGradientKernel< NUM_FACES >::template launch< POLICY >( subRegion.size(), kernel ); + } ); + } +}; + + /******************************** ElementBasedAssemblyKernel ********************************/ /**