Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix: Use aperture table in poromechanics with conforming fractures #3194

Merged
merged 34 commits into from
Jul 10, 2024
Merged

fix: Use aperture table in poromechanics with conforming fractures #3194

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
34 commits
Select commit Hold shift + click to select a range
adbd3c2
decouple debug matrix output from logLevel
Jun 21, 2024
c036976
Merge remote-tracking branch 'origin/develop' into pt/write-matrix
Jun 21, 2024
65fa11a
victor's suggestions
Jun 21, 2024
ffea8de
debug contact
Jun 27, 2024
06c5d84
Merge branch 'develop' into pt/debug-contact
paveltomin Jun 28, 2024
d7df1f0
code style and remove debug
Jun 28, 2024
1ba9393
Update SolidMechanicsLagrangeContact.cpp
paveltomin Jun 28, 2024
194602b
Update SinglePhasePoromechanicsConformingFractures.cpp
paveltomin Jun 28, 2024
2950eaa
Update SolidMechanicsLagrangeContact.cpp
paveltomin Jun 28, 2024
c678634
fix build
Jun 28, 2024
f244566
Merge branch 'develop' into pt/debug-contact
paveltomin Jul 2, 2024
7577b5d
last update
Jul 2, 2024
3da4a7e
fix
Jul 2, 2024
ce6c22f
Merge branch 'develop' into pt/debug-contact
paveltomin Jul 2, 2024
8bf8e02
Merge branch 'develop' into pt/debug-contact
paveltomin Jul 3, 2024
dd88474
fix
Jul 3, 2024
1a0347c
Merge branch 'pt/debug-contact' of https://github.com/GEOS-DEV/GEOS i…
Jul 3, 2024
dcf2308
build fix
Jul 3, 2024
5f22591
fix edfm
Jul 3, 2024
d9b7330
revert and try another way
Jul 4, 2024
d5ddb7f
build fix
Jul 5, 2024
b12f753
this should work
Jul 5, 2024
47e17f1
Update SurfaceElementSubRegion.hpp
paveltomin Jul 5, 2024
8c7a8d6
one last try
Jul 5, 2024
90923b9
Merge branch 'pt/debug-contact' of https://github.com/GEOS-DEV/GEOS i…
Jul 5, 2024
fa97b6a
revert and give up
Jul 5, 2024
3839b4a
Update SinglePhaseBase.cpp
paveltomin Jul 5, 2024
45b2ca0
Update SinglePhaseBase.cpp
paveltomin Jul 5, 2024
17dc67e
revert
Jul 5, 2024
3b5fc93
Merge branch 'pt/debug-contact' of https://github.com/GEOS-DEV/GEOS i…
Jul 5, 2024
cae43ba
Update HydrofractureSolver.cpp
paveltomin Jul 5, 2024
e5ba0ce
Merge branch 'develop' into pt/debug-contact
paveltomin Jul 10, 2024
07a4d1d
Update .integrated_tests.yaml
paveltomin Jul 10, 2024
59f42d2
Update BASELINE_NOTES.md
paveltomin Jul 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 20 additions & 12 deletions src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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];
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is just to avoid gcc warnings

localIndex const f1 = elemsToFaces[k][1];

real64 Nbar[3];
Nbar[0] = faceNormal[f0][0] - faceNormal[f1][0];
Expand Down
93 changes: 75 additions & 18 deletions src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -252,6 +252,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 @@ -301,6 +308,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 @@ -384,15 +398,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 @@ -562,7 +597,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 All @@ -585,10 +620,13 @@ void SolidMechanicsLagrangeContact::assembleContact( DomainPartition & domain,
{
/// assemble Kut
assembleForceResidualDerivativeWrtTraction( mesh, regionNames, dofManager, localMatrix, localRhs );
// std::cout << "assembleForceResidualDerivativeWrtTraction RHS " << localRhs[96609] << std::endl;
/// assemble Ktu, Ktt blocks.
assembleTractionResidualDerivativeWrtDisplacementAndTraction( mesh, regionNames, dofManager, localMatrix, localRhs );
// std::cout << "assembleTractionResidualDerivativeWrtDisplacementAndTraction RHS " << localRhs[96609] << std::endl;
/// assemble stabilization
assembleStabilization( mesh, domain.getNumericalMethodManager(), dofManager, localMatrix, localRhs );
// std::cout << "assembleStabilization RHS " << localRhs[96609] << std::endl;
} );
}

