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

feat: SCEC benchmarks BP6 and BP7 #3523

Draft
wants to merge 11 commits into
base: develop
Choose a base branch
from
+249 −150
Draft

feat: SCEC benchmarks BP6 and BP7 #3523

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
c6e3171
iterative solver + reference solution in plotting script.
CusiniM Jan 28, 2025
3e2796d
remove separate displacement component filter.
CusiniM Jan 29, 2025
40c44dc
Use iterative solvers
VidarStiernstrom Jan 29, 2025
f3467cf
Merge branch 'develop' into cusini/scec-benchmarks
VidarStiernstrom Jan 29, 2025
2006498
Adjust shear modulus in xml to translate from value used in antiplane…
VidarStiernstrom Jan 29, 2025
3fdd2a7
wip: Update xml files and python script to include more mesurement p…
VidarStiernstrom Jan 31, 2025
3af1996
Try resetting variable_n after each step.
CusiniM Jan 31, 2025
d6bd871
Merge branch 'cusini/scec-benchmarks' of github.com:GEOS-DEV/GEOS int…
CusiniM Jan 31, 2025
c5eca39
Restructure reset calls to stress solver to work with the FSAL Runge-…
VidarStiernstrom Feb 1, 2025
8dfaf67
Override resetStateToBeginningOfStep in SpringSlider
VidarStiernstrom Feb 1, 2025
a3e4fca
Add plot of slip rate to python script
VidarStiernstrom Feb 1, 2025
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
9 changes: 5 additions & 4 deletions inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,8 +57,8 @@
lineSearchMaxCuts="2"
maxTimeStepCuts="1"/>
<LinearSolverParameters
solverType="direct"
directParallel="0"
solverType="gmres"
preconditionerType="mgr"
logLevel="0"/>
</SolidMechanicsLagrangeContactBubbleStab>

Expand All @@ -78,8 +78,9 @@
newtonTol = "1.0e-6"
newtonMaxIter ="8"/>
<LinearSolverParameters
solverType="direct"
directParallel="0"/>
solverType="gmres"
preconditionerType="amg"
logLevel="0"/>
</SinglePhaseFVM>
</Solvers>

Expand Down
75 changes: 0 additions & 75 deletions inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,54 +1,5 @@
<?xml version="1.0" ?>
<Problem>
<Solvers
gravityVector="{0.0, 0.0, 0.0}">
<!-- <ImplicitQuasiDynamicEQ
name="QDSolver"
targetRegions="{ Domain, Fault }"
shearImpedance="4.6247113164e6"
initialDt="1000"
logLevel="1"
discretization="FE1"
stressSolverName="stressSolver"
targetSlipIncrement="1.0e-6">
<NonlinearSolverParameters
newtonTol="1.0e-8"
logLevel="2"
newtonMaxIter="20"
maxNumConfigurationAttempts="1"
lineSearchAction="Require"
lineSearchMaxCuts="2"
maxTimeStepCuts="1"/>
</QuasiDynamicEQ> -->

<SolidMechanicsLagrangeContactBubbleStab
name="stressSolver"
timeIntegrationOption="QuasiStatic"
logLevel="2"
writeLinearSystem="0"
discretization="FE1"
targetRegions="{ Domain, Fault }">
<NonlinearSolverParameters
newtonTol="1.0e-8"
logLevel="2"
newtonMaxIter="10"
maxNumConfigurationAttempts="1"
lineSearchAction="Require"
lineSearchMaxCuts="2"
maxTimeStepCuts="1"/>
<LinearSolverParameters
solverType="direct"
directParallel="0"
logLevel="0"/>
</SolidMechanicsLagrangeContactBubbleStab>

<SurfaceGenerator
name="SurfaceGen"
targetRegions="{ Domain }"
initialRockToughness="1.0"
mpiCommOrder="1"
fractureRegion="Fault"/>
</Solvers>

<NumericalMethods>
<FiniteElements>
Expand Down Expand Up @@ -222,30 +173,4 @@
fieldName="pressure"
setNames="{source}"/>
</Tasks> -->

