diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 3d9e32fee7a..084f4bfc9a9 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3849-14514-aaaf0f9 + baseline: integratedTests/baseline_integratedTests-pr3765-14543-498692f allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index e70070b3746..c14da10239c 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,10 +6,16 @@ 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 #3765 (2025-10-30) +===================== +Adding boundary fields to handle Dirichlet conditions in compositional MFD. + PR #3849 (2025-10-23) +===================== Add multiphase contact with wells PR #3880 (2025-10-13) +===================== Fix a bug introduced in #3485: mass that is used in accumulation term was not updated with porosity change after mechanics leading to always converged sequential outer loop. PR #3299 (2025-10-13) diff --git a/inputFiles/compositionalMultiphaseFlow/issue_3497/flowOnly_staircase_co2_3d_mfd.xml b/inputFiles/compositionalMultiphaseFlow/issue_3497/flowOnly_staircase_co2_3d_mfd.xml new file mode 100644 index 00000000000..4954c251290 --- /dev/null +++ b/inputFiles/compositionalMultiphaseFlow/issue_3497/flowOnly_staircase_co2_3d_mfd.xml @@ -0,0 +1,73 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/compositionalMultiphaseFlow/issue_3497/flowOnly_staircase_co2_3d_tpfa.xml b/inputFiles/compositionalMultiphaseFlow/issue_3497/flowOnly_staircase_co2_3d_tpfa.xml new file mode 100644 index 00000000000..a7038800200 --- /dev/null +++ b/inputFiles/compositionalMultiphaseFlow/issue_3497/flowOnly_staircase_co2_3d_tpfa.xml @@ -0,0 +1,73 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/cmake/GeosxOptions.cmake b/src/cmake/GeosxOptions.cmake index fe3adabcc4b..ebc45c3e929 100644 --- a/src/cmake/GeosxOptions.cmake +++ b/src/cmake/GeosxOptions.cmake @@ -147,7 +147,7 @@ message( "CMAKE_CXX_COMPILER_ID = ${CMAKE_CXX_COMPILER_ID}" ) blt_append_custom_compiler_flag( FLAGS_VAR CMAKE_CXX_FLAGS DEFAULT "${OpenMP_CXX_FLAGS}" ) blt_append_custom_compiler_flag( FLAGS_VAR CMAKE_CXX_FLAGS GNU "-Wpedantic -pedantic-errors -Wshadow -Wfloat-equal -Wcast-align -Wcast-qual" - CLANG "-Wpedantic -pedantic-errors -Wshadow -Wfloat-equal -Wno-cast-align -Wcast-qual" + CLANG "-Wpedantic -pedantic-errors -Wshadow -Wfloat-equal -Wno-cast-align -Wcast-qual -Wno-shorten-64-to-32" ) blt_append_custom_compiler_flag( FLAGS_VAR CMAKE_CXX_FLAGS_DEBUG diff --git a/src/coreComponents/finiteVolume/MimeticInnerProductDispatch.hpp b/src/coreComponents/finiteVolume/MimeticInnerProductDispatch.hpp index d9362edd0a3..a418182b715 100644 --- a/src/coreComponents/finiteVolume/MimeticInnerProductDispatch.hpp +++ b/src/coreComponents/finiteVolume/MimeticInnerProductDispatch.hpp @@ -142,10 +142,14 @@ mimeticInnerProductReducedDispatch( MimeticInnerProductBase const & input, { lambda( *ptr1 ); } - else if( auto const * const ptr2 = dynamic_cast< BdVLMInnerProduct const * >(&input) ) + else if( auto const * const ptr2 = dynamic_cast< QuasiTPFAInnerProduct const * >(&input) ) { lambda( *ptr2 ); } + else if( auto const * const ptr3 = dynamic_cast< BdVLMInnerProduct const * >(&input) ) + { + lambda( *ptr3 ); + } else { GEOS_ERROR( "mimeticInnerProductReducedDispatch() is not implemented for input of " << LvArray::system::demangleType( input ) ); @@ -169,10 +173,14 @@ mimeticInnerProductReducedDispatch( MimeticInnerProductBase & input, { lambda( *ptr1 ); } - else if( auto * const ptr2 = dynamic_cast< BdVLMInnerProduct * >(&input) ) + else if( auto * const ptr2 = dynamic_cast< QuasiTPFAInnerProduct * >(&input) ) { lambda( *ptr2 ); } + else if( auto * const ptr3 = dynamic_cast< BdVLMInnerProduct * >(&input) ) + { + lambda( *ptr3 ); + } else { GEOS_ERROR( "mimeticInnerProductReducedDispatch() is not supported for input of " << LvArray::system::demangleType( input ) ); diff --git a/src/coreComponents/integrationTests/fluidFlowTests/CMakeLists.txt b/src/coreComponents/integrationTests/fluidFlowTests/CMakeLists.txt index 129b10676c6..ce446028ce0 100644 --- a/src/coreComponents/integrationTests/fluidFlowTests/CMakeLists.txt +++ b/src/coreComponents/integrationTests/fluidFlowTests/CMakeLists.txt @@ -10,7 +10,8 @@ if( ENABLE_PVTPackage ) list( APPEND gtest_geosx_tests testCompMultiphaseFlow.cpp testCompMultiphaseFlowHybrid.cpp - testReactiveCompositionalMultiphaseOBL.cpp ) + testReactiveCompositionalMultiphaseOBL.cpp + testCompositionalMultiPhaseMFDPolyhedral.cpp ) endif() set( tplDependencyList ${parallelDeps} gtest ) @@ -56,6 +57,12 @@ foreach(test ${gtest_geosx_tests}) ${CMAKE_CURRENT_SOURCE_DIR}/polyhedral_voronoi_complex.vtu $ ) + add_custom_command( + TARGET ${test_name} POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy_if_different + ${CMAKE_CURRENT_SOURCE_DIR}/hex_pyr_tet_nested_mixed.vtu + $ + ) endforeach() diff --git a/src/coreComponents/integrationTests/fluidFlowTests/hex_pyr_tet_nested_mixed.vtu b/src/coreComponents/integrationTests/fluidFlowTests/hex_pyr_tet_nested_mixed.vtu new file mode 100644 index 00000000000..246b49aae24 --- /dev/null +++ b/src/coreComponents/integrationTests/fluidFlowTests/hex_pyr_tet_nested_mixed.vtu @@ -0,0 +1,38 @@ + + + + + + + + + AQAAAACAAAAkBAAAKgAAAA==eJxjZGBgYBwgzEQlTG3zaGHmQIUNsX5hxoIHKsxH0w/lYUNNDACQdQG2 + + + + + AQAAAACAAABYCwAA6wEAAA==eJxtlY1xQyEMg7NGd+gMz0yTObRGO0CnSXZqKafnz+DccSjBCP/IzuPBjy6j9+t14++vj0ibASxg3p/7wotn4cVjmwGs6N/WZTx348mTNgNY0futy3ieG889bQawoo95+jhuHuNlOyJtb/vo8/WfH/AIPAKPwNPlWkeeGCv95ZtnnXJlvczjswGseLQ1zvtZr8WTNgNY0esj3856rd/SZgArem2l31kvc4/I89s+el1mzFmv20fwCDydpjNf1DW1WfVFjbDOOnLMPDFW+nv2UnJwuV7mgX3Qd9aVd4xdL/PAPvoe1vbeiIxrRMZ120ff/9p8FXgEHoGnmx3a4qwzg33P3j3njkqOdi1Tj9TUObN01Ic5Zp4Y6znvdNzbl+tlnn5WaqvrOO67XhlXN2e1aULH264X8tPMaG16oq+24ZzmrOW8VNHiOW84M9j3+3+Dmh6q+aq6pjapLzW1rbmu+WbOGLeaN2ud6hKw52GdlexD1rjeF7A/nI8qPUx91Lf5f+h61RnN/qe2qt+c05y1auZc1WWNmfODM0BN/1VN13xR19SmGl3Ufqi5Zr6ZMzX+1l6qHLmez5/L64/w/v5+fbZnczeeNt2d+Zvx3Dsu8vEez82/v8cz8+9+8A594ZvkIh/v8ZzrFw53A8o= + + + 0 + + + 1.7320508076 + + + + + + + AQAAAACAAAAQMgAAagcAAA==eJyFl1d3kEUURbH33iv2CqiASlAh9i5JlARNghCwd0nQQDQJSrF3KRpASVDA3rvYe+9dUcRe/oEP3v2y17prfDnL5Jw9cyczfPf26vX/f8uFrhK6cugOoTuG7h66m3LLh64qDrmdQvuIQ26F0NXEIbdzaF9xyK0Yuro45HYJ7ScOuZVC1xCH3K6he4jDOVHvWqFrhnJO1LtX6J7KUe/a4pCj3v7ikKPedcQhR70DxCFHveuKQ456B4pDjnNbTxxynNve4nBO1LtB6PqhnBP17hu6j3LUu6E45Kh3kDjkqHcjcchRb4U45Kh3Y3HIUe9gcchxbpuIQ45z208czol6NwvdNJRzot4DQvdXjno3F4cc9Q4Rhxz1biEOOeodKg456t1SHHLUWykOOc5tK3HIcW4HisM5UW/v0K1DOSfqPTj0IOWodxtxyFHvIeKQo95txSFHvYeKQ456txOHXGXoYeKQ49y2F4cc53a4OKXvwxGhR4YOCz1Ouez7QO6o0CpxSt8HckeHVotT+j6QOya0RpzS94HcsaHHi1P6PnBO1Ds89ATlsu/D50m+NrRKPv97iW+4fKwLF06fxAcXzjD5yMFh3f6JL/uefZFw6kKr5fO/7/hq5fP+4PRNfD6/KvnIwWHdAYkv+/5+mXBGhNbI5+8Rvjr5vD84/RKfz69aPnJwWHdg4sv6BfK8oxPFK/ULvBvuYX3oScpl/cJXSb4htFY+fz/x1cvn9+H7bx9cOMPl8/1j3UGJL+tvvpb6PTXo9/y8MVm3Ivk9ee/b78I+n2utfL6XrFuR+LI+7ZuEMzJ0hHzuW/A1yuf9+V3Y5/Ork8/vmnUHJ76sryTP+zpZvFJfyXvifo4OHaVc1ld+m+SbQhvkc5+Fb7R8fje+n/bBhVMvn+8f6w5JfFkf/F3CGRPaKJ/7QnxN8nl/vv/2+fwa5PP9Y92hiS/r279POGNDR8rnPhbfGPm8P99/+3x+jfL5/bJuZeLL5gzyvKNTxCvNGaNCuYenhZ6qXDZn+P2cLk5pzvB9PEOc0pzhv++Z4pDL5gxynNtZ4pTmjLNDzwkdF3qhctmcQe7c0GZxSnMGufNCW8QpzRnkzg8dL05pziB3QehF4pTmDM6JeltDL1bOfTecH5L8hNBm+dwH4WuVj3XhwqlKfHDhjJOPHBzWrU187rvxLUk4E0Nb5ENdb51+X6ffw/O+4bO/JfJ5nWb5yMFhHe/H84HniR8TTlvoePncB+GbKJ/3B6cm8cH1+eMjB4d1RyS+bP4gz/u6RLzS/MF74n62h16qnPt6OD8l+Y7QCfK5D8LXLp/fjd+FfXDhtMqHug7vz/vyup5D/G48Z3QWcp2Fn7OP7B2a7z6c/SxNcpNC2+RzP4OvUz7fc9/jpYkfjvdnHj6/T/bD/pbKl80R5Hknl4lXmiN4F9yzyaGXK+e+Gs7PSX5KaId87qPwTZbP9x9OQ+KDC6ddPr8D1m1KfO6r8S1LOFNDO+Vz34dvinwo6/odLdP/+z357+Pzz/jsZ0zic9+O75eEMy10knzuX/FNlc/787uwz+faKZ/fO+uOTXzZfEGe93WFeKNCs/mC98T9vCr0SuWy+cLv6mpxSvOF7+k14pTmC/99rxWnNF+Q49yuE6c0X1wfekPo9NBblcvmC3I3hs4QpzRfkLspdKY4pfmC3M2hs8QpzRfkbgm9TZzSfME5UW9X6O3Kud+G82uSnx06Qz73Ofi65GNduHCaEx9cONPlIweHdSckPvfV+H5LOHNCZ8rnPgbfbPm8Pzgtic/nN0M+cnBYd2Lic3+O7/eEMzd0lnzuj/DNkc/7gzM+8fn8ZspHDg7rtiW+bI4gzzu6Q7zSHMG74R7OC71TOf5+HeL8keS7Q2fL5z4I3zz5/D58/+2DC6dLPt8/1u2Q70/5+X1n8nt+zt/b/J5Q7gF5c/F1y+d9+11k+/X54/O9ZF3X53fhOeOvhDM/dK587oPw9cjn/cFpS3w+vzny+V2z7qTEl80f5Hlfd4lXmj94T9zPBaF3K9cR6nni7yS/MLRbPvdZ+BbI53fjd2EfXDjz5CMHh3WnJD733fj+STiLQnvkc1+Ib6F83p/vv30+v275fP9Yd2ric/+O79+Ec0/ofPncx+JbJJ/35/tvn8+vRz6/X9adlviyOYM87+he8UpzBu+Ge3h/6H3KZXOG388D4pTmDN/HB8UpzRn++z4kTmnOIMe5PSxOac54JPTR0KdDn1IumzPIPRb6jDilOYPc46HPilOaM8g9EfqcOKU5g9yToc+LU5ozOCfqfSF0sXLuu+GQo94XxXFf7T6cHPW+JI771LnikKPel8Up9ZXkOLdXxCn1lYtDqfe10FeVcz8Fhxz1vi6O+yX3V+So9w1x3H+4XyFHvW+KQy7rF8hxbm+JU+oXOCfqfSf0beX8/YRDjnrfFcffR39PyVHve+L4e+PvEznqfV+c0veBHOf2gTil7wPnRL0fhX6oXPZ9IEe9H4tT+j6Qo95PxCl9H8hR76filL4P5Di3z8T5D0nt7Wg= + + + AQAAAACAAABICAAA2QEAAA==eJwtxVlMDgAAAOBKohjFkBq1YiVMNWcMaci1xdhoNmqGMHOMUQ/uTbaMbGpkq7YcMXK0lc2VB8cYmiHMis1iY2lzMx7+73v5ugcFRDra8U52qsd7qrOd48Ve7tXe4G3e6f0+5HJX+rTr3OCbvuvHfuFWt7vDP/zXIcGBu7mno9zfsY73UKd4lEd7gic7y9me5wXOdb5XeZ03equLvNP7XOwSl7rMFa7ySZ91net91Td8x4/81C1+43du9yd3+rv/ODgkcJh7ONL9HOM4D/Ewp3mCMz3L873EeS7wRm/3Lhf7oI+43Cdc7VM+54uud6Ov+7bv+qGb/dyv3eb3/uxv/u2gLoG6OsK93dcDHOs4JzrJw53qMc7wFGd5pud6oZc6zyu91hu8xYXe4b0udokP+6iPu9I1rvUFX3Gjb/men7jFbf7gL/7p4NDA4Y7yQA92opM90uke50nO9AzPcY4XOdfLvMIFXu/NLvIe73eJS13mCle5xmd83pfd4Gtu8h0/8BM/8yu3ut0d/upf/ufQroHD3ct9He1BTnCSRzjNYz3RUz3ds73Auc73Gm9yoXf7gEt9zNWu9SVfdZPvu9kv/dYf3elfDgkLHOE+jnGCU5zuDE/zf6h1gA4= + + + AQAAAACAAAAJAQAALgAAAA==eJzj4cEL+DAAXIyLC4sgikJMY5D0QAW5IAAqiWQOQilx7kF2DkH3YAEA/7ENdw== + + + + + diff --git a/src/coreComponents/integrationTests/fluidFlowTests/testCompositionalMultiPhaseMFDPolyhedral.cpp b/src/coreComponents/integrationTests/fluidFlowTests/testCompositionalMultiPhaseMFDPolyhedral.cpp new file mode 100644 index 00000000000..8dfd432d91a --- /dev/null +++ b/src/coreComponents/integrationTests/fluidFlowTests/testCompositionalMultiPhaseMFDPolyhedral.cpp @@ -0,0 +1,699 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * ------------------------------------------------------------------------------------------------------------ + */ + +#include +#include + +#include "integrationTests/fluidFlowTests/testCompFlowUtils.hpp" +#include "mainInterface/initialization.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mainInterface/GeosxState.hpp" +#include "mesh/DomainPartition.hpp" +#include "physicsSolvers/PhysicsSolverManager.hpp" + +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" + +using namespace geos; +using namespace geos::dataRepository; +using namespace geos::testing; + +CommandLineOptions g_commandLineOptions; + +// Pressure and saturation L2 error tolerances +static constexpr real64 PRESSURE_L2_TOLERANCE = 1.0e-8; +static constexpr real64 SATURATION_L2_TOLERANCE = 1.0e-8; + +// Single time step +static constexpr real64 TIME_STEP = 1.0e-2; // 0.01 + +static constexpr auto INNER_TPFA = "TPFA"; +static constexpr auto INNER_quasiTPFA = "quasiTPFA"; +static constexpr auto INNER_BDVLM = "beiraoDaVeigaLipnikovManzini"; + +// Generate XML for compositional multiphase FV-TPFA (simplified from the user-provided template) +static std::string generateXmlInputCompTPFA( std::string const & meshFile ) +{ + std::ostringstream oss; + oss << + R"xml( + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +)xml"; + return oss.str(); +} + +// Generate XML for compositional multiphase Hybrid MFD with configurable inner product (default TPFA) +static std::string generateXmlInputCompMFD( std::string const & innerProductType, + std::string const & meshFile ) +{ + std::ostringstream oss; + oss << + R"xml( + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +)xml"; + return oss.str(); +} + +// Helper: copy arrayView1d to std::vector +static inline std::vector< real64 > arrayViewToVector( arrayView1d< real64 const > & arr ) +{ + arr.move( hostMemorySpace, false ); + return std::vector< real64 >( arr.data(), arr.data() + arr.size() ); +} + +// Helper: copy a 2D view to flat vector (row-major by phase index), accepts any 2D view type +template< typename View2D > +static inline std::vector< real64 > arrayView2dToVector( View2D & arr ) +{ + arr.move( hostMemorySpace, false ); + localIndex nCells = arr.size( 0 ); + localIndex nCols = arr.size( 1 ); + std::vector< real64 > out; + out.reserve( static_cast< size_t >( nCells ) * static_cast< size_t >( nCols ) ); + for( localIndex i = 0; i < nCells; ++i ) + { + for( localIndex j = 0; j < nCols; ++j ) + { + out.push_back( static_cast< real64 >( arr( i, j ) ) ); + } + } + return out; +} + +// Helper: compute normalized L2 error for linear pressure profile p(x) = 2*(1-x) + x = 2 - x +static inline real64 computeNormalizedPressureL2Error( CellElementSubRegion & subRegion ) +{ + arrayView2d< real64 const > centers = subRegion.getElementCenter(); + arrayView1d< real64 const > volumes = subRegion.getElementVolume(); + arrayView1d< real64 const > const p_h = subRegion.getField< fields::flow::pressure >(); + + RAJA::ReduceSum< parallelDeviceReduce, real64 > squaredPressureError_ReduceSum( 0.0 ); + RAJA::ReduceSum< parallelDeviceReduce, real64 > squaredExactPressure_ReduceSum( 0.0 ); + localIndex const n_cells = subRegion.size(); + + forAll< geos::parallelDevicePolicy<> >( n_cells, [=] GEOS_HOST_DEVICE ( localIndex const i ) + { + real64 const x = centers[i][0]; + real64 const cellVolume = volumes[i]; + real64 const pNumeric = p_h[i]; + real64 const pExact = 2.0 * (1.0 - x) + 1.0 * x; + squaredPressureError_ReduceSum += (pNumeric - pExact) * (pNumeric - pExact) * cellVolume; + squaredExactPressure_ReduceSum += pExact * pExact * cellVolume; + } ); + + // Gather global reductions + real64 const squaredPressureError = squaredPressureError_ReduceSum.get(); + real64 const squaredExactPressure = squaredExactPressure_ReduceSum.get(); + + // Compute pressure relative L2 error + real64 const normalizedL2Error = std::sqrt( squaredPressureError ) / std::sqrt( squaredExactPressure ); + + return normalizedL2Error; +} + +// TPFA test: check linear pressure profile on regular polyhedral mesh +class CompositionalTPFAIntegrationTest : public ::testing::TestWithParam< const char * > +{ +public: + CompositionalTPFAIntegrationTest() + : state( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ) {} + +protected: + void SetUp() override + { + testBinaryDir = TEST_BINARY_DIR; + std::string meshFile = testBinaryDir + std::string( "/" ) + GetParam(); + std::string xmlInput = generateXmlInputCompTPFA( meshFile ); + setupProblemFromXML( state.getProblemManager(), xmlInput.c_str() ); + } + + GeosxState state; + std::string testBinaryDir; +}; + +INSTANTIATE_TEST_SUITE_P( + MeshFiles, + CompositionalTPFAIntegrationTest, + ::testing::Values( + "polyhedral_voronoi_regular.vtu", + "hex_pyr_tet_nested_mixed.vtu", + "polyhedral_voronoi_lattice.vtu", + "polyhedral_voronoi_complex.vtu" + ) + ); + +TEST_P( CompositionalTPFAIntegrationTest, PressureFieldL2Error ) +{ + ProblemManager & problemManager = state.getProblemManager(); + DomainPartition & domain = problemManager.getDomainPartition(); + + auto & solver = + dynamic_cast< CompositionalMultiphaseFVM & >( + problemManager.getPhysicsSolverManager().getGroup< CompositionalMultiphaseFVM >( "FlowSolverTPFA" ) ); + + // One implicit step with dt = 0.01 + solver.setupSystem( domain, solver.getDofManager(), + solver.getLocalMatrix(), solver.getSystemRhs(), + solver.getSystemSolution() ); + solver.implicitStepSetup( 0.0, TIME_STEP, domain ); + solver.solverStep( 0.0, TIME_STEP, 0, domain ); + solver.updateState( domain ); + solver.implicitStepComplete( 0.0, TIME_STEP, domain ); + + MeshLevel & mesh = domain.getMeshBody( 0 ).getBaseDiscretization(); + CellElementSubRegion & subRegion = mesh.getElemManager().getRegion( 0 ).getSubRegion< CellElementSubRegion >( 0 ); + + real64 l2Error = computeNormalizedPressureL2Error( subRegion ); + + // Expect exact solution only on the k-orthogonal regular mesh; expect non-zero error on mixed mesh + std::string meshFile = GetParam(); + if( meshFile == std::string( "polyhedral_voronoi_regular.vtu" ) ) + { + EXPECT_NEAR( l2Error, 0.0, PRESSURE_L2_TOLERANCE ); + } + else + { + EXPECT_GT( l2Error, PRESSURE_L2_TOLERANCE ); + } +} + +// MFD-TPFA test: check linear pressure profile on polyhedral meshes +class CompositionalMFDTPFAIntegrationTest : public ::testing::TestWithParam< const char * > +{ +public: + CompositionalMFDTPFAIntegrationTest() + : state( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ) {} + +protected: + void SetUp() override + { + testBinaryDir = TEST_BINARY_DIR; + std::string meshFile = testBinaryDir + std::string( "/" ) + GetParam(); + std::string xmlInput = generateXmlInputCompMFD( INNER_TPFA, meshFile ); + setupProblemFromXML( state.getProblemManager(), xmlInput.c_str() ); + } + + GeosxState state; + std::string testBinaryDir; +}; + +INSTANTIATE_TEST_SUITE_P( + MeshFiles, + CompositionalMFDTPFAIntegrationTest, + ::testing::Values( + "polyhedral_voronoi_regular.vtu", + "hex_pyr_tet_nested_mixed.vtu", + "polyhedral_voronoi_lattice.vtu", + "polyhedral_voronoi_complex.vtu" + ) + ); + +TEST_P( CompositionalMFDTPFAIntegrationTest, PressureFieldL2Error ) +{ + ProblemManager & problemManager = state.getProblemManager(); + DomainPartition & domain = problemManager.getDomainPartition(); + + auto & solver = + dynamic_cast< CompositionalMultiphaseHybridFVM & >( + problemManager.getPhysicsSolverManager().getGroup< CompositionalMultiphaseHybridFVM >( "FlowSolverMFD" ) ); + + // One implicit step with dt = 0.01 + solver.setupSystem( domain, solver.getDofManager(), + solver.getLocalMatrix(), solver.getSystemRhs(), + solver.getSystemSolution() ); + solver.implicitStepSetup( 0.0, TIME_STEP, domain ); + solver.solverStep( 0.0, TIME_STEP, 0, domain ); + solver.updateState( domain ); + solver.implicitStepComplete( 0.0, TIME_STEP, domain ); + + MeshLevel & mesh = domain.getMeshBody( 0 ).getBaseDiscretization(); + CellElementSubRegion & subRegion = mesh.getElemManager().getRegion( 0 ).getSubRegion< CellElementSubRegion >( 0 ); + + real64 l2Error = computeNormalizedPressureL2Error( subRegion ); + + // Expect exact solution only on the k-orthogonal regular mesh; expect non-zero error on mixed mesh + std::string meshFile = GetParam(); + if( meshFile == std::string( "polyhedral_voronoi_regular.vtu" ) ) + { + EXPECT_NEAR( l2Error, 0.0, PRESSURE_L2_TOLERANCE ); + } + else + { + EXPECT_GT( l2Error, PRESSURE_L2_TOLERANCE ); + } +} + +// Cross-check: TPFA vs MFD(TPFA): pressure and saturation profiles must match exactly +class CompositionalTPFAvsMFDTPFA : public ::testing::TestWithParam< const char * > {}; + +INSTANTIATE_TEST_SUITE_P( + MeshFiles, + CompositionalTPFAvsMFDTPFA, + ::testing::Values( + "polyhedral_voronoi_regular.vtu", + "hex_pyr_tet_nested_mixed.vtu", + "polyhedral_voronoi_lattice.vtu", + "polyhedral_voronoi_complex.vtu" + ) + ); + +TEST_P( CompositionalTPFAvsMFDTPFA, PressureAndSaturationComparison ) +{ + std::string const testBinaryDir = TEST_BINARY_DIR; + std::string const meshFile = testBinaryDir + std::string( "/" ) + GetParam(); + + std::vector< real64 > p_tpfa, p_mfd; + std::vector< real64 > sat_tpfa, sat_mfd; // phase volume fractions flattened [cell,phase] + localIndex n_cells_tpfa = 0, n_cells_mfd = 0; + localIndex n_phases_tpfa = 2, n_phases_mfd = 2; // From XML phaseNames + + // --- Run TPFA solver --- + { + GeosxState tpfaState( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ); + std::string xml = generateXmlInputCompTPFA( meshFile ); + setupProblemFromXML( tpfaState.getProblemManager(), xml.c_str() ); + + ProblemManager & pm = tpfaState.getProblemManager(); + DomainPartition & domain = pm.getDomainPartition(); + + auto & solver = + dynamic_cast< CompositionalMultiphaseFVM & >( + pm.getPhysicsSolverManager().getGroup< CompositionalMultiphaseFVM >( "FlowSolverTPFA" ) ); + + solver.setupSystem( domain, solver.getDofManager(), + solver.getLocalMatrix(), solver.getSystemRhs(), + solver.getSystemSolution() ); + solver.implicitStepSetup( 0.0, TIME_STEP, domain ); + solver.solverStep( 0.0, TIME_STEP, 0, domain ); + solver.updateState( domain ); + solver.implicitStepComplete( 0.0, TIME_STEP, domain ); + + MeshLevel & mesh = domain.getMeshBody( 0 ).getBaseDiscretization(); + CellElementSubRegion & subRegion = mesh.getElemManager().getRegion( 0 ).getSubRegion< CellElementSubRegion >( 0 ); + + arrayView1d< real64 const > p_h = subRegion.getField< fields::flow::pressure >(); + p_tpfa = arrayViewToVector( p_h ); + + auto sat = subRegion.getField< fields::flow::phaseVolumeFraction >(); + sat_tpfa = arrayView2dToVector( sat ); + } + + // --- Run MFD solver (innerProductType = TPFA) --- + { + GeosxState mfdState( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ); + std::string xml = generateXmlInputCompMFD( INNER_TPFA, meshFile ); + setupProblemFromXML( mfdState.getProblemManager(), xml.c_str() ); + + ProblemManager & pm = mfdState.getProblemManager(); + DomainPartition & domain = pm.getDomainPartition(); + + auto & solver = + dynamic_cast< CompositionalMultiphaseHybridFVM & >( + pm.getPhysicsSolverManager().getGroup< CompositionalMultiphaseHybridFVM >( "FlowSolverMFD" ) ); + + solver.setupSystem( domain, solver.getDofManager(), + solver.getLocalMatrix(), solver.getSystemRhs(), + solver.getSystemSolution() ); + solver.implicitStepSetup( 0.0, TIME_STEP, domain ); + solver.solverStep( 0.0, TIME_STEP, 0, domain ); + solver.updateState( domain ); + solver.implicitStepComplete( 0.0, TIME_STEP, domain ); + + MeshLevel & mesh = domain.getMeshBody( 0 ).getBaseDiscretization(); + CellElementSubRegion & subRegion = mesh.getElemManager().getRegion( 0 ).getSubRegion< CellElementSubRegion >( 0 ); + + arrayView1d< real64 const > p_h = subRegion.getField< fields::flow::pressure >(); + p_mfd = arrayViewToVector( p_h ); + + auto sat = subRegion.getField< fields::flow::phaseVolumeFraction >(); + sat_mfd = arrayView2dToVector( sat ); + } + + ASSERT_EQ( n_cells_tpfa, n_cells_mfd ); + ASSERT_EQ( n_phases_tpfa, n_phases_mfd ); + + // Pressure profile equality + for( localIndex i = 0; i < n_cells_tpfa; ++i ) + { + real64 const diff = std::abs( p_tpfa[i] - p_mfd[i] ); + EXPECT_NEAR( diff, 0.0, PRESSURE_L2_TOLERANCE ); + } + + // Saturation profile (phase volume fraction) equality + size_t const nTot = static_cast< size_t >( n_cells_tpfa ) * static_cast< size_t >( n_phases_tpfa ); + for( size_t k = 0; k < nTot; ++k ) + { + real64 const diff = std::abs( sat_tpfa[k] - sat_mfd[k] ); + EXPECT_NEAR( diff, 0.0, SATURATION_L2_TOLERANCE ); + } +} + +// MFD non-TPFA inner products: check linear pressure profile is exact on all meshes +class CompositionalMFDNonTPFAExactnessTest : public ::testing::TestWithParam< std::tuple< const char *, const char * > > +{ +public: + CompositionalMFDNonTPFAExactnessTest() + : state( std::make_unique< CommandLineOptions >( g_commandLineOptions ) ) {} + +protected: + void SetUp() override + { + testBinaryDir = TEST_BINARY_DIR; + auto const params = GetParam(); + innerProduct = std::get< 1 >( params ); + std::string const meshRel = std::get< 0 >( params ); + std::string const meshFile = testBinaryDir + std::string( "/" ) + meshRel; + std::string xmlInput = generateXmlInputCompMFD( innerProduct, meshFile ); + setupProblemFromXML( state.getProblemManager(), xmlInput.c_str() ); + } + + GeosxState state; + std::string testBinaryDir; + std::string innerProduct; +}; + +INSTANTIATE_TEST_SUITE_P( + MeshAndInnerProducts, + CompositionalMFDNonTPFAExactnessTest, + ::testing::Combine( + ::testing::Values( + "polyhedral_voronoi_regular.vtu", + "hex_pyr_tet_nested_mixed.vtu", + "polyhedral_voronoi_lattice.vtu", + "polyhedral_voronoi_complex.vtu" + ), + ::testing::Values( + INNER_quasiTPFA, + INNER_BDVLM + ) + ) + ); + +TEST_P( CompositionalMFDNonTPFAExactnessTest, PressureFieldL2ErrorExact ) +{ + ProblemManager & problemManager = state.getProblemManager(); + DomainPartition & domain = problemManager.getDomainPartition(); + + auto & solver = + dynamic_cast< CompositionalMultiphaseHybridFVM & >( + problemManager.getPhysicsSolverManager().getGroup< CompositionalMultiphaseHybridFVM >( "FlowSolverMFD" ) ); + + // One implicit step with dt = 0.01 + solver.setupSystem( domain, solver.getDofManager(), + solver.getLocalMatrix(), solver.getSystemRhs(), + solver.getSystemSolution() ); + solver.implicitStepSetup( 0.0, TIME_STEP, domain ); + solver.solverStep( 0.0, TIME_STEP, 0, domain ); + solver.updateState( domain ); + solver.implicitStepComplete( 0.0, TIME_STEP, domain ); + + MeshLevel & mesh = domain.getMeshBody( 0 ).getBaseDiscretization(); + CellElementSubRegion & subRegion = mesh.getElemManager().getRegion( 0 ).getSubRegion< CellElementSubRegion >( 0 ); + + real64 const normalizedL2Error = computeNormalizedPressureL2Error( subRegion ); + + // Expect exact solution for these inner products on both meshes + EXPECT_NEAR( normalizedL2Error, 0.0, PRESSURE_L2_TOLERANCE ); +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + g_commandLineOptions = *geos::basicSetup( argc, argv ); + int result = RUN_ALL_TESTS(); + geos::basicCleanup(); + return result; +} diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/CompositionalMultiphaseHybridFVM.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/CompositionalMultiphaseHybridFVM.hpp index 5f1b8a87f98..7081d021e9e 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/CompositionalMultiphaseHybridFVM.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/CompositionalMultiphaseHybridFVM.hpp @@ -81,10 +81,12 @@ class CompositionalMultiphaseHybridFVM : public MGRStrategyBase< 3 > // Level 1 m_levelFRelaxType[1] = MGRFRelaxationType::none; - m_levelInterpType[1] = MGRInterpolationType::jacobi; + m_levelFRelaxIters[1] = 1; + m_levelInterpType[1] = MGRInterpolationType::injection; m_levelRestrictType[1] = MGRRestrictionType::blockColLumped; // True-IMPES m_levelCoarseGridMethod[1] = MGRCoarseGridMethod::galerkin; - m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::none; + m_levelGlobalSmootherType[1] = MGRGlobalSmootherType::ilu0; + m_levelGlobalSmootherIters[1] = 1; // Level 2 m_levelFRelaxType[2] = MGRFRelaxationType::jacobi; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt index 28ffa57e587..a12f17fa1c2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt @@ -152,26 +152,19 @@ set( fluidFlowSolvers_sources wells/WellControls.cpp wells/WellSolverBase.cpp proppantTransport/ProppantTransport.cpp - proppantTransport/ProppantTransportKernels.cpp ) - -set( dependencyList ${parallelDeps} physicsSolversBase ) - -geos_decorate_link_dependencies( LIST decoratedDependencies DEPENDENCIES ${dependencyList} ) - - - + proppantTransport/ProppantTransportKernels.cpp + # Macro-based explicit instantiations for hybrid FVM kernels + kernels/compositional/CompositionalMultiphaseHybridFVMKernelsInstantiations.cpp + kernels/singlePhase/SinglePhaseHybridFVMKernelsInstantiations.cpp ) +# Remove the old generator-based hybrid flux/dirichlet kernels and switch to macro approach file( READ "${CMAKE_CURRENT_LIST_DIR}/kernelSpecs.json" kernelSpecs ) set( kernelTemplateFileList "" ) +# Keep only the templates that are unrelated to hybrid flux/dirichlet list( APPEND kernelTemplateFileList - SinglePhaseHybridFVMKernels.cpp.template - CompositionalMultiphaseHybridFVMKernels_upwinding.cpp.template - CompositionalMultiphaseHybridFVMKernels_assembly.cpp.template - CompositionalMultiphaseHybridFVMKernels_flux.cpp.template ReactiveCompositionalMultiphaseOBLKernels.cpp.template ) - foreach( kernelTemplateFile ${kernelTemplateFileList} ) get_filename_component( jsonKey ${kernelTemplateFile} NAME_WE ) generateKernels( TEMPLATE ${kernelTemplateFile} @@ -179,9 +172,8 @@ foreach( kernelTemplateFile ${kernelTemplateFileList} ) KEY ${jsonKey} HEADERS headerFiles SOURCES sourceFiles ) - -list(APPEND fluidFlowSolvers_headers ${headerFiles}) -list(APPEND fluidFlowSolvers_sources ${sourceFiles}) + list(APPEND fluidFlowSolvers_headers ${headerFiles}) + list(APPEND fluidFlowSolvers_sources ${sourceFiles}) endforeach() # TODO: The two kernels below have non-matching file names and JSON keys. @@ -215,6 +207,9 @@ blt_add_library( NAME fluidFlowSolvers SHARED ${GEOS_BUILD_SHARED_LIBS} ) +# Always enable the manual instantiation compile flag +target_compile_definitions( fluidFlowSolvers PUBLIC GEOS_ENABLE_MANUAL_HYBRID_FVM_INST ) + target_include_directories( fluidFlowSolvers PUBLIC ${CMAKE_SOURCE_DIR}/coreComponents ) install( TARGETS fluidFlowSolvers LIBRARY DESTINATION ${CMAKE_INSTALL_PREFIX}/lib ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp index 3f520747ae8..31120bcbd40 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp @@ -89,6 +89,14 @@ DECLARE_FIELD( globalCompFraction_n, // NO_WRITE, // "Global component fraction updates at the previous sequential iteration" ); +DECLARE_FIELD( bcGlobalCompFraction, + "bcGlobalCompFraction", + array2dLayoutComp, + 0, + NOPLOT, + WRITE_AND_READ, + "Boundary condition global component fraction" ); + DECLARE_FIELD( faceGlobalCompFraction, "faceGlobalCompFraction", array2dLayoutComp, @@ -137,6 +145,32 @@ DECLARE_FIELD( dPhaseMobility, NO_WRITE, "Derivative of phase volume fraction with respect to pressure, temperature, global component density" ); +// Face-based properties for boundary conditions +// These store constitutive properties evaluated at BC face conditions +DECLARE_FIELD( facePhaseMobility, + "facePhaseMobility", + array2dLayoutPhase, + 0, + NOPLOT, + NO_WRITE, + "Phase mobility at boundary faces evaluated at BC conditions" ); + +DECLARE_FIELD( facePhaseMassDensity, + "facePhaseMassDensity", + array2dLayoutPhase, + 0, + NOPLOT, + NO_WRITE, + "Phase mass density at boundary faces evaluated at BC conditions" ); + +DECLARE_FIELD( facePhaseCompFraction, + "facePhaseCompFraction", + array3dLayoutPhaseComp, + 0, + NOPLOT, + NO_WRITE, + "Phase component fraction at boundary faces evaluated at BC conditions" ); + // this is needed for time step selector DECLARE_FIELD( phaseVolumeFraction_n, "phaseVolumeFraction_n", diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp index a022ad76443..be78e7a85bf 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.cpp @@ -22,11 +22,14 @@ #include "mesh/DomainPartition.hpp" #include "constitutive/ConstitutivePassThru.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "constitutive/fluid/multifluid/MultiFluidSelector.hpp" #include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" #include "fieldSpecification/AquiferBoundaryCondition.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" #include "finiteVolume/HybridMimeticDiscretization.hpp" #include "finiteVolume/MimeticInnerProductDispatch.hpp" +#include "finiteVolume/FluxApproximationBase.hpp" +#include "finiteVolume/BoundaryStencil.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" #include "physicsSolvers/LogLevelsInfo.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" @@ -35,6 +38,7 @@ #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/SolutionCheckKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ResidualNormKernel.hpp" +#include "mesh/CellElementSubRegion.hpp" /** * @namespace the geos namespace that encapsulates the majority of the code @@ -73,8 +77,38 @@ void CompositionalMultiphaseHybridFVM::registerDataOnMesh( Group & meshBodies ) faceManager.registerField< flow::facePressure_n >( getName() ); - // auxiliary data for the buoyancy coefficient - faceManager.registerField< flow::mimGravityCoefficient >( getName() ); + // Register the face data for global component fraction + faceManager.registerField< flow::faceGlobalCompFraction >( getName() ); + + // Register the face data for temperature + faceManager.registerField< flow::faceTemperature >( getName() ); + + // Register the bc face data for pressure + faceManager.registerField< flow::bcPressure >( getName() ); + + // Register the bc face data for global component fraction + faceManager.registerField< flow::bcGlobalCompFraction >( getName() ); + + // Register the bc face data for temperature + faceManager.registerField< flow::bcTemperature >( getName() ); + + // Register face-based constitutive properties for BC faces + // These will store fluid properties evaluated at BC conditions + faceManager.registerField< flow::facePhaseMobility >( getName() ). + reference().resizeDimension< 1 >( m_numPhases ); + + faceManager.registerField< flow::facePhaseMassDensity >( getName() ). + reference().resizeDimension< 1 >( m_numPhases ); + + faceManager.registerField< flow::facePhaseCompFraction >( getName() ). + reference().resizeDimension< 1, 2 >( m_numPhases, m_numComponents ); + + // Register boundary face indicator (1 for boundary faces with Dirichlet BCs, 0 for interior) + faceManager.registerField< flow::isBoundaryFace >( getName() ); + + // precomputed mimetic gravity-driven trans coefficient on faces + faceManager.registerField< flow::mimeticTransGgradZ >( getName() ); + } ); } @@ -109,11 +143,10 @@ void CompositionalMultiphaseHybridFVM::initializePostInitialConditionsPreSubGrou MimeticInnerProductBase const & mimeticInnerProductBase = hmDiscretization.getReference< MimeticInnerProductBase >( HybridMimeticDiscretization::viewKeyStruct::innerProductString() ); if( dynamicCast< QuasiRTInnerProduct const * >( &mimeticInnerProductBase ) || - dynamicCast< QuasiTPFAInnerProduct const * >( &mimeticInnerProductBase ) || dynamicCast< SimpleInnerProduct const * >( &mimeticInnerProductBase ) ) { GEOS_ERROR( getCatalogName() << " " << getDataContext() << - "The QuasiRT, QuasiTPFA, and Simple inner products are only available in SinglePhaseHybridFVM" ); + "The QuasiRT, and Simple inner products are only available in SinglePhaseHybridFVM" ); } m_lengthTolerance = domain.getMeshBody( 0 ).getGlobalLengthScale() * 1e-8; @@ -123,7 +156,7 @@ void CompositionalMultiphaseHybridFVM::initializePostInitialConditionsPreSubGrou string_array const & regionNames ) { ElementRegionManager const & elemManager = mesh.getElemManager(); - FaceManager const & faceManager = mesh.getFaceManager(); + FaceManager & faceManager = mesh.getFaceManager(); CompositionalMultiphaseBase::initializePostInitialConditionsPreSubGroups(); @@ -146,10 +179,52 @@ void CompositionalMultiphaseHybridFVM::initializePostInitialConditionsPreSubGrou GEOS_THROW_IF( minVal.get() <= 0.0, getCatalogName() << " " << getDataContext() << - ": the transmissibility multipliers used in SinglePhaseHybridFVM must strictly larger than 0.0", - std::runtime_error, getDataContext() ); + ": the transmissibility multipliers used in SinglePhaseHybridFVM must be strictly larger than 0.0", + std::runtime_error ); + + // Initialize face-based constitutive property arrays to zero to prevent uninitialized memory usage on GPU + arrayView2d< real64, compflow::USD_PHASE > facePhaseMob = faceManager.getField< flow::facePhaseMobility >(); + arrayView2d< real64, compflow::USD_PHASE > facePhaseMassDens = faceManager.getField< flow::facePhaseMassDensity >(); + arrayView3d< real64, compflow::USD_PHASE_COMP > facePhaseCompFrac = faceManager.getField< flow::facePhaseCompFraction >(); + + localIndex const numFaces = faceManager.size(); + forAll< parallelDevicePolicy<> >( numFaces, [=] GEOS_HOST_DEVICE ( localIndex const iface ) + { + for( integer ip = 0; ip < facePhaseMob.size( 1 ); ++ip ) + { + facePhaseMob[iface][ip] = 0.0; + facePhaseMassDens[iface][ip] = 0.0; + for( integer ic = 0; ic < facePhaseCompFrac.size( 2 ); ++ic ) + { + facePhaseCompFrac[iface][ip][ic] = 0.0; + } + } + } ); + + // Mark boundary faces (faces with Dirichlet BCs) to skip flux continuity constraint + // Initialize all faces as interior (0), then mark boundary faces (1) + arrayView1d< integer > const isBoundaryFaceView = faceManager.getReference< array1d< integer > >( flow::isBoundaryFace::key() ); + // isBoundaryFaceView is default-initialized to zero FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); + fsManager.forSubGroups< FieldSpecificationBase >( [&]( FieldSpecificationBase const & fs ) + { + string const & fieldName = fs.getFieldName(); + if( fieldName == flow::bcPressure::key() || + fieldName == flow::bcGlobalCompFraction::key() || + fieldName == flow::bcTemperature::key() ) + { + for( string const & setName : fs.getSetNames() ) + { + SortedArrayView< localIndex const > const targetSet = faceManager.getSet( setName ).toViewConst(); + forAll< serialPolicy >( targetSet.size(), [=]( localIndex const i ) + { + isBoundaryFaceView[targetSet[i]] = 1; + } ); + } + } + } ); + fsManager.forSubGroups< AquiferBoundaryCondition >( [&] ( AquiferBoundaryCondition const & bc ) { GEOS_LOG_RANK_0( getCatalogName() << " " << getDataContext() << ": An aquifer boundary condition named " << @@ -164,67 +239,57 @@ void CompositionalMultiphaseHybridFVM::precomputeData( MeshLevel & mesh, string_ { FlowSolverBase::precomputeData( mesh, regionNames ); - NodeManager const & nodeManager = mesh.getNodeManager(); - FaceManager & faceManager = mesh.getFaceManager(); - - array1d< RAJA::ReduceSum< serialReduce, real64 > > mimFaceGravCoefNumerator; - array1d< RAJA::ReduceSum< serialReduce, real64 > > mimFaceGravCoefDenominator; - mimFaceGravCoefNumerator.resize( faceManager.size() ); - mimFaceGravCoefDenominator.resize( faceManager.size() ); - - // node data - - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); - - // face data - - arrayView1d< real64 const > const & transMultiplier = - faceManager.getField< flow::transMultiplier >(); - - arrayView1d< real64 > const mimFaceGravCoef = - faceManager.getField< flow::mimGravityCoefficient >(); - - ArrayOfArraysView< localIndex const > const & faceToNodes = faceManager.nodeList().toViewConst(); - - real64 const lengthTolerance = m_lengthTolerance; - - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - CellElementSubRegion & subRegion ) + // Pre-compute and initialize face mimeticTransGgradZ once (only for HybridMimetic discretization) { - arrayView2d< real64 const > const & elemCenter = - subRegion.template getReference< array2d< real64 > >( CellElementSubRegion::viewKeyStruct::elementCenterString() ); - string const & permModelName = subRegion.getReference< string >( viewKeyStruct::permeabilityNamesString() ); - arrayView3d< real64 const > const & elemPerm = - getConstitutiveModel< PermeabilityBase >( subRegion, permModelName ).permeability(); - arrayView1d< real64 const > const elemGravCoef = - subRegion.template getReference< array1d< real64 > >( flow::gravityCoefficient::key() ); - arrayView1d< real64 const > const & elemVolume = subRegion.getElementVolume(); - arrayView2d< localIndex const > const & elemToFaces = subRegion.faceList(); - - // here we precompute some quantities (mimFaceFracCoef) used in the FluxKernel to assemble the one-sided gravity term in the transport - // scheme - // This one-sided gravity term is currently always treated with TPFA, as in MRST. - // In the future, I will change that (here and in the FluxKernel) to have a consistent inner product for the gravity term as well - compositionalMultiphaseHybridFVMKernels:: - simpleKernelLaunchSelector< PrecomputeKernel, - mimeticInnerProduct::TPFAInnerProduct >( subRegion.numFacesPerElement(), - subRegion.size(), - faceManager.size(), - nodePosition, - faceToNodes, - elemCenter, - elemVolume, - elemPerm, - elemGravCoef, - elemToFaces, - transMultiplier, - lengthTolerance, - mimFaceGravCoefNumerator.toView(), - mimFaceGravCoefDenominator.toView(), - mimFaceGravCoef ); + DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); + NumericalMethodsManager const & nm = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = nm.getFiniteVolumeManager(); - } ); + HybridMimeticDiscretization const & hm = fvManager.getHybridMimeticDiscretization( m_discretizationName ); + mimeticInnerProduct::MimeticInnerProductBase const & ip = hm.getReference< mimeticInnerProduct::MimeticInnerProductBase >( HybridMimeticDiscretization::viewKeyStruct::innerProductString() ); + real64 const lengthTolerance = domain.getMeshBody( 0 ).getGlobalLengthScale() * m_lengthTolerance; + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, string_array const & regionNames ) + { + FaceManager & faceManager = mesh.getFaceManager(); + // Temporary accumulators + array1d< real64 > invSum( faceManager.size() ); invSum.setValues< parallelDevicePolicy<> >( 0.0 ); + array1d< integer > count( faceManager.size() ); count.setValues< parallelDevicePolicy<> >( 0 ); + + NodeManager const & nodeManager = mesh.getNodeManager(); + mesh.getElemManager().forElementSubRegionsComplete< CellElementSubRegion >( regionNames, + [&]( localIndex const, localIndex const er, localIndex const esr, ElementRegionBase const &, + CellElementSubRegion const & subRegion ) + { + GEOS_UNUSED_VAR( er ); + GEOS_UNUSED_VAR( esr ); + string const & permName = subRegion.getReference< string >( viewKeyStruct::permeabilityNamesString() ); + PermeabilityBase const & permeability = getConstitutiveModel< PermeabilityBase >( subRegion, permName ); + PrecomputeMimeticTransGgradZKernel::createAndLaunch< parallelDevicePolicy<> >( ip, + nodeManager, + faceManager, + subRegion, + permeability, + lengthTolerance, + invSum.toView(), + count.toView() ); + } ); + + // Reduce to effective value per face and write to field + arrayView1d< real64 > mimeticTransGgradZ = faceManager.getField< flow::mimeticTransGgradZ >(); + arrayView1d< integer const > ghost = faceManager.ghostRank(); + forAll< parallelDevicePolicy<> >( faceManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const kf ) + { + if( ghost[kf] >= 0 ) + { + return; + } + real64 const s = invSum[kf]; + integer const c = count[kf]; + mimeticTransGgradZ[kf] = (c > 0 && s > 0.0) ? static_cast< real64 >( c ) / s : 0.0; + } ); + } ); + } } void CompositionalMultiphaseHybridFVM::implicitStepSetup( real64 const & time_n, @@ -270,10 +335,6 @@ void CompositionalMultiphaseHybridFVM::setupDofs( DomainPartition const & GEOS_U viewKeyStruct::elemDofFieldString(), DofManager::Connector::Face ); - // this call with instruct GEOS to reorder the dof numbers - //dofManager.setLocalReorderingType( viewKeyStruct::elemDofFieldString(), - // DofManager::LocalReorderingType::ReverseCutHillMcKee ); - // for the volume balance equation, disable global coupling // this equation is purely local (not coupled to neighbors or other physics) dofManager.disableGlobalCouplingForEquation( viewKeyStruct::elemDofFieldString(), @@ -318,12 +379,8 @@ void CompositionalMultiphaseHybridFVM::assembleFluxTerms( real64 const dt, NodeManager const & nodeManager = mesh.getNodeManager(); FaceManager const & faceManager = mesh.getFaceManager(); - - // node data (for transmissibility computation) - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); - // face data // get the face-based DOF numbers for the assembly string const faceDofKey = dofManager.getKey( viewKeyStruct::faceDofFieldString() ); @@ -331,6 +388,10 @@ void CompositionalMultiphaseHybridFVM::assembleFluxTerms( real64 const dt, faceManager.getReference< array1d< globalIndex > >( faceDofKey ); arrayView1d< integer const > const & faceGhostRank = faceManager.ghostRank(); + // Get boundary face indicator (initialized during initializePostInitialConditionsPreSubGroups) + arrayView1d< integer const > const isBoundaryFaceView = + faceManager.getField< flow::isBoundaryFace >(); + // get the element dof numbers for the assembly string const & elemDofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > elemDofNumber = @@ -344,8 +405,8 @@ void CompositionalMultiphaseHybridFVM::assembleFluxTerms( real64 const dt, // get the face-centered depth arrayView1d< real64 const > const & faceGravCoef = faceManager.getField< flow::gravityCoefficient >(); - arrayView1d< real64 const > const & mimFaceGravCoef = - faceManager.getField< flow::mimGravityCoefficient >(); + arrayView1d< real64 const > const & mimeticTransGgradZ = + faceManager.getField< flow::mimeticTransGgradZ >(); // get the face-centered transMultiplier arrayView1d< real64 const > const & transMultiplier = @@ -392,9 +453,10 @@ void CompositionalMultiphaseHybridFVM::assembleFluxTerms( real64 const dt, faceToNodes, faceDofNumber, faceGhostRank, + isBoundaryFaceView, facePres, faceGravCoef, - mimFaceGravCoef, + mimeticTransGgradZ, transMultiplier, compFlowAccessors.get( flow::phaseMobility{} ), compFlowAccessors.get( flow::dPhaseMobility{} ), @@ -570,7 +632,375 @@ void CompositionalMultiphaseHybridFVM::applyBoundaryConditions( real64 const tim CompositionalMultiphaseBase::applyBoundaryConditions( time_n, dt, domain, dofManager, localMatrix, localRhs ); - // TODO: implement face boundary conditions here + // Apply face-based Dirichlet boundary conditions + applyFaceDirichletBC( time_n, dt, dofManager, domain, localMatrix, localRhs ); +} + +void CompositionalMultiphaseHybridFVM::applyFaceDirichletBC( real64 const time_n, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + using namespace compositionalMultiphaseHybridFVMKernels; + using namespace mimeticInnerProduct; + + FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance(); + + // Get the inner product type from the discretization + NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); + FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + HybridMimeticDiscretization const & hmDiscretization = fvManager.getHybridMimeticDiscretization( m_discretizationName ); + string const & innerProductType = hmDiscretization.getReference< string >( HybridMimeticDiscretization::viewKeyStruct::innerProductTypeString() ); + + string const elemDofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); + + // Log message for Dirichlet BC application + static char const faceBcLogMessage[] = + "CompositionalMultiphaseHybridFVM {}: at time {}s, " + "the <{}> boundary condition '{}' is applied to the face set '{}' in '{}'. " + "\nThe scale of this boundary condition is {} and multiplies the value of the provided function (if any). " + "\nThe total number of target faces (including ghost faces) is {}." + "\nNote that if this number is equal to zero, the boundary condition will not be applied on this face set."; + + // Apply Dirichlet BC by computing boundary fluxes + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + FaceManager & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + ElementRegionManager & elemManager = mesh.getElemManager(); + string const faceDofKey = dofManager.getKey( viewKeyStruct::faceDofFieldString() ); + arrayView1d< globalIndex const > const faceDofNumber = + faceManager.getReference< array1d< globalIndex > >( faceDofKey ); + arrayView1d< integer const > const faceGhostRank = faceManager.ghostRank(); + globalIndex const rankOffset = dofManager.rankOffset(); + + // Get node positions + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = nodeManager.referencePosition(); + + // Get face data + ArrayOfArraysView< localIndex const > const faceToNodes = faceManager.nodeList().toViewConst(); + arrayView2d< localIndex const > const elemRegionList = faceManager.elementRegionList(); + arrayView2d< localIndex const > const elemSubRegionList = faceManager.elementSubRegionList(); + arrayView2d< localIndex const > const elemList = faceManager.elementList(); + + arrayView1d< real64 const > const transMultiplier = faceManager.getReference< array1d< real64 > >( fields::flow::transMultiplier::key() ); + + // Apply boundary values to face fields first + applyFieldValue< FaceManager >( time_n, dt, mesh, faceBcLogMessage, + flow::bcPressure::key(), flow::facePressure::key() ); + + applyFieldValue< FaceManager >( time_n, dt, mesh, faceBcLogMessage, + flow::bcGlobalCompFraction::key(), flow::faceGlobalCompFraction::key() ); + + applyFieldValue< FaceManager >( time_n, dt, mesh, faceBcLogMessage, + flow::bcTemperature::key(), flow::faceTemperature::key() ); + + // Get face boundary values + arrayView1d< real64 const > const facePres = faceManager.getField< fields::flow::facePressure >(); + arrayView1d< real64 const > const faceTemp = faceManager.getField< fields::flow::faceTemperature >(); + arrayView2d< real64 const, compflow::USD_COMP > const faceCompFrac = faceManager.getField< fields::flow::faceGlobalCompFraction >(); + arrayView1d< real64 const > const faceGravCoef = faceManager.getField< fields::flow::gravityCoefficient >(); + + // Loop over regions and apply Dirichlet flux kernel + elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const erIndex, + CellElementSubRegion & subRegion ) + { + string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); + MultiFluidBase & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); + + string const & permName = subRegion.getReference< string >( viewKeyStruct::permeabilityNamesString() ); + PermeabilityBase const & permeabilityModel = getConstitutiveModel< PermeabilityBase >( subRegion, permName ); + + string const & relPermName = subRegion.getReference< string >( viewKeyStruct::relPermNamesString() ); + RelativePermeabilityBase & relperm = getConstitutiveModel< RelativePermeabilityBase >( subRegion, relPermName ); + + real64 const lengthTolerance = m_lengthTolerance; + + // Get all face sets that have Dirichlet BCs + std::set< string > bcFaceSets; + fsManager.forSubGroups< FieldSpecificationBase >( [&]( FieldSpecificationBase const & fs ) + { + string const & fieldName = fs.getFieldName(); + if( fieldName == flow::bcPressure::key() || + fieldName == flow::bcGlobalCompFraction::key() || + fieldName == flow::bcTemperature::key() ) + { + for( string const & setName : fs.getSetNames() ) + { + bcFaceSets.insert( setName ); + } + } + } ); + + // Get writable face arrays for storing BC face properties + arrayView2d< real64, compflow::USD_PHASE > facePhaseMob = + faceManager.getField< flow::facePhaseMobility >(); + arrayView2d< real64, compflow::USD_PHASE > facePhaseMassDens = + faceManager.getField< flow::facePhaseMassDensity >(); + arrayView3d< real64, compflow::USD_PHASE_COMP > facePhaseCompFrac = + faceManager.getField< flow::facePhaseCompFraction >(); + + // Move arrays to host memory before evaluateBCFaceProperties runs with serialPolicy + facePhaseMob.move( hostMemorySpace, true ); + facePhaseMassDens.move( hostMemorySpace, true ); + facePhaseCompFrac.move( hostMemorySpace, true ); + + // Evaluate constitutive properties at BC face conditions for each face set + for( string const & setName : bcFaceSets ) + { + SortedArrayView< localIndex const > const targetSet = faceManager.getSet( setName ).toViewConst(); + if( targetSet.size() == 0 ) + continue; + + // Call evaluateBCFaceProperties to compute face properties at BC conditions + constitutive::constitutiveComponentUpdatePassThru( fluid, m_numComponents, [&]( auto & fluidWrapper, auto NC ) + { + GEOS_UNUSED_VAR( fluidWrapper ); + integer constexpr NUM_COMP = NC(); + + auto evaluateWithPhases = [&]( auto NP_VALUE ) + { + integer constexpr NUM_PHASES = decltype( NP_VALUE )::value; + + compositionalMultiphaseHybridFVMKernels::evaluateBCFaceProperties< NUM_COMP, NUM_PHASES > + ( m_numPhases, + targetSet, + facePres, + faceTemp, + faceCompFrac, + elemRegionList, + elemSubRegionList, + elemList, + erIndex, + subRegion.getIndexInParent(), + fluid, + relperm, + facePhaseMob, + facePhaseMassDens, + facePhaseCompFrac ); + }; + + if( m_numPhases == 2 ) + { + evaluateWithPhases( std::integral_constant< integer, 2 >() ); + } + else if( m_numPhases == 3 ) + { + evaluateWithPhases( std::integral_constant< integer, 3 >() ); + } + } ); + } + + // CRITICAL: Move the SAME array views that were just modified back to device memory + // Don't get new views - use the views that evaluateBCFaceProperties actually modified + // evaluateBCFaceProperties uses serialPolicy (host), DirichletFluxKernel uses parallelDevicePolicy (device) + facePhaseMob.move( parallelDeviceMemorySpace, false ); + facePhaseMassDens.move( parallelDeviceMemorySpace, false ); + facePhaseCompFrac.move( parallelDeviceMemorySpace, false ); + + // Get const views to the face properties for use in DirichletFluxKernel + arrayView2d< real64 const, compflow::USD_PHASE > const facePhaseMobField = + faceManager.getField< flow::facePhaseMobility >(); + arrayView2d< real64 const, compflow::USD_PHASE > const facePhaseMassDensField = + faceManager.getField< flow::facePhaseMassDensity >(); + arrayView3d< real64 const, compflow::USD_PHASE_COMP > const facePhaseCompFracField = + faceManager.getField< flow::facePhaseCompFraction >(); + + // Apply Dirichlet boundary fluxes for each face set using DirichletFluxKernel + for( string const & setName : bcFaceSets ) + { + SortedArrayView< localIndex const > const targetSet = faceManager.getSet( setName ).toViewConst(); + if( targetSet.size() == 0 ) + continue; + + // Launch the Dirichlet flux kernel with compile-time dispatch + constitutive::constitutiveComponentUpdatePassThru( fluid, m_numComponents, [&]( auto & fluidWrapper, auto NC ) + { + GEOS_UNUSED_VAR( fluidWrapper ); + integer constexpr NUM_COMP = NC(); + + typename DirichletFluxKernel::CompFlowAccessors compFlowAccessors( elemManager, getName() ); + typename DirichletFluxKernel::MultiFluidAccessors multiFluidAccessors( elemManager, getName() ); + ElementRegionManager::ElementViewAccessor< arrayView1d< globalIndex const > > elemDofNumberAccessor = + elemManager.constructArrayViewAccessor< globalIndex, 1 >( elemDofKey ); + + localIndex const numFacesPerElement = subRegion.numFacesPerElement(); + + auto launchWithPhases = [&]( auto NP_VALUE ) + { + integer constexpr NP = decltype( NP_VALUE )::value; + + auto launchKernel = [&]( auto IP_TYPE_WRAPPER, auto NF_VALUE ) + { + using IP_TYPE = decltype( IP_TYPE_WRAPPER ); + integer constexpr NF = decltype( NF_VALUE )::value; + + DirichletFluxKernel::launch< NF, NUM_COMP, NP, IP_TYPE > + ( m_numPhases, + erIndex, + subRegion.getIndexInParent(), + subRegion, + permeabilityModel, + targetSet, + facePres, + faceTemp, + facePhaseMobField, + facePhaseMassDensField, + facePhaseCompFracField, + nodePosition, + faceToNodes, + elemRegionList, + elemSubRegionList, + elemList, + faceDofNumber, + faceGravCoef, + transMultiplier, + lengthTolerance, + dt, + rankOffset, + m_useTotalMassEquation, + compFlowAccessors, + multiFluidAccessors, + elemDofNumberAccessor.toNestedViewConst(), + localMatrix, + localRhs ); + }; + + // Helper lambda to dispatch on number of faces per element + auto launchForFaces = [&]( auto IP_TAG ) + { + if( numFacesPerElement == 4 ) + launchKernel( IP_TAG, std::integral_constant< integer, 4 >{} ); + else if( numFacesPerElement == 5 ) + launchKernel( IP_TAG, std::integral_constant< integer, 5 >{} ); + else if( numFacesPerElement == 6 ) + launchKernel( IP_TAG, std::integral_constant< integer, 6 >{} ); + else if( numFacesPerElement == 7 ) + launchKernel( IP_TAG, std::integral_constant< integer, 7 >{} ); + else if( numFacesPerElement == 8 ) + launchKernel( IP_TAG, std::integral_constant< integer, 8 >{} ); + else if( numFacesPerElement == 9 ) + launchKernel( IP_TAG, std::integral_constant< integer, 9 >{} ); + else if( numFacesPerElement == 10 ) + launchKernel( IP_TAG, std::integral_constant< integer, 10 >{} ); + else if( numFacesPerElement == 11 ) + launchKernel( IP_TAG, std::integral_constant< integer, 11 >{} ); + else if( numFacesPerElement == 12 ) + launchKernel( IP_TAG, std::integral_constant< integer, 12 >{} ); + else if( numFacesPerElement == 13 ) + launchKernel( IP_TAG, std::integral_constant< integer, 13 >{} ); + else + GEOS_ERROR( "Unsupported number of faces per element: " << numFacesPerElement ); + }; + + // Inner-product selection + if( innerProductType == "TPFA" ) + { + launchForFaces( TPFAInnerProduct{} ); + } + else if( innerProductType == "quasiTPFA" ) + { + launchForFaces( QuasiTPFAInnerProduct{} ); + } + else if( innerProductType == "beiraoDaVeigaLipnikovManzini" ) + { + launchForFaces( BdVLMInnerProduct{} ); + } + else + { + GEOS_ERROR( "Unsupported inner product type: " << innerProductType ); + } + }; + + if( m_numPhases == 2 ) + launchWithPhases( std::integral_constant< integer, 2 >{} ); + else if( m_numPhases == 3 ) + launchWithPhases( std::integral_constant< integer, 3 >{} ); + else + GEOS_ERROR( "Unsupported number of phases: " << m_numPhases ); + } ); + } + } ); + + // Keep strong enforcement that the face pressure equals the informed bcPressure value (original behavior) + arrayView1d< real64 const > const presFace = + faceManager.getField< flow::facePressure >(); + arrayView1d< real64 const > const presFaceBC = + faceManager.getField< flow::bcPressure >(); + + fsManager.apply< FaceManager >( time_n + dt, + mesh, + flow::bcPressure::key(), + [&] ( FieldSpecificationBase const & fs, + string const & setName, + SortedArrayView< localIndex const > const & targetSet, + FaceManager & targetGroup, + string const & ) + { + GEOS_UNUSED_VAR( setName ); + // Using the field specification functions to apply the boundary conditions to the system + fs.applyFieldValue< FieldSpecificationEqual, + parallelDevicePolicy<> >( targetSet, + time_n + dt, + targetGroup, + flow::bcPressure::key() ); + + forAll< parallelDevicePolicy<> >( targetSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const a ) + { + + localIndex const kf = targetSet[a]; + if( faceGhostRank[kf] >= 0 ) + { + return; + } + + // Get the dof number of this face + globalIndex const dofIndex = faceDofNumber[kf]; + localIndex const localRow = dofIndex - rankOffset; + + // ENFORCE DIRICHLET CONSTRAINT: x_i = x_spec + // Mathematical procedure to enforce prescribed value in Ax = b: + // 1. For row i: Zero all entries except diagonal, set A[i,i] = 1, set b[i] = x_spec + // 2. For all other rows k: subtract A[k,i] * x_spec from b[k], then set A[k,i] = 0 + // + // Note: Step 2 (removing column influence) should ideally be done, but in practice + // for boundary face DOFs in hybrid FVM, the column entries from interior faces to + // boundary faces are typically zero or minimal because boundary faces don't strongly + // couple to interior equations. The strong coupling is boundary->interior (via fluxes). + // For exact enforcement in non-trivial cases, column zeroing would be needed. + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + arraySlice1d< globalIndex const > const columns = localMatrix.getColumns( localRow ); + arraySlice1d< real64 > const entries = localMatrix.getEntries( localRow ); + localIndex const numEntries = localMatrix.numNonZeros( localRow ); + + // Step 1: Zero out all row entries and set diagonal to 1 + for( localIndex j = 0; j < numEntries; ++j ) + { + if( columns[j] == dofIndex ) + { + entries[j] = 1.0; // Set diagonal to 1 + } + else + { + entries[j] = 0.0; // Zero out off-diagonal entries + } + } + + // Set RHS to the prescribed boundary value (absolute value) + localRhs[localRow] = presFace[kf] - presFaceBC[kf]; + } + } ); + } ); + + } ); } void CompositionalMultiphaseHybridFVM::applyAquiferBC( real64 const time, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp index e0be3309446..38976e66374 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp @@ -168,6 +168,23 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase virtual void initializePreSubGroups() override; + /** + * @brief Function to perform the application of Dirichlet BC data on faces + * @param[in] time_n 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[inout] localMatrix the system matrix + * @param[inout] localRhs the system right-hand side vector + */ + void + applyFaceDirichletBC( real64 const time_n, + real64 const dt, + DofManager const & dofManager, + DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + protected: /// precompute the minGravityCoefficient for the buoyancy term @@ -175,10 +192,10 @@ class CompositionalMultiphaseHybridFVM : public CompositionalMultiphaseBase private: - /// tolerance used in the computation of the transmissibility matrix + /// Tolerance used in the computation of the transmissibility matrix (typically domain_length * 1e-8) real64 m_lengthTolerance; - /// region filter used in flux assembly + /// Set of region indices that are targeted by this solver (used to filter flux assembly) SortedArray< localIndex > m_regionFilter; }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_assembly.cpp.template b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_assembly.cpp.template deleted file mode 100644 index d06fd73408c..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_assembly.cpp.template +++ /dev/null @@ -1,162 +0,0 @@ -#include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels_impl.hpp" - -namespace geos -{ -namespace compositionalMultiphaseHybridFVMKernels -{ - constexpr int NF = @NFACES@; - constexpr int NC = @NCOMPS@; - constexpr int NP = @NPHASES@; - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - AssemblerKernelHelper:: - applyGradient< NF, NC, NP >( arrayView1d< real64 const > const & facePres, - arrayView1d< real64 const > const & faceGravCoef, - arraySlice1d< localIndex const > const & elemToFaces, - real64 const & elemPres, - real64 const & elemGravCoef, - arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE-2 > const & elemPhaseMassDens, - arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC-2 > const & dElemPhaseMassDens_dCompFrac, - arraySlice1d< real64 const, compflow::USD_PHASE-1 > const & elemPhaseMob, - arraySlice2d< real64 const, compflow::USD_PHASE_DC-1 > const & dElemPhaseMob, - arraySlice2d< real64 const, compflow::USD_COMP_DC-1 > const & dElemCompFrac_dCompDens, - arraySlice2d< real64 const > const & transMatrix, - real64 ( &oneSidedVolFlux )[ NF ], - real64 ( &dOneSidedVolFlux_dPres )[ NF ], - real64 ( &dOneSidedVolFlux_dFacePres )[ NF ][ NF ], - real64 ( &dOneSidedVolFlux_dCompDens )[ NF ][ NC ] ); - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - AssemblerKernelHelper:: - assembleFluxDivergence< NF, NC, NP >( localIndex const (&localIds)[ 3 ], - globalIndex const rankOffset, - arrayView2d< localIndex const > const & elemRegionList, - arrayView2d< localIndex const > const & elemSubRegionList, - arrayView2d< localIndex const > const & elemList, - SortedArrayView< localIndex const > const & regionFilter, - arrayView1d< globalIndex const > const & faceDofNumber, - arrayView1d< real64 const > const & mimFaceGravCoef, - arraySlice1d< localIndex const > const & elemToFaces, - real64 const & elemGravCoef, - integer const useTotalMassEquation, - ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseDens, - ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, - ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, - ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, - ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, - arraySlice2d< real64 const > const & transMatrixGrav, - real64 const (&oneSidedVolFlux)[ NF ], - real64 const (&dOneSidedVolFlux_dPres)[ NF ], - real64 const (&dOneSidedVolFlux_dFacePres)[ NF ][ NF ], - real64 const (&dOneSidedVolFlux_dCompDens)[ NF ][ NC ], - real64 const & dt, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - AssemblerKernelHelper:: - assembleViscousFlux< NF, NC, NP >( localIndex const ifaceLoc, - real64 const (&oneSidedVolFlux)[ NF ], - real64 const (&dOneSidedVolFlux_dPres)[ NF ], - real64 const (&dOneSidedVolFlux_dFacePres)[ NF ][ NF ], - real64 const (&dOneSidedVolFlux_dCompDens)[ NF ][ NC ], - real64 const (&upwPhaseViscCoef)[ NP ][ NC ], - real64 const (&dUpwPhaseViscCoef_dPres)[ NP ][ NC ], - real64 const (&dUpwPhaseViscCoef_dCompDens)[ NP ][ NC ][ NC ], - globalIndex const elemDofNumber, - globalIndex const neighborDofNumber, - globalIndex const upwViscDofNumber, - globalIndex const faceDofNumber, - real64 const & dt, - real64 ( &divMassFluxes )[ NC ], - real64 ( &dDivMassFluxes_dElemVars )[ NC ][ (NC+1)*(NF+1) ], - real64 ( &dDivMassFluxes_dFaceVars )[ NC ][ NF ], - globalIndex ( &dofColIndicesElemVars )[ (NC+1)*(NF+1) ], - globalIndex ( &dofColIndicesFaceVars )[ NF ] ); - @EXTERN@ - template - GEOS_HOST_DEVICE - void - AssemblerKernelHelper:: - assembleBuoyancyFlux< NF, NC, NP >( localIndex const ifaceLoc, - real64 const (&phaseGravTerm)[ NP ][ NP-1 ], - real64 const (&dPhaseGravTerm_dPres)[ NP ][ NP-1 ][ 2 ], - real64 const (&dPhaseGravTerm_dCompDens)[ NP ][ NP-1 ][ 2 ][ NC ], - real64 const (&upwPhaseGravCoef)[ NP ][ NP-1 ][ NC ], - real64 const (&dUpwPhaseGravCoef_dPres)[ NP ][ NP-1 ][ NC ][ 2 ], - real64 const (&dUpwPhaseGravCoef_dCompDens)[ NP ][ NP-1 ][ NC ][ 2 ][ NC ], - real64 const & dt, - real64 ( &divMassFluxes )[ NC ], - real64 ( &dDivMassFluxes_dElemVars )[ NC ][ (NC+1)*(NF+1) ] ); - @EXTERN@ - template - GEOS_HOST_DEVICE - void - AssemblerKernelHelper:: - assembleFaceConstraints< NF, NC >( arrayView1d< globalIndex const > const & faceDofNumber, - arrayView1d< integer const > const & faceGhostRank, - arraySlice1d< localIndex const > const & elemToFaces, - globalIndex const elemDofNumber, - globalIndex const rankOffset, - real64 const (&oneSidedVolFlux)[ NF ], - real64 const (&dOneSidedVolFlux_dPres)[ NF ], - real64 const (&dOneSidedVolFlux_dFacePres)[ NF ][ NF ], - real64 const (&dOneSidedVolFlux_dCompDens)[ NF ][ NC ], - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - AssemblerKernel:: - compute< NF, NC, NP >( localIndex const er, localIndex const esr, localIndex const ei, - SortedArrayView< localIndex const > const & regionFilter, - arrayView2d< localIndex const > const & elemRegionList, - arrayView2d< localIndex const > const & elemSubRegionList, - arrayView2d< localIndex const > const & elemList, - arrayView1d< globalIndex const > const & faceDofNumber, - arrayView1d< integer const > const & faceGhostRank, - arrayView1d< real64 const > const & facePres, - arrayView1d< real64 const > const & faceGravCoef, - arrayView1d< real64 const > const & mimFaceGravCoef, - arraySlice1d< localIndex const > const & elemToFaces, - real64 const & elemPres, - real64 const & elemGravCoef, - integer const useTotalMassEquation, - ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseDens, - ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, - ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, - ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, - ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, - integer const elemGhostRank, - globalIndex const rankOffset, - real64 const & dt, - arraySlice2d< real64 const > const & transMatrix, - arraySlice2d< real64 const > const & transMatrixGrav, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - -} -} - diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_flux.cpp.template b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_flux.cpp.template deleted file mode 100644 index 76c8fef1aac..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_flux.cpp.template +++ /dev/null @@ -1,51 +0,0 @@ -#include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels_impl.hpp" - -namespace geos -{ -namespace compositionalMultiphaseHybridFVMKernels -{ - constexpr int NF = @NFACES@; - constexpr int NC = @NCOMPS@; - constexpr int NP = @NPHASES@; - using IP_TYPE = @IP_TYPE@; - - @EXTERN@ - template - void - FluxKernel:: - launch< NF, NC, NP, IP_TYPE >( localIndex er, localIndex esr, - CellElementSubRegion const & subRegion, - constitutive::PermeabilityBase const & permeabilityModel, - SortedArrayView< localIndex const > const & regionFilter, - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, - arrayView2d< localIndex const > const & elemRegionList, - arrayView2d< localIndex const > const & elemSubRegionList, - arrayView2d< localIndex const > const & elemList, - ArrayOfArraysView< localIndex const > const & faceToNodes, - arrayView1d< globalIndex const > const & faceDofNumber, - arrayView1d< integer const > const & faceGhostRank, - arrayView1d< real64 const > const & facePres, - arrayView1d< real64 const > const & faceGravCoef, - arrayView1d< real64 const > const & mimFaceGravCoef, - arrayView1d< real64 const > const & transMultiplier, - ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, - ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, - ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, - ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, - ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseMassDens, - ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac, - ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, - ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, - globalIndex const rankOffset, - real64 const lengthTolerance, - real64 const dt, - integer const useTotalMassEquation, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ); - -} -} - diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_upwinding.cpp.template b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_upwinding.cpp.template deleted file mode 100644 index caa6b2405fe..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVMKernels_upwinding.cpp.template +++ /dev/null @@ -1,98 +0,0 @@ -#include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp" -#include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels_impl.hpp" - -namespace geos -{ -namespace compositionalMultiphaseHybridFVMKernels -{ - constexpr int NC = @NCOMPS@; - constexpr int NP = @NPHASES@; - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - UpwindingHelper:: - upwindViscousCoefficient< NC, NP >( localIndex const (&localIds)[ 3 ], - localIndex const (&neighborIds)[ 3 ], - ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseDens, - ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, - ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, - ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, - real64 const & oneSidedVolFlux, - real64 ( &upwPhaseViscCoef )[ NP ][ NC ], - real64 ( &dUpwPhaseViscCoef_dPres )[ NP ][ NC ], - real64 ( &dUpwPhaseViscCoef_dCompDens )[ NP ][ NC ][ NC ], - globalIndex & upwViscDofNumber ); - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - UpwindingHelper:: - upwindBuoyancyCoefficient< NC, NP >( localIndex const (&localIds)[ 3 ], - localIndex const (&neighborIds)[ 3 ], - real64 const & transGravCoef, - ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseDens, - ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, - ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, - ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, - real64 ( &phaseGravTerm )[ NP ][ NP-1 ], - real64 ( &dPhaseGravTerm_dPres )[ NP ][ NP-1 ][ 2 ], - real64 ( &dPhaseGravTerm_dCompDens )[ NP ][ NP-1 ][ 2 ][ NC ], - real64 ( &upwPhaseGravCoef )[ NP ][ NP-1 ][ NC ], - real64 ( &dUpwPhaseGravCoef_dPres )[ NP ][ NP-1 ][ NC ][ 2 ], - real64 ( &dUpwPhaseGravCoef_dCompDens )[ NP ][ NP-1 ][ NC ][ 2 ][ NC ] ); - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - UpwindingHelper:: - computePhaseGravTerm< NC, NP >( localIndex const (&localIds)[ 3 ], - localIndex const (&neighborIds)[ 3 ], - real64 const & transGravCoef, - ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, - ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, - ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, - real64 ( &phaseGravTerm )[ NP ][ NP-1 ], - real64 ( &dPhaseGravTerm_dPres )[ NP ][ NP-1 ][ 2 ], - real64 ( &dPhaseGravTerm_dCompDens )[ NP ][ NP-1 ][ 2 ][ NC ] ); - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - UpwindingHelper:: - computeUpwindedTotalMobility< NC, NP >( localIndex const (&localIds)[ 3 ], - localIndex const (&neighborIds)[ 3 ], - ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, - ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, - real64 const (&phaseGravTerm)[ NP ][ NP-1 ], - real64 & totalMob, - real64 ( &dTotalMob_dPres )[ 2 ], - real64 ( &dTotalMob_dCompDens )[ 2 ][ NC ] ); - - @EXTERN@ - template - GEOS_HOST_DEVICE - void - UpwindingHelper:: - setIndicesForTotalMobilityUpwinding< NP >( localIndex const (&localIds)[ 3 ], - localIndex const (&neighborIds)[ 3 ], - real64 const (&gravTerm)[ NP ][ NP-1 ], - localIndex ( &totalMobIds )[ NP ][ 3 ], - localIndex ( &totalMobPos )[ NP ] ); - -} -} - diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp index 50e79137a4b..a8023ba1821 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp @@ -89,6 +89,14 @@ DECLARE_FIELD( facePressure_n, NO_WRITE, "Face pressure at the previous converged time step" ); +DECLARE_FIELD( isBoundaryFace, + "isBoundaryFace", + array1d< integer >, + 0, + NOPLOT, + WRITE_AND_READ, + "Boundary face indicator: 1 for faces with Dirichlet BCs, 0 for interior faces" ); + DECLARE_FIELD( pressureGradient, "pressureGradient", array2d< real64 >, @@ -177,13 +185,13 @@ DECLARE_FIELD( gravityCoefficient, WRITE_AND_READ, "Gravity coefficient (dot product of gravity acceleration by gravity vector)" ); -DECLARE_FIELD( mimGravityCoefficient, - "mimGravityCoefficient", +DECLARE_FIELD( mimeticTransGgradZ, + "mimeticTransGgradZ", array1d< real64 >, 0, NOPLOT, WRITE_AND_READ, - "Mimetic gravity coefficient" ); + "Consistent mimetic operator applied to g grad(z) (harmonic average of element contributions)" ); DECLARE_FIELD( macroElementIndex, "macroElementIndex", diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp index 531fd4054ec..9a50003023a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBLFields.hpp @@ -39,14 +39,6 @@ using array2dLayoutComp = array2d< real64, compflow::LAYOUT_COMP >; using array2dLayoutOBLOpVals = array2d< real64, compflow::LAYOUT_OBL_OPERATOR_VALUES >; using array3dLayoutOBLOpDers = array3d< real64, compflow::LAYOUT_OBL_OPERATOR_DERIVATIVES >; -DECLARE_FIELD( bcGlobalCompFraction, - "bcGlobalCompFraction", - array2dLayoutComp, - 0, - NOPLOT, - WRITE_AND_READ, - "Boundary condition global component fraction" ); - DECLARE_FIELD( referencePorosity, "referencePorosity", array1d< real64 >, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.cpp.template b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.cpp.template deleted file mode 100644 index 5a8a3464124..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseHybridFVMKernels.cpp.template +++ /dev/null @@ -1,21 +0,0 @@ -#include "physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernels.hpp" -#include "finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp" -#include "finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp" -#include "finiteVolume/mimeticInnerProducts/QuasiRTInnerProduct.hpp" -#include "finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp" -#include "finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp" - -using ElementBasedAssemblyKernelPolicy = @ElementBasedAssemblyKernelPolicy@; - -namespace geos -{ -namespace singlePhaseHybridFVMKernels -{ - @EXTERN@ template class @NAME@ < @NFACES@, @IP@ >; - @EXTERN@ template void @NAME@ < @NFACES@, @IP@ >::launch< @NAME@Policy, - @NAME@ < @NFACES@, @IP@ > > - ( const localIndex, - @NAME@ < @NFACES@, @IP@ > const & ); -} -} - diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernelSpecs.json b/src/coreComponents/physicsSolvers/fluidFlow/kernelSpecs.json index c398a1e7bd9..f80c23861e4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernelSpecs.json +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernelSpecs.json @@ -24,125 +24,6 @@ "generation call). Duplicate entries are unified internally" ] }, - "SinglePhaseHybridFVMKernels": { - "vars": [ - "NAME", - "NFACES", - "IP" - ], - "constants": [ - [ "ElementBasedAssemblyKernelPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ] - ], - "combinations": { - "NAME": [ - "ElementBasedAssemblyKernel" - ], - "NFACES": [ - 4, - 5, - 6, - 7, - 8, - 9, - 10, - 11, - 12, - 13 - ], - "IP": [ - "mimeticInnerProduct::TPFAInnerProduct", - "mimeticInnerProduct::QuasiTPFAInnerProduct", - "mimeticInnerProduct::QuasiRTInnerProduct", - "mimeticInnerProduct::SimpleInnerProduct", - "mimeticInnerProduct::BdVLMInnerProduct" - ] - }, - "explicit": [] - }, - - "CompositionalMultiphaseHybridFVMKernels_upwinding": { - "vars": [ - "NCOMPS", - "NPHASES" - ], - "constants": [], - "combinations": { - "NCOMPS": [ - 1, - 2, - 3, - 4, - 5 - ], - "NPHASES": [ - 2, - 3 - ] - }, - "explicit": [] - }, - - "CompositionalMultiphaseHybridFVMKernels_assembly": { - "vars": [ - "NFACES", - "NCOMPS", - "NPHASES" - ], - "constants": [], - "combinations": { - "NFACES": [ - 4, - 5, - 6 - ], - "NCOMPS": [ - 1, - 2, - 3, - 4, - 5 - ], - "NPHASES": [ - 2, - 3 - ] - }, - "explicit": [] - }, - - "CompositionalMultiphaseHybridFVMKernels_flux": { - "vars": [ - "NFACES", - "NCOMPS", - "NPHASES", - "IP_TYPE" - ], - "constants": [], - "combinations": { - "NFACES": [ - 4, - 5, - 6 - ], - "NCOMPS": [ - 1, - 2, - 3, - 4, - 5 - ], - "NPHASES": [ - 2, - 3 - ], - "IP_TYPE": [ - "mimeticInnerProduct::TPFAInnerProduct", - "mimeticInnerProduct::BdVLMInnerProduct" - ] - }, - "explicit": [] - }, - "ReactiveCompositionalMultiphaseOBLKernels": { "vars": [ "NAME", diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp index 2d7fcccc661..69259d917d0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp @@ -31,10 +31,12 @@ #include "mesh/ElementRegionManager.hpp" #include "mesh/ObjectManagerBase.hpp" #include "physicsSolvers/PhysicsSolverBaseKernels.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/StencilAccessors.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/PropertyKernelBase.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp" +#include "finiteVolume/MimeticInnerProductDispatch.hpp" namespace geos @@ -43,6 +45,66 @@ namespace geos namespace compositionalMultiphaseHybridFVMKernels { +//namespace internal +/******************************** Kernel switches ********************************/ + +namespace internal +{ + +template< typename T, typename LAMBDA > +void kernelLaunchSelectorFaceSwitch( T value, LAMBDA && lambda ) +{ + static_assert( std::is_integral< T >::value, "KernelLaunchSelectorFaceSwitch: type should be integral" ); + + switch( value ) + { + case 4: + { + return lambda( std::integral_constant< int, 4 >{} ); + } + case 5: + { + return lambda( std::integral_constant< int, 5 >{} ); + } + case 6: + { + return lambda( std::integral_constant< int, 6 >{} ); + } + case 7: + { + return lambda( std::integral_constant< int, 7 >{} ); + } + case 8: + { + return lambda( std::integral_constant< int, 8 >{} ); + } + case 9: + { + return lambda( std::integral_constant< int, 9 >{} ); + } + case 10: + { + return lambda( std::integral_constant< int, 10 >{} ); + } + case 11: + { + return lambda( std::integral_constant< int, 11 >{} ); + } + case 12: + { + return lambda( std::integral_constant< int, 12 >{} ); + } + case 13: + { + return lambda( std::integral_constant< int, 13 >{} ); + } + default: + GEOS_ERROR( "Unknown numFacesInElem value: " << value ); + } +} + +} // namespace internal + // struct to specify local and neighbor derivatives struct Pos { @@ -68,7 +130,7 @@ struct UpwindingHelper /** * @brief At a given one-sided face, compute the upwind viscous transport coefficient * @param[in] localIds region, subRegion, and element indices of the local element - * @param[in] neighborIds region, subRegion, and element indices of the neigbhbor element + * @param[in] neighborIds region, subRegion, and element indices of the neighbor element * @param[in] phaseDens the phase densities in the domain (non-local) * @param[in] dPhaseDens the derivatives of the phase densities in the domain wrt pressure and component fractions (non-local) * @param[in] phaseMob the phase mobilities in the domain (non-local) @@ -106,7 +168,7 @@ struct UpwindingHelper /** * @brief At a given one-sided face, compute the upwind viscous transport coefficient * @param[in] localIds region, subRegion, and element indices of the local element - * @param[in] neighborIds region, subRegion, and element indices of the neigbhbor element + * @param[in] neighborIds region, subRegion, and element indices of the neighbor element * @param[in] transGravCoef * @param[in] phaseDens the phase densities in the domain (non-local) * @param[in] dPhaseDens the derivatives of the phase densities in the domain wrt pressure and component fractions (non-local) @@ -153,7 +215,7 @@ struct UpwindingHelper /** * @brief At a given one-sided face, compute the gravCoef multiplied by the difference in phase densities * @param[in] localIds region, subRegion, and element indices of the local element - * @param[in] neighborIds region, subRegion, and element indices of the neigbhbor element + * @param[in] neighborIds region, subRegion, and element indices of the neighbor element * @param[in] transGravCoef * @param[in] phaseDens the phase densities in the domain (non-local) * @param[in] dPhaseDens the derivatives of the phase densities in the domain wrt pressure and component fractions (non-local) @@ -179,7 +241,7 @@ struct UpwindingHelper /** * @brief At a given one-sided face, compute the upwinded total mobility * @param[in] localIds region, subRegion, and element indices of the local element - * @param[in] neighborIds region, subRegion, and element indices of the neigbhbor element + * @param[in] neighborIds region, subRegion, and element indices of the neighbor element * @param[in] phaseMob the phase mobilities in the domain (non-local) * @param[in] dPhaseMob the derivatives of the phase mobilities in the domain (non-local) * @param[in] phaseGravTerm the gravCoef multiplied by the difference in phase densities @@ -246,6 +308,14 @@ struct UpwindingHelper /******************************** AssemblerKernelHelper ********************************/ +/** + * @struct AssemblerKernelHelper + * @brief Helper structure containing static methods for flux assembly in hybrid FVM + * + * This helper provides methods to compute pressure gradients, flux divergence, and assemble + * face-based continuity equations for the hybrid finite volume method discretization. + * All methods are designed to work on both host and device (GEOS_HOST_DEVICE). + */ struct AssemblerKernelHelper { @@ -312,7 +382,7 @@ struct AssemblerKernelHelper * @param[in] dPhaseMob the derivatives of the phase mobilities in the element * @param[in] dCompFrac_dCompDens the derivatives of the component fractions wrt component density * @param[in] phaseCompFrac the phase component fractions in the domain - * @param[in] dPhaseCompFrac the derivatives of the phase component fractions wrt pressure and component fractions + * @param[in] dPhaseCompFrac the derivatives of the phase component fractions in the domain wrt pressure and component fractions * @param[in] elemDofNumber the dof numbers of the element in the domain * @param[in] oneSidedVolFlux the volumetric fluxes at this element's faces * @param[in] dOneSidedVolFlux_dPres the derivatives of the vol fluxes wrt to this element's cell centered pressure @@ -332,7 +402,9 @@ struct AssemblerKernelHelper arrayView2d< localIndex const > const & elemList, SortedArrayView< localIndex const > const & regionFilter, arrayView1d< globalIndex const > const & faceDofNumber, - arrayView1d< real64 const > const & mimFaceGravCoef, + arrayView1d< integer const > const & isBoundaryFace, + arrayView1d< real64 const > const & mimeticTransGgradZ, + arrayView1d< real64 const > const & faceGravCoef, arraySlice1d< localIndex const > const & elemToFaces, real64 const & elemGravCoef, integer const useTotalMassEquation, @@ -346,7 +418,7 @@ struct AssemblerKernelHelper ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, - arraySlice2d< real64 const > const & transMatrixGrav, + arraySlice2d< real64 const > const & transMatrix, real64 const (&oneSidedVolFlux)[ NF ], real64 const (&dOneSidedVolFlux_dPres)[ NF ], real64 const (&dOneSidedVolFlux_dFacePres)[ NF ][ NF ], @@ -366,13 +438,11 @@ struct AssemblerKernelHelper * @param[in] dUpwPhaseViscCoef_dPres the derivative of the upwinded viscous transport coefficient wrt pressure * @param[in] dUpwPhaseViscCoef_dCompDens the derivative of the upwinded viscous transport coefficient wrt component density * @param[in] upwViscDofNumber degree of freedom number of the upwind element - * @param[in] faceDofNumber degree of freedom number of the face * @param[in] dt the time step * @param[inout] divMassFluxes the divergence of the fluxes in the element * @param[inout] dDivMassFluxes_dElemVars the derivatives of the flux divergence wrt the element centered vars (pres and comp dens) * @param[inout] dDivMassFluxes_dFaceVars the derivatives of the flux divergence wrt the face centered vars * @param[inout] dofColIndicesElemVars degrees of freedom of the cells involved in the flux divergence - * @param[inout] dofColIndicesFaceVars degrees of freedom of the faces involved in the flux divergence */ template< integer NF, integer NC, integer NP > GEOS_HOST_DEVICE @@ -388,13 +458,11 @@ struct AssemblerKernelHelper globalIndex const elemDofNumber, globalIndex const neighborDofNumber, globalIndex const upwViscDofNumber, - globalIndex const faceDofNumber, real64 const & dt, real64 ( &divMassFluxes )[ NC ], real64 ( &dDivMassFluxes_dElemVars )[ NC ][ (NC+1)*(NF+1) ], real64 ( &dDivMassFluxes_dFaceVars )[ NC ][ NF ], - globalIndex ( &dofColIndicesElemVars )[ (NC+1)*(NF+1) ], - globalIndex ( &dofColIndicesFaceVars )[ NF ] ); + globalIndex ( &dofColIndicesElemVars )[ (NC+1)*(NF+1) ] ); /** * @brief In a given element, compute the buoyancy flux divergence, i.e, sum the buoyancy fluxes at this element's faces @@ -445,6 +513,7 @@ struct AssemblerKernelHelper static void assembleFaceConstraints( arrayView1d< globalIndex const > const & faceDofNumber, arrayView1d< integer const > const & faceGhostRank, + arrayView1d< integer const > const & isBoundaryFace, arraySlice1d< localIndex const > const & elemToFaces, globalIndex const elemDofNumber, globalIndex const rankOffset, @@ -517,9 +586,10 @@ struct AssemblerKernel arrayView2d< localIndex const > const & elemList, arrayView1d< globalIndex const > const & faceDofNumber, arrayView1d< integer const > const & faceGhostRank, + arrayView1d< integer const > const & isBoundaryFace, arrayView1d< real64 const > const & facePres, arrayView1d< real64 const > const & faceGravCoef, - arrayView1d< real64 const > const & mimFaceGravCoef, + arrayView1d< real64 const > const & mimeticTransGgradZ, arraySlice1d< localIndex const > const & elemToFaces, real64 const & elemPres, real64 const & elemGravCoef, @@ -538,12 +608,107 @@ struct AssemblerKernel globalIndex const rankOffset, real64 const & dt, arraySlice2d< real64 const > const & transMatrix, - arraySlice2d< real64 const > const & transMatrixGrav, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); }; +/******************************** DirichletFluxKernel ********************************/ + +struct DirichletFluxKernel +{ + + /** + * @brief The type for element-based non-constitutive data parameters. + * Consists entirely of ArrayView's. + * + * Can be converted from ElementRegionManager::ElementViewAccessor + * by calling .toView() or .toViewConst() on an accessor instance + */ + template< typename VIEWTYPE > + using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >; + + using CompFlowAccessors = + StencilAccessors< fields::flow::phaseMobility, + fields::flow::dPhaseMobility, + fields::flow::dGlobalCompFraction_dGlobalCompDensity >; + + using MultiFluidAccessors = + StencilMaterialAccessors< constitutive::MultiFluidBase, + fields::multifluid::phaseDensity, + fields::multifluid::dPhaseDensity, + fields::multifluid::phaseMassDensity, + fields::multifluid::dPhaseMassDensity, + fields::multifluid::phaseCompFraction, + fields::multifluid::dPhaseCompFraction >; + + /** + * @brief Compute the Dirichlet boundary fluxes for compositional multiphase flow + * @tparam NF number of faces per element + * @tparam NC number of components + * @tparam NP number of phases + * @tparam IP_TYPE mimetic inner product type + * @param[in] numPhases number of phases + * @param[in] er index of this element's region + * @param[in] esr index of this element's subregion + * @param[in] subRegion the subRegion + * @param[in] permeabilityModel the permeability model + * @param[in] boundaryFaceSet set of faces on the Dirichlet boundary + * @param[in] facePres prescribed face pressures + * @param[in] faceTemp prescribed face temperatures + * @param[in] facePhaseMob prescribed face phase mobilities + * @param[in] facePhaseMassDens prescribed face phase mass densities + * @param[in] facePhaseCompFrac prescribed face phase component fractions + * @param[in] nodePosition node positions + * @param[in] faceToNodes face to nodes map + * @param[in] elemRegionList face to element region map + * @param[in] elemSubRegionList face to element subregion map + * @param[in] elemList face to element index map + * @param[in] faceDofNumber face DOF numbers + * @param[in] transMultiplier transmissibility multipliers + * @param[in] lengthTolerance geometric tolerance + * @param[in] dt time step size + * @param[in] rankOffset MPI rank offset + * @param[in] useTotalMassEquation flag for total mass equation + * @param[in] compFlowAccessors accessor for comp flow fields + * @param[in] multiFluidAccessors accessor for fluid properties + * @param[in] elemDofNumber element DOF numbers + * @param[inout] localMatrix the local matrix + * @param[inout] localRhs the local RHS vector + */ + template< integer NF, integer NC, integer NP, typename IP_TYPE > + static void + launch( integer const numPhases, + localIndex const er, + localIndex const esr, + CellElementSubRegion const & subRegion, + constitutive::PermeabilityBase const & permeabilityModel, + SortedArrayView< localIndex const > const & boundaryFaceSet, + arrayView1d< real64 const > const & facePres, + arrayView1d< real64 const > const & faceTemp, + arrayView2d< real64 const, compflow::USD_PHASE > const & facePhaseMob, + arrayView2d< real64 const, compflow::USD_PHASE > const & facePhaseMassDens, + arrayView3d< real64 const, compflow::USD_PHASE_COMP > const & facePhaseCompFrac, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arrayView2d< localIndex const > const & elemRegionList, + arrayView2d< localIndex const > const & elemSubRegionList, + arrayView2d< localIndex const > const & elemList, + arrayView1d< globalIndex const > const & faceDofNumber, + arrayView1d< real64 const > const & faceGravCoef, + arrayView1d< real64 const > const & transMultiplier, + real64 const lengthTolerance, + real64 const dt, + globalIndex const rankOffset, + integer const useTotalMassEquation, + CompFlowAccessors const & compFlowAccessors, + MultiFluidAccessors const & multiFluidAccessors, + ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + +}; + /******************************** FluxKernel ********************************/ struct FluxKernel @@ -605,6 +770,7 @@ struct FluxKernel * @param[in] rankOffset the offset of this rank * @param[in] lengthTolerance tolerance used in the transmissibility matrix computation * @param[in] dt time step size + * @param[in] transMatrix the transmissibility matrix in this element * @param[inout] matrix the system matrix * @param[inout] rhs the system right-hand side vector */ @@ -621,9 +787,10 @@ struct FluxKernel ArrayOfArraysView< localIndex const > const & faceToNodes, arrayView1d< globalIndex const > const & faceDofNumber, arrayView1d< integer const > const & faceGhostRank, + arrayView1d< integer const > const & isBoundaryFace, arrayView1d< real64 const > const & facePres, arrayView1d< real64 const > const & faceGravCoef, - arrayView1d< real64 const > const & mimFaceGravCoef, + arrayView1d< real64 const > const & mimeticTransGgradZ, arrayView1d< real64 const > const & transMultiplier, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, @@ -1070,33 +1237,32 @@ struct SolutionCheckKernel }; -/******************************** PrecomputeKernel ********************************/ +/******************************** PrecomputeMimeticTransGgradZKernel ********************************/ -struct PrecomputeKernel +struct PrecomputeMimeticTransGgradZKernel { template< typename IP_TYPE, integer NF > static void launch( localIndex const subRegionSize, - localIndex const faceManagerSize, arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, ArrayOfArraysView< localIndex const > const & faceToNodes, arrayView2d< real64 const > const & elemCenter, arrayView1d< real64 const > const & elemVolume, arrayView3d< real64 const > const & elemPerm, arrayView1d< real64 const > const & elemGravCoef, + arrayView1d< real64 const > const & faceGravCoef, arrayView2d< localIndex const > const & elemToFaces, arrayView1d< real64 const > const & transMultiplier, real64 const & lengthTolerance, - arrayView1d< RAJA::ReduceSum< serialReduce, real64 > > const & mimFaceGravCoefNumerator, - arrayView1d< RAJA::ReduceSum< serialReduce, real64 > > const & mimFaceGravCoefDenominator, - arrayView1d< real64 > const & mimFaceGravCoef ) + arrayView1d< real64 > const & faceInvSum, + arrayView1d< integer > const & faceCount ) { - forAll< serialPolicy >( subRegionSize, [=] ( localIndex const ei ) + forAll< parallelDevicePolicy<> >( subRegionSize, [=] GEOS_HOST_DEVICE ( localIndex const ei ) { - stackArray2d< real64, NF *NF > transMatrix( NF, NF ); + stackArray2d< real64, NF * NF > transMatrix( NF, NF ); - real64 const perm[ 3 ] = { elemPerm[ei][0][0], elemPerm[ei][0][1], elemPerm[ei][0][2] }; + real64 const perm[3] = { elemPerm[ei][0][0], elemPerm[ei][0][1], elemPerm[ei][0][2] }; IP_TYPE::template compute< NF >( nodePosition, transMultiplier, @@ -1108,46 +1274,127 @@ struct PrecomputeKernel lengthTolerance, transMatrix ); - for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + real64 const ccGravCoef = elemGravCoef[ei]; + for( integer i = 0; i < NF; ++i ) { - mimFaceGravCoefNumerator[elemToFaces[ei][ifaceLoc]] += elemGravCoef[ei] * transMatrix[ifaceLoc][ifaceLoc]; - mimFaceGravCoefDenominator[elemToFaces[ei][ifaceLoc]] += transMatrix[ifaceLoc][ifaceLoc]; + localIndex const kf = elemToFaces[ei][i]; + + real64 T_g_delta_z = 0.0; + for( integer j = 0; j < NF; ++j ) + { + real64 const fGravCoef = faceGravCoef[elemToFaces[ei][j]]; + real64 const gravCoefDif = ccGravCoef - fGravCoef; + T_g_delta_z += transMatrix[i][j] * gravCoefDif; + } + RAJA::atomicAdd( parallelDeviceAtomic{}, &faceInvSum[kf], 1.0 / LvArray::math::abs( T_g_delta_z ) ); + RAJA::atomicAdd( parallelDeviceAtomic{}, &faceCount[kf], 1 ); } } ); + } - forAll< serialPolicy >( faceManagerSize, [=] ( localIndex const iface ) + template< typename POLICY > + static void + createAndLaunch( mimeticInnerProduct::MimeticInnerProductBase const & ip, + NodeManager const & nodeManager, + FaceManager const & faceManager, + CellElementSubRegion const & subRegion, + constitutive::PermeabilityBase const & permeability, + real64 const lengthTolerance, + arrayView1d< real64 > const & faceInvSum, + arrayView1d< integer > const & faceCount ) + { + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + ArrayOfArraysView< localIndex const > const & faceToNodes = faceManager.nodeList().toViewConst(); + arrayView1d< real64 const > const & transMultiplier = faceManager.getField< fields::flow::transMultiplier >(); + arrayView1d< real64 const > const & faceGravCoef = faceManager.getField< fields::flow::gravityCoefficient >(); + + arrayView2d< real64 const > const & elemCenter = subRegion.getElementCenter(); + arrayView1d< real64 const > const & elemVolume = subRegion.getElementVolume(); + arrayView3d< real64 const > const & elemPerm = permeability.permeability(); + arrayView1d< real64 const > const & elemGravCoef = subRegion.getField< fields::flow::gravityCoefficient >(); + arrayView2d< localIndex const > const & elemToFaces = subRegion.faceList(); + + mimeticInnerProductDispatch( ip, [&]( auto const & innerProduct ) { - if( !isZero( mimFaceGravCoefDenominator[iface].get() ) ) + using IP_TYPE = std::remove_const_t< TYPEOFREF( innerProduct ) >; + internal::kernelLaunchSelectorFaceSwitch( subRegion.numFacesPerElement(), [&]( auto NF ) { - mimFaceGravCoef[iface] = mimFaceGravCoefNumerator[iface].get() / mimFaceGravCoefDenominator[iface].get(); - } + PrecomputeMimeticTransGgradZKernel::launch< IP_TYPE, NF() >( subRegion.size(), + nodePosition, + faceToNodes, + elemCenter, + elemVolume, + elemPerm, + elemGravCoef, + faceGravCoef, + elemToFaces, + transMultiplier, + lengthTolerance, + faceInvSum, + faceCount ); + } ); } ); } }; /******************************** Kernel switches ********************************/ -namespace internal -{ - -template< typename T, typename LAMBDA > -void kernelLaunchSelectorFaceSwitch( T value, LAMBDA && lambda ) -{ - static_assert( std::is_integral< T >::value, "KernelLaunchSelectorFaceSwitch: type should be integral" ); - - switch( value ) - { - case 4: - { lambda( std::integral_constant< T, 4 >() ); return;} - case 5: - { lambda( std::integral_constant< T, 5 >() ); return;} - case 6: - { lambda( std::integral_constant< T, 6 >() ); return;} - default: GEOS_ERROR( "Unknown numFacesInElem value: " << value ); - } -} - -} // namespace internal +//namespace internal +//{ +// +//template< typename T, typename LAMBDA > +//void kernelLaunchSelectorFaceSwitch( T value, LAMBDA && lambda ) +//{ +// static_assert( std::is_integral< T >::value, "KernelLaunchSelectorFaceSwitch: type should be integral" ); +// +// switch( value ) +// { +// case 4: +// { +// return lambda( std::integral_constant< int, 4 >{} ); +// } +// case 5: +// { +// return lambda( std::integral_constant< int, 5 >{} ); +// } +// case 6: +// { +// return lambda( std::integral_constant< int, 6 >{} ); +// } +// case 7: +// { +// return lambda( std::integral_constant< int, 7 >{} ); +// } +// case 8: +// { +// return lambda( std::integral_constant< int, 8 >{} ); +// } +// case 9: +// { +// return lambda( std::integral_constant< int, 9 >{} ); +// } +// case 10: +// { +// return lambda( std::integral_constant< int, 10 >{} ); +// } +// case 11: +// { +// return lambda( std::integral_constant< int, 11 >{} ); +// } +// case 12: +// { +// return lambda( std::integral_constant< int, 12 >{} ); +// } +// case 13: +// { +// return lambda( std::integral_constant< int, 13 >{} ); +// } +// default: +// GEOS_ERROR( "Unknown numFacesInElem value: " << value ); +// } +//} +// +//} // namespace internal template< typename KERNELWRAPPER, typename INNER_PRODUCT, typename ... ARGS > void simpleKernelLaunchSelector( localIndex numFacesInElem, ARGS && ... args ) @@ -1223,6 +1470,46 @@ void kernelLaunchSelector( integer numFacesInElem, integer numComps, integer num } } +/******************************** EvaluateBCFacePropertiesKernel ********************************/ + +/** + * @brief Evaluate constitutive properties at BC face conditions + * @tparam NC number of components + * @tparam NP number of phases + * @param[in] numPhases number of phases + * @param[in] boundaryFaceSet set of boundary faces + * @param[in] facePres face pressures at BC + * @param[in] faceTemp face temperatures at BC + * @param[in] faceCompFrac face component fractions at BC + * @param[in] elemRegionList face to element region list + * @param[in] elemSubRegionList face to element subregion list + * @param[in] elemList face to element list + * @param[in] er target element region index + * @param[in] esr target element subregion index + * @param[in] fluid multifluid model + * @param[in] relperm relative permeability model + * @param[out] facePhaseMob output: face phase mobility at BC + * @param[out] facePhaseMassDens output: face phase mass density at BC + * @param[out] facePhaseCompFrac output: face phase component fraction at BC + */ +template< integer NC, integer NP > +void +evaluateBCFaceProperties( integer const numPhases, + SortedArrayView< localIndex const > const & boundaryFaceSet, + arrayView1d< real64 const > const & facePres, + arrayView1d< real64 const > const & faceTemp, + arrayView2d< real64 const, compflow::USD_COMP > const & faceCompFrac, + arrayView2d< localIndex const > const & elemRegionList, + arrayView2d< localIndex const > const & elemSubRegionList, + arrayView2d< localIndex const > const & elemList, + localIndex const er, + localIndex const esr, + constitutive::MultiFluidBase & fluid, + constitutive::RelativePermeabilityBase & relperm, + arrayView2d< real64, compflow::USD_PHASE > const & facePhaseMob, + arrayView2d< real64, compflow::USD_PHASE > const & facePhaseMassDens, + arrayView3d< real64, compflow::USD_PHASE_COMP > const & facePhaseCompFrac ); + } // namespace compositionalMultiphaseHybridFVMKernels } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernelsInstantiations.cpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernelsInstantiations.cpp new file mode 100644 index 00000000000..351724f9478 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernelsInstantiations.cpp @@ -0,0 +1,247 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CompositionalMultiphaseHybridFVMKernelsInstantiations.cpp + * @brief Centralized explicit template instantiations for compositional hybrid FVM kernels. + * + * Define GEOS_ENABLE_MANUAL_HYBRID_FVM_INST to activate the explicit instantiations below. + * When enabled, ensure overlapping template-generated sources are not compiled to avoid duplicates. + */ + +#include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels_impl.hpp" + +#include "finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp" +#include "finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp" +#include "finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp" + +namespace geos +{ +namespace compositionalMultiphaseHybridFVMKernels +{ + +#if defined(GEOS_ENABLE_MANUAL_HYBRID_FVM_INST) + +// ----------------------------------------------------------------------------- +// Range driver +// ----------------------------------------------------------------------------- +#define GEOS_NFACES_LIST \ + X_NF( 4 ) \ + X_NF( 5 ) \ + X_NF( 6 ) \ + X_NF( 7 ) \ + X_NF( 8 ) \ + X_NF( 9 ) \ + X_NF( 10 ) \ + X_NF( 11 ) \ + X_NF( 12 ) \ + X_NF( 13 ) + +// ----------------------------------------------------------------------------- +// Common parameter macros (kept in sync with kernel signatures) +// ----------------------------------------------------------------------------- +#define GEOS_DIRICHLET_LAUNCH_PARAMS \ + integer const numPhases, \ + localIndex const er, \ + localIndex const esr, \ + CellElementSubRegion const & subRegion, \ + constitutive::PermeabilityBase const & permeabilityModel, \ + SortedArrayView< localIndex const > const & boundaryFaceSet, \ + arrayView1d< real64 const > const & facePres, \ + arrayView1d< real64 const > const & faceTemp, \ + arrayView2d< real64 const, compflow::USD_PHASE > const & facePhaseMob, \ + arrayView2d< real64 const, compflow::USD_PHASE > const & facePhaseMassDens, \ + arrayView3d< real64 const, compflow::USD_PHASE_COMP > const & facePhaseCompFrac, \ + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, \ + ArrayOfArraysView< localIndex const > const & faceToNodes, \ + arrayView2d< localIndex const > const & elemRegionList, \ + arrayView2d< localIndex const > const & elemSubRegionList, \ + arrayView2d< localIndex const > const & elemList, \ + arrayView1d< globalIndex const > const & faceDofNumber, \ + arrayView1d< real64 const > const & faceGravCoef, \ + arrayView1d< real64 const > const & transMultiplier, \ + real64 const lengthTolerance, \ + real64 const dt, \ + globalIndex const rankOffset, \ + integer const useTotalMassEquation, \ + DirichletFluxKernel::CompFlowAccessors const & compFlowAccessors, \ + DirichletFluxKernel::MultiFluidAccessors const & multiFluidAccessors, \ + DirichletFluxKernel::ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, \ + CRSMatrixView< real64, globalIndex const > const & localMatrix, \ + arrayView1d< real64 > const & localRhs + +#define GEOS_FLUX_LAUNCH_PARAMS \ + localIndex er, localIndex esr, \ + CellElementSubRegion const & subRegion, \ + constitutive::PermeabilityBase const & permeabilityModel, \ + SortedArrayView< localIndex const > const & regionFilter, \ + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, \ + arrayView2d< localIndex const > const & elemRegionList, \ + arrayView2d< localIndex const > const & elemSubRegionList, \ + arrayView2d< localIndex const > const & elemList, \ + ArrayOfArraysView< localIndex const > const & faceToNodes, \ + arrayView1d< globalIndex const > const & faceDofNumber, \ + arrayView1d< integer const > const & faceGhostRank, \ + arrayView1d< integer const > const & isBoundaryFace, \ + arrayView1d< real64 const > const & facePres, \ + arrayView1d< real64 const > const & faceGravCoef, \ + arrayView1d< real64 const > const & mimFaceGravCoef, \ + arrayView1d< real64 const > const & transMultiplier, \ + FluxKernel::ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, \ + FluxKernel::ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, \ + FluxKernel::ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, \ + FluxKernel::ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, \ + FluxKernel::ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, \ + FluxKernel::ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, \ + FluxKernel::ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseMassDens, \ + FluxKernel::ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac, \ + FluxKernel::ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, \ + FluxKernel::ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, \ + globalIndex const rankOffset, \ + real64 const lengthTolerance, \ + real64 const dt, \ + integer const useTotalMassEquation, \ + CRSMatrixView< real64, globalIndex const > const & localMatrix, \ + arrayView1d< real64 > const & localRhs + +// ----------------------------------------------------------------------------- +// IP list helper +// ----------------------------------------------------------------------------- +#define GEOS_FOR_EACH_IP( NF, NC, NP, MACRO ) \ + MACRO( NF, NC, NP, mimeticInnerProduct::TPFAInnerProduct ) \ + MACRO( NF, NC, NP, mimeticInnerProduct::QuasiTPFAInnerProduct ) \ + MACRO( NF, NC, NP, mimeticInnerProduct::BdVLMInnerProduct ) + +// ----------------------------------------------------------------------------- +// DirichletFluxKernel explicit instantiations +// ----------------------------------------------------------------------------- +#define INSTANTIATE_DIRICHLET( NF, NC, NP, IP ) \ + template void DirichletFluxKernel::launch< NF, NC, NP, IP >( GEOS_DIRICHLET_LAUNCH_PARAMS ); + +#define INSTANTIATE_DIRICHLET_FOR_NP( NF, NC, NP ) \ + GEOS_FOR_EACH_IP( NF, NC, NP, INSTANTIATE_DIRICHLET ) + +#define INSTANTIATE_DIRICHLET_FOR_NC( NF, NC ) \ + INSTANTIATE_DIRICHLET_FOR_NP( NF, NC, 2 ) \ + INSTANTIATE_DIRICHLET_FOR_NP( NF, NC, 3 ) + +#define INSTANTIATE_DIRICHLET_FOR_NF( NF ) \ + INSTANTIATE_DIRICHLET_FOR_NC( NF, 1 ) \ + INSTANTIATE_DIRICHLET_FOR_NC( NF, 2 ) \ + INSTANTIATE_DIRICHLET_FOR_NC( NF, 3 ) \ + INSTANTIATE_DIRICHLET_FOR_NC( NF, 4 ) \ + INSTANTIATE_DIRICHLET_FOR_NC( NF, 5 ) + +#define X_NF( NF ) INSTANTIATE_DIRICHLET_FOR_NF( NF ) +GEOS_NFACES_LIST +#undef X_NF + +// ----------------------------------------------------------------------------- +// FluxKernel explicit instantiations +// ----------------------------------------------------------------------------- +#define INSTANTIATE_FLUX( NF, NC, NP, IP ) \ + template void FluxKernel::launch< NF, NC, NP, IP >( GEOS_FLUX_LAUNCH_PARAMS ); + +#define INSTANTIATE_FLUX_FOR_NP( NF, NC, NP ) \ + GEOS_FOR_EACH_IP( NF, NC, NP, INSTANTIATE_FLUX ) + +#define INSTANTIATE_FLUX_FOR_NC( NF, NC ) \ + INSTANTIATE_FLUX_FOR_NP( NF, NC, 2 ) \ + INSTANTIATE_FLUX_FOR_NP( NF, NC, 3 ) + +#define INSTANTIATE_FLUX_FOR_NF( NF ) \ + INSTANTIATE_FLUX_FOR_NC( NF, 1 ) \ + INSTANTIATE_FLUX_FOR_NC( NF, 2 ) \ + INSTANTIATE_FLUX_FOR_NC( NF, 3 ) \ + INSTANTIATE_FLUX_FOR_NC( NF, 4 ) \ + INSTANTIATE_FLUX_FOR_NC( NF, 5 ) + +#define X_NF( NF ) INSTANTIATE_FLUX_FOR_NF( NF ) +GEOS_NFACES_LIST +#undef X_NF + +// Cleanup local helper macros +#undef GEOS_DIRICHLET_LAUNCH_PARAMS +#undef GEOS_FLUX_LAUNCH_PARAMS +#undef GEOS_FOR_EACH_IP +#undef INSTANTIATE_DIRICHLET +#undef INSTANTIATE_DIRICHLET_FOR_NP +#undef INSTANTIATE_DIRICHLET_FOR_NC +#undef INSTANTIATE_DIRICHLET_FOR_NF +#undef INSTANTIATE_FLUX +#undef INSTANTIATE_FLUX_FOR_NP +#undef INSTANTIATE_FLUX_FOR_NC +#undef INSTANTIATE_FLUX_FOR_NF + +// ----------------------------------------------------------------------------- +// evaluateBCFaceProperties explicit instantiations +// ----------------------------------------------------------------------------- +#define INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES( NC, NP ) \ + template void evaluateBCFaceProperties< NC, NP >( \ + integer const numPhases, \ + SortedArrayView< localIndex const > const & boundaryFaceSet, \ + arrayView1d< real64 const > const & facePres, \ + arrayView1d< real64 const > const & faceTemp, \ + arrayView2d< real64 const, compflow::USD_COMP > const & faceCompFrac, \ + arrayView2d< localIndex const > const & elemRegionList, \ + arrayView2d< localIndex const > const & elemSubRegionList, \ + arrayView2d< localIndex const > const & elemList, \ + localIndex const er, \ + localIndex const esr, \ + constitutive::MultiFluidBase & fluid, \ + constitutive::RelativePermeabilityBase & relperm, \ + arrayView2d< real64, compflow::USD_PHASE > const & facePhaseMob, \ + arrayView2d< real64, compflow::USD_PHASE > const & facePhaseMassDens, \ + arrayView3d< real64, compflow::USD_PHASE_COMP > const & facePhaseCompFrac ); + +#define INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NP( NC, NP ) \ + INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES( NC, NP ) + +#define INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NC( NC ) \ + INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NP( NC, 2 ) \ + INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NP( NC, 3 ) + +INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NC( 1 ) +INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NC( 2 ) +INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NC( 3 ) +INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NC( 4 ) +INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NC( 5 ) + +// Cleanup local helper macros +#undef INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES +#undef INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NP +#undef INSTANTIATE_EVALUATE_BC_FACE_PROPERTIES_FOR_NC + +// ----------------------------------------------------------------------------- +// NOTE: Explicit instantiations for AssemblerKernelHelper and AssemblerKernel +// functions are NOT included here because they are declared as 'inline' in +// CompositionalMultiphaseHybridFVMKernels_impl.hpp. Inline template functions +// are automatically instantiated at the point of use and cannot be explicitly +// instantiated. +// +// The following functions are inline and therefore auto-instantiated: +// - AssemblerKernelHelper::applyGradient +// - AssemblerKernelHelper::assembleFluxDivergence +// - AssemblerKernelHelper::assembleViscousFlux +// - AssemblerKernelHelper::assembleBuoyancyFlux +// - AssemblerKernelHelper::assembleFaceConstraints +// - AssemblerKernel::compute +// ----------------------------------------------------------------------------- + +#endif // GEOS_ENABLE_MANUAL_HYBRID_FVM_INST + +} // namespace compositionalMultiphaseHybridFVMKernels +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels_impl.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels_impl.hpp index 91d6475b76d..3075c22a3e9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels_impl.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CompositionalMultiphaseHybridFVMKernels_impl.hpp @@ -20,12 +20,15 @@ #include "CompositionalMultiphaseHybridFVMKernels.hpp" #include "finiteVolume/mimeticInnerProducts/MimeticInnerProductBase.hpp" -#include "finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp" #include "finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp" +#include "finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp" +#include "finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseUtilities.hpp" #include "physicsSolvers/fluidFlow/kernels/HybridFVMHelperKernels.hpp" +#include "constitutive/fluid/multifluid/MultiFluidSelector.hpp" +#include "constitutive/relativePermeability/RelativePermeabilitySelector.hpp" namespace geos { @@ -446,7 +449,7 @@ UpwindingHelper:: localIndex ( & totalMobIds )[ NP ][ 3 ], localIndex ( & totalMobPos )[ NP ] ) { - if( NP == 2 ) + if constexpr ( NP == 2 ) { if( gravTerm[0][0] > 0 ) { @@ -471,7 +474,7 @@ UpwindingHelper:: totalMobPos[1] = Pos::LOCAL; } } - else if( NP == 3 ) + else if constexpr ( NP == 3 ) { // TODO Francois: this should be improved // currently this implements the algorithm proposed by SH Lee @@ -615,21 +618,23 @@ AssemblerKernelHelper:: arrayView2d< localIndex const > const & elemList, SortedArrayView< localIndex const > const & regionFilter, arrayView1d< globalIndex const > const & faceDofNumber, - arrayView1d< real64 const > const & mimFaceGravCoef, + arrayView1d< integer const > const & isBoundaryFace, + arrayView1d< real64 const > const & mimeticTransGgradZ, + arrayView1d< real64 const > const & faceGravCoef, arraySlice1d< localIndex const > const & elemToFaces, real64 const & elemGravCoef, integer const useTotalMassEquation, - ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseDens, - ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseDens, - ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, - ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_DC > > const & dPhaseMassDens, + ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseDens, + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseDens, + ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, - ElementViewConst< arrayView4d< real64 const, multifluid::USD_PHASE_COMP > > const & phaseCompFrac, - ElementViewConst< arrayView5d< real64 const, multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac, + ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac, ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, - arraySlice2d< real64 const > const & transMatrixGrav, + arraySlice2d< real64 const > const & transMatrix, real64 const (&oneSidedVolFlux)[ NF ], real64 const (&dOneSidedVolFlux_dPres)[ NF ], real64 const (&dOneSidedVolFlux_dFacePres)[ NF ][ NF ], @@ -644,6 +649,8 @@ AssemblerKernelHelper:: // dof numbers globalIndex dofColIndicesElemVars[ NDOF*(NF+1) ]{}; globalIndex dofColIndicesFaceVars[ NF ]{}; + integer faceIndexMap[ NF ]{}; // Maps local face index to position in compact DOF array + integer numNonBoundaryFaces = 0; for( integer idof = 0; idof < NDOF; ++idof ) { dofColIndicesElemVars[idof] = elemDofNumber[localIds[0]][localIds[1]][localIds[2]] + idof; @@ -676,95 +683,119 @@ AssemblerKernelHelper:: for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) { - // 1) Find if there is a neighbor, and if there is, grab the indices of the neighbor element - - localIndex neighborIds[ 3 ] = { localIds[0], localIds[1], localIds[2] }; - hybridFVMKernels::CellConnectivity::isNeighborFound( localIds, - ifaceLoc, - elemRegionList, - elemSubRegionList, - elemList, - regionFilter, - elemToFaces, - neighborIds ); - localIndex const neighborDofNumber = elemDofNumber[neighborIds[0]][neighborIds[1]][neighborIds[2]]; - - // 2) *************** Assemble viscous terms ****************** - - // 2.a) Compute the upwinded x_{c, \ell} \rho_{\ell} \frac{\lambda_{\ell}}{\lambda_T} for each phase at this face - UpwindingHelper::upwindViscousCoefficient< NC, NP >( localIds, - neighborIds, - phaseDens, - dPhaseDens, - phaseMob, - dPhaseMob, - dCompFrac_dCompDens, - phaseCompFrac, - dPhaseCompFrac, - elemDofNumber, - oneSidedVolFlux[ifaceLoc], - upwPhaseViscCoef, - dUpwPhaseViscCoef_dPres, - dUpwPhaseViscCoef_dCompDens, - upwViscDofNumber ); - - // 2.b) Add the \x_{c,\ell} \rho_{\ell} \frac{\lambda_{\ell}}{\lambda_T} q_T of this face to the divergence of the flux in this cell - assembleViscousFlux< NF, NC, NP >( ifaceLoc, - oneSidedVolFlux, - dOneSidedVolFlux_dPres, - dOneSidedVolFlux_dFacePres, - dOneSidedVolFlux_dCompDens, - upwPhaseViscCoef, - dUpwPhaseViscCoef_dPres, - dUpwPhaseViscCoef_dCompDens, - elemDofNumber[localIds[0]][localIds[1]][localIds[2]], - neighborDofNumber, - upwViscDofNumber, - faceDofNumber[elemToFaces[ifaceLoc]], - dt, - divMassFluxes, - dDivMassFluxes_dElemVars, - dDivMassFluxes_dFaceVars, - dofColIndicesElemVars, - dofColIndicesFaceVars ); - - // 3) *************** Assemble buoyancy terms ****************** - - real64 const transGravCoef = (localIds[0] != neighborIds[0] || localIds[1] != neighborIds[1] || localIds[2] != neighborIds[2]) - * transMatrixGrav[ifaceLoc][ifaceLoc] * (elemGravCoef - mimFaceGravCoef[elemToFaces[ifaceLoc]]); - - // 3.a) Compute the upwinded x_{c, \ell} \rho_{\ell} \frac{\lambda_{\ell}\lambda_m}{\lambda_T} - // and (\rho_{\ell} - \rho_m) g \Delta z for each phase at this face - UpwindingHelper::upwindBuoyancyCoefficient< NC, NP >( localIds, - neighborIds, - transGravCoef, - phaseDens, - dPhaseDens, - phaseMassDens, - dPhaseMassDens, - phaseMob, - dPhaseMob, - dCompFrac_dCompDens, - phaseCompFrac, - dPhaseCompFrac, - phaseGravTerm, - dPhaseGravTerm_dPres, - dPhaseGravTerm_dCompDens, - upwPhaseGravCoef, - dUpwPhaseGravCoef_dPres, - dUpwPhaseGravCoef_dCompDens ); - - // 3.b) Add the buoyancy term of this face to the divergence of the flux in this cell - assembleBuoyancyFlux< NF, NC, NP >( ifaceLoc, - phaseGravTerm, - dPhaseGravTerm_dPres, - dPhaseGravTerm_dCompDens, - upwPhaseGravCoef, - dUpwPhaseGravCoef_dPres, - dUpwPhaseGravCoef_dCompDens, - dt, - divMassFluxes, - dDivMassFluxes_dElemVars ); + // Contribute and Collect face DOF number only if this is not a boundary face + if( isBoundaryFace[elemToFaces[ifaceLoc]] == 0 ) + { + + faceIndexMap[ifaceLoc] = numNonBoundaryFaces; + dofColIndicesFaceVars[numNonBoundaryFaces] = faceDofNumber[elemToFaces[ifaceLoc]]; + numNonBoundaryFaces++; + + + // 1) Find if there is a neighbor, and if there is, grab the indices of the neighbor element + localIndex neighborIds[ 3 ] = { localIds[0], localIds[1], localIds[2] }; + hybridFVMKernels::CellConnectivity::isNeighborFound( localIds, + ifaceLoc, + elemRegionList, + elemSubRegionList, + elemList, + regionFilter, + elemToFaces, + neighborIds ); + localIndex const neighborDofNumber = elemDofNumber[neighborIds[0]][neighborIds[1]][neighborIds[2]]; + + // 2) *************** Assemble viscous terms ****************** + // 2.a) Compute the upwinded x_{c, \ell} \rho_{\ell} \frac{\lambda_{\ell}}{\lambda_T} for each phase at this face + UpwindingHelper::upwindViscousCoefficient< NC, NP >( localIds, + neighborIds, + phaseDens, + dPhaseDens, + phaseMob, + dPhaseMob, + dCompFrac_dCompDens, + phaseCompFrac, + dPhaseCompFrac, + elemDofNumber, + oneSidedVolFlux[ifaceLoc], + upwPhaseViscCoef, + dUpwPhaseViscCoef_dPres, + dUpwPhaseViscCoef_dCompDens, + upwViscDofNumber ); + + // 2.b) Add the \x_{c,\ell} \rho_{\ell} \frac{\lambda_{\ell}}{\lambda_T} q_T of this face to the divergence of the flux in this cell + assembleViscousFlux< NF, NC, NP >( ifaceLoc, + oneSidedVolFlux, + dOneSidedVolFlux_dPres, + dOneSidedVolFlux_dFacePres, + dOneSidedVolFlux_dCompDens, + upwPhaseViscCoef, + dUpwPhaseViscCoef_dPres, + dUpwPhaseViscCoef_dCompDens, + elemDofNumber[localIds[0]][localIds[1]][localIds[2]], + neighborDofNumber, + upwViscDofNumber, + dt, + divMassFluxes, + dDivMassFluxes_dElemVars, + dDivMassFluxes_dFaceVars, + dofColIndicesElemVars ); + + // 3) *************** Assemble buoyancy terms ****************** + + // The sign of gFlux_i determines direction for buoyancy w.r.t the consistent inner product. + real64 gFlux_i = 0.0; + for( integer jfaceLoc=0; jfaceLoc 0.0 ) - integer( gFlux_i < 0.0 ); + + // Buoyancy transmissibility applied on ifaceLoc face. The boolean comparison acts + // as a 0/1 mask that disables the buoyancy flux for domain boundaries. + real64 const transGravCoef = sign * (localIds[0] != neighborIds[0] || localIds[1] != neighborIds[1] || localIds[2] != neighborIds[2]) * mimeticTransGgradZ[elemToFaces[ifaceLoc]]; + + // 3.a) Compute the upwinded x_{c, \ell} \rho_{\ell} \frac{\lambda_{\ell}\lambda_m}{\lambda_T} + // and (\rho_{\ell} - \rho_m) g \Delta z for each phase at this face + UpwindingHelper::upwindBuoyancyCoefficient< NC, NP >( localIds, + neighborIds, + transGravCoef, + phaseDens, + dPhaseDens, + phaseMassDens, + dPhaseMassDens, + phaseMob, + dPhaseMob, + dCompFrac_dCompDens, + phaseCompFrac, + dPhaseCompFrac, + phaseGravTerm, + dPhaseGravTerm_dPres, + dPhaseGravTerm_dCompDens, + upwPhaseGravCoef, + dUpwPhaseGravCoef_dPres, + dUpwPhaseGravCoef_dCompDens ); + + // 3.b) Add the buoyancy term of this face to the divergence of the flux in this cell + assembleBuoyancyFlux< NF, NC, NP >( ifaceLoc, + phaseGravTerm, + dPhaseGravTerm_dPres, + dPhaseGravTerm_dCompDens, + upwPhaseGravCoef, + dUpwPhaseGravCoef_dPres, + dUpwPhaseGravCoef_dCompDens, + dt, + divMassFluxes, + dDivMassFluxes_dElemVars ); + + } + else + { + faceIndexMap[ifaceLoc] = -1; // Mark boundary faces + } } @@ -792,16 +823,56 @@ AssemblerKernelHelper:: localRhs[eqnRowLocalIndex] = localRhs[eqnRowLocalIndex] + divMassFluxes[ic]; // jacobian -- derivative wrt elem centered vars - localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( eqnRowLocalIndex, - &dofColIndicesElemVars[0], - &dDivMassFluxes_dElemVars[0][0] + ic * NDOF * (NF+1), - NDOF * (NF+1) ); + // Need to compact both DOF indices and derivatives to only include non-boundary neighbors + globalIndex compactElemDofs[ NDOF*(NF+1) ]; + real64 compactElemDerivs[ NDOF*(NF+1) ]; + + // Copy current element DOFs and derivatives + for( integer idof = 0; idof < NDOF; ++idof ) + { + compactElemDofs[idof] = dofColIndicesElemVars[idof]; + compactElemDerivs[idof] = dDivMassFluxes_dElemVars[ic][idof]; + } + + // Copy neighbor DOFs and derivatives only for non-boundary faces + integer compactIdx = NDOF; + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + { + if( faceIndexMap[jfaceLoc] >= 0 ) // non-boundary face + { + integer const neighborOffset = NDOF * (jfaceLoc + 1); + for( integer idof = 0; idof < NDOF; ++idof ) + { + compactElemDofs[compactIdx] = dofColIndicesElemVars[neighborOffset + idof]; + compactElemDerivs[compactIdx] = dDivMassFluxes_dElemVars[ic][neighborOffset + idof]; + compactIdx++; + } + } + } - // jacobian -- derivatives wrt face centered vars localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( eqnRowLocalIndex, - &dofColIndicesFaceVars[0], - &dDivMassFluxes_dFaceVars[0][0] + ic * NF, - NF ); + compactElemDofs, + compactElemDerivs, + compactIdx ); + + // jacobian -- derivatives wrt face centered vars (only non-boundary faces) + // Need to compact the derivatives array to only include non-boundary faces + if( numNonBoundaryFaces > 0 ) + { + real64 compactFaceDerivs[ NF ]{}; + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + { + if( faceIndexMap[jfaceLoc] >= 0 ) + { + compactFaceDerivs[faceIndexMap[jfaceLoc]] = dDivMassFluxes_dFaceVars[ic][jfaceLoc]; + } + } + + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( eqnRowLocalIndex, + &dofColIndicesFaceVars[0], + compactFaceDerivs, + numNonBoundaryFaces ); + } } } @@ -821,13 +892,11 @@ AssemblerKernelHelper:: globalIndex const elemDofNumber, globalIndex const neighborDofNumber, globalIndex const upwViscDofNumber, - globalIndex const faceDofNumber, real64 const & dt, real64 ( & divMassFluxes )[ NC ], real64 ( & dDivMassFluxes_dElemVars )[ NC ][ (NC+1)*(NF+1) ], real64 ( & dDivMassFluxes_dFaceVars )[ NC ][ NF ], - globalIndex ( & dofColIndicesElemVars )[ (NC+1)*(NF+1) ], - globalIndex ( & dofColIndicesFaceVars )[ NF ] ) + globalIndex ( & dofColIndicesElemVars )[ (NC+1)*(NF+1) ] ) { integer constexpr NDOF = NC+1; localIndex const elemVarsOffset = NDOF*(ifaceLoc+1); @@ -884,7 +953,6 @@ AssemblerKernelHelper:: { dofColIndicesElemVars[elemVarsOffset+idof] = neighborDofNumber + idof; } - dofColIndicesFaceVars[ifaceLoc] = faceDofNumber; } template< integer NF, integer NC, integer NP > @@ -956,6 +1024,7 @@ void AssemblerKernelHelper:: assembleFaceConstraints( arrayView1d< globalIndex const > const & faceDofNumber, arrayView1d< integer const > const & faceGhostRank, + arrayView1d< integer const > const & isBoundaryFace, arraySlice1d< localIndex const > const & elemToFaces, globalIndex const elemDofNumber, globalIndex const rankOffset, @@ -983,7 +1052,16 @@ AssemblerKernelHelper:: // for each element, loop over the local (one-sided) faces for( integer ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) { - if( faceGhostRank[elemToFaces[ifaceLoc]] >= 0 ) + localIndex const kf = elemToFaces[ifaceLoc]; + + // Skip ghost faces + if( faceGhostRank[kf] >= 0 ) + { + continue; + } + + // Skip boundary faces - they have Dirichlet constraints, not flux continuity + if( isBoundaryFace[kf] > 0 ) { continue; } @@ -1039,9 +1117,10 @@ AssemblerKernel:: arrayView2d< localIndex const > const & elemList, arrayView1d< globalIndex const > const & faceDofNumber, arrayView1d< integer const > const & faceGhostRank, + arrayView1d< integer const > const & isBoundaryFace, arrayView1d< real64 const > const & facePres, arrayView1d< real64 const > const & faceGravCoef, - arrayView1d< real64 const > const & mimFaceGravCoef, + arrayView1d< real64 const > const & mimeticTransGgradZ, arraySlice1d< localIndex const > const & elemToFaces, real64 const & elemPres, real64 const & elemGravCoef, @@ -1060,7 +1139,6 @@ AssemblerKernel:: globalIndex const rankOffset, real64 const & dt, arraySlice2d< real64 const > const & transMatrix, - arraySlice2d< real64 const > const & transMatrixGrav, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) { @@ -1116,7 +1194,9 @@ AssemblerKernel:: elemList, regionFilter, faceDofNumber, - mimFaceGravCoef, + isBoundaryFace, + mimeticTransGgradZ, + faceGravCoef, elemToFaces, elemGravCoef, useTotalMassEquation, @@ -1130,7 +1210,7 @@ AssemblerKernel:: phaseCompFrac, dPhaseCompFrac, elemDofNumber, - transMatrixGrav, + transMatrix, oneSidedVolFlux, dOneSidedVolFlux_dPres, dOneSidedVolFlux_dFacePres, @@ -1144,6 +1224,7 @@ AssemblerKernel:: // enforcing flux continuity at this element's faces AssemblerKernelHelper::assembleFaceConstraints< NF, NC >( faceDofNumber, faceGhostRank, + isBoundaryFace, elemToFaces, elemDofNumber[er][esr][ei], rankOffset, @@ -1171,9 +1252,10 @@ FluxKernel:: ArrayOfArraysView< localIndex const > const & faceToNodes, arrayView1d< globalIndex const > const & faceDofNumber, arrayView1d< integer const > const & faceGhostRank, + arrayView1d< integer const > const & isBoundaryFace, arrayView1d< real64 const > const & facePres, arrayView1d< real64 const > const & faceGravCoef, - arrayView1d< real64 const > const & mimFaceGravCoef, + arrayView1d< real64 const > const & mimeticTransGgradZ, arrayView1d< real64 const > const & transMultiplier, ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseMob, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob, @@ -1223,8 +1305,6 @@ FluxKernel:: { // transmissibility matrix stackArray2d< real64, NF *NF > transMatrix( NF, NF ); - stackArray2d< real64, NF *NF > transMatrixGrav( NF, NF ); - real64 const perm[ 3 ] = { elemPerm[ei][0][0], elemPerm[ei][0][1], elemPerm[ei][0][2] }; // recompute the local transmissibility matrix at each iteration @@ -1239,19 +1319,6 @@ FluxKernel:: lengthTolerance, transMatrix ); - // currently the gravity term in the transport scheme is treated as in MRST, that is, always with TPFA - // this is why below we have to recompute the TPFA transmissibility in addition to the transmissibility matrix above - // TODO: treat the gravity term with a consistent inner product - mimeticInnerProduct::TPFAInnerProduct::compute< NF >( nodePosition, - transMultiplier, - faceToNodes, - elemToFaces[ei], - elemCenter[ei], - elemVolume[ei], - perm, - lengthTolerance, - transMatrixGrav ); - // perform flux assembly in this element compositionalMultiphaseHybridFVMKernels::AssemblerKernel::compute< NF, NC, NP >( er, esr, ei, regionFilter, @@ -1260,9 +1327,10 @@ FluxKernel:: elemList, faceDofNumber, faceGhostRank, + isBoundaryFace, facePres, faceGravCoef, - mimFaceGravCoef, + mimeticTransGgradZ, elemToFaces[ei], elemPres[ei], elemGravCoef[ei], @@ -1281,12 +1349,461 @@ FluxKernel:: rankOffset, dt, transMatrix, - transMatrixGrav, localMatrix, localRhs ); } ); } +/******************************** DirichletFluxKernel ********************************/ + +template< integer NF, integer NC, integer NP, typename IP_TYPE > +void +DirichletFluxKernel:: + launch( integer const numPhases, + localIndex const er, + localIndex const esr, + CellElementSubRegion const & subRegion, + constitutive::PermeabilityBase const & permeabilityModel, + SortedArrayView< localIndex const > const & boundaryFaceSet, + arrayView1d< real64 const > const & facePres, + arrayView1d< real64 const > const & faceTemp, + arrayView2d< real64 const, compflow::USD_PHASE > const & facePhaseMob, + arrayView2d< real64 const, compflow::USD_PHASE > const & facePhaseMassDens, + arrayView3d< real64 const, compflow::USD_PHASE_COMP > const & facePhaseCompFrac, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arrayView2d< localIndex const > const & elemRegionList, + arrayView2d< localIndex const > const & elemSubRegionList, + arrayView2d< localIndex const > const & elemList, + arrayView1d< globalIndex const > const & faceDofNumber, + arrayView1d< real64 const > const & faceGravCoef, + arrayView1d< real64 const > const & transMultiplier, + real64 const lengthTolerance, + real64 const dt, + globalIndex const rankOffset, + integer const useTotalMassEquation, + CompFlowAccessors const & compFlowAccessors, + MultiFluidAccessors const & multiFluidAccessors, + ElementViewConst< arrayView1d< globalIndex const > > const & elemDofNumber, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_UNUSED_VAR( faceTemp ); + using Deriv = multifluid::DerivativeOffset; + + // Get element data + arrayView2d< localIndex const > const & elemToFaces = subRegion.faceList(); + arrayView1d< real64 const > const & elemPres = subRegion.getField< fields::flow::pressure >(); + arrayView1d< real64 const > const & elemGravCoef = subRegion.getField< fields::flow::gravityCoefficient >(); + arrayView2d< real64 const > const & elemCenter = subRegion.getReference< array2d< real64 > >( CellElementSubRegion::viewKeyStruct::elementCenterString() ); + arrayView1d< real64 const > const & elemVolume = subRegion.getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString() ); + arrayView3d< real64 const > const & elemPerm = permeabilityModel.permeability(); + arrayView1d< integer const > const & elemGhostRank = subRegion.ghostRank(); + + // Get compositional flow fields + auto phaseMob = compFlowAccessors.get( fields::flow::phaseMobility{} ); + auto dPhaseMob = compFlowAccessors.get( fields::flow::dPhaseMobility{} ); + auto dCompFrac_dCompDens = compFlowAccessors.get( fields::flow::dGlobalCompFraction_dGlobalCompDensity{} ); + + // Get multifluid fields + auto phaseMassDens = multiFluidAccessors.get( fields::multifluid::phaseMassDensity{} ); + auto dPhaseMassDens = multiFluidAccessors.get( fields::multifluid::dPhaseMassDensity{} ); + auto phaseCompFrac = multiFluidAccessors.get( fields::multifluid::phaseCompFraction{} ); + auto dPhaseCompFrac = multiFluidAccessors.get( fields::multifluid::dPhaseCompFraction{} ); + + // Loop over boundary faces + forAll< parallelDevicePolicy<> >( boundaryFaceSet.size(), [=] GEOS_DEVICE ( localIndex const iset ) + { + localIndex const kf = boundaryFaceSet[iset]; + + // Find the element adjacent to this boundary face + localIndex erAdj = -1, esrAdj = -1, eiAdj = -1; + for( integer ke = 0; ke < elemRegionList.size( 1 ); ++ke ) + { + if( elemRegionList[kf][ke] == er && elemSubRegionList[kf][ke] == esr ) + { + erAdj = elemRegionList[kf][ke]; + esrAdj = elemSubRegionList[kf][ke]; + eiAdj = elemList[kf][ke]; + break; + } + } + + // Skip if no adjacent element in target region + if( eiAdj < 0 || elemGhostRank[eiAdj] >= 0 ) + { + return; + } + + // Compute one-sided transmissibility + stackArray2d< real64, NF * NF > transMatrix( NF, NF ); + real64 const perm[3] = { elemPerm[eiAdj][0][0], elemPerm[eiAdj][0][1], elemPerm[eiAdj][0][2] }; + + IP_TYPE::template compute< NF >( nodePosition, + transMultiplier, + faceToNodes, + elemToFaces[eiAdj], + elemCenter[eiAdj], + elemVolume[eiAdj], + perm, + lengthTolerance, + transMatrix ); + + // Find local face index + integer ifaceLoc = -1; + for( integer j = 0; j < NF; ++j ) + { + if( elemToFaces[eiAdj][j] == kf ) + { + ifaceLoc = j; + break; + } + } + if( ifaceLoc < 0 ) + { + return; + } + + // Get all global face indices of the cell touching the boundary + stackArray1d< localIndex, NF > cellFaces( NF ); + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + { + cellFaces[jfaceLoc] = elemToFaces[eiAdj][jfaceLoc]; + } + + // Component fluxes + real64 compFlux[NC]{}; + real64 dCompFlux_dP[NC]{}; + real64 dCompFlux_dC[NC][NC]{}; + real64 dCompFlux_dFaceP[NC][NF]{}; // Derivatives w.r.t. face pressures + + // Loop over phases to compute component fluxes + // Use face-based properties for boundary conditions (upwinding) + for( integer ip = 0; ip < numPhases; ++ip ) + { + real64 dDensMean_dC[NC]{}; + real64 dF_dC[NC]{}; + real64 dProp_dC[NC]{}; + + // Use average of element and face phase mass density for gravity term + real64 const elemDens = phaseMassDens[erAdj][esrAdj][eiAdj][0][ip]; + real64 const faceDens = facePhaseMassDens[kf][ip]; + real64 const densMean = 0.5 * (elemDens + faceDens); + + applyChainRule( NC, + dCompFrac_dCompDens[erAdj][esrAdj][eiAdj], + dPhaseMassDens[erAdj][esrAdj][eiAdj][0][ip], + dProp_dC, + Deriv::dC ); + + real64 const dDensMean_dP = 0.5 * dPhaseMassDens[erAdj][esrAdj][eiAdj][0][ip][Deriv::dP]; + for( integer jc = 0; jc < NC; ++jc ) + { + dDensMean_dC[jc] = 0.5 * dProp_dC[jc]; + } + + // Compute flux by looping over all connected faces with consistent transmissibility + real64 f = 0.0; + real64 dF_dP = 0.0; + real64 dF_dFaceP[NF]{}; // Derivatives of flux w.r.t. face pressures + for( integer jc = 0; jc < NC; ++jc ) + { + dF_dC[jc] = 0.0; + } + + // Sum contributions from all connected faces + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + { + // Each stencil face j uses its own face pressure + localIndex const kfj = cellFaces[jfaceLoc]; + real64 const gravTimesDz_j = elemGravCoef[eiAdj] - faceGravCoef[kfj]; + real64 const Tij = transMatrix[ifaceLoc][jfaceLoc]; + real64 const potDif_j = elemPres[eiAdj] - facePres[kfj] - densMean * gravTimesDz_j; + + // Accumulate flux and derivatives + f += Tij * potDif_j; + dF_dP += Tij * ( 1.0 - dDensMean_dP * gravTimesDz_j ); + + // Correct Jacobian: d(potDif_j)/d(facePres[kfj]) = -1, so dF/d(facePres_j) = -T_ij + dF_dFaceP[jfaceLoc] = -Tij; + + for( integer jc = 0; jc < NC; ++jc ) + { + dF_dC[jc] += -Tij * dDensMean_dC[jc] * gravTimesDz_j; + } + } + + // Use element mobility (simplified upwinding for Dirichlet BC) + // Upwind phase component fraction based on flow direction + // Use the total flux f for upwinding decision + real64 const beta = (f > 0.0) ? 1.0 : 0.0; + real64 const facePhaseMobility = facePhaseMob[kf][ip]; + real64 const phaseMobility = beta * phaseMob[erAdj][esrAdj][eiAdj][ip] + (1.0 - beta) * facePhaseMobility; + real64 const phaseFlux = phaseMobility * f; + real64 const dPhaseFlux_dP = phaseMobility * dF_dP + beta * dPhaseMob[erAdj][esrAdj][eiAdj][ip][Deriv::dP] * f; + real64 dPhaseFlux_dC[NC]; + real64 dPhaseFlux_dFaceP[NF]; + for( integer jc = 0; jc < NC; ++jc ) + { + dPhaseFlux_dC[jc] = phaseMobility * dF_dC[jc] + beta * dPhaseMob[erAdj][esrAdj][eiAdj][ip][Deriv::dC+jc] * f; + } + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + { + dPhaseFlux_dFaceP[jfaceLoc] = phaseMobility * dF_dFaceP[jfaceLoc]; + } + + // Use face-based phase composition for boundary conditions when flow is INTO the domain (potDif < 0) + // Use element phase composition when flow is OUT of the domain (potDif > 0) + for( integer ic = 0; ic < NC; ++ic ) + { + real64 const ycpElem = phaseCompFrac[erAdj][esrAdj][eiAdj][0][ip][ic]; + real64 const ycpFace = facePhaseCompFrac[kf][ip][ic]; + + // Upwind phase component fraction based on flow direction + real64 const ycp = beta * ycpElem + (1.0 - beta) * ycpFace; + + compFlux[ic] += ycp * phaseFlux; + + // For derivatives, use element properties (simplified for boundary conditions) + dCompFlux_dP[ic] += dPhaseFlux_dP * ycp + beta * phaseFlux * dPhaseCompFrac[erAdj][esrAdj][eiAdj][0][ip][ic][Deriv::dP]; + + // Compute dycpElem/dC for this component ic + real64 dycpElem_dC[NC]; + applyChainRule( NC, + dCompFrac_dCompDens[erAdj][esrAdj][eiAdj], + dPhaseCompFrac[erAdj][esrAdj][eiAdj][0][ip][ic], + dycpElem_dC, + Deriv::dC ); + for( integer jc = 0; jc < NC; ++jc ) + { + // d(ycp * phaseFlux)/dC = dycp/dC * phaseFlux + ycp * dPhaseFlux/dC + // dycp/dC = beta * dycpElem/dC (ycpFace is BC, independent of C) + dCompFlux_dC[ic][jc] += ycp * dPhaseFlux_dC[jc] + beta * phaseFlux * dycpElem_dC[jc]; + } + + // Add derivatives w.r.t. face pressures + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + { + dCompFlux_dFaceP[ic][jfaceLoc] += dPhaseFlux_dFaceP[jfaceLoc] * ycp; + } + } + } + + // Assemble into residual and Jacobian + real64 localFlux[NC]; + real64 localFluxJacobian[NC][NC+1]; + + for( integer ic = 0; ic < NC; ++ic ) + { + localFlux[ic] = dt * compFlux[ic]; + localFluxJacobian[ic][0] = dt * dCompFlux_dP[ic]; + for( integer jc = 0; jc < NC; ++jc ) + { + localFluxJacobian[ic][jc+1] = dt * dCompFlux_dC[ic][jc]; + } + } + + // Apply total mass equation transformation if needed + if( useTotalMassEquation ) + { + real64 work[NC+1]{}; + compositionalMultiphaseUtilities::shiftRowsAheadByOneAndReplaceFirstRowWithColumnSum( NC, NC+1, localFluxJacobian, work ); + compositionalMultiphaseUtilities::shiftElementsAheadByOneAndReplaceFirstElementWithSum( NC, localFlux ); + } + + // Add to global system + globalIndex const elemDof = elemDofNumber[erAdj][esrAdj][eiAdj]; + localIndex const localRow = LvArray::integerConversion< localIndex >( elemDof - rankOffset ); + + for( integer ic = 0; ic < NC; ++ic ) + { + RAJA::atomicAdd( parallelDeviceAtomic{}, &localRhs[localRow + ic], localFlux[ic] ); + + globalIndex dofColIndices[NC+1]; + dofColIndices[0] = elemDof; + for( integer jc = 0; jc < NC; ++jc ) + { + dofColIndices[jc+1] = elemDof + jc + 1; + } + + localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow + ic, + dofColIndices, + localFluxJacobian[ic], + NC+1 ); + + // Add contributions from face pressure derivatives + for( integer jfaceLoc = 0; jfaceLoc < NF; ++jfaceLoc ) + { + localIndex const kfj = cellFaces[jfaceLoc]; + globalIndex const faceDof = faceDofNumber[kfj]; + + real64 facePressureJacobian = dt * dCompFlux_dFaceP[ic][jfaceLoc]; + + // Apply total mass equation transformation if needed + if( useTotalMassEquation && ic == 0 ) + { + // For total mass equation, sum all component derivatives + facePressureJacobian = 0.0; + for( integer icc = 0; icc < NC; ++icc ) + { + facePressureJacobian += dt * dCompFlux_dFaceP[icc][jfaceLoc]; + } + } + + localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( localRow + ic, + &faceDof, + &facePressureJacobian, + 1 ); + } + + } + + } ); +} + +/******************************** EvaluateBCFacePropertiesKernel ********************************/ + +/** + * @brief Evaluate constitutive properties at BC face conditions using flash calculations + * @tparam NC number of components + * @tparam NP number of phases + * @param[in] numPhases number of phases in the simulation + * @param[in] boundaryFaceSet sorted array view of boundary face indices + * @param[in] facePres pressure values at boundary faces + * @param[in] faceTemp temperature values at boundary faces + * @param[in] faceCompFrac component fraction values at boundary faces + * @param[in] elemRegionList face-to-element region mapping + * @param[in] elemSubRegionList face-to-element subregion mapping + * @param[in] elemList face-to-element index mapping + * @param[in] er target element region index + * @param[in] esr target element subregion index + * @param[in] fluid reference to multifluid constitutive model + * @param[in] relperm reference to relative permeability constitutive model + * @param[out] facePhaseMob computed phase mobility at BC faces + * @param[out] facePhaseMassDens computed phase mass density at BC faces + * @param[out] facePhaseCompFrac computed phase component fractions at BC faces + * + * @note This function uses serialPolicy (host execution) because CUDA does not allow + * extended device lambdas inside generic lambda expressions. The output arrays + * must be explicitly moved to device memory after this function completes. + */ +template< integer NC, integer NP > +void +evaluateBCFaceProperties( integer const numPhases, + SortedArrayView< localIndex const > const & boundaryFaceSet, + arrayView1d< real64 const > const & facePres, + arrayView1d< real64 const > const & faceTemp, + arrayView2d< real64 const, compflow::USD_COMP > const & faceCompFrac, + arrayView2d< localIndex const > const & elemRegionList, + arrayView2d< localIndex const > const & elemSubRegionList, + arrayView2d< localIndex const > const & elemList, + localIndex const er, + localIndex const esr, + constitutive::MultiFluidBase & fluid, + constitutive::RelativePermeabilityBase & relperm, + arrayView2d< real64, compflow::USD_PHASE > const & facePhaseMob, + arrayView2d< real64, compflow::USD_PHASE > const & facePhaseMassDens, + arrayView3d< real64, compflow::USD_PHASE_COMP > const & facePhaseCompFrac ) +{ + // Use constitutiveUpdatePassThru to dispatch to concrete fluid and relperm types + constitutiveUpdatePassThru( fluid, [&] ( auto & castedFluid ) + { + using FluidType = TYPEOFREF( castedFluid ); + typename FluidType::KernelWrapper fluidWrapper = castedFluid.createKernelWrapper(); + + constitutive::constitutiveUpdatePassThru( relperm, [&] ( auto & castedRelperm ) + { + using RelPermType = TYPEOFREF( castedRelperm ); + typename RelPermType::KernelWrapper relPermWrapper = castedRelperm.createKernelWrapper(); + GEOS_UNUSED_VAR( relPermWrapper ); + + // Loop over BC faces and evaluate properties at BC conditions + // Note: Using serialPolicy (host) because extended device lambdas cannot be defined inside generic lambdas + // The output arrays will be moved to device after this completes + forAll< serialPolicy >( boundaryFaceSet.size(), [=, &facePhaseMob, &facePhaseMassDens, &facePhaseCompFrac] ( localIndex const iset ) + { + localIndex const kf = boundaryFaceSet[iset]; + + // Find adjacent element in target region + localIndex eiAdj = -1; + for( integer ke = 0; ke < elemRegionList.size( 1 ); ++ke ) + { + if( elemRegionList[kf][ke] == er && elemSubRegionList[kf][ke] == esr ) + { + eiAdj = elemList[kf][ke]; + break; + } + } + if( eiAdj < 0 ) + return; + + // Allocate temporary storage for face constitutive properties + StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseFrac( 1, 1, numPhases ); + StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseDens( 1, 1, numPhases ); + StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseMassDensLocal( 1, 1, numPhases ); + StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseVisc( 1, 1, numPhases ); + StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseEnthalpy( 1, 1, numPhases ); + StackArray< real64, 3, constitutive::MultiFluidBase::MAX_NUM_PHASES, constitutive::multifluid::LAYOUT_PHASE > facePhaseInternalEnergy( 1, 1, numPhases ); + StackArray< real64, 4, constitutive::MultiFluidBase::MAX_NUM_PHASES * NC, + constitutive::multifluid::LAYOUT_PHASE_COMP > facePhaseCompFracLocal( 1, 1, numPhases, NC ); + + // Initialize all temporary arrays to zero to ensure deterministic behavior on GPU + for( integer ip = 0; ip < numPhases; ++ip ) + { + facePhaseFrac[0][0][ip] = 0.0; + facePhaseDens[0][0][ip] = 0.0; + facePhaseMassDensLocal[0][0][ip] = 0.0; + facePhaseVisc[0][0][ip] = 0.0; + facePhaseEnthalpy[0][0][ip] = 0.0; + facePhaseInternalEnergy[0][0][ip] = 0.0; + for( integer ic = 0; ic < NC; ++ic ) + { + facePhaseCompFracLocal[0][0][ip][ic] = 0.0; + } + } + + real64 faceTotalDens = 0.0; + + // Evaluate fluid properties at BC face conditions using flash calculation + constitutive::MultiFluidBase::KernelWrapper::computeValues( fluidWrapper, + facePres[kf], + faceTemp[kf], + faceCompFrac[kf], + facePhaseFrac[0][0], + facePhaseDens[0][0], + facePhaseMassDensLocal[0][0], + facePhaseVisc[0][0], + facePhaseEnthalpy[0][0], + facePhaseInternalEnergy[0][0], + facePhaseCompFracLocal[0][0], + faceTotalDens ); + + // Evaluate relative permeability at face saturation from flash calculation +// relPermWrapper.compute( facePhaseFrac[0][0], 0, 0 ); + + // Store computed properties in output arrays + for( integer ip = 0; ip < numPhases; ++ip ) + { + // Store phase mass density from flash calculation + facePhaseMassDens[kf][ip] = facePhaseMassDensLocal[0][0][ip]; + + // Compute mobility from relative permeability evaluated at face conditions + real64 const faceKr = facePhaseFrac[0][0][ip]; // phaseRelPerm[eiAdj][0][ip]; + real64 const mu = facePhaseVisc[0][0][ip]; + // Safety check: ensure faceTotalDens and mu are valid before computing mobility + facePhaseMob[kf][ip] = (mu > 0 && faceTotalDens > 0) ? faceTotalDens * faceKr / mu : 0.0; + + // Store phase composition from flash calculation + for( integer ic = 0; ic < NC; ++ic ) + { + facePhaseCompFrac[kf][ip][ic] = facePhaseCompFracLocal[0][0][ip][ic]; + } + } + } ); + } ); // end nested constitutiveUpdatePassThru for relperm + } ); // end constitutiveUpdatePassThru for fluid +} + } // namespace compositionalMultiphaseHybridFVMKernels } // namespace geos diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernelsInstantiations.cpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernelsInstantiations.cpp new file mode 100644 index 00000000000..a573b1db6cb --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernelsInstantiations.cpp @@ -0,0 +1,97 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SinglePhaseHybridFVMKernelsInstantiations.cpp + * @brief Centralized explicit template instantiations for single-phase hybrid FVM kernels. + * + * Define GEOS_ENABLE_MANUAL_HYBRID_FVM_INST to activate the explicit instantiations below. + * When enabled, ensure overlapping template-generated sources are not compiled to avoid duplicates. + */ + +#include "physicsSolvers/fluidFlow/kernels/singlePhase/SinglePhaseHybridFVMKernels.hpp" + +#include "finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp" +#include "finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp" +#include "finiteVolume/mimeticInnerProducts/QuasiRTInnerProduct.hpp" +#include "finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp" +#include "finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp" + +namespace geos +{ +namespace singlePhaseHybridFVMKernels +{ + +#if defined(GEOS_ENABLE_MANUAL_HYBRID_FVM_INST) + +// ----------------------------------------------------------------------------- +// Faces-per-element list (kept in sync with kernelLaunchSelectorFaceSwitch) +// ----------------------------------------------------------------------------- +#define GEOS_NFACES_LIST \ + X_NF( 4 ) \ + X_NF( 5 ) \ + X_NF( 6 ) \ + X_NF( 7 ) \ + X_NF( 8 ) \ + X_NF( 9 ) \ + X_NF( 10 ) \ + X_NF( 11 ) \ + X_NF( 12 ) \ + X_NF( 13 ) + +// ----------------------------------------------------------------------------- +// Inner product list helper +// ----------------------------------------------------------------------------- +#define GEOS_FOR_EACH_IP( NF, MACRO ) \ + MACRO( NF, mimeticInnerProduct::TPFAInnerProduct ) \ + MACRO( NF, mimeticInnerProduct::QuasiTPFAInnerProduct ) \ + MACRO( NF, mimeticInnerProduct::QuasiRTInnerProduct ) \ + MACRO( NF, mimeticInnerProduct::SimpleInnerProduct ) \ + MACRO( NF, mimeticInnerProduct::BdVLMInnerProduct ) + +// ----------------------------------------------------------------------------- +// AveragePressureGradientKernel explicit class instantiations +// ----------------------------------------------------------------------------- +#define INSTANTIATE_AVG_PRES_GRAD( NF ) \ + template class AveragePressureGradientKernel< NF >; + +#define X_NF( NF ) INSTANTIATE_AVG_PRES_GRAD( NF ) +GEOS_NFACES_LIST +#undef X_NF +#undef INSTANTIATE_AVG_PRES_GRAD + +// ----------------------------------------------------------------------------- +// ElementBasedAssemblyKernel explicit class instantiations +// ----------------------------------------------------------------------------- +#define INSTANTIATE_ELEM_KERNEL( NF, IP ) \ + template class ElementBasedAssemblyKernel< NF, IP >; + +#define INSTANTIATE_ELEM_KERNEL_FOR_IPS( NF ) \ + GEOS_FOR_EACH_IP( NF, INSTANTIATE_ELEM_KERNEL ) + +#define X_NF( NF ) INSTANTIATE_ELEM_KERNEL_FOR_IPS( NF ) +GEOS_NFACES_LIST +#undef X_NF +#undef INSTANTIATE_ELEM_KERNEL_FOR_IPS +#undef INSTANTIATE_ELEM_KERNEL + +// Cleanup local helper macros +#undef GEOS_FOR_EACH_IP +#undef GEOS_NFACES_LIST + +#endif // GEOS_ENABLE_MANUAL_HYBRID_FVM_INST + +} // namespace singlePhaseHybridFVMKernels +} // namespace geos