From a599f70d32e7730e5829985a49c9909d126f613d Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 11 Nov 2024 13:18:58 -0600 Subject: [PATCH 01/11] EDFM bugfixes, part of #3348 --- .../constitutive/contact/CoulombFriction.hpp | 6 ++++-- .../SinglePhasePoromechanicsEFEM.hpp | 4 +++- .../SinglePhasePoromechanicsEFEM_impl.hpp | 13 +++++++------ .../ThermalSinglePhasePoromechanicsEFEM.hpp | 2 +- .../ThermalSinglePhasePoromechanicsEFEM_impl.hpp | 4 ++-- 5 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp index d117d8895e5..276bf82fe07 100644 --- a/src/coreComponents/constitutive/contact/CoulombFriction.hpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.hpp @@ -265,6 +265,8 @@ inline void CoulombFrictionUpdates::computeShearTraction( localIndex const k, real64 dLimitTau_dNormalTraction; real64 const limitTau = computeLimitTangentialTractionNorm( tractionVector[0], dLimitTau_dNormalTraction ); + // dLimitTau_dNormalTraction from the function above has a wrong sign, flip it + dLimitTau_dNormalTraction = -dLimitTau_dNormalTraction; real64 const slipNorm = LvArray::tensorOps::l2Norm< 2 >( slip ); @@ -274,10 +276,10 @@ inline void CoulombFrictionUpdates::computeShearTraction( localIndex const k, dTractionVector_dJump[1][0] = dTractionVector_dJump[0][0] * dLimitTau_dNormalTraction * slip[0] / slipNorm; dTractionVector_dJump[1][1] = limitTau * pow( slip[1], 2 ) / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); - dTractionVector_dJump[1][2] = limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); + dTractionVector_dJump[1][2] = -limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); dTractionVector_dJump[2][0] = dTractionVector_dJump[0][0] * dLimitTau_dNormalTraction * slip[1] / slipNorm; - dTractionVector_dJump[2][1] = limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); + dTractionVector_dJump[2][1] = -limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); dTractionVector_dJump[2][2] = limitTau * pow( slip[0], 2 ) / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); // Compute elastic component of the slip for this case diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp index ad5ccf38ab1..e88ce14509f 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp @@ -278,7 +278,9 @@ class SinglePhasePoromechanicsEFEM : arrayView1d< real64 const > const m_surfaceArea; - arrayView1d< real64 const > const m_elementVolume; + arrayView1d< real64 const > const m_elementVolumeCell; + + arrayView1d< real64 const > const m_elementVolumeFrac; arrayView1d< real64 const > const m_deltaVolume; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp index 0aac71f9acb..2a414c29d84 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp @@ -91,8 +91,9 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, m_tVec2( embeddedSurfSubRegion.getTangentVector2() ), m_surfaceCenter( embeddedSurfSubRegion.getElementCenter() ), m_surfaceArea( embeddedSurfSubRegion.getElementArea() ), - m_elementVolume( elementSubRegion.getElementVolume() ), - m_deltaVolume( elementSubRegion.template getField< fields::flow::deltaVolume >() ), + m_elementVolumeCell( elementSubRegion.getElementVolume() ), + m_elementVolumeFrac( embeddedSurfSubRegion.getElementVolume() ), + m_deltaVolume( embeddedSurfSubRegion.template getField< fields::flow::deltaVolume >() ), m_fracturedElems( elementSubRegion.fracturedElementsList() ), m_cellsToEmbeddedSurfaces( elementSubRegion.embeddedSurfacesList().toViewConst() ), m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] }, @@ -148,7 +149,7 @@ setup( localIndex const k, { localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0]; - stack.hInv = m_surfaceArea[embSurfIndex] / m_elementVolume[k]; + stack.hInv = m_surfaceArea[embSurfIndex] / m_elementVolumeCell[k]; for( localIndex a=0; a Date: Mon, 11 Nov 2024 18:38:00 -0600 Subject: [PATCH 02/11] fix sign directly in computeLimitTangentialTractionNorm --- src/coreComponents/constitutive/contact/CoulombFriction.hpp | 4 +--- .../physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp | 4 ++-- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp index 276bf82fe07..ae1147ee4a4 100644 --- a/src/coreComponents/constitutive/contact/CoulombFriction.hpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.hpp @@ -221,7 +221,7 @@ GEOS_HOST_DEVICE real64 CoulombFrictionUpdates::computeLimitTangentialTractionNorm( real64 const & normalTraction, real64 & dLimitTangentialTractionNorm_dTraction ) const { - dLimitTangentialTractionNorm_dTraction = m_frictionCoefficient; + dLimitTangentialTractionNorm_dTraction = - m_frictionCoefficient; return ( m_cohesion - normalTraction * m_frictionCoefficient ); } @@ -265,8 +265,6 @@ inline void CoulombFrictionUpdates::computeShearTraction( localIndex const k, real64 dLimitTau_dNormalTraction; real64 const limitTau = computeLimitTangentialTractionNorm( tractionVector[0], dLimitTau_dNormalTraction ); - // dLimitTau_dNormalTraction from the function above has a wrong sign, flip it - dLimitTau_dNormalTraction = -dLimitTau_dNormalTraction; real64 const slipNorm = LvArray::tensorOps::l2Norm< 2 >( slip ); diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp index 08d5fec172e..3271776a9db 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp @@ -1594,7 +1594,7 @@ void SolidMechanicsLagrangeContact:: { elemRHS[i] = Ja * ( traction[kfe][i] - limitTau * sliding[ i-1 ] / slidingNorm ); - dRdT( i, 0 ) = Ja * dLimitTau_dNormalTraction * sliding[ i-1 ] / slidingNorm; + dRdT( i, 0 ) = - Ja * dLimitTau_dNormalTraction * sliding[ i-1 ] / slidingNorm; dRdT( i, i ) = Ja; } @@ -1642,7 +1642,7 @@ void SolidMechanicsLagrangeContact:: { elemRHS[i] = Ja * traction[kfe][i] * ( 1.0 - limitTau / vauxNorm ); - dRdT( i, 0 ) = Ja * traction[kfe][i] * dLimitTau_dNormalTraction / vauxNorm; + dRdT( i, 0 ) = - Ja * traction[kfe][i] * dLimitTau_dNormalTraction / vauxNorm; dRdT( i, i ) = Ja; } } From 922951b2c3cdc2d27575163dbdafa407f926f83b Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 12 Nov 2024 09:32:19 -0600 Subject: [PATCH 03/11] code style --- src/coreComponents/constitutive/contact/CoulombFriction.hpp | 2 +- .../physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp index ae1147ee4a4..abe165ae4bf 100644 --- a/src/coreComponents/constitutive/contact/CoulombFriction.hpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.hpp @@ -221,7 +221,7 @@ GEOS_HOST_DEVICE real64 CoulombFrictionUpdates::computeLimitTangentialTractionNorm( real64 const & normalTraction, real64 & dLimitTangentialTractionNorm_dTraction ) const { - dLimitTangentialTractionNorm_dTraction = - m_frictionCoefficient; + dLimitTangentialTractionNorm_dTraction = -m_frictionCoefficient; return ( m_cohesion - normalTraction * m_frictionCoefficient ); } diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp index 3271776a9db..c418110f491 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp @@ -1594,7 +1594,7 @@ void SolidMechanicsLagrangeContact:: { elemRHS[i] = Ja * ( traction[kfe][i] - limitTau * sliding[ i-1 ] / slidingNorm ); - dRdT( i, 0 ) = - Ja * dLimitTau_dNormalTraction * sliding[ i-1 ] / slidingNorm; + dRdT( i, 0 ) = -Ja * dLimitTau_dNormalTraction * sliding[ i-1 ] / slidingNorm; dRdT( i, i ) = Ja; } @@ -1642,7 +1642,7 @@ void SolidMechanicsLagrangeContact:: { elemRHS[i] = Ja * traction[kfe][i] * ( 1.0 - limitTau / vauxNorm ); - dRdT( i, 0 ) = - Ja * traction[kfe][i] * dLimitTau_dNormalTraction / vauxNorm; + dRdT( i, 0 ) = -Ja * traction[kfe][i] * dLimitTau_dNormalTraction / vauxNorm; dRdT( i, i ) = Ja; } } From 78a90da439ae2d6f85473df1acf1a6605851800b Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 14 Nov 2024 17:22:22 -0600 Subject: [PATCH 04/11] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 00502a919e4..0c5343a6f52 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3339-8707-7c55c70 + baseline: integratedTests/baseline_integratedTests-pr3439-8742-0bffafb allow_fail: all: '' streak: '' From 4cfc999eb9182ca0c4c72e42b7412a19e93160eb Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 14 Nov 2024 17:23:13 -0600 Subject: [PATCH 05/11] Update BASELINE_NOTES.md --- BASELINE_NOTES.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index f4dcac6a76e..30b91b77430 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,9 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3439 (2024-11-14) +===================== +EDFM bugfixes: derivatives sign, frac/cell element volume. PR #3339 (2024-11-14) ===================== From 0b3cad54d13294aea4e8490bfbe8e8c220d7a924 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 15 Nov 2024 14:43:30 -0600 Subject: [PATCH 06/11] make that case running --- .../ExponentialDecayPermeability_edfm_base.xml | 4 +++- .../constitutive/contact/FrictionBase.hpp | 6 +++++- .../constitutive/contact/FrictionlessContact.hpp | 14 -------------- .../physicsSolvers/fluidFlow/SinglePhaseBase.cpp | 12 ++++++------ 4 files changed, 14 insertions(+), 22 deletions(-) mode change 100755 => 100644 inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml diff --git a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml old mode 100755 new mode 100644 index 64b958cd590..a9f461c4e41 --- a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml +++ b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml @@ -13,7 +13,9 @@ + maxTimeStepCuts="1" + maxAllowedResidualNorm="1e20" + lineSearchAction="None"/> diff --git a/src/coreComponents/constitutive/contact/FrictionBase.hpp b/src/coreComponents/constitutive/contact/FrictionBase.hpp index e05923f2998..84e5aac824d 100644 --- a/src/coreComponents/constitutive/contact/FrictionBase.hpp +++ b/src/coreComponents/constitutive/contact/FrictionBase.hpp @@ -175,7 +175,11 @@ class FrictionBaseUpdates inline virtual real64 computeLimitTangentialTractionNorm( real64 const & normalTraction, real64 & dLimitTangentialTractionNorm_dTraction ) const - { GEOS_UNUSED_VAR( normalTraction, dLimitTangentialTractionNorm_dTraction ); return 0; }; + { + GEOS_UNUSED_VAR( normalTraction ); + dLimitTangentialTractionNorm_dTraction = 0.0; + return 0; + } protected: diff --git a/src/coreComponents/constitutive/contact/FrictionlessContact.hpp b/src/coreComponents/constitutive/contact/FrictionlessContact.hpp index 83e53fd8e14..6b3552acf75 100644 --- a/src/coreComponents/constitutive/contact/FrictionlessContact.hpp +++ b/src/coreComponents/constitutive/contact/FrictionlessContact.hpp @@ -63,20 +63,6 @@ class FrictionlessContactUpdates : public FrictionBaseUpdates arraySlice1d< real64 const > const & tractionVector, integer & fractureState ) const override final; - - /** - * @brief Evaluate the limit tangential traction norm and return the derivative wrt normal traction - * @param[in] normalTraction the normal traction - * @param[out] dLimitTangentialTractionNorm_dTraction the derivative of the limit tangential traction norm wrt normal traction - * @return the limit tangential traction norm - */ - GEOS_HOST_DEVICE - inline - virtual real64 computeLimitTangentialTractionNorm( real64 const & normalTraction, - real64 & dLimitTangentialTractionNorm_dTraction ) const override final - { GEOS_UNUSED_VAR( normalTraction, dLimitTangentialTractionNorm_dTraction ); return 0.0; } - -private: }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 8261e822868..2a9f36b56db 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -471,7 +471,7 @@ void SinglePhaseBase::initializePostInitialConditionsPreSubGroups() [&]( localIndex const, SurfaceElementRegion & region ) { - region.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) + region.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) { subRegion.getWrapper< real64_array >( fields::flow::hydraulicAperture::key() ). setApplyDefaultValue( region.getDefaultAperture() ); @@ -707,8 +707,8 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_ } ); - mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, - FaceElementSubRegion & subRegion ) + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) { arrayView1d< real64 const > const aper = subRegion.getField< fields::flow::hydraulicAperture >(); arrayView1d< real64 > const aper0 = subRegion.getField< fields::flow::aperture0 >(); @@ -728,6 +728,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_ fluid.saveConvergedState(); } ); + } ); } @@ -791,8 +792,8 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, } ); - mesh.getElemManager().forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, - FaceElementSubRegion & subRegion ) + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) { arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank(); arrayView1d< real64 const > const volume = subRegion.getElementVolume(); @@ -817,7 +818,6 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, } } ); } ); - } ); } From 701a258a06fc267710659bd28d9ca9954c627806 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 15 Nov 2024 17:42:44 -0600 Subject: [PATCH 07/11] code style --- src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 2a9f36b56db..429ea336ab7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -793,7 +793,7 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, } ); mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, [&]( localIndex const, - SurfaceElementSubRegion & subRegion ) + SurfaceElementSubRegion & subRegion ) { arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank(); arrayView1d< real64 const > const volume = subRegion.getElementVolume(); From 63aa00ccaddc4a7ce53b0864e41ed7f1a1d7bcf5 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 18 Nov 2024 17:23:27 -0600 Subject: [PATCH 08/11] delta volume fix plus some other stuff --- .../fluidFlow/SinglePhaseBase.cpp | 23 +++--- .../fluidFlow/SinglePhaseBase.hpp | 2 + ...glePhasePoromechanicsEmbeddedFractures.cpp | 5 +- .../SinglePhasePoromechanicsEFEM.hpp | 70 +------------------ .../SinglePhasePoromechanicsEFEM_impl.hpp | 7 +- .../ThermalSinglePhasePoromechanicsEFEM.hpp | 7 +- ...ermalSinglePhasePoromechanicsEFEM_impl.hpp | 8 +-- 7 files changed, 28 insertions(+), 94 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 429ea336ab7..5bf3f88fdfd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -691,8 +691,7 @@ void SinglePhaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM( time_ { saveConvergedState( subRegion ); - arrayView1d< real64 > const & dVol = subRegion.template getField< fields::flow::deltaVolume >(); - dVol.zero(); + applyDeltaVolume( subRegion ); // This should fix NaN density in newly created fracture elements updatePorosityAndPermeability( subRegion ); @@ -756,14 +755,7 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time, singlePhaseBaseKernels::StatisticsKernel:: saveDeltaPressure( subRegion.size(), pres, initPres, deltaPres ); - arrayView1d< real64 > const dVol = subRegion.getField< fields::flow::deltaVolume >(); - arrayView1d< real64 > const vol = subRegion.getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString() ); - - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) - { - vol[ei] += dVol[ei]; - dVol[ei] = 0.0; - } ); + applyDeltaVolume( subRegion ); SingleFluidBase const & fluid = getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) ); @@ -1403,4 +1395,15 @@ void SinglePhaseBase::saveConvergedState( ElementSubRegionBase & subRegion ) con mass_n.setValues< parallelDevicePolicy<> >( mass ); } +void SinglePhaseBase::applyDeltaVolume( ElementSubRegionBase & subRegion ) +{ + arrayView1d< real64 > const dVol = subRegion.template getField< fields::flow::deltaVolume >(); + arrayView1d< real64 > const vol = subRegion.template getReference< array1d< real64 > >( CellElementSubRegion::viewKeyStruct::elementVolumeString()); + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + vol[ei] += dVol[ei]; + dVol[ei] = 0.0; + } ); +} + } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp index bf1cd19eadd..fc0d9de9b91 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp @@ -387,6 +387,8 @@ class SinglePhaseBase : public FlowSolverBase */ virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override; + void applyDeltaVolume( ElementSubRegionBase & subRegion ); + /** * @brief Structure holding views into fluid properties used by the base solver. */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index 394b95d08d3..ca1340ffcb7 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -26,6 +26,7 @@ #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp" +#include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" #include "finiteVolume/FluxApproximationBase.hpp" @@ -510,10 +511,10 @@ void SinglePhasePoromechanicsEmbeddedFractures::updateState( DomainPartition & d using HydraulicApertureModelType = TYPEOFREF( castedHydraulicApertureModel ); typename HydraulicApertureModelType::KernelWrapper hydraulicApertureModelWrapper = castedHydraulicApertureModel.createKernelWrapper(); - poromechanicsEFEMKernels::StateUpdateKernel:: + poromechanicsFracturesKernels::StateUpdateKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), - hydraulicApertureModelWrapper, porousMaterialWrapper, + hydraulicApertureModelWrapper, dispJump, pressure, area, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp index e88ce14509f..387cf5c2090 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp @@ -250,7 +250,6 @@ class SinglePhasePoromechanicsEFEM : /// The rank global densities arrayView2d< real64 const > const m_solidDensity; arrayView2d< real64 const > const m_fluidDensity; - arrayView2d< real64 const > const m_fluidDensity_n; arrayView2d< real64 const > const m_dFluidDensity_dPressure; /// The rank-global fluid pressure array. @@ -259,9 +258,6 @@ class SinglePhasePoromechanicsEFEM : /// The rank-global fluid pressure array. arrayView1d< real64 const > const m_fracturePressure; - /// The rank-global delta-fluid pressure array. - arrayView2d< real64 const > const m_porosity_n; - arrayView2d< real64 const > const m_tractionVec; arrayView3d< real64 const > const m_dTraction_dJump; @@ -284,6 +280,8 @@ class SinglePhasePoromechanicsEFEM : arrayView1d< real64 const > const m_deltaVolume; + arrayView1d< real64 const > const m_mass_n; + SortedArrayView< localIndex const > const m_fracturedElems; ArrayOfArraysView< localIndex const > const m_cellsToEmbeddedSurfaces; @@ -307,70 +305,6 @@ using SinglePhaseKernelFactory = finiteElement::KernelFactory< SinglePhasePorome real64 const (&)[3], string const >; -/** - * @brief A struct to perform volume, aperture and fracture traction updates - */ -struct StateUpdateKernel -{ - - /** - * @brief Launch the kernel function doing volume, aperture and fracture traction updates - * @tparam POLICY the type of policy used in the kernel launch - * @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates - * @param[in] size the size of the subregion - * @param[in] contactWrapper the wrapper implementing the contact relationship - * @param[in] dispJump the displacement jump - * @param[in] pressure the pressure - * @param[in] area the area - * @param[in] volume the volume - * @param[out] deltaVolume the change in volume - * @param[out] aperture the aperture - * @param[out] hydraulicAperture the effecture aperture - * @param[out] fractureContactTraction the fracture contact traction - */ - template< typename POLICY, typename POROUS_WRAPPER, typename CONTACT_WRAPPER > - static void - launch( localIndex const size, - CONTACT_WRAPPER const & contactWrapper, - POROUS_WRAPPER const & porousMaterialWrapper, - arrayView2d< real64 const > const & dispJump, - arrayView1d< real64 const > const & pressure, - arrayView1d< real64 const > const & area, - arrayView1d< real64 const > const & volume, - arrayView1d< real64 > const & deltaVolume, - arrayView1d< real64 > const & aperture, - arrayView1d< real64 const > const & oldHydraulicAperture, - arrayView1d< real64 > const & hydraulicAperture, - arrayView2d< real64 > const & fractureEffectiveTraction ) - { - forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - // update aperture to be equal to the normal displacement jump - aperture[k] = dispJump[k][0]; // the first component of the jump is the normal one. - - real64 dHydraulicAperture_dNormalJump = 0.0; - real64 dHydraulicAperture_dNormalTraction = 0.0; - hydraulicAperture[k] = contactWrapper.computeHydraulicAperture( aperture[k], - fractureEffectiveTraction[k][0], - dHydraulicAperture_dNormalJump, - dHydraulicAperture_dNormalTraction ); - - deltaVolume[k] = hydraulicAperture[k] * area[k] - volume[k]; - - real64 const jump[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( dispJump[k] ); - real64 const effectiveTraction[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( fractureEffectiveTraction[k] ); - - // all perm update models below should need effective traction instead of total traction - // (total traction is combined forces of fluid pressure and effective traction) - porousMaterialWrapper.updateStateFromPressureApertureJumpAndTraction( k, 0, pressure[k], - oldHydraulicAperture[k], hydraulicAperture[k], - dHydraulicAperture_dNormalJump, - jump, effectiveTraction ); - - } ); - } -}; - } // namespace poromechanicsEFEMKernels } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp index 2a414c29d84..774f09d44bf 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM_impl.hpp @@ -77,12 +77,10 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, m_wDofNumber( jumpDofNumber ), m_solidDensity( inputConstitutiveType.getDensity() ), m_fluidDensity( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ), - m_fluidDensity_n( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density_n() ), m_dFluidDensity_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dPressure() ), m_matrixPressure( elementSubRegion.template getField< fields::flow::pressure >() ), m_fracturePressure( embeddedSurfSubRegion.template getField< fields::flow::pressure >() ), - m_porosity_n( inputConstitutiveType.getPorosity_n() ), m_tractionVec( embeddedSurfSubRegion.getField< fields::contact::traction >() ), m_dTraction_dJump( embeddedSurfSubRegion.getField< fields::contact::dTraction_dJump >() ), m_dTraction_dPressure( embeddedSurfSubRegion.getField< fields::contact::dTraction_dPressure >() ), @@ -94,6 +92,7 @@ SinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, m_elementVolumeCell( elementSubRegion.getElementVolume() ), m_elementVolumeFrac( embeddedSurfSubRegion.getElementVolume() ), m_deltaVolume( embeddedSurfSubRegion.template getField< fields::flow::deltaVolume >() ), + m_mass_n( embeddedSurfSubRegion.template getField< fields::flow::mass_n >() ), m_fracturedElems( elementSubRegion.fracturedElementsList() ), m_cellsToEmbeddedSurfaces( elementSubRegion.embeddedSurfacesList().toViewConst() ), m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] }, @@ -322,9 +321,7 @@ complete( localIndex const k, // Mass balance accumulation real64 const newVolume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex ); - real64 const newMass = m_fluidDensity( embSurfIndex, 0 ) * newVolume; - real64 const oldMass = m_fluidDensity_n( embSurfIndex, 0 ) * m_elementVolumeFrac( embSurfIndex ); - real64 const localFlowResidual = ( newMass - oldMass ); + real64 const localFlowResidual = m_fluidDensity( embSurfIndex, 0 ) * newVolume - m_mass_n[embSurfIndex]; real64 const localFlowJumpJacobian = m_fluidDensity( embSurfIndex, 0 ) * m_surfaceArea[ embSurfIndex ]; real64 const localFlowFlowJacobian = m_dFluidDensity_dPressure( embSurfIndex, 0 ) * newVolume; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp index b90d9006683..24ed5dc4496 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM.hpp @@ -57,9 +57,7 @@ class ThermalSinglePhasePoromechanicsEFEM : using Base::m_matrixPresDofNumber; using Base::m_wDofNumber; using Base::m_fluidDensity; - using Base::m_fluidDensity_n; using Base::m_dFluidDensity_dPressure; - using Base::m_porosity_n; using Base::m_surfaceArea; using Base::m_elementVolumeFrac; using Base::m_deltaVolume; @@ -164,17 +162,18 @@ class ThermalSinglePhasePoromechanicsEFEM : arrayView2d< real64 const > const m_dFluidDensity_dTemperature; /// Views on fluid internal energy - arrayView2d< real64 const > const m_fluidInternalEnergy_n; arrayView2d< real64 const > const m_fluidInternalEnergy; arrayView2d< real64 const > const m_dFluidInternalEnergy_dPressure; arrayView2d< real64 const > const m_dFluidInternalEnergy_dTemperature; /// Views on temperature - arrayView1d< real64 const > const m_temperature_n; arrayView1d< real64 const > const m_temperature; /// The rank-global fluid pressure array. arrayView1d< real64 const > const m_matrixTemperature; + + /// Views on energy + arrayView1d< real64 const > const m_energy_n; }; diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp index b838e72f7b2..fd1a33dce86 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp @@ -68,15 +68,14 @@ ThermalSinglePhasePoromechanicsEFEM( NodeManager const & nodeManager, fluidModelKey ), m_dFluidDensity_dTemperature( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dTemperature() ), - m_fluidInternalEnergy_n( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy_n() ), m_fluidInternalEnergy( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).internalEnergy() ), m_dFluidInternalEnergy_dPressure( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dInternalEnergy_dPressure() ), m_dFluidInternalEnergy_dTemperature( embeddedSurfSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dInternalEnergy_dTemperature() ), - m_temperature_n( embeddedSurfSubRegion.template getField< fields::flow::temperature_n >() ), m_temperature( embeddedSurfSubRegion.template getField< fields::flow::temperature >() ), - m_matrixTemperature( elementSubRegion.template getField< fields::flow::temperature >() ) + m_matrixTemperature( elementSubRegion.template getField< fields::flow::temperature >() ), + m_energy_n( embeddedSurfSubRegion.template getField< fields::flow::energy_n >() ) {} @@ -176,9 +175,8 @@ complete( localIndex const k, localIndex const embSurfIndex = m_cellsToEmbeddedSurfaces[k][0]; // Energy balance accumulation real64 const volume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex ); - real64 const volume_n = m_elementVolumeFrac( embSurfIndex ); real64 const fluidEnergy = m_fluidDensity( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) * volume; - real64 const fluidEnergy_n = m_fluidDensity_n( embSurfIndex, 0 ) * m_fluidInternalEnergy_n( embSurfIndex, 0 ) * volume_n; + real64 const fluidEnergy_n = m_energy_n[embSurfIndex]; // TODO where is solid energy? stack.dFluidMassIncrement_dTemperature = m_dFluidDensity_dTemperature( embSurfIndex, 0 ) * volume; From 029bd0740d36f54343bd0893aef51e0a6afff229 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 19 Nov 2024 10:12:50 -0600 Subject: [PATCH 09/11] make default aperture consistent with table --- .../ExponentialDecayPermeability_conformingFracture_base.xml | 2 +- .../ExponentialDecayPermeability_edfm_base.xml | 2 +- .../poromechanicsFractures/SlipPermeability_embeddedFrac.xml | 2 +- .../poromechanicsFractures/SlipPermeability_pEDFM_base.xml | 2 +- .../WillisRichardsPermeability_efem-edfm_base.xml | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) mode change 100755 => 100644 inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml mode change 100755 => 100644 inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml mode change 100755 => 100644 inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml diff --git a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_conformingFracture_base.xml b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_conformingFracture_base.xml index 133e70876ff..71594ab2961 100644 --- a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_conformingFracture_base.xml +++ b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_conformingFracture_base.xml @@ -89,7 +89,7 @@ name="Fracture" faceBlock="FractureSubRegion" materialList="{ water, fractureFilling, fractureContact, rock, hApertureModel}" - defaultAperture="1e-3"/> + defaultAperture="1e-4"/> diff --git a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml index a9f461c4e41..d3dc7017974 100644 --- a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml +++ b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml @@ -67,7 +67,7 @@ faceBlock="embeddedSurfaceSubRegion" subRegionType="embeddedElement" materialList="{ water, fractureFilling, frictionLaw, hApertureModel }" - defaultAperture="1e-3"/> + defaultAperture="1e-4"/> diff --git a/inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml b/inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml old mode 100755 new mode 100644 index 27a00a1005d..214490fcb6d --- a/inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml +++ b/inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml @@ -173,7 +173,7 @@ faceBlock="embeddedSurfaceSubRegion" subRegionType="embeddedElement" materialList="{ water, fractureFilling, fractureContact, hApertureModel }" - defaultAperture="1e-3"/> + defaultAperture="1e-4"/> diff --git a/inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml b/inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml old mode 100755 new mode 100644 index 0d8af6099a1..081c7191e2b --- a/inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml +++ b/inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml @@ -65,7 +65,7 @@ faceBlock="embeddedSurfaceSubRegion" subRegionType="embeddedElement" materialList="{ water, fractureFilling, fractureContact, hApertureModel}" - defaultAperture="1e-3"/> + defaultAperture="1e-4"/> diff --git a/inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml b/inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml old mode 100755 new mode 100644 index 5227b229979..953d63239df --- a/inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml +++ b/inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml @@ -68,7 +68,7 @@ faceBlock="embeddedSurfaceSubRegion" subRegionType="embeddedElement" materialList="{ water, fractureFilling, fractureContact, hApertureModel }" - defaultAperture="1e-3"/> + defaultAperture="1e-4"/> From c48028f09d0cb85551988681fb369c85e4899425 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 21 Nov 2024 11:00:04 -0600 Subject: [PATCH 10/11] Update ExponentialDecayPermeability_edfm_base.xml --- .../ExponentialDecayPermeability_edfm_base.xml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml index e14dccc51ed..2d6d1b26f68 100644 --- a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml +++ b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml @@ -13,9 +13,7 @@ + maxTimeStepCuts="1"/> From 805e889f1b988aad7280f1aadc5e46a66f56a1b9 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 12 Dec 2024 10:44:27 -0600 Subject: [PATCH 11/11] Update ThermalSinglePhasePoromechanicsEFEM_impl.hpp --- .../ThermalSinglePhasePoromechanicsEFEM_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp index fd1a33dce86..ac4475339d2 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanicsEFEM_impl.hpp @@ -176,7 +176,7 @@ complete( localIndex const k, // Energy balance accumulation real64 const volume = m_elementVolumeFrac( embSurfIndex ) + m_deltaVolume( embSurfIndex ); real64 const fluidEnergy = m_fluidDensity( embSurfIndex, 0 ) * m_fluidInternalEnergy( embSurfIndex, 0 ) * volume; - real64 const fluidEnergy_n = m_energy_n[embSurfIndex]; // TODO where is solid energy? + real64 const fluidEnergy_n = m_energy_n[embSurfIndex]; stack.dFluidMassIncrement_dTemperature = m_dFluidDensity_dTemperature( embSurfIndex, 0 ) * volume;