From 52439ee968a48bc1ed7d04c8601f19eee95254b5 Mon Sep 17 00:00:00 2001 From: "Yoojin S. Cha" Date: Thu, 12 Mar 2026 13:20:52 -0700 Subject: [PATCH 01/11] implement mfd for mixed method and update unit tests --- .../BdVLMInnerProduct.hpp | 142 +++++++- .../QuasiTPFAInnerProduct.hpp | 137 +++++++- .../SimpleInnerProduct.hpp | 131 +++++++- .../mimeticInnerProducts/TPFAInnerProduct.hpp | 75 ++++- .../unitTests/testMimeticInnerProducts.cpp | 309 ++++++++++++++++++ 5 files changed, 790 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp index 5bde9c3a3c3..62ed1c5d1cf 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp @@ -40,7 +40,7 @@ class BdVLMInnerProduct : public MimeticInnerProductBase /** * @brief In a given element, recompute the transmissibility matrix in a cell using the inner product of Beirao da Veiga, Lipnikov, - * Manzini + * Manzini (page 113) * @param[in] nodePosition the position of the nodes * @param[in] transMultiplier the transmissibility multipliers at the mesh faces * @param[in] faceToNodes the map from the face to their nodes @@ -65,6 +65,18 @@ class BdVLMInnerProduct : public MimeticInnerProductBase real64 const (&elemPerm)[ 3 ], real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); + + template< localIndex NF > + GEOS_HOST_DEVICE + static void + computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ); }; @@ -207,6 +219,134 @@ BdVLMInnerProduct::compute( arrayView2d< real64 const, nodes::REFERENCE_POSITION } +template< localIndex NF > +GEOS_HOST_DEVICE +void +BdVLMInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ) +{ + GEOS_UNUSED_VAR( elemVolume ); + real64 const areaTolerance = lengthTolerance * lengthTolerance; + + real64 cellToFaceMat[ NF ][ 3 ] = {{ 0 }}; + real64 normalsMat[ NF ][ 3 ] = {{ 0 }}; + real64 permMat[ 3 ][ 3 ] = {{ 0 }}; + real64 faceAreaMat[ NF ][ NF ] = {{ 0 }}; + real64 faceArea[ NF ] = { 0.0 }; // store diag(A) + + real64 work_dimByDim[ 3 ][ 3 ] = {{ 0 }}; + real64 work_numFacesByDim[ NF ][ 3 ] = {{ 0 }}; + real64 work_dimByNumFaces[ 3 ][ NF ] = {{ 0 }}; + real64 work_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; + real64 tmp_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; + + // 0) assemble full coefficient tensor from principal axis/components + MimeticInnerProductHelpers::makeFullTensor( elemPerm, permMat ); + + // 1) fill R (= cellToFaceMat) and normalsMat + for ( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + { + real64 faceCenter[ 3 ], faceNormal[ 3 ], cellToFaceVec[ 3 ]; + + faceAreaMat[ ifaceLoc ][ ifaceLoc ] = + computationalGeometry::centroid_3DPolygon( faceToNodes[ elemToFaces[ ifaceLoc ] ], + nodePosition, + faceCenter, + faceNormal, + areaTolerance ); + + faceArea[ ifaceLoc ] = faceAreaMat[ ifaceLoc ][ ifaceLoc ]; + + LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); + LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); + + // R row: A_f * (x_f - x_c) + cellToFaceMat[ ifaceLoc ][ 0 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 0 ]; + cellToFaceMat[ ifaceLoc ][ 1 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 1 ]; + cellToFaceMat[ ifaceLoc ][ 2 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 2 ]; + + // orient normal outward + if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) + { + LvArray::tensorOps::scale< 3 >( faceNormal, -1.0 ); + } + + // normalsMat row + normalsMat[ ifaceLoc ][ 0 ] = faceNormal[ 0 ]; + normalsMat[ ifaceLoc ][ 1 ] = faceNormal[ 1 ]; + normalsMat[ ifaceLoc ][ 2 ] = faceNormal[ 2 ]; + } + + // 2) compute N + LvArray::tensorOps::Rij_eq_AikBkj< NF, 3, 3 >( work_numFacesByDim, + normalsMat, + permMat ); + + // BdVLM inner product matrix M = M0 + M1 by Beirao da Veiga, Lipnikov, Manzini (pg.89) + // M0 = R (R^T N)^(-1) R^T + + // 3) compute (R^T N)^-1 -> (3 X 3), work_dimByDim = R^T N + LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, NF >( work_dimByDim, + cellToFaceMat, + work_numFacesByDim ); + LvArray::tensorOps::invert< 3 >( work_dimByDim ); + + // 4) compute R (R^T N)^(-1) R^T into M + // work_dimByNumFaces = (R^T N)^-1 R^T -> (3 X NF) + LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work_dimByNumFaces, + work_dimByDim, + cellToFaceMat ); + + // M = R * work_dimByNumFaces (NF x NF) + LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( M, + cellToFaceMat, + work_dimByNumFaces ); + + // P_N = I - N (N^T N)^(-1) N^T + // 5) compute (N^T N)^-1 + LvArray::tensorOps::Rij_eq_AkiAkj< 3, NF >( work_dimByDim, + work_numFacesByDim ); // N + LvArray::tensorOps::invert< 3 >( work_dimByDim ); + + // 6) tmp = (N^T N)^-1 N^T + LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work_dimByNumFaces, + work_dimByDim, + work_numFacesByDim ); + + // 7) work_numFacesByNumFaces = -I + N * tmp + LvArray::tensorOps::addIdentity< NF >( work_numFacesByNumFaces, -1.0 ); + LvArray::tensorOps::Rij_add_AikBkj< NF, NF, 3 >( work_numFacesByNumFaces, + work_numFacesByDim, // N + work_dimByNumFaces // tmp + ); + + // 8) M = M0 + gamma * P_N + real64 const gamma = 1.0 / static_cast< real64 >( NF ); + LvArray::tensorOps::scaledAdd< NF, NF >( M, work_numFacesByNumFaces, -gamma ); + + // convert to flux-based inner product: M_flux = A^{-1} M_vel A^{-1} + // here, A is diag(faceArea) + real64 invA[ NF ]; + for( localIndex i = 0; i < NF; ++i ) + { + invA[ i ] = 1.0 / faceArea[ i ]; + } + + for( localIndex i = 0; i < NF; ++i ) + { + for( localIndex j = 0; j < NF; ++j ) + { + M[ i ][ j ] *= invA[ i ] * invA[ j ]; + } + } +} + } // end namespace mimeticInnerProduct } // end namespace geos diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp index 4f481e7181b..4082941a677 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp @@ -62,7 +62,18 @@ class QuasiTPFAInnerProduct : public MimeticInnerProductBase real64 const (&elemPerm)[ 3 ], real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); - + + template< localIndex NF > + GEOS_HOST_DEVICE + static void + computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ); }; template< localIndex NF > @@ -90,6 +101,130 @@ QuasiTPFAInnerProduct::compute( arrayView2d< real64 const, nodes::REFERENCE_POSI transMatrix ); } +template< localIndex NF > +GEOS_HOST_DEVICE +void +QuasiTPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ) +{ + real64 const areaTolerance = lengthTolerance * lengthTolerance; + + // 0) stabilization parameter for quasi-TPFA + constexpr real64 tParam = 2.0; + + // 1) assemble permeability tensor and its inverse + real64 K[3][3] = {{0}}; + MimeticInnerProductHelpers::makeFullTensor( elemPerm, K ); + + real64 Kinv[3][3] = {{0}}; + for( int i = 0; i < 3; ++i ) + { + Kinv[i][i] = 1.0 / elemPerm[i]; + } + + // 2) compute C and N + real64 C[ NF ][ 3 ] = {{ 0 }}; + real64 N[ NF ][ 3 ] = {{ 0 }}; + + for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + { + real64 faceCenter[3], faceNormal[3]; + + real64 const faceArea = + computationalGeometry::centroid_3DPolygon( + faceToNodes[elemToFaces[ifaceLoc]], + nodePosition, + faceCenter, + faceNormal, + areaTolerance ); + + real64 cellToFaceVec[3]; + LvArray::tensorOps::copy<3>( cellToFaceVec, faceCenter ); + LvArray::tensorOps::subtract<3>( cellToFaceVec, elemCenter ); + + if( LvArray::tensorOps::AiBi<3>( cellToFaceVec, faceNormal ) < 0.0 ) + { + LvArray::tensorOps::scale<3>( faceNormal, -1.0 ); + } + + for( int d = 0; d < 3; ++d ) + { + C[ifaceLoc][d] = cellToFaceVec[d]; + N[ifaceLoc][d] = faceArea * faceNormal[d]; + } + } + + // 3) compute consistency part C K^{-1} C^T / volume + real64 CKCt[ NF ][ NF ] = {{ 0 }}; + real64 work[ 3 ][ NF ] = {{ 0 }}; + + LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work, Kinv, C ); + LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( CKCt, C, work ); + LvArray::tensorOps::scale< NF, NF >( CKCt, 1.0 / elemVolume ); + + // 4) compute W = N K N' + real64 W[ NF ][ NF ] = {{ 0 }}; + real64 workNK[ 3 ][ NF ] = {{ 0 }}; + + LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( workNK, K, N ); + LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( W, N, workNK ); + + // 5) build Q from C (orthonormal basis of consistency space) + real64 q0[ NF ], q1[ NF ], q2[ NF ]; + real64 Qmat[ NF ][ 3 ]; + + for( localIndex i = 0; i < NF; ++i ) + { + q0[i] = C[i][0]; + q1[i] = C[i][1]; + q2[i] = C[i][2]; + } + + MimeticInnerProductHelpers::orthonormalize< NF >( q0, q1, q2, Qmat ); + + // 6) compute P = I - Q Q^T + real64 P[ NF ][ NF ] = {{ 0 }}; + LvArray::tensorOps::addIdentity< NF >( P, -1.0 ); + LvArray::tensorOps::Rij_add_AikAjk< NF, 3 >( P, Qmat ); + LvArray::tensorOps::scale< NF, NF >( P, -1.0 ); + + // 7) compute stabilization term (v/t) P diag(W)^{-1} P + real64 stab[ NF ][ NF ] = {{ 0 }}; + + for( localIndex i = 0; i < NF; ++i ) + { + for( localIndex j = 0; j < NF; ++j ) + { + real64 val = 0.0; + for( localIndex k = 0; k < NF; ++k ) + { + if( LvArray::math::abs( W[k][k] ) > 0 ) + { + val += P[i][k] * ( 1.0 / W[k][k] ) * P[k][j]; + } + } + stab[i][j] = val; + } + } + + real64 const scale = elemVolume / tParam; + + // 8) assemble final M + for( localIndex i = 0; i < NF; ++i ) + { + for( localIndex j = 0; j < NF; ++j ) + { + M[i][j] = CKCt[i][j] + scale * stab[i][j]; + } + } +} + } // end namespace mimeticInnerProduct } // end namespace geos diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp index da7f4b2b330..5f243b8a589 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp @@ -23,6 +23,9 @@ #include "finiteVolume/mimeticInnerProducts/MimeticInnerProductBase.hpp" #include "finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp" #include "mesh/utilities/ComputationalGeometry.hpp" +#include +#include + namespace geos { @@ -64,7 +67,18 @@ class SimpleInnerProduct : public MimeticInnerProductBase real64 const (&elemPerm)[ 3 ], real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); - + + template< localIndex NF > + GEOS_HOST_DEVICE + static void + computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ); }; template< localIndex NF > @@ -195,6 +209,121 @@ SimpleInnerProduct::compute( arrayView2d< real64 const, nodes::REFERENCE_POSITIO } +template< localIndex NF > +GEOS_HOST_DEVICE +void +SimpleInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ) +{ + real64 const areaTolerance = lengthTolerance * lengthTolerance; + // 0) Helper: diagonal inverse of [k1, k2, k3] + real64 Kinv[3][3] = {{ 0 }}; + for( int i =0; i < 3; ++i) + { + Kinv[i][i] = 1.0 / elemPerm[i]; + } + + // 1) Compute C, N, A + real64 C[ NF ][ 3 ] = {{ 0 }}; + real64 N[ NF ][ 3 ] = {{ 0 }}; + real64 A[ NF ] = { 0.0 }; + + for( localIndex i = 0; i < NF; ++i) + { + real64 fCenter[3], fNormal[3]; + real64 area = computationalGeometry::centroid_3DPolygon( + faceToNodes[elemToFaces[i]], nodePosition, fCenter, fNormal, areaTolerance ); + + real64 vec_cf[3]; + LvArray::tensorOps::copy<3>(vec_cf, fCenter); + LvArray::tensorOps::subtract<3>(vec_cf, elemCenter); + if ( LvArray::tensorOps::AiBi<3>(vec_cf, fNormal) < 0.0 ) + { + LvArray::tensorOps::scale<3>(fNormal, -1.0); + } + + for ( int d = 0; d < 3; ++d) + { + C[i][d] = vec_cf[d]; + N[i][d] = fNormal[d] * area; + } + + A[i] = area; + } + + // 2) Compute C K^{-1} C^T / volume + real64 CKCt[ NF ][ NF ] = {{0}}; + real64 work3xNF[ 3 ][ NF ] = {{0}}; + + LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work3xNF, Kinv, C ); + LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( CKCt, C, work3xNF ); + LvArray::tensorOps::scale< NF, NF >( CKCt, 1.0 / elemVolume ); + + // 3) Q = orth(N / A) + real64 q0[ NF ], q1[ NF ], q2[ NF ]; + real64 Qmat[ NF ][ 3 ]; + for( localIndex i = 0; i < NF; ++i) + { + q0[i] = N[i][0] / A[i]; + q1[i] = N[i][1] / A[i]; + q2[i] = N[i][2] / A[i]; + } + + MimeticInnerProductHelpers::orthonormalize< NF >( q0, q1, q2, Qmat ); + + // 4) U = I - Q Q^T + real64 U[ NF ][ NF ] = {{0}}; + LvArray::tensorOps::addIdentity< NF >( U, -1.0 ); + LvArray::tensorOps::Rij_add_AikAjk< NF, 3 >( U, Qmat ); + LvArray::tensorOps::scale< NF, NF >( U, -1.0 ); + + // 5) A^{-1} * U * A^{-1} + real64 invA[ NF ]; + for( localIndex i = 0; i < NF; ++i) + { + invA[i] = 1.0 / A[i]; + } + + // temp = A^{-1} * U + real64 temp[ NF ][ NF ]; + for ( localIndex i = 0; i < NF; ++i ) + { + for ( localIndex j = 0; j < NF; ++j ) + { + temp[i][j] = invA[i] * U[i][j]; + } + } + + // Usc = A^{-1} * U * A^{-1} + real64 Usc[ NF ][ NF ]; + for ( localIndex i = 0; i < NF; ++i ) + { + for ( localIndex j = 0; j < NF; ++j ) + { + Usc[i][j] = temp[i][j] * invA[j]; + } + } + + // 6) M = CKCt + (v / t) * Usc + real64 tParam = 2.0 * ( elemPerm[0] + elemPerm[1] + elemPerm[2] ); + real64 scale = elemVolume / tParam; + + for( localIndex i = 0; i < NF; ++i ) + { + for ( localIndex j = 0; j < NF; ++j) + { + M[i][j] = CKCt[i][j] + scale * Usc[i][j]; + } + } +} + + } // end namespace mimeticInnerProduct } // end namespace geos diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp index dde1042e41b..e9bb2a5a159 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp @@ -60,7 +60,18 @@ class TPFAInnerProduct : public MimeticInnerProductBase real64 const (&elemPerm)[ 3 ], real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); - + + template< localIndex NF> + GEOS_HOST_DEVICE + static void + computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ); }; template< localIndex NF > @@ -129,6 +140,68 @@ TPFAInnerProduct::compute( arrayView2d< real64 const, nodes::REFERENCE_POSITION_ } } +template< localIndex NF > +GEOS_HOST_DEVICE +void +TPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ) +{ + GEOS_UNUSED_VAR( elemVolume ); + + real64 const areaTolerance = lengthTolerance * lengthTolerance; + real64 const weightTolerance = 1e-30 * lengthTolerance; + + // initialize M to zero + for( localIndex i = 0; i < NF; ++i ) + { + for ( localIndex j = 0; j < NF; ++j) + { + M[i][j] = 0.0; + } + } + + for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + { + real64 faceCenter[3], faceNormal[3], faceConormal[3], cellToFaceVec[3]; + + // 1) face geometry + real64 const faceArea = + computationalGeometry::centroid_3DPolygon( faceToNodes[elemToFaces[ifaceLoc]], + nodePosition, + faceCenter, + faceNormal, + areaTolerance ); + + LvArray::tensorOps::copy<3>(cellToFaceVec, faceCenter); + LvArray::tensorOps::subtract<3>(cellToFaceVec, elemCenter); + + if( LvArray::tensorOps::AiBi<3>(cellToFaceVec, faceNormal) < 0.0 ) + { + LvArray::tensorOps::scale< 3 >( faceNormal, -1 ); + } + + real64 const c2fDistance = LvArray::tensorOps::normalize<3>( cellToFaceVec ); + + // 2) K * n (TPFA assumes diagonal K) + LvArray::tensorOps::hadamardProduct<3>(faceConormal, elemPerm, faceNormal); + + // 3) compute T_ii + real64 Tii = LvArray::tensorOps::AiBi<3>( cellToFaceVec, faceConormal ) * faceArea / c2fDistance; + + Tii = LvArray::math::max( Tii, weightTolerance ); + + // 4) M = |T|^{-1} + M[ifaceLoc][ifaceLoc] = 1.0 / LvArray::math::abs( Tii ); + } +} + + } // end namespace mimeticInnerProduct } // end namespace geos diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp index 49c1df7a4c5..da4ae05498f 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp @@ -24,6 +24,8 @@ #include "mainInterface/initialization.hpp" #include "mesh/FaceManager.hpp" #include "mesh/utilities/ComputationalGeometry.hpp" +#include +#include // TPL includes #include @@ -60,6 +62,20 @@ void compareTransmissibilityMatrices( arraySlice2d< real64 const > const & trans } } +void compareInnerProductMatrices( arraySlice2d< real64 const > const & M, + arraySlice2d< real64 const > const & Mref ) +{ + for( localIndex ifaceLoc = 0; ifaceLoc < M.size( 0 ); ++ifaceLoc ) + { + for( localIndex jfaceLoc = 0; jfaceLoc < M.size( 1 ); ++jfaceLoc ) + { + checkRelativeError( M( ifaceLoc, jfaceLoc ), + Mref( ifaceLoc, jfaceLoc ), + 1e-15 ); + } + } +} + void computeVolumeAndCenter( array2d< real64, nodes::REFERENCE_POSITION_PERM > const & nodePosition, array1d< localIndex > const & toNodes, real64 ( & elemCenter )[3], @@ -527,6 +543,110 @@ void makeTetra( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodePosition } } +void makeCube( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodePosition, + FaceManager::NodeMapType & faceToNodes, + array1d< localIndex > & elemToFaces, + real64 ( & elemCenter )[3], + real64 & elemVolume, + real64 ( & elemPerm )[3], + real64 & lengthTolerance, + integer const ipType, + arraySlice2d< real64 > const & transMatrixRef ) +{ + localIndex constexpr numNodes = 8; + localIndex constexpr numFaces = 6; + localIndex constexpr numNodesPerFace = 4; + + lengthTolerance = 1e-12; + + // permeability + elemPerm[0] = 1.0; + elemPerm[1] = 1.0; + elemPerm[2] = 1.0; + + // elem-to-faces + elemToFaces.resize( numFaces ); + for( localIndex i = 0; i < numFaces; ++i ) + elemToFaces(i) = i; + + // face-to-nodes + faceToNodes.resize( numFaces ); + for( localIndex f = 0; f < numFaces; ++f ) + faceToNodes.resizeArray( f, numNodesPerFace ); + + // vertices + nodePosition.resize( numNodes, 3 ); + + nodePosition(0,0) = 0; nodePosition(0,1) = 0; nodePosition(0,2) = 0; + nodePosition(1,0) = 1; nodePosition(1,1) = 0; nodePosition(1,2) = 0; + nodePosition(2,0) = 1; nodePosition(2,1) = 1; nodePosition(2,2) = 0; + nodePosition(3,0) = 0; nodePosition(3,1) = 1; nodePosition(3,2) = 0; + nodePosition(4,0) = 0; nodePosition(4,1) = 0; nodePosition(4,2) = 1; + nodePosition(5,0) = 1; nodePosition(5,1) = 0; nodePosition(5,2) = 1; + nodePosition(6,0) = 1; nodePosition(6,1) = 1; nodePosition(6,2) = 1; + nodePosition(7,0) = 0; nodePosition(7,1) = 1; nodePosition(7,2) = 1; + + // x = 0 + faceToNodes(0,0) = 0; faceToNodes(0,1) = 3; faceToNodes(0,2) = 7; faceToNodes(0,3) = 4; + + // x = 1 + faceToNodes(1,0) = 1; faceToNodes(1,1) = 2; faceToNodes(1,2) = 6; faceToNodes(1,3) = 5; + + // y = 0 + faceToNodes(2,0) = 0; faceToNodes(2,1) = 1; faceToNodes(2,2) = 5; faceToNodes(2,3) = 4; + + // y = 1 + faceToNodes(3,0) = 3; faceToNodes(3,1) = 2; faceToNodes(3,2) = 6; faceToNodes(3,3) = 7; + + // z = 0 + faceToNodes(4,0) = 0; faceToNodes(4,1) = 1; faceToNodes(4,2) = 2; faceToNodes(4,3) = 3; + + // z = 1 + faceToNodes(5,0) = 4; faceToNodes(5,1) = 5; faceToNodes(5,2) = 6; faceToNodes(5,3) = 7; + + // center & volume + array1d< localIndex > toNodes; + toNodes.resize( numNodes ); + for( localIndex i = 0; i < numNodes; ++i ) + toNodes(i) = i; + + computeVolumeAndCenter( nodePosition, + toNodes, + elemCenter, + elemVolume ); + + // reference matrix (MRST / MATLAB) + if( ipType == InnerProductType::TPFA || + ipType == InnerProductType::QUASI_TPFA ) + { + // diagonal + transMatrixRef(0,0) = 0.5; + transMatrixRef(1,1) = 0.5; + transMatrixRef(2,2) = 0.5; + transMatrixRef(3,3) = 0.5; + transMatrixRef(4,4) = 0.5; + transMatrixRef(5,5) = 0.5; + } + else if( ipType == InnerProductType::SIMPLE || + ipType == InnerProductType::BDVLM ) + { + // SIMPLE = TPFA + stabilization + transMatrixRef(0,0) = 0.333333333333; + transMatrixRef(1,1) = 0.333333333333; + transMatrixRef(2,2) = 0.333333333333; + transMatrixRef(3,3) = 0.333333333333; + transMatrixRef(4,4) = 0.333333333333; + transMatrixRef(5,5) = 0.333333333333; + + transMatrixRef(0,1) = -0.166666666666; + transMatrixRef(1,0) = -0.166666666666; + transMatrixRef(2,3) = -0.166666666666; + transMatrixRef(3,2) = -0.166666666666; + transMatrixRef(4,5) = -0.166666666666; + transMatrixRef(5,4) = -0.166666666666; + } +} + template< localIndex NF, typename ARRAY_VIEW_T > static void runConsistencyTest( array2d< real64, nodes::REFERENCE_POSITION_PERM > const & nodePosition, FaceManager::NodeMapType const & faceToNodes, @@ -1003,6 +1123,195 @@ TEST( testMimeticInnerProducts, BdVLMtetra ) transMatrix.toViewConst() ); } +TEST( testMimeticInnerProducts, TPFAM_cube ) +{ + localIndex constexpr NF = 6; + + array2d< real64, nodes::REFERENCE_POSITION_PERM > nodePosition; + FaceManager::NodeMapType faceToNodes; + array1d< localIndex > elemToFaces; + real64 elemCenter[3] = { 0.0 }; + real64 elemPerm[3] = { 0.0 }; + real64 elemVolume = 0; + real64 lengthTolerance = 0; + + stackArray2d< real64, NF * NF > Mref( NF, NF ); + Mref.setValues< parallelHostPolicy >( 0.0 ); + + makeCube( nodePosition, + faceToNodes, + elemToFaces, + elemCenter, + elemVolume, + elemPerm, + lengthTolerance, + InnerProductType::TPFA, + Mref.toSlice() ); + + stackArray2d< real64, NF * NF > M( NF, NF ); + + stackArray1d< real64, 3 > center( 3 ); + center[0] = elemCenter[0]; + center[1] = elemCenter[1]; + center[2] = elemCenter[2]; + + real64 const perm[ 3 ] = { elemPerm[0], elemPerm[1], elemPerm[2] }; + + TPFAInnerProduct::computeM< NF >( nodePosition.toViewConst(), + faceToNodes.toViewConst(), + elemToFaces.toSliceConst(), + center, + 1, + perm, + lengthTolerance, + M.toSlice() ); + + compareInnerProductMatrices( M.toSliceConst(), + Mref.toSliceConst() ); +} + +TEST( testMimeticInnerProducts, QuasiTPFAM_cube ) +{ + localIndex constexpr NF = 6; + + array2d< real64, nodes::REFERENCE_POSITION_PERM > nodePosition; + FaceManager::NodeMapType faceToNodes; + array1d< localIndex > elemToFaces; + real64 elemCenter[3] = { 0.0 }; + real64 elemPerm[3] = { 0.0 }; + real64 elemVolume = 0; + real64 lengthTolerance = 0; + + stackArray2d< real64, NF * NF > Mref( NF, NF ); + Mref.setValues< parallelHostPolicy >( 0.0 ); + + makeCube( nodePosition, + faceToNodes, + elemToFaces, + elemCenter, + elemVolume, + elemPerm, + lengthTolerance, + InnerProductType::QUASI_TPFA, + Mref.toSlice() ); + + stackArray2d< real64, NF * NF > M( NF, NF ); + + stackArray1d< real64, 3 > center( 3 ); + center[0] = elemCenter[0]; + center[1] = elemCenter[1]; + center[2] = elemCenter[2]; + + real64 const perm[ 3 ] = { elemPerm[0], elemPerm[1], elemPerm[2] }; + + QuasiTPFAInnerProduct::computeM< NF >( nodePosition.toViewConst(), + faceToNodes.toViewConst(), + elemToFaces.toSliceConst(), + center, + 1, + perm, + lengthTolerance, + M.toSlice() ); + + compareInnerProductMatrices( M.toSliceConst(), + Mref.toSliceConst() ); +} + +TEST( testMimeticInnerProducts, SimpleM_cube ) +{ + localIndex constexpr NF = 6; + + array2d< real64, nodes::REFERENCE_POSITION_PERM > nodePosition; + FaceManager::NodeMapType faceToNodes; + array1d< localIndex > elemToFaces; + real64 elemCenter[3] = { 0.0 }; + real64 elemPerm[3] = { 0.0 }; + real64 elemVolume = 0; + real64 lengthTolerance = 0; + + stackArray2d< real64, NF * NF > Mref( NF, NF ); + Mref.setValues< parallelHostPolicy >( 0.0 ); + + makeCube( nodePosition, + faceToNodes, + elemToFaces, + elemCenter, + elemVolume, + elemPerm, + lengthTolerance, + InnerProductType::SIMPLE, + Mref.toSlice() ); + + stackArray2d< real64, NF * NF > M( NF, NF ); + + stackArray1d< real64, 3 > center( 3 ); + center[0] = elemCenter[0]; + center[1] = elemCenter[1]; + center[2] = elemCenter[2]; + + real64 const perm[ 3 ] = { elemPerm[0], elemPerm[1], elemPerm[2] }; + + SimpleInnerProduct::computeM< NF >( nodePosition.toViewConst(), + faceToNodes.toViewConst(), + elemToFaces.toSliceConst(), + center, + 1, + perm, + lengthTolerance, + M.toSlice() ); + + compareInnerProductMatrices( M.toSliceConst(), + Mref.toSliceConst() ); +} + +TEST( testMimeticInnerProducts, BdVLMM_cube ) +{ + localIndex constexpr NF = 6; + + array2d< real64, nodes::REFERENCE_POSITION_PERM > nodePosition; + FaceManager::NodeMapType faceToNodes; + array1d< localIndex > elemToFaces; + real64 elemCenter[3] = { 0.0 }; + real64 elemPerm[3] = { 0.0 }; + real64 elemVolume = 0; + real64 lengthTolerance = 0; + + stackArray2d< real64, NF * NF > M( NF, NF ); + M.setValues< parallelHostPolicy >( 0.0 ); + + stackArray2d< real64, NF * NF > Mref( NF, NF ); + Mref.setValues< parallelHostPolicy >( 0.0 ); + + makeCube( nodePosition, + faceToNodes, + elemToFaces, + elemCenter, + elemVolume, + elemPerm, + lengthTolerance, + InnerProductType::BDVLM, + Mref.toSlice() ); + + stackArray1d< real64, 3 > center( 3 ); + center[0] = elemCenter[0]; + center[1] = elemCenter[1]; + center[2] = elemCenter[2]; + + real64 const perm[ 3 ] = { elemPerm[0], elemPerm[1], elemPerm[2] }; + + BdVLMInnerProduct::computeM< NF >( nodePosition.toViewConst(), + faceToNodes.toViewConst(), + elemToFaces.toSliceConst(), + center, + elemVolume, + perm, + lengthTolerance, + M.toSlice() ); + + compareInnerProductMatrices( M.toSliceConst(), + Mref.toSliceConst() ); +} + //======================== Linear Pressure Recovery Test ============================= // Three cases: without distortion (unit cube cell), with distortion: (1) planar (2) nonplanar From f7d1850881c8032a3855b8645f9d19383694d0b1 Mon Sep 17 00:00:00 2001 From: "Yoojin S. Cha" Date: Thu, 12 Mar 2026 14:26:45 -0700 Subject: [PATCH 02/11] apply uncrustify --- .../BdVLMInnerProduct.hpp | 218 +++++++++--------- .../QuasiTPFAInnerProduct.hpp | 24 +- .../SimpleInnerProduct.hpp | 182 +++++++-------- .../mimeticInnerProducts/TPFAInnerProduct.hpp | 108 ++++----- .../unitTests/testMimeticInnerProducts.cpp | 102 ++++---- 5 files changed, 317 insertions(+), 317 deletions(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp index 62ed1c5d1cf..aa96cae10de 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp @@ -65,7 +65,7 @@ class BdVLMInnerProduct : public MimeticInnerProductBase real64 const (&elemPerm)[ 3 ], real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); - + template< localIndex NF > GEOS_HOST_DEVICE static void @@ -231,120 +231,120 @@ BdVLMInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITIO real64 const & lengthTolerance, arraySlice2d< real64 > const & M ) { - GEOS_UNUSED_VAR( elemVolume ); - real64 const areaTolerance = lengthTolerance * lengthTolerance; - - real64 cellToFaceMat[ NF ][ 3 ] = {{ 0 }}; - real64 normalsMat[ NF ][ 3 ] = {{ 0 }}; - real64 permMat[ 3 ][ 3 ] = {{ 0 }}; - real64 faceAreaMat[ NF ][ NF ] = {{ 0 }}; - real64 faceArea[ NF ] = { 0.0 }; // store diag(A) - - real64 work_dimByDim[ 3 ][ 3 ] = {{ 0 }}; - real64 work_numFacesByDim[ NF ][ 3 ] = {{ 0 }}; - real64 work_dimByNumFaces[ 3 ][ NF ] = {{ 0 }}; - real64 work_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; - real64 tmp_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; - - // 0) assemble full coefficient tensor from principal axis/components - MimeticInnerProductHelpers::makeFullTensor( elemPerm, permMat ); - - // 1) fill R (= cellToFaceMat) and normalsMat - for ( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + GEOS_UNUSED_VAR( elemVolume ); + real64 const areaTolerance = lengthTolerance * lengthTolerance; + + real64 cellToFaceMat[ NF ][ 3 ] = {{ 0 }}; + real64 normalsMat[ NF ][ 3 ] = {{ 0 }}; + real64 permMat[ 3 ][ 3 ] = {{ 0 }}; + real64 faceAreaMat[ NF ][ NF ] = {{ 0 }}; + real64 faceArea[ NF ] = { 0.0 }; // store diag(A) + + real64 work_dimByDim[ 3 ][ 3 ] = {{ 0 }}; + real64 work_numFacesByDim[ NF ][ 3 ] = {{ 0 }}; + real64 work_dimByNumFaces[ 3 ][ NF ] = {{ 0 }}; + real64 work_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; + real64 tmp_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; + + // 0) assemble full coefficient tensor from principal axis/components + MimeticInnerProductHelpers::makeFullTensor( elemPerm, permMat ); + + // 1) fill R (= cellToFaceMat) and normalsMat + for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + { + real64 faceCenter[ 3 ], faceNormal[ 3 ], cellToFaceVec[ 3 ]; + + faceAreaMat[ ifaceLoc ][ ifaceLoc ] = + computationalGeometry::centroid_3DPolygon( faceToNodes[ elemToFaces[ ifaceLoc ] ], + nodePosition, + faceCenter, + faceNormal, + areaTolerance ); + + faceArea[ ifaceLoc ] = faceAreaMat[ ifaceLoc ][ ifaceLoc ]; + + LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); + LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); + + // R row: A_f * (x_f - x_c) + cellToFaceMat[ ifaceLoc ][ 0 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 0 ]; + cellToFaceMat[ ifaceLoc ][ 1 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 1 ]; + cellToFaceMat[ ifaceLoc ][ 2 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 2 ]; + + // orient normal outward + if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) { - real64 faceCenter[ 3 ], faceNormal[ 3 ], cellToFaceVec[ 3 ]; - - faceAreaMat[ ifaceLoc ][ ifaceLoc ] = - computationalGeometry::centroid_3DPolygon( faceToNodes[ elemToFaces[ ifaceLoc ] ], - nodePosition, - faceCenter, - faceNormal, - areaTolerance ); - - faceArea[ ifaceLoc ] = faceAreaMat[ ifaceLoc ][ ifaceLoc ]; - - LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); - LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); - - // R row: A_f * (x_f - x_c) - cellToFaceMat[ ifaceLoc ][ 0 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 0 ]; - cellToFaceMat[ ifaceLoc ][ 1 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 1 ]; - cellToFaceMat[ ifaceLoc ][ 2 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 2 ]; - - // orient normal outward - if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( faceNormal, -1.0 ); - } - - // normalsMat row - normalsMat[ ifaceLoc ][ 0 ] = faceNormal[ 0 ]; - normalsMat[ ifaceLoc ][ 1 ] = faceNormal[ 1 ]; - normalsMat[ ifaceLoc ][ 2 ] = faceNormal[ 2 ]; + LvArray::tensorOps::scale< 3 >( faceNormal, -1.0 ); } - - // 2) compute N - LvArray::tensorOps::Rij_eq_AikBkj< NF, 3, 3 >( work_numFacesByDim, - normalsMat, - permMat ); - - // BdVLM inner product matrix M = M0 + M1 by Beirao da Veiga, Lipnikov, Manzini (pg.89) - // M0 = R (R^T N)^(-1) R^T - - // 3) compute (R^T N)^-1 -> (3 X 3), work_dimByDim = R^T N - LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, NF >( work_dimByDim, - cellToFaceMat, - work_numFacesByDim ); - LvArray::tensorOps::invert< 3 >( work_dimByDim ); - - // 4) compute R (R^T N)^(-1) R^T into M - // work_dimByNumFaces = (R^T N)^-1 R^T -> (3 X NF) - LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work_dimByNumFaces, - work_dimByDim, - cellToFaceMat ); - - // M = R * work_dimByNumFaces (NF x NF) - LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( M, - cellToFaceMat, - work_dimByNumFaces ); - - // P_N = I - N (N^T N)^(-1) N^T - // 5) compute (N^T N)^-1 - LvArray::tensorOps::Rij_eq_AkiAkj< 3, NF >( work_dimByDim, - work_numFacesByDim ); // N - LvArray::tensorOps::invert< 3 >( work_dimByDim ); - - // 6) tmp = (N^T N)^-1 N^T - LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work_dimByNumFaces, - work_dimByDim, - work_numFacesByDim ); - - // 7) work_numFacesByNumFaces = -I + N * tmp - LvArray::tensorOps::addIdentity< NF >( work_numFacesByNumFaces, -1.0 ); - LvArray::tensorOps::Rij_add_AikBkj< NF, NF, 3 >( work_numFacesByNumFaces, - work_numFacesByDim, // N - work_dimByNumFaces // tmp + + // normalsMat row + normalsMat[ ifaceLoc ][ 0 ] = faceNormal[ 0 ]; + normalsMat[ ifaceLoc ][ 1 ] = faceNormal[ 1 ]; + normalsMat[ ifaceLoc ][ 2 ] = faceNormal[ 2 ]; + } + + // 2) compute N + LvArray::tensorOps::Rij_eq_AikBkj< NF, 3, 3 >( work_numFacesByDim, + normalsMat, + permMat ); + + // BdVLM inner product matrix M = M0 + M1 by Beirao da Veiga, Lipnikov, Manzini (pg.89) + // M0 = R (R^T N)^(-1) R^T + + // 3) compute (R^T N)^-1 -> (3 X 3), work_dimByDim = R^T N + LvArray::tensorOps::Rij_eq_AkiBkj< 3, 3, NF >( work_dimByDim, + cellToFaceMat, + work_numFacesByDim ); + LvArray::tensorOps::invert< 3 >( work_dimByDim ); + + // 4) compute R (R^T N)^(-1) R^T into M + // work_dimByNumFaces = (R^T N)^-1 R^T -> (3 X NF) + LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work_dimByNumFaces, + work_dimByDim, + cellToFaceMat ); + + // M = R * work_dimByNumFaces (NF x NF) + LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( M, + cellToFaceMat, + work_dimByNumFaces ); + + // P_N = I - N (N^T N)^(-1) N^T + // 5) compute (N^T N)^-1 + LvArray::tensorOps::Rij_eq_AkiAkj< 3, NF >( work_dimByDim, + work_numFacesByDim ); // N + LvArray::tensorOps::invert< 3 >( work_dimByDim ); + + // 6) tmp = (N^T N)^-1 N^T + LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work_dimByNumFaces, + work_dimByDim, + work_numFacesByDim ); + + // 7) work_numFacesByNumFaces = -I + N * tmp + LvArray::tensorOps::addIdentity< NF >( work_numFacesByNumFaces, -1.0 ); + LvArray::tensorOps::Rij_add_AikBkj< NF, NF, 3 >( work_numFacesByNumFaces, + work_numFacesByDim, // N + work_dimByNumFaces // tmp ); - - // 8) M = M0 + gamma * P_N - real64 const gamma = 1.0 / static_cast< real64 >( NF ); - LvArray::tensorOps::scaledAdd< NF, NF >( M, work_numFacesByNumFaces, -gamma ); - - // convert to flux-based inner product: M_flux = A^{-1} M_vel A^{-1} - // here, A is diag(faceArea) - real64 invA[ NF ]; - for( localIndex i = 0; i < NF; ++i ) - { - invA[ i ] = 1.0 / faceArea[ i ]; - } - for( localIndex i = 0; i < NF; ++i ) + // 8) M = M0 + gamma * P_N + real64 const gamma = 1.0 / static_cast< real64 >( NF ); + LvArray::tensorOps::scaledAdd< NF, NF >( M, work_numFacesByNumFaces, -gamma ); + + // convert to flux-based inner product: M_flux = A^{-1} M_vel A^{-1} + // here, A is diag(faceArea) + real64 invA[ NF ]; + for( localIndex i = 0; i < NF; ++i ) + { + invA[ i ] = 1.0 / faceArea[ i ]; + } + + for( localIndex i = 0; i < NF; ++i ) + { + for( localIndex j = 0; j < NF; ++j ) { - for( localIndex j = 0; j < NF; ++j ) - { - M[ i ][ j ] *= invA[ i ] * invA[ j ]; - } + M[ i ][ j ] *= invA[ i ] * invA[ j ]; } + } } } // end namespace mimeticInnerProduct diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp index 4082941a677..0f6b1b3fe52 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp @@ -62,7 +62,7 @@ class QuasiTPFAInnerProduct : public MimeticInnerProductBase real64 const (&elemPerm)[ 3 ], real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); - + template< localIndex NF > GEOS_HOST_DEVICE static void @@ -105,13 +105,13 @@ template< localIndex NF > GEOS_HOST_DEVICE void QuasiTPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, - ArrayOfArraysView< localIndex const > const & faceToNodes, - arraySlice1d< localIndex const > const & elemToFaces, - arraySlice1d< real64 const > const & elemCenter, - real64 const & elemVolume, - real64 const (&elemPerm)[ 3 ], - real64 const & lengthTolerance, - arraySlice2d< real64 > const & M ) + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ) { real64 const areaTolerance = lengthTolerance * lengthTolerance; @@ -145,12 +145,12 @@ QuasiTPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POS areaTolerance ); real64 cellToFaceVec[3]; - LvArray::tensorOps::copy<3>( cellToFaceVec, faceCenter ); - LvArray::tensorOps::subtract<3>( cellToFaceVec, elemCenter ); + LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); + LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); - if( LvArray::tensorOps::AiBi<3>( cellToFaceVec, faceNormal ) < 0.0 ) + if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) { - LvArray::tensorOps::scale<3>( faceNormal, -1.0 ); + LvArray::tensorOps::scale< 3 >( faceNormal, -1.0 ); } for( int d = 0; d < 3; ++d ) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp index 5f243b8a589..e44e03aafee 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp @@ -67,7 +67,7 @@ class SimpleInnerProduct : public MimeticInnerProductBase real64 const (&elemPerm)[ 3 ], real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); - + template< localIndex NF > GEOS_HOST_DEVICE static void @@ -221,106 +221,106 @@ SimpleInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITI real64 const & lengthTolerance, arraySlice2d< real64 > const & M ) { - real64 const areaTolerance = lengthTolerance * lengthTolerance; - // 0) Helper: diagonal inverse of [k1, k2, k3] - real64 Kinv[3][3] = {{ 0 }}; - for( int i =0; i < 3; ++i) - { - Kinv[i][i] = 1.0 / elemPerm[i]; - } - - // 1) Compute C, N, A - real64 C[ NF ][ 3 ] = {{ 0 }}; - real64 N[ NF ][ 3 ] = {{ 0 }}; - real64 A[ NF ] = { 0.0 }; - - for( localIndex i = 0; i < NF; ++i) - { - real64 fCenter[3], fNormal[3]; - real64 area = computationalGeometry::centroid_3DPolygon( - faceToNodes[elemToFaces[i]], nodePosition, fCenter, fNormal, areaTolerance ); - - real64 vec_cf[3]; - LvArray::tensorOps::copy<3>(vec_cf, fCenter); - LvArray::tensorOps::subtract<3>(vec_cf, elemCenter); - if ( LvArray::tensorOps::AiBi<3>(vec_cf, fNormal) < 0.0 ) - { - LvArray::tensorOps::scale<3>(fNormal, -1.0); - } - - for ( int d = 0; d < 3; ++d) - { - C[i][d] = vec_cf[d]; - N[i][d] = fNormal[d] * area; - } - - A[i] = area; - } - - // 2) Compute C K^{-1} C^T / volume - real64 CKCt[ NF ][ NF ] = {{0}}; - real64 work3xNF[ 3 ][ NF ] = {{0}}; - - LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work3xNF, Kinv, C ); - LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( CKCt, C, work3xNF ); - LvArray::tensorOps::scale< NF, NF >( CKCt, 1.0 / elemVolume ); - - // 3) Q = orth(N / A) - real64 q0[ NF ], q1[ NF ], q2[ NF ]; - real64 Qmat[ NF ][ 3 ]; - for( localIndex i = 0; i < NF; ++i) + real64 const areaTolerance = lengthTolerance * lengthTolerance; + // 0) Helper: diagonal inverse of [k1, k2, k3] + real64 Kinv[3][3] = {{ 0 }}; + for( int i =0; i < 3; ++i ) + { + Kinv[i][i] = 1.0 / elemPerm[i]; + } + + // 1) Compute C, N, A + real64 C[ NF ][ 3 ] = {{ 0 }}; + real64 N[ NF ][ 3 ] = {{ 0 }}; + real64 A[ NF ] = { 0.0 }; + + for( localIndex i = 0; i < NF; ++i ) + { + real64 fCenter[3], fNormal[3]; + real64 area = computationalGeometry::centroid_3DPolygon( + faceToNodes[elemToFaces[i]], nodePosition, fCenter, fNormal, areaTolerance ); + + real64 vec_cf[3]; + LvArray::tensorOps::copy< 3 >( vec_cf, fCenter ); + LvArray::tensorOps::subtract< 3 >( vec_cf, elemCenter ); + if( LvArray::tensorOps::AiBi< 3 >( vec_cf, fNormal ) < 0.0 ) { - q0[i] = N[i][0] / A[i]; - q1[i] = N[i][1] / A[i]; - q2[i] = N[i][2] / A[i]; + LvArray::tensorOps::scale< 3 >( fNormal, -1.0 ); } - - MimeticInnerProductHelpers::orthonormalize< NF >( q0, q1, q2, Qmat ); - - // 4) U = I - Q Q^T - real64 U[ NF ][ NF ] = {{0}}; - LvArray::tensorOps::addIdentity< NF >( U, -1.0 ); - LvArray::tensorOps::Rij_add_AikAjk< NF, 3 >( U, Qmat ); - LvArray::tensorOps::scale< NF, NF >( U, -1.0 ); - - // 5) A^{-1} * U * A^{-1} - real64 invA[ NF ]; - for( localIndex i = 0; i < NF; ++i) + + for( int d = 0; d < 3; ++d ) { - invA[i] = 1.0 / A[i]; + C[i][d] = vec_cf[d]; + N[i][d] = fNormal[d] * area; } - - // temp = A^{-1} * U - real64 temp[ NF ][ NF ]; - for ( localIndex i = 0; i < NF; ++i ) + + A[i] = area; + } + + // 2) Compute C K^{-1} C^T / volume + real64 CKCt[ NF ][ NF ] = {{0}}; + real64 work3xNF[ 3 ][ NF ] = {{0}}; + + LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( work3xNF, Kinv, C ); + LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( CKCt, C, work3xNF ); + LvArray::tensorOps::scale< NF, NF >( CKCt, 1.0 / elemVolume ); + + // 3) Q = orth(N / A) + real64 q0[ NF ], q1[ NF ], q2[ NF ]; + real64 Qmat[ NF ][ 3 ]; + for( localIndex i = 0; i < NF; ++i ) + { + q0[i] = N[i][0] / A[i]; + q1[i] = N[i][1] / A[i]; + q2[i] = N[i][2] / A[i]; + } + + MimeticInnerProductHelpers::orthonormalize< NF >( q0, q1, q2, Qmat ); + + // 4) U = I - Q Q^T + real64 U[ NF ][ NF ] = {{0}}; + LvArray::tensorOps::addIdentity< NF >( U, -1.0 ); + LvArray::tensorOps::Rij_add_AikAjk< NF, 3 >( U, Qmat ); + LvArray::tensorOps::scale< NF, NF >( U, -1.0 ); + + // 5) A^{-1} * U * A^{-1} + real64 invA[ NF ]; + for( localIndex i = 0; i < NF; ++i ) + { + invA[i] = 1.0 / A[i]; + } + + // temp = A^{-1} * U + real64 temp[ NF ][ NF ]; + for( localIndex i = 0; i < NF; ++i ) + { + for( localIndex j = 0; j < NF; ++j ) { - for ( localIndex j = 0; j < NF; ++j ) - { - temp[i][j] = invA[i] * U[i][j]; - } + temp[i][j] = invA[i] * U[i][j]; } - - // Usc = A^{-1} * U * A^{-1} - real64 Usc[ NF ][ NF ]; - for ( localIndex i = 0; i < NF; ++i ) + } + + // Usc = A^{-1} * U * A^{-1} + real64 Usc[ NF ][ NF ]; + for( localIndex i = 0; i < NF; ++i ) + { + for( localIndex j = 0; j < NF; ++j ) { - for ( localIndex j = 0; j < NF; ++j ) - { - Usc[i][j] = temp[i][j] * invA[j]; - } + Usc[i][j] = temp[i][j] * invA[j]; } - - // 6) M = CKCt + (v / t) * Usc - real64 tParam = 2.0 * ( elemPerm[0] + elemPerm[1] + elemPerm[2] ); - real64 scale = elemVolume / tParam; - - for( localIndex i = 0; i < NF; ++i ) + } + + // 6) M = CKCt + (v / t) * Usc + real64 tParam = 2.0 * ( elemPerm[0] + elemPerm[1] + elemPerm[2] ); + real64 scale = elemVolume / tParam; + + for( localIndex i = 0; i < NF; ++i ) + { + for( localIndex j = 0; j < NF; ++j ) { - for ( localIndex j = 0; j < NF; ++j) - { - M[i][j] = CKCt[i][j] + scale * Usc[i][j]; - } + M[i][j] = CKCt[i][j] + scale * Usc[i][j]; } + } } diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp index e9bb2a5a159..9e23dc1649b 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp @@ -60,18 +60,18 @@ class TPFAInnerProduct : public MimeticInnerProductBase real64 const (&elemPerm)[ 3 ], real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); - - template< localIndex NF> + + template< localIndex NF > GEOS_HOST_DEVICE static void computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, - ArrayOfArraysView< localIndex const > const & faceToNodes, - arraySlice1d< localIndex const > const & elemToFaces, - arraySlice1d< real64 const > const & elemCenter, - real64 const & elemVolume, - real64 const (&elemPerm)[ 3 ], - real64 const & lengthTolerance, - arraySlice2d< real64 > const & M ); + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ); }; template< localIndex NF > @@ -144,61 +144,61 @@ template< localIndex NF > GEOS_HOST_DEVICE void TPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, - ArrayOfArraysView< localIndex const > const & faceToNodes, - arraySlice1d< localIndex const > const & elemToFaces, - arraySlice1d< real64 const > const & elemCenter, - real64 const & elemVolume, - real64 const (&elemPerm)[ 3 ], - real64 const & lengthTolerance, - arraySlice2d< real64 > const & M ) + ArrayOfArraysView< localIndex const > const & faceToNodes, + arraySlice1d< localIndex const > const & elemToFaces, + arraySlice1d< real64 const > const & elemCenter, + real64 const & elemVolume, + real64 const (&elemPerm)[ 3 ], + real64 const & lengthTolerance, + arraySlice2d< real64 > const & M ) { GEOS_UNUSED_VAR( elemVolume ); real64 const areaTolerance = lengthTolerance * lengthTolerance; real64 const weightTolerance = 1e-30 * lengthTolerance; - + // initialize M to zero - for( localIndex i = 0; i < NF; ++i ) + for( localIndex i = 0; i < NF; ++i ) + { + for( localIndex j = 0; j < NF; ++j ) { - for ( localIndex j = 0; j < NF; ++j) - { - M[i][j] = 0.0; - } + M[i][j] = 0.0; } - - for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + } + + for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) + { + real64 faceCenter[3], faceNormal[3], faceConormal[3], cellToFaceVec[3]; + + // 1) face geometry + real64 const faceArea = + computationalGeometry::centroid_3DPolygon( faceToNodes[elemToFaces[ifaceLoc]], + nodePosition, + faceCenter, + faceNormal, + areaTolerance ); + + LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); + LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); + + if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) { - real64 faceCenter[3], faceNormal[3], faceConormal[3], cellToFaceVec[3]; - - // 1) face geometry - real64 const faceArea = - computationalGeometry::centroid_3DPolygon( faceToNodes[elemToFaces[ifaceLoc]], - nodePosition, - faceCenter, - faceNormal, - areaTolerance ); - - LvArray::tensorOps::copy<3>(cellToFaceVec, faceCenter); - LvArray::tensorOps::subtract<3>(cellToFaceVec, elemCenter); - - if( LvArray::tensorOps::AiBi<3>(cellToFaceVec, faceNormal) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( faceNormal, -1 ); - } - - real64 const c2fDistance = LvArray::tensorOps::normalize<3>( cellToFaceVec ); - - // 2) K * n (TPFA assumes diagonal K) - LvArray::tensorOps::hadamardProduct<3>(faceConormal, elemPerm, faceNormal); - - // 3) compute T_ii - real64 Tii = LvArray::tensorOps::AiBi<3>( cellToFaceVec, faceConormal ) * faceArea / c2fDistance; - - Tii = LvArray::math::max( Tii, weightTolerance ); - - // 4) M = |T|^{-1} - M[ifaceLoc][ifaceLoc] = 1.0 / LvArray::math::abs( Tii ); + LvArray::tensorOps::scale< 3 >( faceNormal, -1 ); } + + real64 const c2fDistance = LvArray::tensorOps::normalize< 3 >( cellToFaceVec ); + + // 2) K * n (TPFA assumes diagonal K) + LvArray::tensorOps::hadamardProduct< 3 >( faceConormal, elemPerm, faceNormal ); + + // 3) compute T_ii + real64 Tii = LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceConormal ) * faceArea / c2fDistance; + + Tii = LvArray::math::max( Tii, weightTolerance ); + + // 4) M = |T|^{-1} + M[ifaceLoc][ifaceLoc] = 1.0 / LvArray::math::abs( Tii ); + } } diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp index da4ae05498f..ea4e7b0ebb7 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp @@ -63,7 +63,7 @@ void compareTransmissibilityMatrices( arraySlice2d< real64 const > const & trans } void compareInnerProductMatrices( arraySlice2d< real64 const > const & M, - arraySlice2d< real64 const > const & Mref ) + arraySlice2d< real64 const > const & Mref ) { for( localIndex ifaceLoc = 0; ifaceLoc < M.size( 0 ); ++ifaceLoc ) { @@ -567,7 +567,7 @@ void makeCube( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodePosition, // elem-to-faces elemToFaces.resize( numFaces ); for( localIndex i = 0; i < numFaces; ++i ) - elemToFaces(i) = i; + elemToFaces( i ) = i; // face-to-nodes faceToNodes.resize( numFaces ); @@ -577,38 +577,38 @@ void makeCube( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodePosition, // vertices nodePosition.resize( numNodes, 3 ); - nodePosition(0,0) = 0; nodePosition(0,1) = 0; nodePosition(0,2) = 0; - nodePosition(1,0) = 1; nodePosition(1,1) = 0; nodePosition(1,2) = 0; - nodePosition(2,0) = 1; nodePosition(2,1) = 1; nodePosition(2,2) = 0; - nodePosition(3,0) = 0; nodePosition(3,1) = 1; nodePosition(3,2) = 0; - nodePosition(4,0) = 0; nodePosition(4,1) = 0; nodePosition(4,2) = 1; - nodePosition(5,0) = 1; nodePosition(5,1) = 0; nodePosition(5,2) = 1; - nodePosition(6,0) = 1; nodePosition(6,1) = 1; nodePosition(6,2) = 1; - nodePosition(7,0) = 0; nodePosition(7,1) = 1; nodePosition(7,2) = 1; + nodePosition( 0, 0 ) = 0; nodePosition( 0, 1 ) = 0; nodePosition( 0, 2 ) = 0; + nodePosition( 1, 0 ) = 1; nodePosition( 1, 1 ) = 0; nodePosition( 1, 2 ) = 0; + nodePosition( 2, 0 ) = 1; nodePosition( 2, 1 ) = 1; nodePosition( 2, 2 ) = 0; + nodePosition( 3, 0 ) = 0; nodePosition( 3, 1 ) = 1; nodePosition( 3, 2 ) = 0; + nodePosition( 4, 0 ) = 0; nodePosition( 4, 1 ) = 0; nodePosition( 4, 2 ) = 1; + nodePosition( 5, 0 ) = 1; nodePosition( 5, 1 ) = 0; nodePosition( 5, 2 ) = 1; + nodePosition( 6, 0 ) = 1; nodePosition( 6, 1 ) = 1; nodePosition( 6, 2 ) = 1; + nodePosition( 7, 0 ) = 0; nodePosition( 7, 1 ) = 1; nodePosition( 7, 2 ) = 1; // x = 0 - faceToNodes(0,0) = 0; faceToNodes(0,1) = 3; faceToNodes(0,2) = 7; faceToNodes(0,3) = 4; + faceToNodes( 0, 0 ) = 0; faceToNodes( 0, 1 ) = 3; faceToNodes( 0, 2 ) = 7; faceToNodes( 0, 3 ) = 4; // x = 1 - faceToNodes(1,0) = 1; faceToNodes(1,1) = 2; faceToNodes(1,2) = 6; faceToNodes(1,3) = 5; + faceToNodes( 1, 0 ) = 1; faceToNodes( 1, 1 ) = 2; faceToNodes( 1, 2 ) = 6; faceToNodes( 1, 3 ) = 5; // y = 0 - faceToNodes(2,0) = 0; faceToNodes(2,1) = 1; faceToNodes(2,2) = 5; faceToNodes(2,3) = 4; + faceToNodes( 2, 0 ) = 0; faceToNodes( 2, 1 ) = 1; faceToNodes( 2, 2 ) = 5; faceToNodes( 2, 3 ) = 4; // y = 1 - faceToNodes(3,0) = 3; faceToNodes(3,1) = 2; faceToNodes(3,2) = 6; faceToNodes(3,3) = 7; + faceToNodes( 3, 0 ) = 3; faceToNodes( 3, 1 ) = 2; faceToNodes( 3, 2 ) = 6; faceToNodes( 3, 3 ) = 7; // z = 0 - faceToNodes(4,0) = 0; faceToNodes(4,1) = 1; faceToNodes(4,2) = 2; faceToNodes(4,3) = 3; + faceToNodes( 4, 0 ) = 0; faceToNodes( 4, 1 ) = 1; faceToNodes( 4, 2 ) = 2; faceToNodes( 4, 3 ) = 3; // z = 1 - faceToNodes(5,0) = 4; faceToNodes(5,1) = 5; faceToNodes(5,2) = 6; faceToNodes(5,3) = 7; + faceToNodes( 5, 0 ) = 4; faceToNodes( 5, 1 ) = 5; faceToNodes( 5, 2 ) = 6; faceToNodes( 5, 3 ) = 7; // center & volume array1d< localIndex > toNodes; toNodes.resize( numNodes ); for( localIndex i = 0; i < numNodes; ++i ) - toNodes(i) = i; + toNodes( i ) = i; computeVolumeAndCenter( nodePosition, toNodes, @@ -620,30 +620,30 @@ void makeCube( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodePosition, ipType == InnerProductType::QUASI_TPFA ) { // diagonal - transMatrixRef(0,0) = 0.5; - transMatrixRef(1,1) = 0.5; - transMatrixRef(2,2) = 0.5; - transMatrixRef(3,3) = 0.5; - transMatrixRef(4,4) = 0.5; - transMatrixRef(5,5) = 0.5; + transMatrixRef( 0, 0 ) = 0.5; + transMatrixRef( 1, 1 ) = 0.5; + transMatrixRef( 2, 2 ) = 0.5; + transMatrixRef( 3, 3 ) = 0.5; + transMatrixRef( 4, 4 ) = 0.5; + transMatrixRef( 5, 5 ) = 0.5; } else if( ipType == InnerProductType::SIMPLE || ipType == InnerProductType::BDVLM ) { // SIMPLE = TPFA + stabilization - transMatrixRef(0,0) = 0.333333333333; - transMatrixRef(1,1) = 0.333333333333; - transMatrixRef(2,2) = 0.333333333333; - transMatrixRef(3,3) = 0.333333333333; - transMatrixRef(4,4) = 0.333333333333; - transMatrixRef(5,5) = 0.333333333333; - - transMatrixRef(0,1) = -0.166666666666; - transMatrixRef(1,0) = -0.166666666666; - transMatrixRef(2,3) = -0.166666666666; - transMatrixRef(3,2) = -0.166666666666; - transMatrixRef(4,5) = -0.166666666666; - transMatrixRef(5,4) = -0.166666666666; + transMatrixRef( 0, 0 ) = 0.333333333333; + transMatrixRef( 1, 1 ) = 0.333333333333; + transMatrixRef( 2, 2 ) = 0.333333333333; + transMatrixRef( 3, 3 ) = 0.333333333333; + transMatrixRef( 4, 4 ) = 0.333333333333; + transMatrixRef( 5, 5 ) = 0.333333333333; + + transMatrixRef( 0, 1 ) = -0.166666666666; + transMatrixRef( 1, 0 ) = -0.166666666666; + transMatrixRef( 2, 3 ) = -0.166666666666; + transMatrixRef( 3, 2 ) = -0.166666666666; + transMatrixRef( 4, 5 ) = -0.166666666666; + transMatrixRef( 5, 4 ) = -0.166666666666; } } @@ -1158,13 +1158,13 @@ TEST( testMimeticInnerProducts, TPFAM_cube ) real64 const perm[ 3 ] = { elemPerm[0], elemPerm[1], elemPerm[2] }; TPFAInnerProduct::computeM< NF >( nodePosition.toViewConst(), - faceToNodes.toViewConst(), - elemToFaces.toSliceConst(), - center, - 1, - perm, - lengthTolerance, - M.toSlice() ); + faceToNodes.toViewConst(), + elemToFaces.toSliceConst(), + center, + 1, + perm, + lengthTolerance, + M.toSlice() ); compareInnerProductMatrices( M.toSliceConst(), Mref.toSliceConst() ); @@ -1205,13 +1205,13 @@ TEST( testMimeticInnerProducts, QuasiTPFAM_cube ) real64 const perm[ 3 ] = { elemPerm[0], elemPerm[1], elemPerm[2] }; QuasiTPFAInnerProduct::computeM< NF >( nodePosition.toViewConst(), - faceToNodes.toViewConst(), - elemToFaces.toSliceConst(), - center, - 1, - perm, - lengthTolerance, - M.toSlice() ); + faceToNodes.toViewConst(), + elemToFaces.toSliceConst(), + center, + 1, + perm, + lengthTolerance, + M.toSlice() ); compareInnerProductMatrices( M.toSliceConst(), Mref.toSliceConst() ); @@ -1307,7 +1307,7 @@ TEST( testMimeticInnerProducts, BdVLMM_cube ) perm, lengthTolerance, M.toSlice() ); - + compareInnerProductMatrices( M.toSliceConst(), Mref.toSliceConst() ); } From 817ed6a5fb5b280ae1c04fc7a12597df40fca14e Mon Sep 17 00:00:00 2001 From: Nicola Castelletto Date: Sat, 14 Mar 2026 12:35:15 -0700 Subject: [PATCH 03/11] Added helper functions --- .../BdVLMInnerProduct.hpp | 21 ++--- .../MimeticInnerProductHelpers.hpp | 32 ++++++++ .../QuasiTPFAInnerProduct.hpp | 10 +-- .../SimpleInnerProduct.hpp | 79 +++++++------------ .../mimeticInnerProducts/TPFAInnerProduct.hpp | 26 ++---- 5 files changed, 74 insertions(+), 94 deletions(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp index aa96cae10de..12fa949b753 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp @@ -124,18 +124,13 @@ BdVLMInnerProduct::compute( arrayView2d< real64 const, nodes::REFERENCE_POSITION faceNormal, areaTolerance ); - LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); - LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); + MimeticInnerProductHelpers::computeCellToFacetVector( cellToFaceVec, faceCenter, elemCenter ); + MimeticInnerProductHelpers::orientNormalOutward( cellToFaceVec, faceNormal ); cellToFaceMat[ ifaceLoc ][0] = faceAreaMat[ ifaceLoc ][ ifaceLoc ] * cellToFaceVec[ 0 ]; cellToFaceMat[ ifaceLoc ][1] = faceAreaMat[ ifaceLoc ][ ifaceLoc ] * cellToFaceVec[ 1 ]; cellToFaceMat[ ifaceLoc ][2] = faceAreaMat[ ifaceLoc ][ ifaceLoc ] * cellToFaceVec[ 2 ]; - if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( faceNormal, -1 ); - } - // the two-point transmissibility is computed to computed here because it is needed // in the implementation of the transmissibility multiplier (see below) // TODO: see what it would take to bring the (harmonically averaged) two-point trans here @@ -244,7 +239,7 @@ BdVLMInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITIO real64 work_numFacesByDim[ NF ][ 3 ] = {{ 0 }}; real64 work_dimByNumFaces[ 3 ][ NF ] = {{ 0 }}; real64 work_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; - real64 tmp_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; + // real64 tmp_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; // 0) assemble full coefficient tensor from principal axis/components MimeticInnerProductHelpers::makeFullTensor( elemPerm, permMat ); @@ -263,20 +258,14 @@ BdVLMInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITIO faceArea[ ifaceLoc ] = faceAreaMat[ ifaceLoc ][ ifaceLoc ]; - LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); - LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); + MimeticInnerProductHelpers::computeCellToFacetVector( cellToFaceVec, faceCenter, elemCenter ); + MimeticInnerProductHelpers::orientNormalOutward( cellToFaceVec, faceNormal ); // R row: A_f * (x_f - x_c) cellToFaceMat[ ifaceLoc ][ 0 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 0 ]; cellToFaceMat[ ifaceLoc ][ 1 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 1 ]; cellToFaceMat[ ifaceLoc ][ 2 ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 2 ]; - // orient normal outward - if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( faceNormal, -1.0 ); - } - // normalsMat row normalsMat[ ifaceLoc ][ 0 ] = faceNormal[ 0 ]; normalsMat[ ifaceLoc ][ 1 ] = faceNormal[ 1 ]; diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp index 7eb8f5558c6..aa240a1ec65 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp @@ -48,6 +48,38 @@ struct MimeticInnerProductHelpers result[ 2 ][ 2 ] = values[ 2 ]; } + /** + * @brief Compute the vector from cell center to facet center. + * @param[out] cellToFacetVec output vector + * @param[in] facetCenter facet center coordinate + * @param[in] cellCenter cell center coordinate + */ + GEOS_HOST_DEVICE + static + void computeCellToFacetVector( real64 (& cellToFacetVec)[ 3 ], + real64 const (& facetCenter)[ 3 ], + arraySlice1d< real64 const > const & cellCenter ) + { + LvArray::tensorOps::copy< 3 >( cellToFacetVec, facetCenter ); + LvArray::tensorOps::subtract< 3 >( cellToFacetVec, cellCenter ); + } + + /** + * @brief Ensure the facet normal points outward from the cell. + * @param[in] cellToFacetVec vector from cell center to face center + * @param[in,out] faceNormal face normal (may be flipped) + */ + GEOS_HOST_DEVICE + static + void orientNormalOutward( real64 const (& cellToFacetVec)[ 3 ], + real64 (& faceNormal)[ 3 ] ) + { + if( LvArray::tensorOps::AiBi< 3 >( cellToFacetVec, faceNormal ) < 0.0 ) + { + LvArray::tensorOps::scale< 3 >( faceNormal, -1.0 ); + } + } + /** * @brief Orthonormalize a set of three vectors * @tparam NF number of faces in the element diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp index 0f6b1b3fe52..853b4af58d1 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp @@ -21,6 +21,7 @@ #define GEOS_FINITEVOLUME_MIMETICINNERPRODUCTS_QUASITPFAINNERPRODUCT_HPP_ #include "finiteVolume/mimeticInnerProducts/MimeticInnerProductBase.hpp" +#include "finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp" namespace geos { @@ -145,13 +146,8 @@ QuasiTPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POS areaTolerance ); real64 cellToFaceVec[3]; - LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); - LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); - - if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( faceNormal, -1.0 ); - } + MimeticInnerProductHelpers::computeCellToFacetVector( cellToFaceVec, faceCenter, elemCenter ); + MimeticInnerProductHelpers::orientNormalOutward( cellToFaceVec, faceNormal ); for( int d = 0; d < 3; ++d ) { diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp index e44e03aafee..57178fb773b 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp @@ -127,18 +127,13 @@ SimpleInnerProduct::compute( arrayView2d< real64 const, nodes::REFERENCE_POSITIO faceNormal, areaTolerance ); - LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); - LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); + MimeticInnerProductHelpers::computeCellToFacetVector( cellToFaceVec, faceCenter, elemCenter ); + MimeticInnerProductHelpers::orientNormalOutward( cellToFaceVec, faceNormal ); q0[ ifaceLoc ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 0 ]; q1[ ifaceLoc ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 1 ]; q2[ ifaceLoc ] = faceArea[ ifaceLoc ] * cellToFaceVec[ 2 ]; - if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( faceNormal, -1 ); - } - // the two-point transmissibility is computed to computed here because it is needed // in the implementation of the transmissibility multiplier (see below) // TODO: see what it would take to bring the (harmonically averaged) two-point trans here @@ -241,12 +236,8 @@ SimpleInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITI faceToNodes[elemToFaces[i]], nodePosition, fCenter, fNormal, areaTolerance ); real64 vec_cf[3]; - LvArray::tensorOps::copy< 3 >( vec_cf, fCenter ); - LvArray::tensorOps::subtract< 3 >( vec_cf, elemCenter ); - if( LvArray::tensorOps::AiBi< 3 >( vec_cf, fNormal ) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( fNormal, -1.0 ); - } + MimeticInnerProductHelpers::computeCellToFacetVector( vec_cf, fCenter, elemCenter ); + MimeticInnerProductHelpers::orientNormalOutward( vec_cf, fNormal ); for( int d = 0; d < 3; ++d ) { @@ -277,48 +268,36 @@ SimpleInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITI MimeticInnerProductHelpers::orthonormalize< NF >( q0, q1, q2, Qmat ); - // 4) U = I - Q Q^T - real64 U[ NF ][ NF ] = {{0}}; - LvArray::tensorOps::addIdentity< NF >( U, -1.0 ); - LvArray::tensorOps::Rij_add_AikAjk< NF, 3 >( U, Qmat ); - LvArray::tensorOps::scale< NF, NF >( U, -1.0 ); - - // 5) A^{-1} * U * A^{-1} - real64 invA[ NF ]; - for( localIndex i = 0; i < NF; ++i ) + // 4) M = CKCt + (v / t) * A^{-1} * ( I - Q Q^T ) * A^{-1} + real64 invA[NF]; + for (localIndex i = 0; i < NF; ++i) { - invA[i] = 1.0 / A[i]; + invA[i] = real64(1) / A[i]; } - - // temp = A^{-1} * U - real64 temp[ NF ][ NF ]; - for( localIndex i = 0; i < NF; ++i ) + + // scale = elemVolume / tParam, where tParam = 2 * trace(K) + real64 const tParam = real64(2) * (elemPerm[0] + elemPerm[1] + elemPerm[2]); + real64 const scale = elemVolume / tParam; + + for (localIndex i = 0; i < NF; ++i) { - for( localIndex j = 0; j < NF; ++j ) - { - temp[i][j] = invA[i] * U[i][j]; - } - } - - // Usc = A^{-1} * U * A^{-1} - real64 Usc[ NF ][ NF ]; - for( localIndex i = 0; i < NF; ++i ) - { - for( localIndex j = 0; j < NF; ++j ) - { - Usc[i][j] = temp[i][j] * invA[j]; - } - } - - // 6) M = CKCt + (v / t) * Usc - real64 tParam = 2.0 * ( elemPerm[0] + elemPerm[1] + elemPerm[2] ); - real64 scale = elemVolume / tParam; + real64 const invAiScaled = scale * invA[i]; - for( localIndex i = 0; i < NF; ++i ) - { - for( localIndex j = 0; j < NF; ++j ) + real64 const qi0 = Qmat[i][0]; + real64 const qi1 = Qmat[i][1]; + real64 const qi2 = Qmat[i][2]; + + for (localIndex j = i; j < NF; ++j) { - M[i][j] = CKCt[i][j] + scale * Usc[i][j]; + // qdot = (Q Q^T)_{ij} = sum_k Qik * Qjk, with k in [0,2] + real64 const qdot = qi0*Qmat[j][0] + qi1*Qmat[j][1] + qi2*Qmat[j][2]; + + real64 const u = (i == j ? real64(1) : real64(0)) - qdot; + + real64 const mij = CKCt[i][j] + (invAiScaled * invA[j]) * u; + + M[i][j] = mij; + M[j][i] = mij; // symmetry fill } } } diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp index 9e23dc1649b..2567544c271 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp @@ -115,13 +115,8 @@ TPFAInnerProduct::compute( arrayView2d< real64 const, nodes::REFERENCE_POSITION_ faceNormal, areaTolerance ); - LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); - LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); - - if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( faceNormal, -1 ); - } + MimeticInnerProductHelpers::computeCellToFacetVector( cellToFaceVec, faceCenter, elemCenter ); + MimeticInnerProductHelpers::orientNormalOutward( cellToFaceVec, faceNormal ); real64 const c2fDistance = LvArray::tensorOps::normalize< 3 >( cellToFaceVec ); @@ -158,13 +153,7 @@ TPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION real64 const weightTolerance = 1e-30 * lengthTolerance; // initialize M to zero - for( localIndex i = 0; i < NF; ++i ) - { - for( localIndex j = 0; j < NF; ++j ) - { - M[i][j] = 0.0; - } - } + LvArray::tensorOps::fill< NF, NF >( M, 0.0 ); for( localIndex ifaceLoc = 0; ifaceLoc < NF; ++ifaceLoc ) { @@ -178,13 +167,8 @@ TPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITION faceNormal, areaTolerance ); - LvArray::tensorOps::copy< 3 >( cellToFaceVec, faceCenter ); - LvArray::tensorOps::subtract< 3 >( cellToFaceVec, elemCenter ); - - if( LvArray::tensorOps::AiBi< 3 >( cellToFaceVec, faceNormal ) < 0.0 ) - { - LvArray::tensorOps::scale< 3 >( faceNormal, -1 ); - } + MimeticInnerProductHelpers::computeCellToFacetVector( cellToFaceVec, faceCenter, elemCenter ); + MimeticInnerProductHelpers::orientNormalOutward( cellToFaceVec, faceNormal ); real64 const c2fDistance = LvArray::tensorOps::normalize< 3 >( cellToFaceVec ); From 5560168c0d6bba2264833ae9d89cda958ce973e7 Mon Sep 17 00:00:00 2001 From: "Yoojin S. Cha" Date: Mon, 16 Mar 2026 16:20:03 -0700 Subject: [PATCH 04/11] unify geometry variable naming --- .../mimeticInnerProducts/SimpleInnerProduct.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp index 57178fb773b..ce1037a9259 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp @@ -231,18 +231,18 @@ SimpleInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITI for( localIndex i = 0; i < NF; ++i ) { - real64 fCenter[3], fNormal[3]; + real64 faceCenter[3], faceNormal[3]; real64 area = computationalGeometry::centroid_3DPolygon( - faceToNodes[elemToFaces[i]], nodePosition, fCenter, fNormal, areaTolerance ); + faceToNodes[elemToFaces[i]], nodePosition, faceCenter, faceNormal, areaTolerance ); - real64 vec_cf[3]; - MimeticInnerProductHelpers::computeCellToFacetVector( vec_cf, fCenter, elemCenter ); - MimeticInnerProductHelpers::orientNormalOutward( vec_cf, fNormal ); + real64 cellToFaceVec[3]; + MimeticInnerProductHelpers::computeCellToFacetVector( cellToFaceVec, faceCenter, elemCenter ); + MimeticInnerProductHelpers::orientNormalOutward( cellToFaceVec, faceNormal ); for( int d = 0; d < 3; ++d ) { - C[i][d] = vec_cf[d]; - N[i][d] = fNormal[d] * area; + C[i][d] = cellToFaceVec[d]; + N[i][d] = faceNormal[d] * area; } A[i] = area; From 3caf737c55ba0b8c41a87c151a1ef6633619592a Mon Sep 17 00:00:00 2001 From: "Yoojin S. Cha" Date: Mon, 16 Mar 2026 16:43:14 -0700 Subject: [PATCH 05/11] apply uncrustify formatting --- .../MimeticInnerProductHelpers.hpp | 4 ++-- .../SimpleInnerProduct.hpp | 24 +++++++++---------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp index aa240a1ec65..54afd79bdcf 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp @@ -57,7 +57,7 @@ struct MimeticInnerProductHelpers GEOS_HOST_DEVICE static void computeCellToFacetVector( real64 (& cellToFacetVec)[ 3 ], - real64 const (& facetCenter)[ 3 ], + real64 const (&facetCenter)[ 3 ], arraySlice1d< real64 const > const & cellCenter ) { LvArray::tensorOps::copy< 3 >( cellToFacetVec, facetCenter ); @@ -71,7 +71,7 @@ struct MimeticInnerProductHelpers */ GEOS_HOST_DEVICE static - void orientNormalOutward( real64 const (& cellToFacetVec)[ 3 ], + void orientNormalOutward( real64 const (&cellToFacetVec)[ 3 ], real64 (& faceNormal)[ 3 ] ) { if( LvArray::tensorOps::AiBi< 3 >( cellToFacetVec, faceNormal ) < 0.0 ) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp index ce1037a9259..e8123bbf8c9 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp @@ -270,32 +270,32 @@ SimpleInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITI // 4) M = CKCt + (v / t) * A^{-1} * ( I - Q Q^T ) * A^{-1} real64 invA[NF]; - for (localIndex i = 0; i < NF; ++i) + for( localIndex i = 0; i < NF; ++i ) { - invA[i] = real64(1) / A[i]; + invA[i] = real64( 1 ) / A[i]; } - + // scale = elemVolume / tParam, where tParam = 2 * trace(K) - real64 const tParam = real64(2) * (elemPerm[0] + elemPerm[1] + elemPerm[2]); + real64 const tParam = real64( 2 ) * (elemPerm[0] + elemPerm[1] + elemPerm[2]); real64 const scale = elemVolume / tParam; - - for (localIndex i = 0; i < NF; ++i) + + for( localIndex i = 0; i < NF; ++i ) { real64 const invAiScaled = scale * invA[i]; real64 const qi0 = Qmat[i][0]; real64 const qi1 = Qmat[i][1]; real64 const qi2 = Qmat[i][2]; - - for (localIndex j = i; j < NF; ++j) + + for( localIndex j = i; j < NF; ++j ) { // qdot = (Q Q^T)_{ij} = sum_k Qik * Qjk, with k in [0,2] real64 const qdot = qi0*Qmat[j][0] + qi1*Qmat[j][1] + qi2*Qmat[j][2]; - - real64 const u = (i == j ? real64(1) : real64(0)) - qdot; - + + real64 const u = (i == j ? real64( 1 ) : real64( 0 )) - qdot; + real64 const mij = CKCt[i][j] + (invAiScaled * invA[j]) * u; - + M[i][j] = mij; M[j][i] = mij; // symmetry fill } From 56659293e11e75af209ca2d6a39be770846b9da4 Mon Sep 17 00:00:00 2001 From: "Yoojin S. Cha" Date: Mon, 16 Mar 2026 17:23:55 -0700 Subject: [PATCH 06/11] add doxygen documentation for computeM functions --- .../mimeticInnerProducts/BdVLMInnerProduct.hpp | 16 +++++++++++++++- .../QuasiTPFAInnerProduct.hpp | 13 +++++++++++++ .../mimeticInnerProducts/SimpleInnerProduct.hpp | 13 +++++++++++++ .../mimeticInnerProducts/TPFAInnerProduct.hpp | 11 +++++++++++ 4 files changed, 52 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp index 12fa949b753..22524666ee5 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp @@ -40,7 +40,7 @@ class BdVLMInnerProduct : public MimeticInnerProductBase /** * @brief In a given element, recompute the transmissibility matrix in a cell using the inner product of Beirao da Veiga, Lipnikov, - * Manzini (page 113) + *Manzini (page 113) * @param[in] nodePosition the position of the nodes * @param[in] transMultiplier the transmissibility multipliers at the mesh faces * @param[in] faceToNodes the map from the face to their nodes @@ -66,6 +66,20 @@ class BdVLMInnerProduct : public MimeticInnerProductBase real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); + /** + * @brief Compute the mimetic inner product matrix M in a given element using the inner product of Beirao da Veiga, Lipnikov, Manzini + *(page 89) + * @param[in] nodePosition the position of the nodes + * @param[in] faceToNodes the map from the face to their nodes + * @param[in] elemToFaces the maps from the one-sided face to the corresponding face + * @param[in] elemCenter the center of the element + * @param[in] elemVolume the volume of the element + * @param[in] elemPerm the permeability in the element + * @param[in] lengthTolerance the tolerance used in the trans calculations + * @param[inout] M the output inner product matrix + * + * @details Reference: Beirao da Veiga, Lipnikov, Manzini, "The mimetic finite-difference method for elliptic problems" + */ template< localIndex NF > GEOS_HOST_DEVICE static void diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp index 853b4af58d1..d966871f120 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp @@ -64,6 +64,19 @@ class QuasiTPFAInnerProduct : public MimeticInnerProductBase real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); + /** + * @brief Compute the mimetic inner product matrix M in a given element using the quasi TPFA inner product. + * @param[in] nodePosition the position of the nodes + * @param[in] faceToNodes the map from the face to their nodes + * @param[in] elemToFaces the maps from the one-sided face to the corresponding face + * @param[in] elemCenter the center of the element + * @param[in] elemVolume the volume of the element + * @param[in] elemPerm the permeability in the element + * @param[in] lengthTolerance the tolerance used in the trans calculations + * @param[inout] M the output inner product matrix + * + * @details Reference: K-A Lie, An Introduction to Reservoir Simulation Using MATLAB/GNU Octave (2019) + */ template< localIndex NF > GEOS_HOST_DEVICE static void diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp index e8123bbf8c9..e3bd12e52a1 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp @@ -68,6 +68,19 @@ class SimpleInnerProduct : public MimeticInnerProductBase real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); + /** + * @brief Compute the mimetic inner product matrix M in a given element using the Simple inner product. + * @param[in] nodePosition the position of the nodes + * @param[in] faceToNodes the map from the face to their nodes + * @param[in] elemToFaces the maps from the one-sided face to the corresponding face + * @param[in] elemCenter the center of the element + * @param[in] elemVolume the volume of the element + * @param[in] elemPerm the permeability in the element + * @param[in] lengthTolerance the tolerance used in the trans calculations + * @param[inout] M the output inner product matrix + * + * @details Reference: K-A Lie, An Introduction to Reservoir Simulation Using MATLAB/GNU Octave (2019) + */ template< localIndex NF > GEOS_HOST_DEVICE static void diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp index 2567544c271..e2c97d905d4 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/TPFAInnerProduct.hpp @@ -61,6 +61,17 @@ class TPFAInnerProduct : public MimeticInnerProductBase real64 const & lengthTolerance, arraySlice2d< real64 > const & transMatrix ); + /** + * @brief Compute the mimetic inner product matrix M in a given element using TPFA. + * @param[in] nodePosition the position of the nodes + * @param[in] faceToNodes the map from the face to their nodes + * @param[in] elemToFaces the maps from the one-sided face to the corresponding face + * @param[in] elemCenter the center of the element + * @param[in] elemVolume the volume of the element + * @param[in] elemPerm the permeability in the element + * @param[in] lengthTolerance the tolerance used in the trans calculations + * @param[inout] M the output inner product matrix + */ template< localIndex NF > GEOS_HOST_DEVICE static void From 41f75ee344fd39ae80968a37b133d8c4a5ed5aac Mon Sep 17 00:00:00 2001 From: "Yoojin S. Cha" Date: Mon, 16 Mar 2026 17:28:50 -0700 Subject: [PATCH 07/11] apply uncrustify formatting --- .../finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp index 22524666ee5..f33449bd780 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp @@ -40,7 +40,7 @@ class BdVLMInnerProduct : public MimeticInnerProductBase /** * @brief In a given element, recompute the transmissibility matrix in a cell using the inner product of Beirao da Veiga, Lipnikov, - *Manzini (page 113) + * Manzini (page 113) * @param[in] nodePosition the position of the nodes * @param[in] transMultiplier the transmissibility multipliers at the mesh faces * @param[in] faceToNodes the map from the face to their nodes From d3c5c1c1653501ea91a90ceb2c9f8c436689ce1a Mon Sep 17 00:00:00 2001 From: Omar Duran Date: Mon, 3 Aug 2026 18:10:15 -0700 Subject: [PATCH 08/11] wip: fix qTPFA and BdVLM --- .../BdVLMInnerProduct.hpp | 5 +- .../QuasiTPFAInnerProduct.hpp | 8 +- .../SimpleInnerProduct.hpp | 3 - .../unitTests/testMimeticInnerProducts.cpp | 525 +++++++++++++++++- 4 files changed, 511 insertions(+), 30 deletions(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp index f33449bd780..25114e7b705 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/BdVLMInnerProduct.hpp @@ -253,7 +253,6 @@ BdVLMInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITIO real64 work_numFacesByDim[ NF ][ 3 ] = {{ 0 }}; real64 work_dimByNumFaces[ 3 ][ NF ] = {{ 0 }}; real64 work_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; - // real64 tmp_numFacesByNumFaces[ NF ][ NF ] = {{ 0 }}; // 0) assemble full coefficient tensor from principal axis/components MimeticInnerProductHelpers::makeFullTensor( elemPerm, permMat ); @@ -329,8 +328,8 @@ BdVLMInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POSITIO work_dimByNumFaces // tmp ); - // 8) M = M0 + gamma * P_N - real64 const gamma = 1.0 / static_cast< real64 >( NF ); + // 8) M = M0 + gamma * P_N, gamma = 2/NF * trace(M0) (mirrors stabCoef in compute()) + real64 const gamma = 2.0 / NF * LvArray::tensorOps::trace< NF >( M ); LvArray::tensorOps::scaledAdd< NF, NF >( M, work_numFacesByNumFaces, -gamma ); // convert to flux-based inner product: M_flux = A^{-1} M_vel A^{-1} diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp index d966871f120..a3b873e1e50 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/QuasiTPFAInnerProduct.hpp @@ -184,15 +184,15 @@ QuasiTPFAInnerProduct::computeM( arrayView2d< real64 const, nodes::REFERENCE_POS LvArray::tensorOps::Rij_eq_AikBjk< 3, NF, 3 >( workNK, K, N ); LvArray::tensorOps::Rij_eq_AikBkj< NF, NF, 3 >( W, N, workNK ); - // 5) build Q from C (orthonormal basis of consistency space) + // 5) build Q from N so the stabilization annihilates range(N) (=> M N K = C) real64 q0[ NF ], q1[ NF ], q2[ NF ]; real64 Qmat[ NF ][ 3 ]; for( localIndex i = 0; i < NF; ++i ) { - q0[i] = C[i][0]; - q1[i] = C[i][1]; - q2[i] = C[i][2]; + q0[i] = N[i][0]; + q1[i] = N[i][1]; + q2[i] = N[i][2]; } MimeticInnerProductHelpers::orthonormalize< NF >( q0, q1, q2, Qmat ); diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp index e3bd12e52a1..2283dafc213 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/SimpleInnerProduct.hpp @@ -23,9 +23,6 @@ #include "finiteVolume/mimeticInnerProducts/MimeticInnerProductBase.hpp" #include "finiteVolume/mimeticInnerProducts/MimeticInnerProductHelpers.hpp" #include "mesh/utilities/ComputationalGeometry.hpp" -#include -#include - namespace geos { diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp index cf64a4a869f..c403693464d 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp @@ -24,8 +24,6 @@ #include "mainInterface/initialization.hpp" #include "mesh/FaceManager.hpp" #include "mesh/utilities/ComputationalGeometry.hpp" -#include -#include // TPL includes #include @@ -617,9 +615,10 @@ void makeCube( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodePosition, // reference matrix (MRST / MATLAB) if( ipType == InnerProductType::TPFA || - ipType == InnerProductType::QUASI_TPFA ) + ipType == InnerProductType::QUASI_TPFA || + ipType == InnerProductType::BDVLM ) { - // diagonal + // these inner products reduce to TPFA on the unit cube: M = inv(T) = 0.5 I transMatrixRef( 0, 0 ) = 0.5; transMatrixRef( 1, 1 ) = 0.5; transMatrixRef( 2, 2 ) = 0.5; @@ -627,23 +626,22 @@ void makeCube( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodePosition, transMatrixRef( 4, 4 ) = 0.5; transMatrixRef( 5, 5 ) = 0.5; } - else if( ipType == InnerProductType::SIMPLE || - ipType == InnerProductType::BDVLM ) + else if( ipType == InnerProductType::SIMPLE ) { // SIMPLE = TPFA + stabilization - transMatrixRef( 0, 0 ) = 0.333333333333; - transMatrixRef( 1, 1 ) = 0.333333333333; - transMatrixRef( 2, 2 ) = 0.333333333333; - transMatrixRef( 3, 3 ) = 0.333333333333; - transMatrixRef( 4, 4 ) = 0.333333333333; - transMatrixRef( 5, 5 ) = 0.333333333333; - - transMatrixRef( 0, 1 ) = -0.166666666666; - transMatrixRef( 1, 0 ) = -0.166666666666; - transMatrixRef( 2, 3 ) = -0.166666666666; - transMatrixRef( 3, 2 ) = -0.166666666666; - transMatrixRef( 4, 5 ) = -0.166666666666; - transMatrixRef( 5, 4 ) = -0.166666666666; + transMatrixRef( 0, 0 ) = 1.0 / 3.0; + transMatrixRef( 1, 1 ) = 1.0 / 3.0; + transMatrixRef( 2, 2 ) = 1.0 / 3.0; + transMatrixRef( 3, 3 ) = 1.0 / 3.0; + transMatrixRef( 4, 4 ) = 1.0 / 3.0; + transMatrixRef( 5, 5 ) = 1.0 / 3.0; + + transMatrixRef( 0, 1 ) = -1.0 / 6.0; + transMatrixRef( 1, 0 ) = -1.0 / 6.0; + transMatrixRef( 2, 3 ) = -1.0 / 6.0; + transMatrixRef( 3, 2 ) = -1.0 / 6.0; + transMatrixRef( 4, 5 ) = -1.0 / 6.0; + transMatrixRef( 5, 4 ) = -1.0 / 6.0; } } @@ -1459,7 +1457,8 @@ static inline void computeDistortedVolumeAndCenter( array2d< real64, nodes::REFE center[1] += X[a][1]; center[2] += X[a][2]; } - vol = hexahedronVolume( X ); + // node ordering is axis-swapped w.r.t. hexahedronVolume's convention -> negative sign + vol = std::fabs( hexahedronVolume( X ) ); center[0] /= 8.0; center[1] /= 8.0; center[2] /= 8.0; } @@ -1754,6 +1753,492 @@ TEST( MimeticIP_Linear, Distortion_NonPlanar_LinearPressure ) } } +//======================== Linear Pressure Recovery Test (mixed form, computeM) ============================= +// Same two-cell setup as the hybrid test above, but solving the mixed saddle-point system +// (M q = p_C 1 - pi per cell, mass conservation, interface flux continuity) built from computeM. + +static constexpr real64 mixed_consistency_tol = 1e-13; + +template< int NF > +static void computeM_dispatch( int ipKind, + array2d< real64, nodes::REFERENCE_POSITION_PERM > const & node, + FaceManager::NodeMapType const & faceTonode, + array1d< localIndex > const & elemToface, + arraySlice1d< real64 const > const & center, + real64 const vol, + real64 const (&Kvec)[3], + real64 const ltol, + arraySlice2d< real64 > const & M ) +{ + if( ipKind == InnerProductType::TPFA ) + { + TPFAInnerProduct::computeM< NF >( node.toViewConst(), faceTonode.toViewConst(), elemToface.toSliceConst(), + center, vol, Kvec, ltol, M ); + } + else if( ipKind == InnerProductType::QUASI_TPFA ) + { + QuasiTPFAInnerProduct::computeM< NF >( node.toViewConst(), faceTonode.toViewConst(), elemToface.toSliceConst(), + center, vol, Kvec, ltol, M ); + } + else if( ipKind == InnerProductType::SIMPLE ) + { + SimpleInnerProduct::computeM< NF >( node.toViewConst(), faceTonode.toViewConst(), elemToface.toSliceConst(), + center, vol, Kvec, ltol, M ); + } + else if( ipKind == InnerProductType::BDVLM ) + { + BdVLMInnerProduct::computeM< NF >( node.toViewConst(), faceTonode.toViewConst(), elemToface.toSliceConst(), + center, vol, Kvec, ltol, M ); + } +} + +// solve a small dense linear system A x = b by Gaussian elimination with partial pivoting +template< int N > +static void solveDense( real64 (& A)[N][N], real64 (& b)[N], real64 (& x)[N] ) +{ + for( int k = 0; k < N; ++k ) + { + int piv = k; + for( int i = k + 1; i < N; ++i ) + { + if( std::abs( A[i][k] ) > std::abs( A[piv][k] ) ) + piv = i; + } + if( piv != k ) + { + for( int j = 0; j < N; ++j ) + std::swap( A[k][j], A[piv][j] ); + std::swap( b[k], b[piv] ); + } + for( int i = k + 1; i < N; ++i ) + { + real64 const f = A[i][k] / A[k][k]; + for( int j = k; j < N; ++j ) + A[i][j] -= f * A[k][j]; + b[i] -= f * b[k]; + } + } + for( int i = N - 1; i >= 0; --i ) + { + real64 s = b[i]; + for( int j = i + 1; j < N; ++j ) + s -= A[i][j] * x[j]; + x[i] = s / A[i][i]; + } +} + +// recover the two cell pressures from the mixed (saddle-point) system based on M and +// compare against the analytical values +template< int NF > +static double computeLinearPressureMixed_error( int ipKind, + DistortionMode mode = DistortionMode::None, + real64 eps = 0.0 ) +{ + array2d< real64, nodes::REFERENCE_POSITION_PERM > node_L; + array2d< real64, nodes::REFERENCE_POSITION_PERM > node_R; + FaceManager::NodeMapType faceTonode_L; + FaceManager::NodeMapType faceTonode_R; + array1d< localIndex > elemToface_L; + array1d< localIndex > elemToface_R; + real64 center_L[3], center_R[3], vol_L = 0, vol_R = 0; + + makeUnitCube( 0.0, node_L, faceTonode_L, elemToface_L, center_L, vol_L ); + makeUnitCube( 1.0, node_R, faceTonode_R, elemToface_R, center_R, vol_R ); + + if( mode == DistortionMode::Planar ) + { + makeDistortedPlanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } + else if( mode == DistortionMode::NonPlanar ) + { + makeDistortedNonplanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } + + computeDistortedVolumeAndCenter( node_L, center_L, vol_L ); + computeDistortedVolumeAndCenter( node_R, center_R, vol_R ); + + // indicate the shared face + localIndex const fL_int = 3; + localIndex const fR_int = 2; + + real64 const alpha_lin = 0.3; // pressure = x + alpha_lin * y + constexpr real64 ltol = 1e-12; + real64 Kvec[3] = {4.0, 1.0, 0.5}; + + stackArray1d< real64, 3 > cLc( 3 ), cRc( 3 ); + for( int d = 0; d < 3; ++d ) + { + cLc[d] = center_L[d]; + cRc[d] = center_R[d]; + } + + stackArray2d< real64, NF * NF > ML( NF, NF ), MR( NF, NF ); + ML.template setValues< parallelHostPolicy >( 0.0 ); + MR.template setValues< parallelHostPolicy >( 0.0 ); + + computeM_dispatch< NF >( ipKind, node_L, faceTonode_L, elemToface_L, cLc.toSliceConst(), vol_L, Kvec, ltol, ML.toSlice() ); + computeM_dispatch< NF >( ipKind, node_R, faceTonode_R, elemToface_R, cRc.toSliceConst(), vol_R, Kvec, ltol, MR.toSlice() ); + + // unknown ordering: [ qL(0..NF-1), qR(0..NF-1), pL, pR, piInt ] + constexpr int N = 2 * NF + 3; + int const ipL = 2 * NF, ipR = 2 * NF + 1, iPi = 2 * NF + 2; + real64 A[N][N] = {{ 0 }}; + real64 b[N] = { 0 }; + real64 x[N] = { 0 }; + + // flux equations + for( int i = 0; i < NF; ++i ) + { + for( int j = 0; j < NF; ++j ) + { + A[i][j] = ML( i, j ); + A[NF + i][NF + j] = MR( i, j ); + } + A[i][ipL] = -1.0; + A[NF + i][ipR] = -1.0; + + if( i == fL_int ) + A[i][iPi] = 1.0; + else + b[i] = -computeBoundaryPressure( faceTonode_L, node_L, i, alpha_lin ); + + if( i == fR_int ) + A[NF + i][iPi] = 1.0; + else + b[NF + i] = -computeBoundaryPressure( faceTonode_R, node_R, i, alpha_lin ); + } + + // mass conservation + for( int j = 0; j < NF; ++j ) + { + A[ipL][j] = 1.0; + A[ipR][NF + j] = 1.0; + } + + // flux continuity at the interface + A[iPi][fL_int] = 1.0; + A[iPi][NF + fR_int] = 1.0; + + solveDense< N >( A, b, x ); + + real64 const pL = x[ipL]; + real64 const pR = x[ipR]; + + // compute analytical pressure values + double pL_exact = center_L[0] + alpha_lin * center_L[1]; + double pR_exact = center_R[0] + alpha_lin * center_R[1]; + double err = std::max( std::abs( pL - pL_exact ) / std::abs( pL_exact ), std::abs( pR - pR_exact ) / std::abs( pR_exact ) ); + + return err; +} + +// ================= mixed form, case 0: without distortion =========================== +TEST( MimeticIP_MixedLinear, UnitCube_LinearPressure_TPFA ) +{ + double err = computeLinearPressureMixed_error< 6 >( InnerProductType::TPFA ); + EXPECT_LT( err, mixed_consistency_tol ); +} + +TEST( MimeticIP_MixedLinear, UnitCube_LinearPressure_QuasiTPFA ) +{ + double err = computeLinearPressureMixed_error< 6 >( InnerProductType::QUASI_TPFA ); + EXPECT_LT( err, mixed_consistency_tol ); +} + +TEST( MimeticIP_MixedLinear, UnitCube_LinearPressure_Simple ) +{ + double err = computeLinearPressureMixed_error< 6 >( InnerProductType::SIMPLE ); + EXPECT_LT( err, mixed_consistency_tol ); +} + +TEST( MimeticIP_MixedLinear, UnitCube_LinearPressure_BdVLM ) +{ + double err = computeLinearPressureMixed_error< 6 >( InnerProductType::BDVLM ); + EXPECT_LT( err, mixed_consistency_tol ); +} + +// =================== mixed form, case 1: with distortion (planar) =========================== +TEST( MimeticIP_MixedLinear, Distortion_Planar_LinearPressure ) +{ + int neps = 3; + stdVector< double > eps_values( 3 ); + eps_values[0] = 0.0; // no distortion + eps_values[1] = 0.2; // moderate distortion + eps_values[2] = 0.9; // severe distortion + + // all mimetic inner products must remain consistent on distorted (planar-face) cells + for( int i = 0; i < neps; ++i ) + { + double eps = eps_values[i]; + + double errQTPFA = computeLinearPressureMixed_error< 6 >( InnerProductType::QUASI_TPFA, DistortionMode::Planar, eps ); + EXPECT_LT( errQTPFA, mixed_consistency_tol ); + + double errSIMPLE = computeLinearPressureMixed_error< 6 >( InnerProductType::SIMPLE, DistortionMode::Planar, eps ); + EXPECT_LT( errSIMPLE, mixed_consistency_tol ); + + double errBDVLM = computeLinearPressureMixed_error< 6 >( InnerProductType::BDVLM, DistortionMode::Planar, eps ); + EXPECT_LT( errBDVLM, mixed_consistency_tol ); + } + + for( int i = 0; i < neps; ++i ) + { + double eps = eps_values[i]; + + double errTPFA = computeLinearPressureMixed_error< 6 >( InnerProductType::TPFA, DistortionMode::Planar, eps ); + + if( i == 0 ) + { + EXPECT_LT( errTPFA, mixed_consistency_tol ); // TPFA is consistent on K-orthogonal grids + } + else + { + EXPECT_GT( errTPFA, mixed_consistency_tol ); // TPFA is inconsistent on non K-orthogonal grids + } + } +} + +// =================== mixed form, case 2: with distortion (nonplanar) =========================== +TEST( MimeticIP_MixedLinear, Distortion_NonPlanar_LinearPressure ) +{ + int neps = 2; + stdVector< double > eps_values( 2 ); + eps_values[0] = 0.2; // moderate distortion + eps_values[1] = 0.9; // severe distortion + + // all schemes are inconsistent with nonplanar faces + for( int i = 0; i < neps; ++i ) + { + double eps = eps_values[i]; + + double errTPFA = computeLinearPressureMixed_error< 6 >( InnerProductType::TPFA, DistortionMode::NonPlanar, eps ); + EXPECT_GT( errTPFA, mixed_consistency_tol ); + + double errQTPFA = computeLinearPressureMixed_error< 6 >( InnerProductType::QUASI_TPFA, DistortionMode::NonPlanar, eps ); + EXPECT_GT( errQTPFA, mixed_consistency_tol ); + + double errSIMPLE = computeLinearPressureMixed_error< 6 >( InnerProductType::SIMPLE, DistortionMode::NonPlanar, eps ); + EXPECT_GT( errSIMPLE, mixed_consistency_tol ); + + double errBDVLM = computeLinearPressureMixed_error< 6 >( InnerProductType::BDVLM, DistortionMode::NonPlanar, eps ); + EXPECT_GT( errBDVLM, mixed_consistency_tol ); + } +} + +//======================== Mixed Inner Product Algebraic Consistency Test ============================= +// Direct check of the mimetic consistency condition M (N K) = C on the element matrix +// (stronger than the pressure recovery above, where a defect can hide in the fluxes). +// Valid only on cells where C^T N = vol * I holds with vertex-averaged face centers. + +template< int NF > +static double computeMixedIPConsistency_error( int ipKind, + DistortionMode mode = DistortionMode::None, + real64 eps = 0.0 ) +{ + array2d< real64, nodes::REFERENCE_POSITION_PERM > node_L; + array2d< real64, nodes::REFERENCE_POSITION_PERM > node_R; + FaceManager::NodeMapType faceTonode_L; + FaceManager::NodeMapType faceTonode_R; + array1d< localIndex > elemToface_L; + array1d< localIndex > elemToface_R; + real64 center_L[3], center_R[3], vol_L = 0, vol_R = 0; + + makeUnitCube( 0.0, node_L, faceTonode_L, elemToface_L, center_L, vol_L ); + makeUnitCube( 1.0, node_R, faceTonode_R, elemToface_R, center_R, vol_R ); + + if( mode == DistortionMode::Planar ) + { + makeDistortedPlanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } + else if( mode == DistortionMode::NonPlanar ) + { + makeDistortedNonplanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } + + computeDistortedVolumeAndCenter( node_L, center_L, vol_L ); + + constexpr real64 ltol = 1e-12; + real64 Kvec[3] = {4.0, 1.0, 0.5}; + + stackArray1d< real64, 3 > cLc( 3 ); + for( int d = 0; d < 3; ++d ) + { + cLc[d] = center_L[d]; + } + + stackArray2d< real64, NF * NF > M( NF, NF ); + M.template setValues< parallelHostPolicy >( 0.0 ); + computeM_dispatch< NF >( ipKind, node_L, faceTonode_L, elemToface_L, cLc.toSliceConst(), vol_L, Kvec, ltol, M.toSlice() ); + + // assemble N (area-weighted outward normals, times K) and C (cell-to-face vectors), + // using the same geometry convention as computeM (centroid_3DPolygon) + real64 NK[NF][3]; + real64 C[NF][3]; + for( int f = 0; f < NF; ++f ) + { + real64 fc[3], fn[3]; + real64 const area = centroid_3DPolygon( faceTonode_L[f], node_L.toViewConst(), fc, fn ); + real64 c2f[3] = { fc[0] - center_L[0], fc[1] - center_L[1], fc[2] - center_L[2] }; + if( LvArray::tensorOps::AiBi< 3 >( c2f, fn ) < 0.0 ) + { + LvArray::tensorOps::scale< 3 >( fn, -1.0 ); + } + for( int d = 0; d < 3; ++d ) + { + NK[f][d] = area * fn[d] * Kvec[d]; // (|f| n_f) K for diagonal K + C[f][d] = c2f[d]; + } + } + + // err = max | ( M NK - C )_fd | + double err = 0.0; + for( int f = 0; f < NF; ++f ) + { + for( int d = 0; d < 3; ++d ) + { + real64 s = 0.0; + for( int k = 0; k < NF; ++k ) + { + s += M( f, k ) * NK[k][d]; + } + err = std::max( err, std::abs( s - C[f][d] ) ); + } + } + return err; +} + +TEST( MimeticIP_MixedConsistency, UnitCube ) +{ + for( int ip : { InnerProductType::TPFA, InnerProductType::QUASI_TPFA, + InnerProductType::SIMPLE, InnerProductType::BDVLM } ) + { + EXPECT_LT( computeMixedIPConsistency_error< 6 >( ip ), mixed_consistency_tol ) << "ipKind = " << ip; + } +} + +TEST( MimeticIP_MixedConsistency, Distortion_Planar ) +{ + for( double eps : { 0.2, 0.9 } ) + { + // all mimetic inner products must satisfy M ( N K ) = C on planar-face cells + EXPECT_LT( computeMixedIPConsistency_error< 6 >( InnerProductType::QUASI_TPFA, DistortionMode::Planar, eps ), + mixed_consistency_tol ) << "eps = " << eps; + EXPECT_LT( computeMixedIPConsistency_error< 6 >( InnerProductType::SIMPLE, DistortionMode::Planar, eps ), + mixed_consistency_tol ) << "eps = " << eps; + EXPECT_LT( computeMixedIPConsistency_error< 6 >( InnerProductType::BDVLM, DistortionMode::Planar, eps ), + mixed_consistency_tol ) << "eps = " << eps; + + // TPFA is inconsistent on non K-orthogonal cells + EXPECT_GT( computeMixedIPConsistency_error< 6 >( InnerProductType::TPFA, DistortionMode::Planar, eps ), + mixed_consistency_tol ) << "eps = " << eps; + } +} + +//======================== Hybrid/Mixed Duality Test ============================= +// M = T^{-1} must hold; exact on parallelepipeds with isotropic K at any cell size h. + +template< int NF > +static double computeHybridMixedDuality_error( int ipKind, real64 h ) +{ + array2d< real64, nodes::REFERENCE_POSITION_PERM > node; + FaceManager::NodeMapType faceTonode; + array1d< localIndex > elemToface; + real64 center[3], vol = 0; + + makeUnitCube( 0.0, node, faceTonode, elemToface, center, vol ); + + // scale the cell by h (stays a parallelepiped) + for( int a = 0; a < 8; ++a ) + { + for( int d = 0; d < 3; ++d ) + { + node( a, d ) *= h; + } + } + computeDistortedVolumeAndCenter( node, center, vol ); + + constexpr real64 ltol = 1e-12; + real64 Kvec[3] = {1.0, 1.0, 1.0}; + + stackArray1d< real64, 3 > cc( 3 ); + for( int d = 0; d < 3; ++d ) + { + cc[d] = center[d]; + } + + array1d< real64 > mult( NF ); + mult.setValues< parallelHostPolicy >( 1.0 ); + + stackArray2d< real64, NF * NF > T( NF, NF ), M( NF, NF ); + T.template setValues< parallelHostPolicy >( 0.0 ); + M.template setValues< parallelHostPolicy >( 0.0 ); + + if( ipKind == InnerProductType::TPFA ) + { + TPFAInnerProduct::compute< NF >( node.toViewConst(), mult.toViewConst(), faceTonode.toViewConst(), + elemToface.toSliceConst(), cc, vol, Kvec, ltol, T.toSlice() ); + } + else if( ipKind == InnerProductType::QUASI_TPFA ) + { + QuasiTPFAInnerProduct::compute< NF >( node.toViewConst(), mult.toViewConst(), faceTonode.toViewConst(), + elemToface.toSliceConst(), cc, vol, Kvec, ltol, T.toSlice() ); + } + else if( ipKind == InnerProductType::SIMPLE ) + { + SimpleInnerProduct::compute< NF >( node.toViewConst(), mult.toViewConst(), faceTonode.toViewConst(), + elemToface.toSliceConst(), cc, vol, Kvec, ltol, T.toSlice() ); + } + else if( ipKind == InnerProductType::BDVLM ) + { + BdVLMInnerProduct::compute< NF >( node.toViewConst(), mult.toViewConst(), faceTonode.toViewConst(), + elemToface.toSliceConst(), cc, vol, Kvec, ltol, T.toSlice() ); + } + + computeM_dispatch< NF >( ipKind, node, faceTonode, elemToface, cc.toSliceConst(), vol, Kvec, ltol, M.toSlice() ); + + // err = max | (M T)_ij - delta_ij | + double err = 0.0; + for( int i = 0; i < NF; ++i ) + { + for( int j = 0; j < NF; ++j ) + { + real64 s = 0.0; + for( int k = 0; k < NF; ++k ) + { + s += M( i, k ) * T( k, j ); + } + err = std::max( err, std::abs( s - ( i == j ? 1.0 : 0.0 ) ) ); + } + } + return err; +} + +static constexpr real64 duality_tol = 1e-10; + +TEST( MimeticIP_MixedDuality, TPFA ) +{ + EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::TPFA, 1.0 ), duality_tol ); + EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::TPFA, 0.2 ), duality_tol ); +} + +TEST( MimeticIP_MixedDuality, QuasiTPFA ) +{ + EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::QUASI_TPFA, 1.0 ), duality_tol ); + EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::QUASI_TPFA, 0.2 ), duality_tol ); +} + +TEST( MimeticIP_MixedDuality, Simple ) +{ + EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::SIMPLE, 1.0 ), duality_tol ); + EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::SIMPLE, 0.2 ), duality_tol ); +} + +TEST( MimeticIP_MixedDuality, BdVLM ) +{ + EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::BDVLM, 1.0 ), duality_tol ); + EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::BDVLM, 0.2 ), duality_tol ); +} + //======================== Hydrostatic Equilibrium Consistency Test ============================= static inline void makeDistortedPlanar_gravity( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeL, array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeU, From 5b450de13f219ac89003bd213675a9f2e7d66cfb Mon Sep 17 00:00:00 2001 From: Omar Duran Date: Mon, 3 Aug 2026 18:36:27 -0700 Subject: [PATCH 09/11] fix: centroid computation exact for planar polygons in R3 --- .../unitTests/testMimeticInnerProducts.cpp | 193 ++++++++++++++++-- .../mesh/utilities/ComputationalGeometry.hpp | 18 ++ 2 files changed, 195 insertions(+), 16 deletions(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp index c403693464d..abcf66f7d03 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp @@ -1313,7 +1313,7 @@ TEST( testMimeticInnerProducts, BdVLMM_cube ) //======================== Linear Pressure Recovery Test ============================= // Three cases: without distortion (unit cube cell), with distortion: (1) planar (2) nonplanar -enum class DistortionMode { None, Planar, NonPlanar }; +enum class DistortionMode { None, Planar, PlanarFunnel, NonPlanar }; static inline void makeUnitCube( double x_start, array2d< real64, nodes::REFERENCE_POSITION_PERM > & node, @@ -1401,6 +1401,25 @@ static inline void makeDistortedPlanar( array2d< real64, nodes::REFERENCE_POSITI nodeL( Lvert2, 0 ) += eps; nodeR( Rvert2, 0 ) += eps; } +// shrink the shared face toward its center: two funnel-shaped cells, all faces planar +static inline void makeDistortedPlanarFunnel( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeL, + array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeR, + FaceManager::NodeMapType const & faceL, + FaceManager::NodeMapType const & faceR, + real64 eps ) +{ + localIndex const fL = 3, fR = 2; + for( int a = 0; a < 4; ++a ) + { + localIndex const vL = faceL( fL, a ); + localIndex const vR = faceR( fR, a ); + nodeL( vL, 1 ) += eps * ( 0.5 - nodeL( vL, 1 ) ); + nodeL( vL, 2 ) += eps * ( 0.5 - nodeL( vL, 2 ) ); + nodeR( vR, 1 ) += eps * ( 0.5 - nodeR( vR, 1 ) ); + nodeR( vR, 2 ) += eps * ( 0.5 - nodeR( vR, 2 ) ); + } +} + // create nonplanar case for distortion test static inline void makeDistortedNonplanar( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeL, array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeR, @@ -1462,19 +1481,6 @@ static inline void computeDistortedVolumeAndCenter( array2d< real64, nodes::REFE center[0] /= 8.0; center[1] /= 8.0; center[2] /= 8.0; } -static inline bool isBoundaryFace( FaceManager::NodeMapType const & faceTonode, - array2d< real64, nodes::REFERENCE_POSITION_PERM > const & node, - localIndex f ) -{ - real64 fc[3], fn[3]; - centroid_3DPolygon( faceTonode[f], node.toViewConst(), fc, fn ); - constexpr real64 tol = 1e-12; - - return (std::abs( fc[0] - 0.0 ) < tol) || (std::abs( fc[0] - 2.0 ) < tol) || - (std::abs( fc[1] - 0.0 ) < tol) || (std::abs( fc[1] - 1.0 ) < tol) || - (std::abs( fc[2] - 0.0 ) < tol) || (std::abs( fc[2] - 1.0 ) < tol); -} - // calculate pressures on cell boundary (for Dirichlet boundary condition) static inline real64 computeBoundaryPressure( FaceManager::NodeMapType const & faceTonode, array2d< real64, nodes::REFERENCE_POSITION_PERM > const & node, @@ -1507,6 +1513,10 @@ static double computeLinearPressure_error( int ipKind, { makeDistortedPlanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); } + else if( mode == DistortionMode::PlanarFunnel ) + { + makeDistortedPlanarFunnel( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } else if( mode == DistortionMode::NonPlanar ) { makeDistortedNonplanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); @@ -1520,8 +1530,9 @@ static double computeLinearPressure_error( int ipKind, localIndex const fR_int = 2; real64 alpha_lin = 0.3; // pressure = x * alpha_lin * y - auto isDirL = [&]( localIndex f ){ return isBoundaryFace( faceTonode_L, node_L, f ); }; - auto isDirR = [&]( localIndex f ){ return isBoundaryFace( faceTonode_R, node_R, f ); }; + // in this two-cell mesh, every face except the shared one is a Dirichlet boundary face + auto isDirL = [&]( localIndex f ){ return f != fL_int; }; + auto isDirR = [&]( localIndex f ){ return f != fR_int; }; // compute the transmissibility matrix T constexpr real64 ltol = 1e-12; @@ -1726,6 +1737,23 @@ TEST( MimeticIP_Linear, Distortion_Planar_LinearPressure ) } } +// =================== case 1b: with distortion (funnel) =========================== +TEST( MimeticIP_Linear, Distortion_PlanarFunnel_LinearPressure ) +{ + for( double eps : { 0.2, 0.4 } ) + { + EXPECT_LT( computeLinearPressure_error< 6 >( InnerProductType::QUASI_TPFA, DistortionMode::PlanarFunnel, eps ), + consistency_tol ) << "eps = " << eps; + EXPECT_LT( computeLinearPressure_error< 6 >( InnerProductType::SIMPLE, DistortionMode::PlanarFunnel, eps ), + consistency_tol ) << "eps = " << eps; + EXPECT_LT( computeLinearPressure_error< 6 >( InnerProductType::BDVLM, DistortionMode::PlanarFunnel, eps ), + consistency_tol ) << "eps = " << eps; + + EXPECT_GT( computeLinearPressure_error< 6 >( InnerProductType::TPFA, DistortionMode::PlanarFunnel, eps ), + consistency_tol ) << "eps = " << eps; + } +} + // =================== case 2: with distortion (nonplanar) =========================== TEST( MimeticIP_Linear, Distortion_NonPlanar_LinearPressure ) { @@ -1849,6 +1877,10 @@ static double computeLinearPressureMixed_error( int ipKind, { makeDistortedPlanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); } + else if( mode == DistortionMode::PlanarFunnel ) + { + makeDistortedPlanarFunnel( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } else if( mode == DistortionMode::NonPlanar ) { makeDistortedNonplanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); @@ -1998,6 +2030,23 @@ TEST( MimeticIP_MixedLinear, Distortion_Planar_LinearPressure ) } } +// =================== mixed form: with distortion (funnel) =========================== +TEST( MimeticIP_MixedLinear, Distortion_PlanarFunnel_LinearPressure ) +{ + for( double eps : { 0.2, 0.4 } ) + { + EXPECT_LT( computeLinearPressureMixed_error< 6 >( InnerProductType::QUASI_TPFA, DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + EXPECT_LT( computeLinearPressureMixed_error< 6 >( InnerProductType::SIMPLE, DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + EXPECT_LT( computeLinearPressureMixed_error< 6 >( InnerProductType::BDVLM, DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + + EXPECT_GT( computeLinearPressureMixed_error< 6 >( InnerProductType::TPFA, DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + } +} + // =================== mixed form, case 2: with distortion (nonplanar) =========================== TEST( MimeticIP_MixedLinear, Distortion_NonPlanar_LinearPressure ) { @@ -2050,6 +2099,10 @@ static double computeMixedIPConsistency_error( int ipKind, { makeDistortedPlanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); } + else if( mode == DistortionMode::PlanarFunnel ) + { + makeDistortedPlanarFunnel( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } else if( mode == DistortionMode::NonPlanar ) { makeDistortedNonplanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); @@ -2134,6 +2187,114 @@ TEST( MimeticIP_MixedConsistency, Distortion_Planar ) } } +TEST( MimeticIP_MixedConsistency, Distortion_PlanarFunnel ) +{ + for( double eps : { 0.2, 0.4 } ) + { + EXPECT_LT( computeMixedIPConsistency_error< 6 >( InnerProductType::QUASI_TPFA, DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + EXPECT_LT( computeMixedIPConsistency_error< 6 >( InnerProductType::SIMPLE, DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + EXPECT_LT( computeMixedIPConsistency_error< 6 >( InnerProductType::BDVLM, DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + + EXPECT_GT( computeMixedIPConsistency_error< 6 >( InnerProductType::TPFA, DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + } +} + +//======================== Geometric Identity Test ============================= +// C^T N = vol * I, with N_f = |f| n_f and C_f = x_f - x_c (any reference point x_c); +// holds iff faces are planar and x_f is the area centroid. This identity is what the +// consistency term of every mimetic inner product relies on. + +template< int NF > +static double computeGeometricIdentity_error( DistortionMode mode = DistortionMode::None, + real64 eps = 0.0 ) +{ + array2d< real64, nodes::REFERENCE_POSITION_PERM > node_L; + array2d< real64, nodes::REFERENCE_POSITION_PERM > node_R; + FaceManager::NodeMapType faceTonode_L; + FaceManager::NodeMapType faceTonode_R; + array1d< localIndex > elemToface_L; + array1d< localIndex > elemToface_R; + real64 center_L[3], center_R[3], vol_L = 0, vol_R = 0; + + makeUnitCube( 0.0, node_L, faceTonode_L, elemToface_L, center_L, vol_L ); + makeUnitCube( 1.0, node_R, faceTonode_R, elemToface_R, center_R, vol_R ); + + if( mode == DistortionMode::Planar ) + { + makeDistortedPlanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } + else if( mode == DistortionMode::PlanarFunnel ) + { + makeDistortedPlanarFunnel( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } + else if( mode == DistortionMode::NonPlanar ) + { + makeDistortedNonplanar( node_L, node_R, faceTonode_L, faceTonode_R, eps ); + } + + computeDistortedVolumeAndCenter( node_L, center_L, vol_L ); + + real64 E[3][3] = {{ 0 }}; + for( int d = 0; d < 3; ++d ) + { + E[d][d] = -vol_L; + } + + for( int f = 0; f < NF; ++f ) + { + real64 fc[3], fn[3]; + real64 const area = centroid_3DPolygon( faceTonode_L[f], node_L.toViewConst(), fc, fn ); + real64 c2f[3] = { fc[0] - center_L[0], fc[1] - center_L[1], fc[2] - center_L[2] }; + if( LvArray::tensorOps::AiBi< 3 >( c2f, fn ) < 0.0 ) + { + LvArray::tensorOps::scale< 3 >( fn, -1.0 ); + } + for( int i = 0; i < 3; ++i ) + { + for( int j = 0; j < 3; ++j ) + { + E[i][j] += c2f[i] * area * fn[j]; + } + } + } + + double err = 0.0; + for( int i = 0; i < 3; ++i ) + { + for( int j = 0; j < 3; ++j ) + { + err = std::max( err, std::abs( E[i][j] ) / vol_L ); + } + } + return err; +} + +TEST( MimeticIP_GeometricIdentity, PlanarFaceCells ) +{ + EXPECT_LT( computeGeometricIdentity_error< 6 >(), mixed_consistency_tol ); + for( double eps : { 0.2, 0.9 } ) + { + EXPECT_LT( computeGeometricIdentity_error< 6 >( DistortionMode::Planar, eps ), + mixed_consistency_tol ) << "eps = " << eps; + } + for( double eps : { 0.2, 0.4 } ) + { + EXPECT_LT( computeGeometricIdentity_error< 6 >( DistortionMode::PlanarFunnel, eps ), + mixed_consistency_tol ) << "eps = " << eps; + } +} + +TEST( MimeticIP_GeometricIdentity, NonPlanarCells ) +{ + // non-planar faces break the identity + EXPECT_GT( computeGeometricIdentity_error< 6 >( DistortionMode::NonPlanar, 0.2 ), + mixed_consistency_tol ); +} + //======================== Hybrid/Mixed Duality Test ============================= // M = T^{-1} must hold; exact on parallelepipeds with isotropic K at any cell size h. diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index b42090eb5e6..4cd4cffc448 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -271,6 +271,24 @@ real64 centroid_3DPolygon( arraySlice1d< localIndex const > const pointsIndices, { LvArray::tensorOps::normalize< 3 >( normal ); area *= 0.5; + + // replace the vertex average by the area centroid (exact for planar polygons) + real64 areaCentroid[ 3 ] = { 0.0 }; + for( localIndex a=0; a( current, points[ pointsIndices[ a++ ]] ); + LvArray::tensorOps::scaledAdd< 3 >( current, origin, -1. ); + LvArray::tensorOps::copy< 3 >( next, points[ pointsIndices[ a % numberOfPoints ] ] ); + LvArray::tensorOps::scaledAdd< 3 >( next, origin, -1. ); + + LvArray::tensorOps::crossProduct( crossProduct, current, next ); + real64 const signedTriArea = 0.5 * LvArray::tensorOps::AiBi< 3 >( crossProduct, normal ); + + LvArray::tensorOps::scaledAdd< 3 >( areaCentroid, current, signedTriArea / 3.0 ); + LvArray::tensorOps::scaledAdd< 3 >( areaCentroid, next, signedTriArea / 3.0 ); + } + LvArray::tensorOps::copy< 3 >( center, origin ); + LvArray::tensorOps::scaledAdd< 3 >( center, areaCentroid, 1.0 / area ); } else if( area < -areaTolerance ) { From b5a7c8fffbc5e7a89135489b0e69aab3138f2e6c Mon Sep 17 00:00:00 2001 From: Omar Duran Date: Mon, 3 Aug 2026 18:44:18 -0700 Subject: [PATCH 10/11] wip --- .../unitTests/testMimeticInnerProducts.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp index abcf66f7d03..c9df233b9f9 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp @@ -1403,10 +1403,10 @@ static inline void makeDistortedPlanar( array2d< real64, nodes::REFERENCE_POSITI // shrink the shared face toward its center: two funnel-shaped cells, all faces planar static inline void makeDistortedPlanarFunnel( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeL, - array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeR, - FaceManager::NodeMapType const & faceL, - FaceManager::NodeMapType const & faceR, - real64 eps ) + array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeR, + FaceManager::NodeMapType const & faceL, + FaceManager::NodeMapType const & faceR, + real64 eps ) { localIndex const fL = 3, fR = 2; for( int a = 0; a < 4; ++a ) From 82fb8b49d8c8493571a0509eab0261f88bfe509c Mon Sep 17 00:00:00 2001 From: Omar Duran Date: Tue, 4 Aug 2026 11:48:31 -0700 Subject: [PATCH 11/11] test: BdVLM inner product reduce to TPFA on hexahedral cells --- .../unitTests/testMimeticInnerProducts.cpp | 112 ++++++++++++++++++ 1 file changed, 112 insertions(+) diff --git a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp index c9df233b9f9..b0876394994 100644 --- a/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp +++ b/src/coreComponents/finiteVolume/mimeticInnerProducts/unitTests/testMimeticInnerProducts.cpp @@ -2400,6 +2400,118 @@ TEST( MimeticIP_MixedDuality, BdVLM ) EXPECT_LT( computeHybridMixedDuality_error< 6 >( InnerProductType::BDVLM, 0.2 ), duality_tol ); } +//======================== TPFA Reduction Test ============================= +// On K-orthogonal grids: T = diag( 2 k_d ), M = T^{-1} on the unit cube. +// BdVLM reduces iff gamma = 2/NF * trace (isotropic K only); Eq. (3.67) gamma = +// trace/NF and SIMPLE do not reduce. This test fixes the adopted convention. + +template< int NF > +static double computeTpfaReduction_error( int ipKind, bool mixedForm, + real64 kx = 1.0, real64 ky = 1.0, real64 kz = 1.0 ) +{ + array2d< real64, nodes::REFERENCE_POSITION_PERM > node; + FaceManager::NodeMapType faceTonode; + array1d< localIndex > elemToface; + real64 center[3], vol = 0; + + makeUnitCube( 0.0, node, faceTonode, elemToface, center, vol ); + computeDistortedVolumeAndCenter( node, center, vol ); + + constexpr real64 ltol = 1e-12; + real64 Kvec[3] = {kx, ky, kz}; + + // makeUnitCube face order: z=0, y=0, x=0, x=1, y=1, z=1 + int const dirOfFace[6] = {2, 1, 0, 0, 1, 2}; + + stackArray1d< real64, 3 > cc( 3 ); + for( int d = 0; d < 3; ++d ) + { + cc[d] = center[d]; + } + + array1d< real64 > mult( NF ); + mult.setValues< parallelHostPolicy >( 1.0 ); + + stackArray2d< real64, NF * NF > A( NF, NF ); + A.template setValues< parallelHostPolicy >( 0.0 ); + + bool refIsM = false; + if( mixedForm ) + { + computeM_dispatch< NF >( ipKind, node, faceTonode, elemToface, cc.toSliceConst(), vol, Kvec, ltol, A.toSlice() ); + refIsM = true; // M_tpfa(f) = 1 / (2 k_dir) + } + else + { + if( ipKind == InnerProductType::TPFA ) + { + TPFAInnerProduct::compute< NF >( node.toViewConst(), mult.toViewConst(), faceTonode.toViewConst(), + elemToface.toSliceConst(), cc, vol, Kvec, ltol, A.toSlice() ); + } + else if( ipKind == InnerProductType::QUASI_TPFA ) + { + QuasiTPFAInnerProduct::compute< NF >( node.toViewConst(), mult.toViewConst(), faceTonode.toViewConst(), + elemToface.toSliceConst(), cc, vol, Kvec, ltol, A.toSlice() ); + } + else if( ipKind == InnerProductType::SIMPLE ) + { + SimpleInnerProduct::compute< NF >( node.toViewConst(), mult.toViewConst(), faceTonode.toViewConst(), + elemToface.toSliceConst(), cc, vol, Kvec, ltol, A.toSlice() ); + } + else if( ipKind == InnerProductType::BDVLM ) + { + BdVLMInnerProduct::compute< NF >( node.toViewConst(), mult.toViewConst(), faceTonode.toViewConst(), + elemToface.toSliceConst(), cc, vol, Kvec, ltol, A.toSlice() ); + } + } + + double err = 0.0; + for( int i = 0; i < NF; ++i ) + { + // T_tpfa(f) = k_dir A / d = 2 k_dir per face + real64 const refDiag = refIsM ? 1.0 / ( 2.0 * Kvec[ dirOfFace[i] ] ) : 2.0 * Kvec[ dirOfFace[i] ]; + for( int j = 0; j < NF; ++j ) + { + err = std::max( err, std::abs( A( i, j ) - ( i == j ? refDiag : 0.0 ) ) ); + } + } + return err; +} + +TEST( MimeticIP_TpfaReduction, Hybrid ) +{ + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::TPFA, false ), mixed_consistency_tol ); + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::QUASI_TPFA, false ), mixed_consistency_tol ); + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::BDVLM, false ), mixed_consistency_tol ); + EXPECT_GT( computeTpfaReduction_error< 6 >( InnerProductType::SIMPLE, false ), mixed_consistency_tol ); +} + +TEST( MimeticIP_TpfaReduction, Mixed ) +{ + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::TPFA, true ), mixed_consistency_tol ); + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::QUASI_TPFA, true ), mixed_consistency_tol ); + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::BDVLM, true ), mixed_consistency_tol ); + EXPECT_GT( computeTpfaReduction_error< 6 >( InnerProductType::SIMPLE, true ), mixed_consistency_tol ); +} + +// with anisotropic (diagonal) K the grid is still K-orthogonal: TPFA and quasi-TPFA +// reduce per direction, but BdVLM's scalar trace-based stabilization cannot +TEST( MimeticIP_TpfaReduction, Hybrid_AnisotropicK ) +{ + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::TPFA, false, 4.0, 1.0, 0.5 ), mixed_consistency_tol ); + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::QUASI_TPFA, false, 4.0, 1.0, 0.5 ), mixed_consistency_tol ); + EXPECT_GT( computeTpfaReduction_error< 6 >( InnerProductType::BDVLM, false, 4.0, 1.0, 0.5 ), mixed_consistency_tol ); + EXPECT_GT( computeTpfaReduction_error< 6 >( InnerProductType::SIMPLE, false, 4.0, 1.0, 0.5 ), mixed_consistency_tol ); +} + +TEST( MimeticIP_TpfaReduction, Mixed_AnisotropicK ) +{ + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::TPFA, true, 4.0, 1.0, 0.5 ), mixed_consistency_tol ); + EXPECT_LT( computeTpfaReduction_error< 6 >( InnerProductType::QUASI_TPFA, true, 4.0, 1.0, 0.5 ), mixed_consistency_tol ); + EXPECT_GT( computeTpfaReduction_error< 6 >( InnerProductType::BDVLM, true, 4.0, 1.0, 0.5 ), mixed_consistency_tol ); + EXPECT_GT( computeTpfaReduction_error< 6 >( InnerProductType::SIMPLE, true, 4.0, 1.0, 0.5 ), mixed_consistency_tol ); +} + //======================== Hydrostatic Equilibrium Consistency Test ============================= static inline void makeDistortedPlanar_gravity( array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeL, array2d< real64, nodes::REFERENCE_POSITION_PERM > & nodeU,