<Events
maxTime="4.5e4"
maxCycle="1">

<SoloEvent
name="generateFault"
target="/Solvers/SurfaceGen"/>

<PeriodicEvent
name="vtkOutput"
cycleFrequency="1"
targetExactTimestep="0"
target="/Outputs/vtkOutput"/>

<PeriodicEvent
name="solverApplications"
maxEventDt="1e4"
target="/Solvers/stressSolver"/>

<!-- <PeriodicEvent
name="resarts"
timeFrequency="2e4"
targetExactTimestep="0"
target="/Outputs/restart"/> -->
</Events>
</Problem>
84 changes: 84 additions & 0 deletions inputFiles/inducedSeismicity/SCEC_BP7_QD_implicit.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
<?xml version="1.0" ?>
<Problem>

<Included>
<File
name="./SCEC_BP7_QD_base.xml"/>
</Included>

<Solvers
gravityVector="{0.0, 0.0, 0.0}">
<ImplicitQuasiDynamicEQ
name="QDSolver"
targetRegions="{ Domain, Fault }"
shearImpedance="4.6247113164e6"
initialDt="1000"
logLevel="1"
discretization="FE1"
stressSolverName="stressSolver"
targetSlipIncrement="1.0e-6">
<NonlinearSolverParameters
newtonTol="1.0e-8"
logLevel="2"
newtonMaxIter="20"
maxNumConfigurationAttempts="1"
lineSearchAction="Require"
lineSearchMaxCuts="2"
maxTimeStepCuts="1"/>
</ImplicitQuasiDynamicEQ>

<SolidMechanicsLagrangeContactBubbleStab
name="stressSolver"
timeIntegrationOption="QuasiStatic"
logLevel="2"
writeLinearSystem="0"
discretization="FE1"
targetRegions="{ Domain, Fault }">
<NonlinearSolverParameters
newtonTol="1.0e-8"
logLevel="2"
newtonMaxIter="10"
maxNumConfigurationAttempts="1"
lineSearchAction="Require"
lineSearchMaxCuts="2"
maxTimeStepCuts="1"/>
<LinearSolverParameters
solverType="direct"
directParallel="0"
logLevel="0"/>
</SolidMechanicsLagrangeContactBubbleStab>

<SurfaceGenerator
name="SurfaceGen"
targetRegions="{ Domain }"
initialRockToughness="1.0"
mpiCommOrder="1"
fractureRegion="Fault"/>
</Solvers>

<Events
maxTime="4.5e4"
maxCycle="1">

<SoloEvent
name="generateFault"
target="/Solvers/SurfaceGen"/>

<PeriodicEvent
name="vtkOutput"
cycleFrequency="1"
targetExactTimestep="0"
target="/Outputs/vtkOutput"/>

<PeriodicEvent
name="solverApplications"
maxEventDt="1e4"
target="/Solvers/stressSolver"/>

<!-- <PeriodicEvent
name="resarts"
timeFrequency="2e4"
targetExactTimestep="0"
target="/Outputs/restart"/> -->
</Events>
</Problem>
25 changes: 23 additions & 2 deletions inputFiles/inducedSeismicity/scripts/plotBP6_results.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
from geos.hdf5_wrapper import hdf5_wrapper
import matplotlib.pyplot as plt
import os
from scipy.special import erfc

def remove_padding(data):
"""Removes trailing zeros from a NumPy array."""
Expand All @@ -22,9 +23,23 @@ def getDataFromHDF5( hdf5FilePath, var_name ):
# Remove padding
var = remove_padding(var)
time = time[:len(var)] # Match the length of the time array to the variable
return time, var
return time, var

def analytical_pressure( time, x ):
alpha = 0.1 # m^2 / s
phi = 0.1
beta = 1e-8 # Pa^-1
t_off = 100 * 24 * 3600
q0 = 1.25 * 1e-6 #m/s

