Skip to content

Commit

Permalink
fix: Use aperture table in poromechanics with conforming fractures (#…
Browse files Browse the repository at this point in the history 
…3194)

* decouple debug matrix output from logLevel

* Update SolidMechanicsLagrangeContact.cpp

* Update SinglePhasePoromechanicsConformingFractures.cpp

* Update SolidMechanicsLagrangeContact.cpp

* Update SurfaceElementSubRegion.hpp

* Update SinglePhaseBase.cpp

* Update HydrofractureSolver.cpp

* Update .integrated_tests.yaml

* Update BASELINE_NOTES.md

---------

Co-authored-by: Pavel Tomin <“[email protected]”>
  • Loading branch information
@paveltomin @Algiane
2 people authored and Algiane committed Jul 30, 2024
1 parent 0661702 commit dc6b2d1
Show file tree
Hide file tree
Showing 13 changed files with 160 additions and 102 deletions.
2 changes: 1 addition & 1 deletion .integrated_tests.yaml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
baselines:
bucket: geosx
baseline: integratedTests/baseline_integratedTests-pr3006-5860-947a907
baseline: integratedTests/baseline_integratedTests-pr3194-6060-e5ba0ce

allow_fail:
all: ''
Expand Down
8 changes: 8 additions & 0 deletions BASELINE_NOTES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,18 +6,26 @@ 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 #3194 (2024-07-10)
======================
Use aperture table in poromechanics with conforming fractures. Rebaseline the corresponding cases.


PR #3006 (2024-07-01)
======================
Added baselines for new tests. Relaxing tolerances for singlePhasePoromechanics_FaultModel_smoke.


PR #3196 (2024-06-28)
======================
Added isLaggingFractureStencilWeightsUpdate to hydrofracture solve. Rebaseline because of the new input.


PR #3177 (2024-06-28)
======================
Added logLevel to TimeHistoryOutput. Rebaseline because of the new input flag.


PR #3181 (2024-06-25)
======================
Decouple debug matrix output from logLevel. Rebaseline because of the new input flag.
Expand Down
32 changes: 20 additions & 12 deletions src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,19 +122,20 @@ void ContactSolverBase::setFractureRegions( dataRepository::Group const & meshBo

void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh,
globalIndex & numStick,
globalIndex & numNewSlip,
globalIndex & numSlip,
globalIndex & numOpen ) const
{
ElementRegionManager const & elemManager = mesh.getElemManager();

array1d< globalIndex > localCounter( 3 );
array1d< globalIndex > localCounter( 4 );

elemManager.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion const & subRegion )
{
arrayView1d< integer const > const & ghostRank = subRegion.ghostRank();
arrayView1d< integer const > const & fractureState = subRegion.getField< fields::contact::fractureState >();

RAJA::ReduceSum< parallelHostReduce, localIndex > stickCount( 0 ), slipCount( 0 ), openCount( 0 );
RAJA::ReduceSum< parallelHostReduce, localIndex > stickCount( 0 ), newSlipCount( 0 ), slipCount( 0 ), openCount( 0 );
forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe )
{
if( ghostRank[kfe] < 0 )
Expand All @@ -147,6 +148,10 @@ void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh,
break;
}
case FractureState::NewSlip:
{
newSlipCount += 1;
break;
}
case FractureState::Slip:
{
slipCount += 1;
Expand All @@ -162,40 +167,43 @@ void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh,
} );

localCounter[0] += stickCount.get();
localCounter[1] += slipCount.get();
localCounter[2] += openCount.get();
localCounter[1] += newSlipCount.get();
localCounter[2] += slipCount.get();
localCounter[3] += openCount.get();
} );

array1d< globalIndex > totalCounter( 3 );
array1d< globalIndex > totalCounter( 4 );

MpiWrapper::allReduce( localCounter.data(),
totalCounter.data(),
3,
4,
MPI_SUM,
MPI_COMM_GEOSX );

numStick = totalCounter[0];
numSlip = totalCounter[1];
numOpen = totalCounter[2];
numStick = totalCounter[0];
numNewSlip = totalCounter[1];
numSlip = totalCounter[2];
numOpen = totalCounter[3];
}