Expand Down Expand Up @@ -715,6 +753,12 @@ real64 SolidMechanicsLagrangeContact::calculateContactResidualNorm( DomainPartit
arrayView1d< integer const > const & fractureState = subRegion.getField< contact::fractureState >();
arrayView1d< real64 const > const & area = subRegion.getElementArea();

// string const & contactRelationName = subRegion.template getReference< string >( viewKeyStruct::contactRelationNameString() );
// ContactBase const & contact = getConstitutiveModel< ContactBase >( subRegion, contactRelationName );
//
// constitutiveUpdatePassThru( contact, [&] ( auto & castedContact )
// {

RAJA::ReduceSum< parallelHostReduce, real64 > stickSum( 0.0 );
RAJA::ReduceSum< parallelHostReduce, real64 > slipSum( 0.0 );
RAJA::ReduceMax< parallelHostReduce, real64 > slipMax( 0.0 );
Expand Down Expand Up @@ -764,6 +808,7 @@ real64 SolidMechanicsLagrangeContact::calculateContactResidualNorm( DomainPartit
slipNormalizer = LvArray::math::max( slipNormalizer, slipMax.get());
openResidual += openSum.get();
openNormalizer = LvArray::math::max( openNormalizer, openMax.get());
// } );
} );
} );

Expand Down Expand Up @@ -1763,6 +1808,8 @@ void SolidMechanicsLagrangeContact::assembleStabilization( MeshLevel const & mes
nDof[1 - kf] );
}

//if(l ocalRow==96609)
//std::cout << "stab " <<localRhs[localRow + idof] << " " << rhs[idof] << std::endl;
// residual
RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow + idof], rhs[idof] );
}
Expand Down Expand Up @@ -1896,20 +1943,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] )
{
fractureState[kfe] = contact::FractureState::Open;
}
else
if( dispJump[kfe][0] <= -normalDisplacementTolerance[kfe] )
{
fractureState[kfe] = contact::FractureState::Stick;
fractureState[kfe] = FractureState::Stick;
if( getLogLevel() >= 3 )
GEOS_LOG( GEOS_FMT( "{}: {} -> {}: dispJump = {}, normalDisplacementTolerance = {}",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it would be better to avoid this logging inside the cell loop.

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() >= 3 )
GEOS_LOG( GEOS_FMT( "{}: {} -> {}: traction = {}, normalTractionTolerance = {}",
kfe, originalFractureState, fractureState[kfe],
traction[kfe][0], normalTractionTolerance[kfe] ) );
}
else
{
Expand All @@ -1920,29 +1971,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() >= 3 && fractureState[kfe] != originalFractureState )
GEOS_LOG( GEOS_FMT( "{}: {} -> {}: currentTau = {}, limitTau = {}",
kfe, originalFractureState, fractureState[kfe],
currentTau, limitTau ) );
}

changed += faceArea[kfe] * !compareFractureStates( originalFractureState, fractureState[kfe] );
Expand All @@ -1964,21 +2019,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
5 changes: 4 additions & 1 deletion src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
View file
Open in desktop
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 All @@ -773,6 +774,8 @@ void SinglePhaseBase::implicitStepComplete( real64 const & time,
porousSolid.saveConvergedState(); // porosity_n <- porosity
}

updateMass( subRegion );
paveltomin marked this conversation as resolved.
Show resolved Hide resolved

if( m_isThermal )
{
arrayView2d< real64 const > const porosity = porousSolid.getPorosity();
Expand Down
Loading