pressure = q0 / (beta * phi * np.sqrt(alpha)) * (G(time, x, alpha) - np.where(time > t_off, G(time - t_off, x, alpha), 0))
return pressure

def G( t, x, alpha ):
A = np.abs(x) / np.sqrt( 4*alpha*t )
B = -x**2 / (4 * alpha * t)
return np.sqrt(t) * ( np.exp(B) / np.sqrt(np.pi) - A * erfc( A ) )

if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('-d', '--dir', type=str, help='Path to hdf5 file')
Expand All @@ -43,20 +58,26 @@ def getDataFromHDF5( hdf5FilePath, var_name ):
time, pressure = getDataFromHDF5( filePath, "pressure" )
time, slipRate = getDataFromHDF5( filePath, "slipRate" )

pressure_analytical = analytical_pressure( time, 0.0 )
print(pressure_analytical)

# Convert time to years
time_in_years = time / (365 * 24 * 3600) # Assuming time is in seconds

# Plotting
fig, ax1 = plt.subplots()

print(len(time_in_years))
print(len(pressure))
print(len(pressure))
print(len(pressure_analytical))

# Plot pressure on the left y-axis
ax1.set_xlabel('Time (years)')
ax1.set_ylabel('Pressure (Pa)', color='tab:blue')
ax1.plot(time_in_years, pressure, label="Pressure", color='tab:blue')
ax1.plot(time_in_years, pressure_analytical, label="Pressure Analytical", color='black', linestyle='--')
ax1.tick_params(axis='y', labelcolor='tab:blue')
ax1.legend()


# Plot slipRate on the right y-axis (log scale)
Expand Down
7 changes: 4 additions & 3 deletions inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_smoke.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,10 @@
lineSearchMaxCuts="2"
maxTimeStepCuts="1"/>
<LinearSolverParameters
solverType="direct"
solverType="gmres"
preconditionerType="mgr"
directParallel="0"
logLevel="0"/>
logLevel="1"/>
</SolidMechanicsLagrangeContactBubbleStab>

<SurfaceGenerator
Expand All @@ -43,7 +44,7 @@
yCoords="{ -5, 5 }"
zCoords="{ -0.05, 0.05 }"
nx="{ 400 }"
ny="{ 100}"
ny="{ 100 }"
nz="{ 1 }"
cellBlockNames="{ cb1 }"/>
</Mesh>
Expand Down
3 changes: 2 additions & 1 deletion src/coreComponents/linearAlgebra/interfaces/hypre/HypreUtils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -541,7 +541,8 @@ enum class MGRFRelaxationType : HYPRE_Int
l1backwardGaussSeidel = 14, //!< \f$\ell_1\f$ Gauss-Seidel, backward solve
l1jacobi = 18, //!< \f$\ell_1\f$-scaled Jacobi
gsElimWPivoting = 99, //!< Gaussian Elimination with pivoting direct solver (for small systems)
gsElimWInverse = 199 //!< Direct Inversion with Gaussian Elimination (OK for larger systems)
gsElimWInverse = 199, //!< Direct Inversion with Gaussian Elimination (OK for larger systems)
blockJacobi = 0
};