void ContactSolverBase::outputConfigurationStatistics( DomainPartition const & domain ) const
{
if( getLogLevel() >=1 )
{
globalIndex numStick = 0;
globalIndex numNewSlip = 0;
globalIndex numSlip = 0;
globalIndex numOpen = 0;

forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
MeshLevel const & mesh,
arrayView1d< string const > const & )
{
computeFractureStateStatistics( mesh, numStick, numSlip, numOpen );
computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen );

GEOS_LOG_RANK_0( GEOS_FMT( " Number of element for each fracture state:"
" stick: {:12} | slip: {:12} | open: {:12}",
numStick, numSlip, numOpen ) );
" stick: {:12} | new slip: {:12} | slip: {:12} | open: {:12}",
numStick, numNewSlip, numSlip, numOpen ) );
} );
}
}
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,7 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM

void computeFractureStateStatistics( MeshLevel const & mesh,
globalIndex & numStick,
globalIndex & numNewSlip,
globalIndex & numSlip,
globalIndex & numOpen ) const;

Expand Down
4 changes: 2 additions & 2 deletions src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,8 +406,8 @@ struct ComputeRotationMatricesKernel
forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
{

localIndex const & f0 = elemsToFaces[k][0];
localIndex const & f1 = elemsToFaces[k][1];
localIndex const f0 = elemsToFaces[k][0];
localIndex const f1 = elemsToFaces[k][1];

real64 Nbar[3];
Nbar[0] = faceNormal[f0][0] - faceNormal[f1][0];
Expand Down
81 changes: 63 additions & 18 deletions src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,6 +253,13 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain
{
GEOS_MARK_FUNCTION;

real64 minNormalTractionTolerance( 1e10 );
real64 maxNormalTractionTolerance( -1e10 );
real64 minNormalDisplacementTolerance( 1e10 );
real64 maxNormalDisplacementTolerance( -1e10 );
real64 minSlidingTolerance( 1e10 );
real64 maxSlidingTolerance( -1e10 );

forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
arrayView1d< string const > const & )
Expand Down Expand Up @@ -302,6 +309,13 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain
arrayView1d< real64 > const & slidingTolerance =
subRegion.getReference< array1d< real64 > >( viewKeyStruct::slidingToleranceString() );

RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionNormalTractionTolerance( 1e10 );
RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionNormalTractionTolerance( -1e10 );
RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionNormalDisplacementTolerance( 1e10 );
RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionNormalDisplacementTolerance( -1e10 );
RAJA::ReduceMin< ReducePolicy< parallelHostPolicy >, real64 > minSubRegionSlidingTolerance( 1e10 );
RAJA::ReduceMax< ReducePolicy< parallelHostPolicy >, real64 > maxSubRegionSlidingTolerance( -1e10 );

forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe )
{
if( ghostRank[kfe] < 0 )
Expand Down Expand Up @@ -385,15 +399,36 @@ void SolidMechanicsLagrangeContact::computeTolerances( DomainPartition & domain
LvArray::tensorOps::scale< 3, 3 >( rotatedInvStiffApprox, area );

// Finally, compute tolerances for the given fracture element

normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus / 2.e+7;
minSubRegionNormalDisplacementTolerance.min( normalDisplacementTolerance[kfe] );
maxSubRegionNormalDisplacementTolerance.max( normalDisplacementTolerance[kfe] );

slidingTolerance[kfe] = sqrt( rotatedInvStiffApprox[ 1 ][ 1 ] * rotatedInvStiffApprox[ 1 ][ 1 ] +
rotatedInvStiffApprox[ 2 ][ 2 ] * rotatedInvStiffApprox[ 2 ][ 2 ] ) * averageYoungModulus / 2.e+7;
minSubRegionSlidingTolerance.min( slidingTolerance[kfe] );
maxSubRegionSlidingTolerance.max( slidingTolerance[kfe] );

normalTractionTolerance[kfe] = 1.0 / 2.0 * averageConstrainedModulus / averageBoxSize0 * normalDisplacementTolerance[kfe];
minSubRegionNormalTractionTolerance.min( normalTractionTolerance[kfe] );
maxSubRegionNormalTractionTolerance.max( normalTractionTolerance[kfe] );
}
} );