/**
Expand Down
8 changes: 6 additions & 2 deletions src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ SolidMechanicsLagrangeContactBubbleStab::SolidMechanicsLagrangeContactBubbleStab

LinearSolverParameters & linSolParams = m_linearSolverParameters.get();
linSolParams.mgr.strategy = LinearSolverParameters::MGR::StrategyType::lagrangianContactMechanicsBubbleStab;
linSolParams.mgr.separateComponents = true;
linSolParams.mgr.separateComponents = false;
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@victorapm for this physics solver, which uses this MGR strategy, the boomer AMG setup crashes in parallel during setup. Without the filter everything seems to work fine.

The case that I have been using to debug it is: https://github.com/GEOS-DEV/GEOS/blob/develop/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_smoke.xml and you can reduce the size to :
nx = 4 and ny = 2 and still reproduce the issue. A partition of -x 2 is sufficient to trigger the crash.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

this is the stack trace of the crash:

Number of element for each fracture state: stick:          380 | new slip:            0 | slip:             0 | open:             0
    Attempt:  0, ConfigurationIter:  0, NewtonIter:  0
        ( Rsolid ) = ( 0.00e+00 )        ( Rt  ) = (    1.549432e+01  )        ( R ) = ( 1.55e+01 )
        MGR preconditioner: numComponentsPerField = [3, 3, 3]
Received signal 11: Segmentation fault: 11

** StackTrace of 0 frames **
Frame 0: _sigtramp 
=====

Received signal 15: Terminated: 15

** StackTrace of 24 frames **
Frame 0: _sigtramp 
Frame 1: opal_progress 
Frame 2: mca_pml_ob1_iprobe 
Frame 3: MPI_Iprobe 
Frame 4: hypre_MPI_Iprobe 
Frame 5: hypre_DataExchangeList 
Frame 6: hypre_ParCSRCommPkgCreateApart_core 
Frame 7: hypre_ParCSRCommPkgCreateApart 
Frame 8: hypre_MatvecCommPkgCreate 
Frame 9: hypre_BoomerAMGBuildModExtPEInterpHost 
Frame 10: hypre_BoomerAMGBuildModExtPEInterp 
Frame 11: hypre_BoomerAMGSetup 
Frame 12: hypre_MGRSetup 
Frame 13: geos::HyprePreconditioner::setup(geos::HypreMatrix const&) 
Frame 14: geos::HypreSolver::setup(geos::HypreMatrix const&) 
Frame 15: geos::PhysicsSolverBase::solveLinearSystem(geos::DofManager const&, geos::HypreMatrix&, geos::HypreVector&, geos::HypreVector&) 
Frame 16: geos::PhysicsSolverBase::solveNonlinearSystem(double const&, double const&, int, geos::DomainPartition&) 
Frame 17: geos::PhysicsSolverBase::nonlinearImplicitStep(double const&, double const&, int, geos::DomainPartition&) 
--------------------------------------------------------------------------
prterun noticed that process rank 1 with PID 5091 on node onix exited on
signal 6 (Abort trap: 6).

Copy link
Contributor

Choose a reason for hiding this comment

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

Thanks for reporting. Let me take a look

linSolParams.dofsPerNode = 3;
}

Expand Down Expand Up @@ -896,6 +896,8 @@ void SolidMechanicsLagrangeContactBubbleStab::updateStickSlipList( DomainPartiti

arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >();

// std::cout << fractureState << std::endl;

forFiniteElementOnFractureSubRegions( meshName, [&] ( string const & finiteElementName,
finiteElement::FiniteElementBase const &,
arrayView1d< localIndex const > const & faceElementList )
Expand All @@ -910,6 +912,8 @@ void SolidMechanicsLagrangeContactBubbleStab::updateStickSlipList( DomainPartiti

arrayView1d< localIndex > const keys_v = keys.toView();
arrayView1d< localIndex > const vals_v = vals.toView();

// std::cout<< faceElementList.size() << std::endl;
forAll< parallelDevicePolicy<> >( faceElementList.size(),
[ = ]
GEOS_HOST_DEVICE ( localIndex const kfe )
Expand Down Expand Up @@ -965,7 +969,7 @@ void SolidMechanicsLagrangeContactBubbleStab::updateStickSlipList( DomainPartiti
this->m_faceTypesToFaceElementsStick[meshName][finiteElementName] = stickList;
this->m_faceTypesToFaceElementsSlip[meshName][finiteElementName] = slipList;

GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::Configuration, GEOS_FMT( "# stick elements: {}, # slip elements: {}", nStick, nSlip ))
GEOS_LOG_LEVEL_INFO_BY_RANK( logInfo::Configuration, GEOS_FMT( "# stick elements: {}, # slip elements: {}", nStick, nSlip ))
} );
} );

Expand Down
Loading