minNormalDisplacementTolerance = std::min( minNormalDisplacementTolerance, minSubRegionNormalDisplacementTolerance.get() );
maxNormalDisplacementTolerance = std::max( maxNormalDisplacementTolerance, maxSubRegionNormalDisplacementTolerance.get() );
minSlidingTolerance = std::min( minSlidingTolerance, minSubRegionSlidingTolerance.get() );
maxSlidingTolerance = std::max( maxSlidingTolerance, maxSubRegionSlidingTolerance.get() );
minNormalTractionTolerance = std::min( minNormalTractionTolerance, minSubRegionNormalTractionTolerance.get() );
maxNormalTractionTolerance = std::max( maxNormalTractionTolerance, maxSubRegionNormalTractionTolerance.get() );
}
} );
} );

GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{}: normal displacement tolerance = [{}, {}], sliding tolerance = [{}, {}], normal traction tolerance = [{}, {}]",
this->getName(), minNormalDisplacementTolerance, maxNormalDisplacementTolerance,
minSlidingTolerance, maxSlidingTolerance,
minNormalTractionTolerance, maxNormalTractionTolerance ) );
}

void SolidMechanicsLagrangeContact::resetStateToBeginningOfStep( DomainPartition & domain )
Expand Down Expand Up @@ -586,7 +621,7 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time,

assembleContact( domain, dofManager, localMatrix, localRhs );

// for sequenatial: add (fixed) pressure force contribution into residual (no derivatives)
// for sequential: add (fixed) pressure force contribution into residual (no derivatives)
if( m_isFixedStressPoromechanicsUpdate )
{
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
Expand Down Expand Up @@ -2221,20 +2256,24 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
if( ghostRank[kfe] < 0 )
{
integer const originalFractureState = fractureState[kfe];
if( originalFractureState == contact::FractureState::Open )
if( originalFractureState == FractureState::Open )
{
if( dispJump[kfe][0] > -normalDisplacementTolerance[kfe] )
if( dispJump[kfe][0] <= -normalDisplacementTolerance[kfe] )
{
fractureState[kfe] = contact::FractureState::Open;
}
else
{
fractureState[kfe] = contact::FractureState::Stick;
fractureState[kfe] = FractureState::Stick;
if( getLogLevel() >= 10 )
GEOS_LOG( GEOS_FMT( "{}: {} -> {}: dispJump = {}, normalDisplacementTolerance = {}",
kfe, originalFractureState, fractureState[kfe],
dispJump[kfe][0], normalDisplacementTolerance[kfe] ) );
}
}
else if( traction[kfe][0] > normalTractionTolerance[kfe] )
{
fractureState[kfe] = contact::FractureState::Open;
fractureState[kfe] = FractureState::Open;
if( getLogLevel() >= 10 )
GEOS_LOG( GEOS_FMT( "{}: {} -> {}: traction = {}, normalTractionTolerance = {}",
kfe, originalFractureState, fractureState[kfe],
traction[kfe][0], normalTractionTolerance[kfe] ) );
}
else
{
Expand All @@ -2245,29 +2284,33 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
contactWrapper.computeLimitTangentialTractionNorm( traction[kfe][0],
dLimitTangentialTractionNorm_dTraction );

if( originalFractureState == contact::FractureState::Stick && currentTau >= limitTau )
if( originalFractureState == FractureState::Stick && currentTau >= limitTau )
{
currentTau *= (1.0 - m_slidingCheckTolerance);
}
else if( originalFractureState != contact::FractureState::Stick && currentTau <= limitTau )
else if( originalFractureState != FractureState::Stick && currentTau <= limitTau )
{
currentTau *= (1.0 + m_slidingCheckTolerance);
}
if( currentTau > limitTau )
{
if( originalFractureState == contact::FractureState::Stick )
if( originalFractureState == FractureState::Stick )
{
fractureState[kfe] = contact::FractureState::NewSlip;
fractureState[kfe] = FractureState::NewSlip;
}
else
{
fractureState[kfe] = contact::FractureState::Slip;
fractureState[kfe] = FractureState::Slip;
}
}
else
{
fractureState[kfe] = contact::FractureState::Stick;
fractureState[kfe] = FractureState::Stick;
}
if( getLogLevel() >= 10 && fractureState[kfe] != originalFractureState )
GEOS_LOG( GEOS_FMT( "{}: {} -> {}: currentTau = {}, limitTau = {}",
kfe, originalFractureState, fractureState[kfe],
currentTau, limitTau ) );
}

changed += faceArea[kfe] * !compareFractureStates( originalFractureState, fractureState[kfe] );
Expand All @@ -2289,21 +2332,23 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
// and total area of fracture elements
totalArea = MpiWrapper::sum( totalArea );

GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( " {}: changed area {} out of {}", getName(), changedArea, totalArea ) );

// Assume converged if changed area is below certain fraction of total area
return changedArea <= m_nonlinearSolverParameters.m_configurationTolerance * totalArea;
}

bool SolidMechanicsLagrangeContact::isFractureAllInStickCondition( DomainPartition const & domain ) const
{
globalIndex numStick, numSlip, numOpen;
globalIndex numStick, numNewSlip, numSlip, numOpen;
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel const & mesh,
arrayView1d< string const > const & )
{
computeFractureStateStatistics( mesh, numStick, numSlip, numOpen );
computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen );
} );

return ( ( numSlip + numOpen ) == 0 );
return ( ( numNewSlip + numSlip + numOpen ) == 0 );
}

real64 SolidMechanicsLagrangeContact::setNextDt( real64 const & currentDt,
Expand Down
3 changes: 0 additions & 3 deletions src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2585,9 +2585,6 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain )
{
GEOS_MARK_FUNCTION;

if( m_keepFlowVariablesConstantDuringInitStep )
return;

real64 maxDeltaPhaseVolFrac = 0.0;
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
MeshLevel & mesh,
Expand Down
3 changes: 0 additions & 3 deletions src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,9 +185,6 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies )
subRegion.registerField< fields::flow::hydraulicAperture >( getName() ).
setApplyDefaultValue( faceRegion.getDefaultAperture() );

subRegion.registerField< fields::flow::minimumHydraulicAperture >( getName() ).
setApplyDefaultValue( faceRegion.getDefaultAperture() );

} );

FaceManager & faceManager = mesh.getFaceManager();
Expand Down
8 changes: 0 additions & 8 deletions src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,14 +168,6 @@ DECLARE_FIELD( hydraulicAperture,
WRITE_AND_READ,
"Hydraulic aperture" );

DECLARE_FIELD( minimumHydraulicAperture,
"minimumHydraulicAperture",
array1d< real64 >,
0,
LEVEL_0,
WRITE_AND_READ,
"minimum value of the hydraulic aperture" );

DECLARE_FIELD( gravityCoefficient,
"gravityCoefficient",
array1d< real64 >,
Expand Down
8 changes: 3 additions & 5 deletions src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -750,12 +750,13 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,
singlePhaseBaseKernels::StatisticsKernel::
saveDeltaPressure( subRegion.size(), pres, initPres, deltaPres );

arrayView1d< real64 const > const dVol = subRegion.getField< fields::flow::deltaVolume >();
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;
} );

SingleFluidBase const & fluid =
Expand Down Expand Up @@ -803,7 +804,7 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,
if( volume[ei] * density_n[ei][0] > 1.1 * creationMass[ei] )
{
creationMass[ei] *= 0.75;
if( creationMass[ei]<1.0e-20 )
if( creationMass[ei] < 1.0e-20 )
{
creationMass[ei] = 0.0;
}
Expand Down Expand Up @@ -1255,9 +1256,6 @@ void SinglePhaseBase::updateState( DomainPartition & domain )
{
GEOS_MARK_FUNCTION;

if( m_keepFlowVariablesConstantDuringInitStep )
return;

forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
MeshLevel & mesh,
arrayView1d< string const > const & regionNames )
Expand Down
Loading

0 comments on commit dc6b2d1

Please sign in to comment.