diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml new file mode 100644 index 00000000000..1257f1ae3ec --- /dev/null +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_base.xml @@ -0,0 +1,217 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml new file mode 100644 index 00000000000..fbbfccd5b38 --- /dev/null +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml @@ -0,0 +1,139 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml new file mode 100644 index 00000000000..1b74521b1ad --- /dev/null +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml @@ -0,0 +1,140 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml new file mode 100644 index 00000000000..08fc7cb9930 --- /dev/null +++ b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml @@ -0,0 +1,251 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml deleted file mode 100644 index a0c0381fb23..00000000000 --- a/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml +++ /dev/null @@ -1,167 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml index 2de09b79a03..a1a67094f7c 100644 --- a/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml +++ b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml @@ -2,8 +2,19 @@ + name="./SpringSlider_base.xml"/> + + + + + @@ -23,6 +34,32 @@ maxEventDt="1e4" target="/Solvers/SpringSlider"/> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SpringSlider_base.xml b/inputFiles/inducedSeismicity/SpringSlider_base.xml index 5f6aefc94ad..181db4397bb 100644 --- a/inputFiles/inducedSeismicity/SpringSlider_base.xml +++ b/inputFiles/inducedSeismicity/SpringSlider_base.xml @@ -1,15 +1,6 @@ - - - - + diff --git a/inputFiles/inducedSeismicity/SpringSlider_smoke.xml b/inputFiles/inducedSeismicity/SpringSlider_smoke.xml deleted file mode 100644 index b10580131d8..00000000000 --- a/inputFiles/inducedSeismicity/SpringSlider_smoke.xml +++ /dev/null @@ -1,32 +0,0 @@ - - - - - - - - - - - - - - - - - \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/inducedSeismicity.ats b/inputFiles/inducedSeismicity/inducedSeismicity.ats index 813bdf08b8b..56aea6b5820 100644 --- a/inputFiles/inducedSeismicity/inducedSeismicity.ats +++ b/inputFiles/inducedSeismicity/inducedSeismicity.ats @@ -23,7 +23,7 @@ decks = [ check_step=100, curvecheck_params=CurveCheckParameters(**curvecheck_params)), TestDeck( - name="SpringSlider_smoke", + name="SpringSliderImplicit_smoke", description="Spring slider 0D system", partitions=((1, 1, 1), ), restart_step=0, @@ -34,7 +34,7 @@ decks = [ description="Spring slider 0D system", partitions=((1, 1, 1), ), restart_step=0, - check_step=532, + check_step=343, restartcheck_params=RestartcheckParameters(atol=1e-4, rtol=1e-3)) ] generate_geos_tests(decks) diff --git a/inputFiles/inducedSeismicity/scripts/plotBP6_results.py b/inputFiles/inducedSeismicity/scripts/plotBP6_results.py new file mode 100644 index 00000000000..e869c3d5b5a --- /dev/null +++ b/inputFiles/inducedSeismicity/scripts/plotBP6_results.py @@ -0,0 +1,78 @@ +import argparse +import numpy as np +from geos.hdf5_wrapper import hdf5_wrapper +import matplotlib.pyplot as plt +import os + +def remove_padding(data): + """Removes trailing zeros from a NumPy array.""" + nonzero_indices = np.nonzero(data)[0] + if nonzero_indices.size == 0: # If all elements are zero + return np.array([]), np.array([]) + last_nonzero = nonzero_indices[-1] + return data[:last_nonzero + 1] + +def getDataFromHDF5( hdf5FilePath, var_name ): + # Read HDF5 + data = hdf5_wrapper(f'{hdf5FilePath}').get_copy() + var = data[f'{var_name} source'] + var = np.asarray(var) + time = data[f'{var_name} Time'] + + # Remove padding + var = remove_padding(var) + time = time[:len(var)] # Match the length of the time array to the variable + return time, var + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument('-d', '--dir', type=str, help='Path to hdf5 file') + args = parser.parse_args() + + normalized_dir = os.path.abspath(args.dir) + if not os.path.isdir(normalized_dir): + print(f"Error: {normalized_dir} is not a valid directory.") + exit(1) + + filePath = f"{normalized_dir}/BP6_DQ_S.hdf5" + if not os.path.exists(filePath): + print(f"Error: {filePath} not found.") + exit(1) + + time, pressure = getDataFromHDF5( filePath, "pressure" ) + time, slipRate = getDataFromHDF5( filePath, "slipRate" ) + + # 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)) + + # 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.tick_params(axis='y', labelcolor='tab:blue') + + + # Plot slipRate on the right y-axis (log scale) + ax2 = ax1.twinx() + ax2.set_ylabel('Slip Rate (m/s)', color='tab:red') + ax2.plot(time_in_years, slipRate, label="Slip Rate", color='tab:red') + ax2.set_yscale('log') + ax2.tick_params(axis='y', labelcolor='tab:red') + + # Set x-axis limits to 0 to 2 years + ax1.set_xlim(0, np.max(time_in_years)) + + + # Add grid and title + plt.title("Pressure and Slip Rate vs Time") + plt.grid() + + # Show plot + plt.show() \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/scripts/plot_springSlider_results.py b/inputFiles/inducedSeismicity/scripts/plot_springSlider_results.py new file mode 100644 index 00000000000..ac56838015b --- /dev/null +++ b/inputFiles/inducedSeismicity/scripts/plot_springSlider_results.py @@ -0,0 +1,74 @@ +import argparse +import numpy as np +from geos.hdf5_wrapper import hdf5_wrapper +import matplotlib.pyplot as plt +import os + +def remove_padding(data): + """Removes trailing zeros from a NumPy array.""" + nonzero_indices = np.nonzero(data)[0] + if nonzero_indices.size == 0: # If all elements are zero + return np.array([]), np.array([]) + last_nonzero = nonzero_indices[-1] + return data[:last_nonzero + 1] + +def getDataFromHDF5( hdf5FilePath, var_name ): + # Read HDF5 + data = hdf5_wrapper(f'{hdf5FilePath}').get_copy() + var = data[f'{var_name}'] + var = np.asarray(var) + time = data[f'{var_name} Time'] + + # Remove padding + var = remove_padding(var) + time = time[:len(var)] # Match the length of the time array to the variable + return time, var + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument('-d', '--dir', type=str, help='Path to hdf5 file') + args = parser.parse_args() + + normalized_dir = os.path.abspath(args.dir) + if not os.path.isdir(normalized_dir): + print(f"Error: {normalized_dir} is not a valid directory.") + exit(1) + + filePath = f"{normalized_dir}/springSlider.hdf5" + if not os.path.exists(filePath): + print(f"Error: {filePath} not found.") + exit(1) + + time, slipRate = getDataFromHDF5( filePath, "slipRate" ) + time, stateVariable = getDataFromHDF5( filePath, "stateVariable" ) + + # Convert time to years + time_in_years = time / (365 * 24 * 3600) # Assuming time is in seconds + + # Plotting + fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8, 8), sharex=True) + + # Plot pressure on the left y-axis + ax1.set_xlabel('Time (years)') + ax1.set_ylabel('slip rate (Pa)', color='tab:blue') + ax1.plot(time_in_years, slipRate, label="slip rate", color='tab:blue') + ax1.set_yscale('log') + ax1.tick_params(axis='y', labelcolor='tab:blue') + + # Set x-axis limits to 0 to 2 years + ax1.set_xlim(0, np.max(time_in_years)) + + # Plot stateVariable on the right y-axis (log scale) + ax2.set_ylabel('state variable', color='tab:red') + ax2.plot(time_in_years, stateVariable, label="state variable", color='tab:red') + ax2.tick_params(axis='y', labelcolor='tab:red') + ax2.set_xlim(0, np.max(time_in_years)) + + # Add grid and title + fig.suptitle('Spring slider solution', fontsize=16) + plt.tight_layout() + plt.grid() + + # Show plot + plt.show() \ No newline at end of file diff --git a/inputFiles/lagrangianContactMechanics/SimpleCubes_base.xml b/inputFiles/lagrangianContactMechanics/SimpleCubes_base.xml index 32e64f498de..e7820ebbcad 100644 --- a/inputFiles/lagrangianContactMechanics/SimpleCubes_base.xml +++ b/inputFiles/lagrangianContactMechanics/SimpleCubes_base.xml @@ -134,7 +134,7 @@ - diff --git a/inputFiles/lagrangianContactMechanics/scripts/generateTables.py b/inputFiles/lagrangianContactMechanics/scripts/generateTables.py new file mode 100644 index 00000000000..9c83614b2f4 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/scripts/generateTables.py @@ -0,0 +1,179 @@ +import numpy as np +import os +import sys +import xml.etree.ElementTree as ElementTree +import matplotlib +import matplotlib.pyplot as plt + +class SingularCrackSlip: + + def __init__(self, mechanicalParameters, length ): + K = mechanicalParameters["bulkModulus"] + G = mechanicalParameters["shearModulus"] + poisson_ratio= (3 * K - 2 * G) / (2 * (3 * K + G)) + + mu_star = G /( 1 - poisson_ratio) + self.tau_0 = 0.0 + self.tau_r = -1.0 + + self.scaling = 2*(self.tau_0 - self.tau_r)/mu_star + self.halfLength = length + + def computeSlip(self, x): + return self.scaling * np.sqrt(self.halfLength**2 - x**2) + + def computeTraction(self, x): + if x < -self.halfLength or x > self.halfLength: + return self.tau_0 + (self.tau_0-self.tau_r) * ( np.abs(x)/np.sqrt(x**2 - self.halfLength**2) - 1 ) + else: + return self.tau_r +class GaussianSlip: + + def __init__(self, peakStrength, length ): + self.scaling = peakStrength + self.halfLength = length + + def computeSlip(self, x): + denom = 1 / (self.halfLength/2) + return self.scaling*np.exp(-0.5*((x)/denom)**2) + +def getMechanicalParametersFromXML(xmlFilePath): + tree = ElementTree.parse(xmlFilePath) + + param = tree.find('Constitutive/ElasticIsotropic') + + mechanicalParameters = dict.fromkeys(["bulkModulus", "shearModulus"]) + mechanicalParameters["bulkModulus"] = float(param.get("defaultBulkModulus")) + mechanicalParameters["shearModulus"] = float(param.get("defaultShearModulus")) + return mechanicalParameters + +def getFractureLengthFromXML(xmlFilePath): + tree = ElementTree.parse(xmlFilePath) + + rectangle = tree.find('Geometry/Box') + xmin = rectangle.get("xMin") + xmax = rectangle.get("xMax") + xmin = [float(i) for i in xmin[1:-1].split(",")] + xmax = [float(i) for i in xmax[1:-1].split(",")] + length = ( xmax[0] - xmin[0] ) / 2 + origin = 0.0 + + return length, origin + +def plot_traction_solution(): + # Read HDF5 + import hdf5_wrapper + hdf5File1Path = "Output/traction.hdf5" + + # Read HDF5 + data = hdf5_wrapper.hdf5_wrapper(hdf5File1Path).get_copy() + traction = data['traction'] + traction = np.asarray(traction) + traction_geos = traction[0, :, 1] + x = data['traction elementCenter'] + x_geos = x[0, :, 0] + + #-------- Extract info from XML + xmlFilePath = "/Users/cusini1/geosx/GEOSX/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_base.xml" + + mechanicalParameters = getMechanicalParametersFromXML(xmlFilePath) + + # Get length of the fracture + xmlFilePath = "/Users/cusini1/geosx/GEOSX/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_smoke.xml" + totalHalfLength, originShift = getFractureLengthFromXML(xmlFilePath) + halfLength = 2.0 + + singularCrackSlipSolution = SingularCrackSlip(mechanicalParameters, halfLength) + x = np.linspace(-totalHalfLength, totalHalfLength, 10000, endpoint=True) + traction_analytical = np.zeros(len(x)) + i = 0 + for xCell in x: + traction_analytical[i] = singularCrackSlipSolution.computeTraction(xCell) + i += 1 + + fsize = 30 + msize = 15 + lw = 2 + fig, ax = plt.subplots(1, figsize=(16, 12)) + cmap = plt.get_cmap("tab10") + + # Plot analytical (continuous line) and numerical (markers) aperture solution + ax.plot(x, traction_analytical, color='r', label='Traction analytical', lw=lw) + ax.plot(x_geos, traction_geos, color='k', label='geos', marker="o", lw=lw) + + ax.set_xlabel('Fault coordinate [m]', size=fsize, weight="bold") + ax.set_ylabel('Shear traction', size=fsize, weight="bold") + ax.legend(bbox_to_anchor=(0.75, 0.9), loc='center', borderaxespad=0., fontsize=fsize) + ax.xaxis.set_tick_params(labelsize=fsize) + ax.yaxis.set_tick_params(labelsize=fsize) + plt.savefig("traction.png") + +def output_tables(x, slip, name): + # Save x to x.csv with one value per row + np.savetxt('x.csv', x, fmt='%f') + + # Save aperture_analytical to jump.csv with one value per row + np.savetxt(f'{name}.csv', slip, fmt='%f') + + +def debug(): + #-------- Extract info from XML + xmlFilePath = "/Users/cusini1/geosx/GEOSX/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_base.xml" + + mechanicalParameters = getMechanicalParametersFromXML(xmlFilePath) + appliedPressure = 1.0 + + # Get length of the fracture + xmlFilePath = "/Users/cusini1/geosx/GEOSX/inputFiles/lagrangianContactMechanics/LagrangeContactBubbleStab_FixedSlip_smoke.xml" + totalHalfLength, originShift = getFractureLengthFromXML(xmlFilePath) + halfLength = 2.0 + + # Initialize Sneddon's analytical solution + singularCrackSlipSolution = SingularCrackSlip(mechanicalParameters, halfLength ) + peakStrength = 3.0 + gaussianSlipSolution = GaussianSlip( peakStrength, halfLength) + + # Plot analytical (continuous line) and numerical (markers) aperture solution + x = np.linspace(-totalHalfLength, totalHalfLength, 100000, endpoint=True) + singularCrackSlip = np.zeros(len(x)) + gaussianSlip = np.zeros(len(x)) + i = 0 + for xCell in x: + if xCell > -halfLength and xCell < halfLength: + singularCrackSlip[i] = singularCrackSlipSolution.computeSlip(xCell) + gaussianSlip[i] = gaussianSlipSolution.computeSlip(xCell) + i += 1 + + fsize = 24 + msize = 15 + lw = 6 + fig, ax = plt.subplots(1, figsize=(16, 12)) + cmap = plt.get_cmap("tab10") + + ax.plot(x, singularCrackSlip , color='k', label='Analytical Solution', lw=lw) + ax.grid() + ax.set_xlabel('Fault coordinate [m]', size=fsize, weight="bold") + ax.set_ylabel('slip [m]', size=fsize, weight="bold") + ax.legend(bbox_to_anchor=(0.7, 1), loc='center', borderaxespad=0., fontsize=fsize) + ax.xaxis.set_tick_params(labelsize=fsize) + ax.yaxis.set_tick_params(labelsize=fsize) + plt.savefig("singularCrackSlip.png") + + fig, ax = plt.subplots(1, figsize=(16, 12)) + cmap = plt.get_cmap("tab10") + + ax.plot(x, gaussianSlip , color='k', label='Analytical Solution', lw=lw) + ax.grid() + ax.set_xlabel('Fault coordinate [m]', size=fsize, weight="bold") + ax.set_ylabel('slip [m]', size=fsize, weight="bold") + ax.legend(bbox_to_anchor=(0.75, 0.9), loc='center', borderaxespad=0., fontsize=fsize) + ax.xaxis.set_tick_params(labelsize=fsize) + ax.yaxis.set_tick_params(labelsize=fsize) + plt.savefig("gaussianSlip.png") + + output_tables(x, singularCrackSlip, "singularCrackSlip") + output_tables(x, gaussianSlip, "gaussianSlip") + +if __name__ == "__main__": + debug() + plot_traction_solution() diff --git a/src/coreComponents/constitutive/ConstitutivePassThru.hpp b/src/coreComponents/constitutive/ConstitutivePassThru.hpp index ec72812cdb1..7b5b1533134 100644 --- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp +++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp @@ -101,7 +101,6 @@ struct ConstitutivePassThru< FrictionBase > } }; - /** * Specialization for models that derive from CoulombFriction. */ diff --git a/src/coreComponents/constitutive/contact/RateAndStateFriction.cpp b/src/coreComponents/constitutive/contact/RateAndStateFriction.cpp index 30d43c9ae7f..6376553ec85 100644 --- a/src/coreComponents/constitutive/contact/RateAndStateFriction.cpp +++ b/src/coreComponents/constitutive/contact/RateAndStateFriction.cpp @@ -31,12 +31,15 @@ RateAndStateFriction::RateAndStateFriction( string const & name, Group * const p FrictionBase( name, parent ) { registerWrapper( viewKeyStruct::aCoefficientString(), &m_a ). + setPlotLevel( PlotLevel::LEVEL_0 ). setDescription( "Rate- and State-dependent friction coefficient a." ); registerWrapper( viewKeyStruct::bCoefficientString(), &m_b ). + setPlotLevel( PlotLevel::LEVEL_0 ). setDescription( "Rate- and State-dependent friction coefficient b." ); registerWrapper( viewKeyStruct::DcCoefficientString(), &m_Dc ). + setPlotLevel( PlotLevel::LEVEL_0 ). setDescription( "Rate- and State-dependent friction characteristic length." ); registerWrapper( viewKeyStruct::referenceVelocityString(), &m_V0 ). diff --git a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp index df37934c0cd..2c61762010c 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp +++ b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp @@ -75,7 +75,7 @@ DECLARE_FIELD( dispJump, "Displacement jump vector in the local reference system" ); DECLARE_FIELD( dispJump_n, - "displacementJump", + "displacementJump_n", array2d< real64 >, 0, NOPLOT, @@ -98,6 +98,14 @@ DECLARE_FIELD( deltaSlip, WRITE_AND_READ, "Slip increment" ); +DECLARE_FIELD( deltaSlip_n, + "deltaSlip_n", + array2d< real64 >, + 0.0, + NOPLOT, + WRITE_AND_READ, + "Initial slip increment at this time step" ); + DECLARE_FIELD( deltaDispJump, "deltaDisplacementJump", array2d< real64 >, diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp index 6fecd60e93c..c1c886aa3dd 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp @@ -26,6 +26,7 @@ #include "physicsSolvers/contact/LogLevelsInfo.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "common/GEOS_RAJA_Interface.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" namespace geos { @@ -66,6 +67,7 @@ void ContactSolverBase::registerDataOnMesh( dataRepository::Group & meshBodies ) forFractureRegionOnMeshTargets( meshBodies, [&] ( SurfaceElementRegion & fractureRegion ) { string const labels[3] = { "normal", "tangent1", "tangent2" }; + string const labelsTangent[2] = { "tangent1", "tangent2" }; fractureRegion.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) { @@ -83,6 +85,10 @@ void ContactSolverBase::registerDataOnMesh( dataRepository::Group & meshBodies ) setDimLabels( 1, labels ). reference().resizeDimension< 1 >( 3 ); + subRegion.registerField< fields::contact::dispJump_n >( getName() ). + setDimLabels( 1, labels ). + reference().resizeDimension< 1 >( 3 ); + subRegion.registerField< fields::contact::traction >( getName() ). setDimLabels( 1, labels ). reference().resizeDimension< 1 >( 3 ); @@ -93,7 +99,11 @@ void ContactSolverBase::registerDataOnMesh( dataRepository::Group & meshBodies ) subRegion.registerField< fields::contact::slip >( getName() ); - subRegion.registerField< fields::contact::deltaSlip >( getName() ); + subRegion.registerField< fields::contact::deltaSlip >( getName() ). + setDimLabels( 1, labelsTangent ).reference().resizeDimension< 1 >( 2 ); + + subRegion.registerField< fields::contact::deltaSlip_n >( this->getName() ). + setDimLabels( 1, labelsTangent ).reference().resizeDimension< 1 >( 2 ); } ); } ); diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp index 2d1f139362b..3602616e148 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.cpp @@ -324,10 +324,6 @@ void SolidMechanicsLagrangeContactBubbleStab::assembleSystem( real64 const time, assembleStabilization( dt, domain, dofManager, localMatrix, localRhs ); assembleContact( dt, domain, dofManager, localMatrix, localRhs ); - - // parallel_matrix.create( localMatrix.toViewConst(), dofManager.numLocalDofs(), MPI_COMM_GEOS ); - // parallel_matrix.write("newMatrix.mtx"); - // std::cout << localRhs << std::endl; } void SolidMechanicsLagrangeContactBubbleStab::assembleStabilization( real64 const dt, diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt b/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt index c8f23c55323..df05497f6f2 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt @@ -1,19 +1,27 @@ # Specify solver headers set( physicsSolvers_headers ${physicsSolvers_headers} + inducedSeismicity/ExplicitQDRateAndState.hpp + inducedSeismicity/ImplicitQDRateAndState.hpp inducedSeismicity/inducedSeismicityFields.hpp + inducedSeismicity/QDRateAndStateBase.hpp + inducedSeismicity/QuasiDynamicEarthQuake.hpp inducedSeismicity/rateAndStateFields.hpp - inducedSeismicity/QuasiDynamicEQ.hpp - inducedSeismicity/QuasiDynamicEQRK32.hpp inducedSeismicity/SeismicityRate.hpp - inducedSeismicity/kernels/RateAndStateKernels.hpp - inducedSeismicity/kernels/SeismicityRateKernels.hpp + inducedSeismicity/SpringSlider.hpp + inducedSeismicity/kernels/ExplicitRateAndStateKernels.hpp + inducedSeismicity/kernels/ImplicitRateAndStateKernels.hpp + inducedSeismicity/kernels/RateAndStateKernelsBase.hpp + inducedSeismicity/kernels/SeismicityRateKernels.hpp PARENT_SCOPE ) # Specify solver sources set( physicsSolvers_sources ${physicsSolvers_sources} - inducedSeismicity/QuasiDynamicEQ.cpp - inducedSeismicity/QuasiDynamicEQRK32.cpp - inducedSeismicity/SeismicityRate.cpp + inducedSeismicity/ExplicitQDRateAndState.cpp + inducedSeismicity/ImplicitQDRateAndState.cpp + inducedSeismicity/QDRateAndStateBase.cpp + inducedSeismicity/QuasiDynamicEarthQuake.cpp + inducedSeismicity/SeismicityRate.cpp + inducedSeismicity/SpringSlider.cpp PARENT_SCOPE ) diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp similarity index 51% rename from src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp rename to src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp index 28ff8f5ff14..d6fc3c74871 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.cpp @@ -14,14 +14,14 @@ */ /** - * @file QuasiDynamicEQRK32.cpp + * @file ExplicitQDRateAndState.cpp */ -#include "QuasiDynamicEQRK32.hpp" +#include "ExplicitQDRateAndState.hpp" #include "dataRepository/InputFlags.hpp" #include "mesh/DomainPartition.hpp" -#include "kernels/RateAndStateKernels.hpp" +#include "kernels/ExplicitRateAndStateKernels.hpp" #include "rateAndStateFields.hpp" #include "physicsSolvers/contact/ContactFields.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" @@ -34,47 +34,24 @@ using namespace fields; using namespace constitutive; using namespace rateAndStateKernels; -QuasiDynamicEQRK32::QuasiDynamicEQRK32( const string & name, - Group * const parent ): - PhysicsSolverBase( name, parent ), - m_stressSolver( nullptr ), - m_stressSolverName( "SpringSlider" ), - m_shearImpedance( 0.0 ), +ExplicitQDRateAndState::ExplicitQDRateAndState( const string & name, + Group * const parent ): + QDRateAndStateBase( name, parent ), m_butcherTable( BogackiShampine32Table()), // TODO: The butcher table should be specified in the XML input. m_successfulStep( false ), - m_controller( PIDController( { 1.0/18.0, 1.0/9.0, 1.0/18.0 }, + m_stepUpdateFactor( 1.0 ), + m_controller( PIDController( { 0.6, -0.2, 0.0 }, 1.0e-6, 1.0e-6, 0.81 )) // TODO: The control parameters should be specified in the XML input -{ - this->registerWrapper( viewKeyStruct::shearImpedanceString(), &m_shearImpedance ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "Shear impedance." ); - - this->registerWrapper( viewKeyStruct::stressSolverNameString(), &m_stressSolverName ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Name of solver for computing stress. If empty, the spring-slider model is run." ); -} - -void QuasiDynamicEQRK32::postInputInitialization() -{ - - // Initialize member stress solver as specified in XML input - if( !m_stressSolverName.empty() ) - { - m_stressSolver = &this->getParent().getGroup< PhysicsSolverBase >( m_stressSolverName ); - } - - PhysicsSolverBase::postInputInitialization(); -} +{} -QuasiDynamicEQRK32::~QuasiDynamicEQRK32() +ExplicitQDRateAndState::~ExplicitQDRateAndState() { // TODO Auto-generated destructor stub } - -void QuasiDynamicEQRK32::registerDataOnMesh( Group & meshBodies ) +void ExplicitQDRateAndState::registerDataOnMesh( Group & meshBodies ) { - PhysicsSolverBase::registerDataOnMesh( meshBodies ); + QDRateAndStateBase::registerDataOnMesh( meshBodies ); forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, MeshLevel & mesh, @@ -86,79 +63,21 @@ void QuasiDynamicEQRK32::registerDataOnMesh( Group & meshBodies ) [&]( localIndex const, SurfaceElementSubRegion & subRegion ) { - // Scalar functions on fault - subRegion.registerField< rateAndState::stateVariable >( getName() ); - subRegion.registerField< rateAndState::stateVariable_n >( getName() ); - subRegion.registerField< rateAndState::slipRate >( getName() ); - - // Tangent (2-component) functions on fault - string const labels2Comp[2] = {"tangent1", "tangent2" }; - subRegion.registerField< rateAndState::slipVelocity >( getName() ). - setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); - subRegion.registerField< rateAndState::slipVelocity_n >( getName() ). - setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); - subRegion.registerField< rateAndState::deltaSlip >( getName() ). - setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); - subRegion.registerField< rateAndState::deltaSlip_n >( getName() ). - setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); - // Runge-Kutta stage rates and error integer const numRKComponents = 3; subRegion.registerField< rateAndState::rungeKuttaStageRates >( getName() ).reference().resizeDimension< 1, 2 >( m_butcherTable.numStages, numRKComponents ); subRegion.registerField< rateAndState::error >( getName() ).reference().resizeDimension< 1 >( numRKComponents ); - - - if( !subRegion.hasWrapper( contact::dispJump::key() )) - { - // 3-component functions on fault - string const labels3Comp[3] = { "normal", "tangent1", "tangent2" }; - subRegion.registerField< contact::dispJump >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - subRegion.registerField< contact::dispJump_n >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - subRegion.registerField< contact::traction >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - subRegion.registerField< contact::traction_n >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - - subRegion.registerWrapper< string >( viewKeyStruct::frictionLawNameString() ). - setPlotLevel( PlotLevel::NOPLOT ). - setRestartFlags( RestartFlags::NO_WRITE ). - setSizedFromParent( 0 ); - - string & frictionLawName = subRegion.getReference< string >( viewKeyStruct::frictionLawNameString() ); - frictionLawName = PhysicsSolverBase::getConstitutiveName< FrictionBase >( subRegion ); - GEOS_ERROR_IF( frictionLawName.empty(), GEOS_FMT( "{}: FrictionBase model not found on subregion {}", - getDataContext(), subRegion.getDataContext() ) ); - } } ); } ); } -real64 QuasiDynamicEQRK32::solverStep( real64 const & time_n, - real64 const & dt, - int const cycleNumber, - DomainPartition & domain ) +real64 ExplicitQDRateAndState::solverStep( real64 const & time_n, + real64 const & dt, + int const cycleNumber, + DomainPartition & domain ) { - if( cycleNumber == 0 ) - { - /// Apply initial conditions to the Fault - FieldSpecificationManager & fieldSpecificationManager = FieldSpecificationManager::getInstance(); - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & ) - - { - fieldSpecificationManager.applyInitialConditions( mesh ); - - } ); - saveState( domain ); - } + applyInitialConditionsToFault( cycleNumber, domain ); + saveState( domain ); real64 dtAdaptive = dt; @@ -184,31 +103,29 @@ real64 QuasiDynamicEQRK32::solverStep( real64 const & time_n, stepRateStateODEAndComputeError( dtAdaptive, domain ); // Update timestep based on the time step error - real64 const dtNext = setNextDt( dtAdaptive, domain ); - if( m_successfulStep ) // set in setNextDt + evalTimestep( domain ); + if( m_successfulStep ) // set in evalTimestep { - // Compute stresses, and slip velocity and save results at updated time, + // Compute stresses, and slip velocity and save results at time_n + dtAdapitve if( !m_butcherTable.FSAL ) { dtStress = updateStresses( time_n, dtAdaptive, cycleNumber, domain ); updateSlipVelocity( time_n, dtAdaptive, domain ); } saveState( domain ); - // update the time step and exit the adaptive time step loop - dtAdaptive = dtNext; break; } else { // Retry with updated time step - dtAdaptive = dtNext; + dtAdaptive = setNextDt( dtAdaptive, domain ); } } - // return time step size achieved by stress solver + // return last successful adaptive time step (passed along to setNextDt) return dtAdaptive; } -void QuasiDynamicEQRK32::stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const +void ExplicitQDRateAndState::stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, @@ -247,9 +164,9 @@ void QuasiDynamicEQRK32::stepRateStateODEInitialSubstage( real64 const dt, Domai } ); } -void QuasiDynamicEQRK32::stepRateStateODESubstage( integer const stageIndex, - real64 const dt, - DomainPartition & domain ) const +void ExplicitQDRateAndState::stepRateStateODESubstage( integer const stageIndex, + real64 const dt, + DomainPartition & domain ) const { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, @@ -276,7 +193,7 @@ void QuasiDynamicEQRK32::stepRateStateODESubstage( integer const stageIndex, } ); } -void QuasiDynamicEQRK32::stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const +void ExplicitQDRateAndState::stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, @@ -316,64 +233,9 @@ void QuasiDynamicEQRK32::stepRateStateODEAndComputeError( real64 const dt, Domai } ); } -real64 QuasiDynamicEQRK32::updateStresses( real64 const & time_n, - real64 const & dt, - const int cycleNumber, - DomainPartition & domain ) const -{ - GEOS_LOG_LEVEL_RANK_0( 1, "Stress solver" ); - // Call member variable stress solver to update the stress state - if( m_stressSolver ) - { - // 1. Solve the momentum balance - real64 const dtStress = m_stressSolver->solverStep( time_n, dt, cycleNumber, domain ); - - return dtStress; - } - else - { - // Spring-slider shear traction computation - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - - { - mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, - [&]( localIndex const, - SurfaceElementSubRegion & subRegion ) - { - - arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >(); - arrayView2d< real64 > const traction = subRegion.getField< fields::contact::traction >(); - arrayView2d< real64 const > const traction_n = subRegion.getField< fields::contact::traction_n >(); - - string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); - RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); - - RateAndStateFriction::KernelWrapper frictionKernelWrapper = frictionLaw.createKernelUpdates(); - - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - SpringSliderParameters springSliderParameters = SpringSliderParameters( traction[k][0], - frictionKernelWrapper.getACoefficient( k ), - frictionKernelWrapper.getBCoefficient( k ), - frictionKernelWrapper.getDcCoefficient( k ) ); - - - traction[k][1] = traction_n[k][1] + springSliderParameters.tauRate * dt - - springSliderParameters.springStiffness * deltaSlip[k][0]; - traction[k][2] = traction_n[k][2] + springSliderParameters.tauRate * dt - - springSliderParameters.springStiffness * deltaSlip[k][1]; - } ); - } ); - } ); - return dt; - } -} - -void QuasiDynamicEQRK32::updateSlipVelocity( real64 const & time_n, - real64 const & dt, - DomainPartition & domain ) const +void ExplicitQDRateAndState::updateSlipVelocity( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) const { GEOS_LOG_LEVEL_RANK_0( 1, "Rate and State solver" ); integer const maxIterNewton = m_nonlinearSolverParameters.m_maxIterNewton; @@ -394,45 +256,8 @@ void QuasiDynamicEQRK32::updateSlipVelocity( real64 const & time_n, } ); } -void QuasiDynamicEQRK32::saveState( DomainPartition & domain ) const +void ExplicitQDRateAndState::evalTimestep( DomainPartition & domain ) { - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - - { - mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, - [&]( localIndex const, - SurfaceElementSubRegion & subRegion ) - { - arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); - arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); - arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >(); - arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); - arrayView2d< real64 const > const traction = subRegion.getField< contact::traction >(); - - arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >(); - arrayView2d< real64 > const slipVelocity_n = subRegion.getField< rateAndState::slipVelocity_n >(); - arrayView2d< real64 > const deltaSlip_n = subRegion.getField< rateAndState::deltaSlip >(); - arrayView2d< real64 > const dispJump_n = subRegion.getField< contact::dispJump_n >(); - arrayView2d< real64 > const traction_n = subRegion.getField< contact::traction_n >(); - - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - stateVariable_n[k] = stateVariable[k]; - LvArray::tensorOps::copy< 2 >( deltaSlip_n[k], deltaSlip[k] ); - LvArray::tensorOps::copy< 2 >( slipVelocity_n[k], slipVelocity[k] ); - LvArray::tensorOps::copy< 3 >( dispJump_n[k], dispJump[k] ); - LvArray::tensorOps::copy< 3 >( traction_n[k], traction[k] ); - } ); - } ); - } ); -} - -real64 QuasiDynamicEQRK32::setNextDt( real64 const & currentDt, DomainPartition & domain ) -{ - - // Spring-slider shear traction computation forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, arrayView1d< string const > const & regionNames ) @@ -455,24 +280,29 @@ real64 QuasiDynamicEQRK32::setNextDt( real64 const & currentDt, DomainPartition } ); // Compute update factor to currentDt using PID error controller + limiter - real64 const dtFactor = m_controller.computeUpdateFactor( m_butcherTable.algHighOrder, m_butcherTable.algLowOrder ); - real64 const nextDt = dtFactor*currentDt; + m_stepUpdateFactor = m_controller.computeUpdateFactor( m_butcherTable.algHighOrder, m_butcherTable.algLowOrder ); // Check if step was acceptable - m_successfulStep = (dtFactor >= m_controller.acceptSafety) ? true : false; + m_successfulStep = (m_stepUpdateFactor >= m_controller.acceptSafety) ? true : false; if( m_successfulStep ) { m_controller.errors[2] = m_controller.errors[1]; m_controller.errors[1] = m_controller.errors[0]; + } +} + +real64 ExplicitQDRateAndState::setNextDt( real64 const & currentDt, DomainPartition & domain ) +{ + GEOS_UNUSED_VAR( domain ); + real64 const nextDt = m_stepUpdateFactor*currentDt; + if( m_successfulStep ) + { GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "Adaptive time step successful. The next dt will be {:.2e} s", nextDt )); } else { - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "Adaptive time step failed. The next dt will be {:.2e} s", nextDt )); + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "Adaptive time step failed. Retry step with dt {:.2e} s", nextDt )); } - return nextDt; } -REGISTER_CATALOG_ENTRY( PhysicsSolverBase, QuasiDynamicEQRK32, string const &, dataRepository::Group * const ) - } // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp similarity index 67% rename from src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp rename to src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp index bcb3dd8a72a..69769af9381 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ExplicitQDRateAndState.hpp @@ -13,46 +13,35 @@ * ------------------------------------------------------------------------------------------------------------ */ -#ifndef GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP -#define GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP +#ifndef GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_QUASIDYNAMICEQRK32_HPP +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_QUASIDYNAMICEQRK32_HPP -#include "physicsSolvers/PhysicsSolverBase.hpp" -#include "kernels/RateAndStateKernels.hpp" +#include "QDRateAndStateBase.hpp" +#include "kernels/ExplicitRateAndStateKernels.hpp" namespace geos { -class QuasiDynamicEQRK32 : public PhysicsSolverBase +class ExplicitQDRateAndState : public QDRateAndStateBase { public: /// The default nullary constructor is disabled to avoid compiler auto-generation: - QuasiDynamicEQRK32() = delete; + ExplicitQDRateAndState() = delete; /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) - QuasiDynamicEQRK32( const string & name, - Group * const parent ); + ExplicitQDRateAndState( const string & name, + Group * const parent ); /// Destructor - virtual ~QuasiDynamicEQRK32() override; + virtual ~ExplicitQDRateAndState() override; - static string catalogName() { return "QuasiDynamicEQRK32"; } - - /** - * @return Get the final class Catalog name - */ - virtual string getCatalogName() const override { return catalogName(); } + static string derivedSolverPrefix() { return "Explicit";}; /// This method ties properties with their supporting mesh virtual void registerDataOnMesh( Group & meshBodies ) override; - struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct + struct viewKeyStruct : public QDRateAndStateBase::viewKeyStruct { - /// stress solver name - constexpr static char const * stressSolverNameString() { return "stressSolverName"; } - /// Friction law name string - constexpr static char const * frictionLawNameString() { return "frictionLawName"; } - /// Friction law name string - constexpr static char const * shearImpedanceString() { return "shearImpedance"; } /// target slip increment constexpr static char const * timeStepTol() { return "timeStepTol"; } }; @@ -66,6 +55,12 @@ class QuasiDynamicEQRK32 : public PhysicsSolverBase virtual real64 setNextDt( real64 const & currentDt, DomainPartition & domain ) override final; + /** + * @brief Evaluates whether an adaptive time step was successful + * @param domain + */ + void evalTimestep( DomainPartition & domain ); + /** * @brief Computes stage rates for the initial Runge-Kutta substage and updates slip and state * @param dt @@ -90,11 +85,6 @@ class QuasiDynamicEQRK32 : public PhysicsSolverBase */ void stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const; - real64 updateStresses( real64 const & time_n, - real64 const & dt, - const int cycleNumber, - DomainPartition & domain ) const; - /** * @brief Updates rate-and-state slip velocity * @param domain @@ -103,25 +93,7 @@ class QuasiDynamicEQRK32 : public PhysicsSolverBase real64 const & dt, DomainPartition & domain ) const; - /** - * @brief save the current state - * @param domain - */ - void saveState( DomainPartition & domain ) const; - -private: - - virtual void postInputInitialization() override; - - - /// pointer to stress solver - PhysicsSolverBase * m_stressSolver; - - /// stress solver name - string m_stressSolverName; - - /// shear impedance - real64 m_shearImpedance; +protected: /// Runge-Kutta Butcher table (specifies the embedded RK method) // TODO: The specific type should not be hardcoded! @@ -130,6 +102,8 @@ class QuasiDynamicEQRK32 : public PhysicsSolverBase bool m_successfulStep; // Flag indicating if the adative time step was accepted + real64 m_stepUpdateFactor; // Factor to update timestep with + /** * @brief Proportional-integral-derivative controller used for updating time step * based error estimate in the current and previous time steps. @@ -195,41 +169,8 @@ class QuasiDynamicEQRK32 : public PhysicsSolverBase PIDController m_controller; - - class SpringSliderParameters - { -public: - - GEOS_HOST_DEVICE - SpringSliderParameters( real64 const normalTraction, real64 const a, real64 const b, real64 const Dc ): - tauRate( 1e-4 ), - springStiffness( 0.0 ) - { - real64 const criticalStiffness = normalTraction * (b - a) / Dc; - springStiffness = 0.9 * criticalStiffness; - } - - /// Default copy constructor - SpringSliderParameters( SpringSliderParameters const & ) = default; - - /// Default move constructor - SpringSliderParameters( SpringSliderParameters && ) = default; - - /// Deleted default constructor - SpringSliderParameters() = delete; - - /// Deleted copy assignment operator - SpringSliderParameters & operator=( SpringSliderParameters const & ) = delete; - - /// Deleted move assignment operator - SpringSliderParameters & operator=( SpringSliderParameters && ) = delete; - - real64 tauRate; - - real64 springStiffness; - }; }; } /* namespace geos */ -#endif /* GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP */ +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_QUASIDYNAMICEQRK32_HPP */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp new file mode 100644 index 00000000000..733fbe4ab7a --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.cpp @@ -0,0 +1,142 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ImplicitQDRateAndState.cpp + */ + +#include "ImplicitQDRateAndState.hpp" + +#include "dataRepository/InputFlags.hpp" +#include "mesh/DomainPartition.hpp" +#include "kernels/ImplicitRateAndStateKernels.hpp" +#include "rateAndStateFields.hpp" +#include "physicsSolvers/contact/ContactFields.hpp" + + +namespace geos +{ + +using namespace dataRepository; +using namespace fields; +using namespace constitutive; + +ImplicitQDRateAndState::ImplicitQDRateAndState( const string & name, + Group * const parent ): + QDRateAndStateBase( name, parent ), + m_targetSlipIncrement( 1.0e-7 ) +{ + this->registerWrapper( viewKeyStruct::targetSlipIncrementString(), &m_targetSlipIncrement ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1.0e-7 ). + setDescription( "Target slip incrmeent for timestep size selction" ); +} + +ImplicitQDRateAndState::~ImplicitQDRateAndState() +{ + // TODO Auto-generated destructor stub +} + +void ImplicitQDRateAndState::solveRateAndStateEquations( real64 const time_n, + real64 const dt, + DomainPartition & domain ) const +{ + integer const maxNewtonIter = m_nonlinearSolverParameters.m_maxIterNewton; + real64 const newtonTol = m_nonlinearSolverParameters.m_newtonTol; + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + // solve rate and state equations. + rateAndStateKernels::createAndLaunch< rateAndStateKernels::ImplicitFixedStressRateAndStateKernel, parallelDevicePolicy<> >( subRegion, + viewKeyStruct::frictionLawNameString(), + m_shearImpedance, + maxNewtonIter, newtonTol, + time_n, + dt ); // save + // old + // state + updateSlip( subRegion, dt ); + } ); + } ); +} + +real64 ImplicitQDRateAndState::solverStep( real64 const & time_n, + real64 const & dt, + int const cycleNumber, + DomainPartition & domain ) +{ + applyInitialConditionsToFault( cycleNumber, domain ); + GEOS_LOG_LEVEL_RANK_0( 1, "Stress solver" ); + updateStresses( time_n, dt, cycleNumber, domain ); + GEOS_LOG_LEVEL_RANK_0( 1, "Rate and state solver" ); + solveRateAndStateEquations( time_n, dt, domain ); + saveState( domain ); + return dt; +} + +void ImplicitQDRateAndState::updateSlip( ElementSubRegionBase & subRegion, real64 const dt ) const +{ + arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); + arrayView2d< real64 > const deltaSlip = subRegion.getField< contact::deltaSlip >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + deltaSlip[k][0] = slipVelocity[k][0] * dt; + deltaSlip[k][1] = slipVelocity[k][1] * dt; + } ); +} + +real64 ImplicitQDRateAndState::setNextDt( real64 const & currentDt, DomainPartition & domain ) +{ + GEOS_UNUSED_VAR( currentDt ); + + real64 maxSlipRate = 0.0; + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & regionNames ) + + { + real64 maxSlipRateOnThisRank = 0.0; + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion const & subRegion ) + { + arrayView1d< real64 const > const slipRate = subRegion.getField< rateAndState::slipRate >(); + + RAJA::ReduceMax< parallelDeviceReduce, real64 > maximumSlipRateOnThisRegion( 0.0 ); + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + maximumSlipRateOnThisRegion.max( slipRate[k] ); + } ); + if( maximumSlipRateOnThisRegion.get() > maxSlipRateOnThisRank ) + maxSlipRateOnThisRank = maximumSlipRateOnThisRegion.get(); + } ); + maxSlipRate = MpiWrapper::max( maxSlipRateOnThisRank ); + } ); + + real64 const nextDt = m_targetSlipIncrement / maxSlipRate; + + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "The next dt will be {:.2e} s", nextDt )); + + return nextDt; +} + +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp new file mode 100644 index 00000000000..907612d0aef --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp @@ -0,0 +1,70 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_IMPLICITQDRATEANDSTATE_HPP +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_IMPLICITQDRATEANDSTATE_HPP + +#include "QDRateAndStateBase.hpp" + +namespace geos +{ + +class ImplicitQDRateAndState : public QDRateAndStateBase +{ +public: + /// The default nullary constructor is disabled to avoid compiler auto-generation: + ImplicitQDRateAndState() = delete; + + /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) + ImplicitQDRateAndState( const string & name, + Group * const parent ); + + /// Destructor + virtual ~ImplicitQDRateAndState() override; + + static string derivedSolverPrefix() { return "Implicit";}; + + struct viewKeyStruct : public QDRateAndStateBase::viewKeyStruct + { + /// target slip increment + constexpr static char const * targetSlipIncrementString() { return "targetSlipIncrement"; } + }; + + virtual real64 setNextDt( real64 const & currentDt, + DomainPartition & domain ) override final; + + /** + * @brief save the old state + * @param subRegion + */ + void updateSlip( ElementSubRegionBase & subRegion, real64 const dt ) const; + + virtual real64 solverStep( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) override final; +protected: + + void solveRateAndStateEquations( real64 const time_n, + real64 const dt, + DomainPartition & domain ) const; + + /// target slip rate + real64 m_targetSlipIncrement; +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_IMPLICITQDRATEANDSTATE_HPP */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QDRateAndStateBase.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QDRateAndStateBase.cpp new file mode 100644 index 00000000000..a515685cab9 --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QDRateAndStateBase.cpp @@ -0,0 +1,237 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file QDRateAndStateBase.cpp + */ + +#include "QDRateAndStateBase.hpp" + +#include "dataRepository/InputFlags.hpp" +#include "mesh/DomainPartition.hpp" +#include "rateAndStateFields.hpp" +#include "physicsSolvers/contact/ContactFields.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" +#include "constitutive/contact/RateAndStateFriction.hpp" +#include "kernels/RateAndStateKernelsBase.hpp" + +namespace geos +{ + +using namespace dataRepository; +using namespace fields; +using namespace constitutive; + +QDRateAndStateBase::QDRateAndStateBase( const string & name, + Group * const parent ): + PhysicsSolverBase( name, parent ), + m_shearImpedance( 0.0 ) +{ + this->registerWrapper( viewKeyStruct::shearImpedanceString(), &m_shearImpedance ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Shear impedance." ); +} + +QDRateAndStateBase::~QDRateAndStateBase() +{ + // TODO Auto-generated destructor stub +} + +void QDRateAndStateBase::registerDataOnMesh( Group & meshBodies ) +{ + PhysicsSolverBase::registerDataOnMesh( meshBodies ); + + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + // Scalar functions on fault + subRegion.registerField< rateAndState::stateVariable >( getName() ); + subRegion.registerField< rateAndState::stateVariable_n >( getName() ); + subRegion.registerField< rateAndState::slipRate >( getName() ); + subRegion.registerField< rateAndState::slipRate_n >( getName() ); + + // Tangent (2-component) functions on fault + string const labels2Comp[2] = {"tangent1", "tangent2" }; + + subRegion.registerField< rateAndState::slipVelocity >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + subRegion.registerField< rateAndState::slipVelocity_n >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + + subRegion.registerField< rateAndState::shearTraction >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + subRegion.registerField< rateAndState::normalTraction >( getName() ); + + subRegion.registerField< rateAndState::shearTraction_n >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + subRegion.registerField< rateAndState::normalTraction_n >( getName() ); + + subRegion.registerField< rateAndState::backgroundShearStress >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + subRegion.registerField< rateAndState::backgroundNormalStress >( getName() ); + } ); + } ); +} + +void QDRateAndStateBase::enforceRateAndVelocityConsistency( SurfaceElementSubRegion & subRegion ) const +{ + arrayView2d< real64 > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); + arrayView1d< real64 > const slipRate = subRegion.getField< rateAndState::slipRate >(); + arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); + + arrayView2d< real64 > const backgroundShearStress = subRegion.getField< rateAndState::backgroundShearStress >(); + arrayView1d< real64 > const backgroundNormalStress = subRegion.getField< rateAndState::backgroundNormalStress >(); + + real64 const shearImpedance = m_shearImpedance; + + string const & frictionaLawName = subRegion.getReference< string >( viewKeyStruct::frictionLawNameString() ); + constitutive::RateAndStateFriction const & frictionLaw = subRegion.getConstitutiveModel< constitutive::RateAndStateFriction >( frictionaLawName ); + + constitutive::RateAndStateFriction::KernelWrapper frictionLawKernelWrapper = frictionLaw.createKernelUpdates(); + + RAJA::ReduceMax< parallelDeviceReduce, int > negativeSlipRate( 0 ); + RAJA::ReduceMax< parallelDeviceReduce, int > bothNonZero( 0 ); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + if( slipRate[k] < 0.0 ) + { + negativeSlipRate.max( 1 ); + } + else if( LvArray::tensorOps::l2Norm< 2 >( slipVelocity[k] ) > 0.0 && slipRate[k] > 0.0 ) + { + bothNonZero.max( 1 ); + } + else if( LvArray::tensorOps::l2Norm< 2 >( slipVelocity[k] ) > 0.0 ) + { + slipRate[k] = LvArray::tensorOps::l2Norm< 2 >( slipVelocity[k] ); + } + else if( slipRate[k] > 0.0 ) + { + real64 const frictionCoefficient = frictionLawKernelWrapper.frictionCoefficient( k, slipRate[k], stateVariable[k] ); + rateAndStateKernels::projectSlipRateBase( k, + frictionCoefficient, + shearImpedance, + backgroundNormalStress, + backgroundShearStress, + slipRate, + slipVelocity ); + } + } ); + + + GEOS_ERROR_IF( negativeSlipRate.get() > 0, "SlipRate cannot be negative." ); + GEOS_ERROR_IF( bothNonZero.get() > 0, "Only one between slipRate and slipVelocity can be specified as i.c." ); + +} + +void QDRateAndStateBase::applyInitialConditionsToFault( int const cycleNumber, + DomainPartition & domain ) const +{ + if( cycleNumber == 0 ) + { + /// Apply initial conditions to the Fault + FieldSpecificationManager & fieldSpecificationManager = FieldSpecificationManager::getInstance(); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + fieldSpecificationManager.applyInitialConditions( mesh ); + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + enforceRateAndVelocityConsistency( subRegion ); + + arrayView1d< real64 const > const slipRate = subRegion.getField< rateAndState::slipRate >(); + arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); + + arrayView1d< real64 > const normalTraction = subRegion.getField< rateAndState::normalTraction >(); + arrayView2d< real64 > const shearTraction = subRegion.getField< rateAndState::shearTraction >(); + + arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >(); + arrayView1d< real64 > const slipRate_n = subRegion.getField< rateAndState::slipRate_n >(); + arrayView1d< real64 > const normalTraction_n = subRegion.getField< rateAndState::normalTraction_n >(); + arrayView2d< real64 > const shearTraction_n = subRegion.getField< rateAndState::shearTraction_n >(); + + arrayView2d< real64 const > const backgroundShearStress = subRegion.getField< rateAndState::backgroundShearStress >(); + arrayView1d< real64 const > const backgroundNormalStress = subRegion.getField< rateAndState::backgroundNormalStress >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + slipRate_n [k] = slipRate[k]; + stateVariable_n[k] = stateVariable[k]; + + normalTraction[k] = backgroundNormalStress[k]; + LvArray::tensorOps::copy< 2 >( shearTraction[k], backgroundShearStress[k] ); + + normalTraction_n[k] = normalTraction[k]; + LvArray::tensorOps::copy< 2 >( shearTraction_n[k], shearTraction[k] ); + } ); + } ); + } ); + } +} + +void QDRateAndStateBase::saveState( DomainPartition & domain ) const +{ + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); + arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); + arrayView2d< real64 const > const deltaSlip = subRegion.getField< contact::deltaSlip >(); + arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + arrayView2d< real64 const > const shearTraction = subRegion.getField< rateAndState::shearTraction >(); + arrayView1d< real64 const > const normalTraction = subRegion.getField< rateAndState::normalTraction >(); + + + arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >(); + arrayView2d< real64 > const slipVelocity_n = subRegion.getField< rateAndState::slipVelocity_n >(); + arrayView2d< real64 > const deltaSlip_n = subRegion.getField< contact::deltaSlip >(); + arrayView2d< real64 > const dispJump_n = subRegion.getField< contact::dispJump_n >(); + arrayView2d< real64 > const shearTraction_n = subRegion.getField< rateAndState::shearTraction_n >(); + arrayView1d< real64 > const normalTraction_n = subRegion.getField< rateAndState::normalTraction_n >(); + + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + stateVariable_n[k] = stateVariable[k]; + normalTraction_n[k] = normalTraction[k]; + LvArray::tensorOps::copy< 2 >( deltaSlip_n[k], deltaSlip[k] ); + LvArray::tensorOps::copy< 2 >( slipVelocity_n[k], slipVelocity[k] ); + LvArray::tensorOps::copy< 3 >( dispJump_n[k], dispJump[k] ); + LvArray::tensorOps::copy< 2 >( shearTraction_n[k], shearTraction[k] ); + } ); + } ); + } ); +} + +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QDRateAndStateBase.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QDRateAndStateBase.hpp new file mode 100644 index 00000000000..17b4ecb0685 --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QDRateAndStateBase.hpp @@ -0,0 +1,89 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_QDRATEANDSTATEBASE_HPP +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_QDRATEANDSTATEBASE_HPP + +#include "physicsSolvers/PhysicsSolverBase.hpp" + +namespace geos +{ + +class QDRateAndStateBase : public PhysicsSolverBase +{ +public: + /// The default nullary constructor is disabled to avoid compiler auto-generation: + QDRateAndStateBase() = delete; + + /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) + QDRateAndStateBase( const string & name, + Group * const parent ); + + /// Destructor + virtual ~QDRateAndStateBase() override; + + /// This method ties properties with their supporting mesh + virtual void registerDataOnMesh( Group & meshBodies ) override; + + struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct + { + /// Friction law name string + constexpr static char const * frictionLawNameString() { return "frictionLawName"; } + /// Friction law name string + constexpr static char const * shearImpedanceString() { return "shearImpedance"; } + }; + + /** + * @brief Save the current state of the solver fields + * @param domain the domain object + */ + void saveState( DomainPartition & domain ) const; + + /** + * @brief Check that only one of slip rate or slip velocity are specified as initial conditions + * and initialize the unspecified field + * @param subRegion the element subregion + */ + void enforceRateAndVelocityConsistency( SurfaceElementSubRegion & subRegion ) const; + + /** + * @brief Compute stresses and update tractions on the fault + * @param time_n the current time + * @param dt the time step + * @param cycleNumber the current cycle number + * @param domain the domain object + */ + virtual real64 updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const = 0; + + /** + * @brief Apply initial conditions to fields on the fault + * @param cycleNumber the current cycle number + * @param domain the domain object + */ + virtual void applyInitialConditionsToFault( int const cycleNumber, + DomainPartition & domain ) const; + +protected: + + /// shear impedance + real64 m_shearImpedance; +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICTY_QDRATEANDSTATEBASE_HPP */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp deleted file mode 100644 index 3ddf5fc1494..00000000000 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp +++ /dev/null @@ -1,287 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2024 TotalEnergies - * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2023-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file QuasiDynamicEQ.cpp - */ - -#include "QuasiDynamicEQ.hpp" - -#include "dataRepository/InputFlags.hpp" -#include "mesh/DomainPartition.hpp" -#include "kernels/RateAndStateKernels.hpp" -#include "rateAndStateFields.hpp" -#include "physicsSolvers/contact/ContactFields.hpp" -#include "fieldSpecification/FieldSpecificationManager.hpp" - -namespace geos -{ - -using namespace dataRepository; -using namespace fields; -using namespace constitutive; -using namespace rateAndStateKernels; - -QuasiDynamicEQ::QuasiDynamicEQ( const string & name, - Group * const parent ): - PhysicsSolverBase( name, parent ), - m_stressSolver( nullptr ), - m_stressSolverName( "SpringSlider" ), - m_shearImpedance( 0.0 ), - m_targetSlipIncrement( 1.0e-7 ) -{ - this->registerWrapper( viewKeyStruct::shearImpedanceString(), &m_shearImpedance ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "Shear impedance." ); - - this->registerWrapper( viewKeyStruct::stressSolverNameString(), &m_stressSolverName ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Name of solver for computing stress. If empty, the spring-slider model is run." ); - - this->registerWrapper( viewKeyStruct::targetSlipIncrementString(), &m_targetSlipIncrement ). - setInputFlag( InputFlags::OPTIONAL ). - setApplyDefaultValue( 1.0e-7 ). - setDescription( "Target slip incrmeent for timestep size selction" ); -} - -void QuasiDynamicEQ::postInputInitialization() -{ - - // Initialize member stress solver as specified in XML input - if( !m_stressSolverName.empty() ) - { - m_stressSolver = &this->getParent().getGroup< PhysicsSolverBase >( m_stressSolverName ); - } - - PhysicsSolverBase::postInputInitialization(); -} - -QuasiDynamicEQ::~QuasiDynamicEQ() -{ - // TODO Auto-generated destructor stub -} - -void QuasiDynamicEQ::registerDataOnMesh( Group & meshBodies ) -{ - PhysicsSolverBase::registerDataOnMesh( meshBodies ); - - forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - ElementRegionManager & elemManager = mesh.getElemManager(); - - elemManager.forElementSubRegions< SurfaceElementSubRegion >( regionNames, - [&]( localIndex const, - SurfaceElementSubRegion & subRegion ) - { - // Scalar functions on fault - subRegion.registerField< rateAndState::stateVariable >( getName() ); - subRegion.registerField< rateAndState::stateVariable_n >( getName() ); - subRegion.registerField< rateAndState::slipRate >( getName() ); - - // Tangent (2-component) functions on fault - string const labels2Comp[2] = {"tangent1", "tangent2" }; - subRegion.registerField< rateAndState::slipVelocity >( getName() ). - setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); - subRegion.registerField< rateAndState::deltaSlip >( getName() ). - setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); - - if( !subRegion.hasWrapper( contact::dispJump::key() )) - { - // 3-component functions on fault - string const labels3Comp[3] = { "normal", "tangent1", "tangent2" }; - subRegion.registerField< contact::dispJump >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - subRegion.registerField< contact::traction >( getName() ). - setDimLabels( 1, labels3Comp ). - reference().resizeDimension< 1 >( 3 ); - - subRegion.registerWrapper< string >( viewKeyStruct::frictionLawNameString() ). - setPlotLevel( PlotLevel::NOPLOT ). - setRestartFlags( RestartFlags::NO_WRITE ). - setSizedFromParent( 0 ); - - string & frictionLawName = subRegion.getReference< string >( viewKeyStruct::frictionLawNameString() ); - frictionLawName =PhysicsSolverBase::getConstitutiveName< FrictionBase >( subRegion ); - GEOS_ERROR_IF( frictionLawName.empty(), GEOS_FMT( "{}: FrictionBase model not found on subregion {}", - getDataContext(), subRegion.getDataContext() ) ); - } - } ); - } ); -} - -real64 QuasiDynamicEQ::solverStep( real64 const & time_n, - real64 const & dt, - int const cycleNumber, - DomainPartition & domain ) -{ - if( cycleNumber == 0 ) - { - /// Apply initial conditions to the Fault - FieldSpecificationManager & fieldSpecificationManager = FieldSpecificationManager::getInstance(); - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & ) - - { - fieldSpecificationManager.applyInitialConditions( mesh ); - } ); - } - - /// 1. Compute shear and normal tractions - GEOS_LOG_LEVEL_RANK_0( 1, "Stress solver" ); - - real64 const dtStress = updateStresses( time_n, dt, cycleNumber, domain ); - - /// 2. Solve for slip rate and state variable and, compute slip - GEOS_LOG_LEVEL_RANK_0( 1, "Rate and State solver" ); - integer const maxIterNewton = m_nonlinearSolverParameters.m_maxIterNewton; - real64 const newtonTol = m_nonlinearSolverParameters.m_newtonTol; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - - { - mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, - [&]( localIndex const, - SurfaceElementSubRegion & subRegion ) - { - // solve rate and state equations. - createAndLaunch< ImplicitFixedStressRateAndStateKernel, parallelDevicePolicy<> >( subRegion, viewKeyStruct::frictionLawNameString(), m_shearImpedance, maxIterNewton, newtonTol, time_n, - dtStress ); - // save old state - saveOldStateAndUpdateSlip( subRegion, dtStress ); - } ); - } ); - - // return time step size achieved by stress solver - return dtStress; -} - -real64 QuasiDynamicEQ::updateStresses( real64 const & time_n, - real64 const & dt, - const int cycleNumber, - DomainPartition & domain ) const -{ - // Call member variable stress solver to update the stress state - if( m_stressSolver ) - { - // 1. Solve the momentum balance - real64 const dtStress = m_stressSolver->solverStep( time_n, dt, cycleNumber, domain ); - - return dtStress; - } - else - { - // Spring-slider shear traction computation - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - - { - mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, - [&]( localIndex const, - SurfaceElementSubRegion & subRegion ) - { - - arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >(); - arrayView2d< real64 > const traction = subRegion.getField< fields::contact::traction >(); - - string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); - RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); - - RateAndStateFriction::KernelWrapper frictionKernelWrapper = frictionLaw.createKernelUpdates(); - - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - SpringSliderParameters springSliderParameters = SpringSliderParameters( traction[k][0], - frictionKernelWrapper.getACoefficient( k ), - frictionKernelWrapper.getBCoefficient( k ), - frictionKernelWrapper.getDcCoefficient( k ) ); - - - traction[k][1] = traction[k][1] + springSliderParameters.tauRate * dt - - springSliderParameters.springStiffness * deltaSlip[k][0]; - traction[k][2] = traction[k][2] + springSliderParameters.tauRate * dt - - springSliderParameters.springStiffness * deltaSlip[k][1]; - } ); - } ); - } ); - return dt; - } -} - -void QuasiDynamicEQ::saveOldStateAndUpdateSlip( ElementSubRegionBase & subRegion, real64 const dt ) const -{ - arrayView1d< real64 > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); - arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >(); - arrayView2d< real64 > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); - arrayView2d< real64 > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >(); - - arrayView2d< real64 > const dispJump = subRegion.getField< contact::dispJump >(); - - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - stateVariable_n[k] = stateVariable[k]; - deltaSlip[k][0] = slipVelocity[k][0] * dt; - deltaSlip[k][1] = slipVelocity[k][1] * dt; - // Update tangential components of the displacement jump - dispJump[k][1] = dispJump[k][1] + slipVelocity[k][0] * dt; - dispJump[k][2] = dispJump[k][2] + slipVelocity[k][1] * dt; - } ); -} - -real64 QuasiDynamicEQ::setNextDt( real64 const & currentDt, DomainPartition & domain ) -{ - GEOS_UNUSED_VAR( currentDt ); - - real64 maxSlipRate = 0.0; - // Spring-slider shear traction computation - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel const & mesh, - arrayView1d< string const > const & regionNames ) - - { - real64 maxSlipRateOnThisRank = 0.0; - mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, - [&]( localIndex const, - SurfaceElementSubRegion const & subRegion ) - { - arrayView1d< real64 const > const slipRate = subRegion.getField< rateAndState::slipRate >(); - - RAJA::ReduceMax< parallelDeviceReduce, real64 > maximumSlipRateOnThisRegion( 0.0 ); - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - maximumSlipRateOnThisRegion.max( slipRate[k] ); - } ); - if( maximumSlipRateOnThisRegion.get() > maxSlipRateOnThisRank ) - maxSlipRateOnThisRank = maximumSlipRateOnThisRegion.get(); - } ); - maxSlipRate = MpiWrapper::max( maxSlipRateOnThisRank ); - } ); - - real64 const nextDt = m_targetSlipIncrement / maxSlipRate; - - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "The next dt will be {:.2e} s", nextDt )); - - return nextDt; -} - -REGISTER_CATALOG_ENTRY( PhysicsSolverBase, QuasiDynamicEQ, string const &, dataRepository::Group * const ) - -} // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp deleted file mode 100644 index 6d85175d3c7..00000000000 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp +++ /dev/null @@ -1,135 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2024 TotalEnergies - * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2023-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -#ifndef GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQ_HPP -#define GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQ_HPP - -#include "physicsSolvers/PhysicsSolverBase.hpp" - -namespace geos -{ - -class QuasiDynamicEQ : public PhysicsSolverBase -{ -public: - /// The default nullary constructor is disabled to avoid compiler auto-generation: - QuasiDynamicEQ() = delete; - - /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) - QuasiDynamicEQ( const string & name, - Group * const parent ); - - /// Destructor - virtual ~QuasiDynamicEQ() override; - - static string catalogName() { return "QuasiDynamicEQ"; } - - /** - * @return Get the final class Catalog name - */ - virtual string getCatalogName() const override { return catalogName(); } - - /// This method ties properties with their supporting mesh - virtual void registerDataOnMesh( Group & meshBodies ) override; - - struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct - { - /// stress solver name - constexpr static char const * stressSolverNameString() { return "stressSolverName"; } - /// Friction law name string - constexpr static char const * frictionLawNameString() { return "frictionLawName"; } - /// Friction law name string - constexpr static char const * shearImpedanceString() { return "shearImpedance"; } - /// target slip increment - constexpr static char const * targetSlipIncrementString() { return "targetSlipIncrement"; } - }; - - virtual real64 solverStep( real64 const & time_n, - real64 const & dt, - integer const cycleNumber, - DomainPartition & domain ) override final; - - virtual real64 setNextDt( real64 const & currentDt, - DomainPartition & domain ) override final; - - real64 updateStresses( real64 const & time_n, - real64 const & dt, - const int cycleNumber, - DomainPartition & domain ) const; - - /** - * @brief save the old state - * @param subRegion - */ - void saveOldStateAndUpdateSlip( ElementSubRegionBase & subRegion, real64 const dt ) const; - - -private: - - - - virtual void postInputInitialization() override; - - - - /// pointer to stress solver - PhysicsSolverBase * m_stressSolver; - - /// stress solver name - string m_stressSolverName; - - /// shear impedance - real64 m_shearImpedance; - - /// target slip rate - real64 m_targetSlipIncrement; - - class SpringSliderParameters - { -public: - - GEOS_HOST_DEVICE - SpringSliderParameters( real64 const normalTraction, real64 const a, real64 const b, real64 const Dc ): - tauRate( 1e-4 ), - springStiffness( 0.0 ) - { - real64 const criticalStiffness = normalTraction * (b - a) / Dc; - springStiffness = 0.9 * criticalStiffness; - } - - /// Default copy constructor - SpringSliderParameters( SpringSliderParameters const & ) = default; - - /// Default move constructor - SpringSliderParameters( SpringSliderParameters && ) = default; - - /// Deleted default constructor - SpringSliderParameters() = delete; - - /// Deleted copy assignment operator - SpringSliderParameters & operator=( SpringSliderParameters const & ) = delete; - - /// Deleted move assignment operator - SpringSliderParameters & operator=( SpringSliderParameters && ) = delete; - - real64 tauRate; - - real64 springStiffness; - }; -}; - -} /* namespace geos */ - -#endif /* GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQ_HPP */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp new file mode 100644 index 00000000000..7c45b81a7a1 --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.cpp @@ -0,0 +1,145 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file QuasiDynamicEarthQuake.cpp + */ + +#include "QuasiDynamicEarthQuake.hpp" + +#include "dataRepository/InputFlags.hpp" +#include "mesh/DomainPartition.hpp" +#include "rateAndStateFields.hpp" +#include "physicsSolvers/contact/ContactFields.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" + +#include "ExplicitQDRateAndState.hpp" + +namespace geos +{ + +using namespace dataRepository; +using namespace fields; +using namespace constitutive; +using namespace rateAndStateKernels; + +template< typename RSSOLVER_TYPE > +QuasiDynamicEarthQuake< RSSOLVER_TYPE >::QuasiDynamicEarthQuake( const string & name, + Group * const parent ): + RSSOLVER_TYPE( name, parent ), + m_stressSolverName(), + m_stressSolver( nullptr ) +{ + this->registerWrapper( viewKeyStruct::stressSolverNameString(), &m_stressSolverName ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of solver for computing stress." ); +} + +template< typename RSSOLVER_TYPE > +void QuasiDynamicEarthQuake< RSSOLVER_TYPE >::postInputInitialization() +{ + + // Initialize member stress solver as specified in XML input + m_stressSolver = &this->getParent().template getGroup< PhysicsSolverBase >( m_stressSolverName ); + + PhysicsSolverBase::postInputInitialization(); +} + +template< typename RSSOLVER_TYPE > +QuasiDynamicEarthQuake< RSSOLVER_TYPE >::~QuasiDynamicEarthQuake() +{ + // TODO Auto-generated destructor stub +} + +template< typename RSSOLVER_TYPE > +real64 QuasiDynamicEarthQuake< RSSOLVER_TYPE >::updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const +{ + // 1. Setup variables for stress solver + setTargetDispJump( domain ); + + // 2. Solve the momentum balance + real64 const dtAccepted = m_stressSolver->solverStep( time_n, dt, cycleNumber, domain ); + + // 3. Add background stress and possible forcing. + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + arrayView2d< real64 const > const traction = subRegion.getField< contact::traction >(); + + arrayView2d< real64 > const shearTraction = subRegion.getField< rateAndState::shearTraction >(); + arrayView1d< real64 > const normalTraction = subRegion.getField< rateAndState::normalTraction >(); + + arrayView2d< real64 const > const backgroundShearStress = subRegion.getField< rateAndState::backgroundShearStress >(); + arrayView1d< real64 const > const backgroundNormalStress = subRegion.getField< rateAndState::backgroundNormalStress >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + normalTraction[k] = backgroundNormalStress[k] - traction[k][0]; // compressive traction is negative in geos + for( int i = 0; i < 2; ++i ) + { + shearTraction( k, i ) = backgroundShearStress( k, i ) + traction( k, i+1 ); + } + } ); + } ); + } ); + + return dtAccepted; +} + +template< typename RSSOLVER_TYPE > +void QuasiDynamicEarthQuake< RSSOLVER_TYPE >::setTargetDispJump( DomainPartition & domain ) const +{ + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + arrayView2d< real64 > const deltaSlip = subRegion.getField< contact::deltaSlip >(); + arrayView2d< real64 > const targetDispJump = subRegion.getField< contact::targetIncrementalJump >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + targetDispJump( k, 0 ) = 0.0; + targetDispJump( k, 1 ) = deltaSlip( k, 0 ); + targetDispJump( k, 2 ) = deltaSlip( k, 1 ); + } ); + } ); + } ); +} + +template class QuasiDynamicEarthQuake< ImplicitQDRateAndState >; +template class QuasiDynamicEarthQuake< ExplicitQDRateAndState >; + +namespace +{ +typedef QuasiDynamicEarthQuake< ImplicitQDRateAndState > ImplicitQuasiDynamicEarthQuake; +typedef QuasiDynamicEarthQuake< ExplicitQDRateAndState > ExplicitQuasiDynamicEarthQuake; +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, ImplicitQuasiDynamicEarthQuake, string const &, dataRepository::Group * const ) +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, ExplicitQuasiDynamicEarthQuake, string const &, dataRepository::Group * const ) +} + +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp new file mode 100644 index 00000000000..728ec97ca3d --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEarthQuake.hpp @@ -0,0 +1,72 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +#ifndef GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEARTHQUAKE_HPP +#define GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEARTHQUAKE_HPP + +#include "ImplicitQDRateAndState.hpp" + +namespace geos +{ + +template< typename RSSOLVER_TYPE = ImplicitQDRateAndState > +class QuasiDynamicEarthQuake : public RSSOLVER_TYPE +{ +public: + + /// The default nullary constructor is disabled to avoid compiler auto-generation: + QuasiDynamicEarthQuake() = delete; + + /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) + QuasiDynamicEarthQuake( const string & name, + dataRepository::Group * const parent ); + + /// Destructor + virtual ~QuasiDynamicEarthQuake() override; + + static string catalogName() { return RSSOLVER_TYPE::derivedSolverPrefix() + "QuasiDynamicEQ"; } + + /** + * @return Get the final class Catalog name + */ + virtual string getCatalogName() const override { return catalogName(); } + + struct viewKeyStruct : public RSSOLVER_TYPE::viewKeyStruct + { + /// stress solver name + static constexpr char const * stressSolverNameString() { return "stressSolverName"; } + }; + + void postInputInitialization() override final; + + void setTargetDispJump( DomainPartition & domain ) const; + + virtual real64 updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const override final; + +private: + + string m_stressSolverName; + + PhysicsSolverBase * m_stressSolver; + +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQBASE_HPP */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.cpp new file mode 100644 index 00000000000..0c427304f9d --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.cpp @@ -0,0 +1,158 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SpringSlider.cpp + */ + +#include "SpringSlider.hpp" + +#include "dataRepository/InputFlags.hpp" +#include "mesh/DomainPartition.hpp" +#include "rateAndStateFields.hpp" +#include "physicsSolvers/contact/ContactFields.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" +#include "constitutive/contact/RateAndStateFriction.hpp" +#include "ExplicitQDRateAndState.hpp" + +namespace geos +{ + +using namespace dataRepository; +using namespace fields; +using namespace constitutive; + +template< typename RSSOLVER_TYPE > +SpringSlider< RSSOLVER_TYPE >::SpringSlider( const string & name, + Group * const parent ): + RSSOLVER_TYPE( name, parent ) +{} + +template< typename RSSOLVER_TYPE > +SpringSlider< RSSOLVER_TYPE >::~SpringSlider() +{ + // TODO Auto-generated destructor stub +} + +template< typename RSSOLVER_TYPE > +void SpringSlider< RSSOLVER_TYPE >::registerDataOnMesh( Group & meshBodies ) +{ + RSSOLVER_TYPE::registerDataOnMesh( meshBodies ); + + this->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + // 3-component functions on fault + string const labels3Comp[3] = { "normal", "tangent1", "tangent2" }; + subRegion.registerField< contact::dispJump >( this->getName() ). + setDimLabels( 1, labels3Comp ). + reference().template resizeDimension< 1 >( 3 ); + + subRegion.registerField< contact::dispJump_n >( this->getName() ). + setDimLabels( 1, labels3Comp ). + reference().template resizeDimension< 1 >( 3 ); + + string const labels2Comp[2] = { "tangent1", "tangent2" }; + + subRegion.registerField< contact::deltaSlip >( this->getName() ). + setDimLabels( 1, labels2Comp ). + reference().template resizeDimension< 1 >( 2 ); + + subRegion.registerField< contact::deltaSlip_n >( this->getName() ). + setDimLabels( 1, labels2Comp ). + reference().template resizeDimension< 1 >( 2 ); + + subRegion.registerWrapper< string >( viewKeyStruct::frictionLawNameString() ). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setSizedFromParent( 0 ); + + string & frictionLawName = subRegion.getReference< string >( viewKeyStruct::frictionLawNameString() ); + frictionLawName =PhysicsSolverBase::getConstitutiveName< FrictionBase >( subRegion ); + GEOS_ERROR_IF( frictionLawName.empty(), GEOS_FMT( "{}: FrictionBase model not found on subregion {}", + this->getDataContext(), subRegion.getDataContext() ) ); + } ); + } ); +} + +template< typename RSSOLVER_TYPE > +real64 SpringSlider< RSSOLVER_TYPE >::updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const + +{ + GEOS_UNUSED_VAR( cycleNumber, time_n ); + // Spring-slider shear traction computation + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + + arrayView2d< real64 const > const deltaSlip = subRegion.getField< fields::contact::deltaSlip >(); + arrayView2d< real64 > const shearTraction = subRegion.getField< fields::rateAndState::shearTraction >(); + arrayView2d< real64 > const shearTraction_n = subRegion.getField< fields::rateAndState::shearTraction_n >(); + + arrayView1d< real64 > const normalTraction = subRegion.getField< fields::rateAndState::normalTraction >(); + + + string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + RateAndStateFriction const & frictionLaw = this->template getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); + + RateAndStateFriction::KernelWrapper frictionKernelWrapper = frictionLaw.createKernelUpdates(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + SpringSliderParameters springSliderParameters = SpringSliderParameters( normalTraction[k], + frictionKernelWrapper.getACoefficient( k ), + frictionKernelWrapper.getBCoefficient( k ), + frictionKernelWrapper.getDcCoefficient( k ) ); + + + + shearTraction[k][0] = shearTraction_n[k][0] + springSliderParameters.tauRate * dt + - springSliderParameters.springStiffness * deltaSlip[k][0]; + shearTraction[k][1] = shearTraction_n[k][1] + springSliderParameters.tauRate * dt + - springSliderParameters.springStiffness * deltaSlip[k][1]; + } ); + } ); + } ); + return dt; +} + +template class SpringSlider< ImplicitQDRateAndState >; +template class SpringSlider< ExplicitQDRateAndState >; + +namespace +{ +typedef SpringSlider< ImplicitQDRateAndState > ImplicitSpringSlider; +typedef SpringSlider< ExplicitQDRateAndState > ExplicitSpringSlider; +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, ImplicitSpringSlider, string const &, dataRepository::Group * const ) +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, ExplicitSpringSlider, string const &, dataRepository::Group * const ) +} + +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp new file mode 100644 index 00000000000..08c9a829c7e --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/SpringSlider.hpp @@ -0,0 +1,91 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + + +#ifndef GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_SPRINGSLIDER_HPP +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_SPRINGSLIDER_HPP + +#include "physicsSolvers/inducedSeismicity/ImplicitQDRateAndState.hpp" + +namespace geos +{ + +template< typename RSSOLVER_TYPE = ImplicitQDRateAndState > +class SpringSlider : public RSSOLVER_TYPE +{ +public: + + SpringSlider() = delete; + + SpringSlider( const string & name, + dataRepository::Group * const parent ); + + /// Destructor + virtual ~SpringSlider() override; + + static string catalogName() { return RSSOLVER_TYPE::derivedSolverPrefix() + "SpringSlider"; } + + virtual string getCatalogName() const override { return catalogName(); } + + virtual void registerDataOnMesh( dataRepository::Group & meshBodies ) override; + + struct viewKeyStruct : public RSSOLVER_TYPE::viewKeyStruct + {}; + + virtual real64 updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const override final; + +private: + + class SpringSliderParameters + { +public: + + GEOS_HOST_DEVICE + SpringSliderParameters( real64 const normalTraction, real64 const a, real64 const b, real64 const Dc ): + tauRate( 1e-4 ), + springStiffness( 0.0 ) + { + real64 const criticalStiffness = normalTraction * (b - a) / Dc; + springStiffness = 0.9 * criticalStiffness; + } + + /// Default copy constructor + SpringSliderParameters( SpringSliderParameters const & ) = default; + + /// Default move constructor + SpringSliderParameters( SpringSliderParameters && ) = default; + + /// Deleted default constructor + SpringSliderParameters() = delete; + + /// Deleted copy assignment operator + SpringSliderParameters & operator=( SpringSliderParameters const & ) = delete; + + /// Deleted move assignment operator + SpringSliderParameters & operator=( SpringSliderParameters && ) = delete; + + real64 tauRate; + + real64 springStiffness; + }; +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_SPRINGSLIDER_HPP */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ExplicitRateAndStateKernels.hpp similarity index 65% rename from src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp rename to src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ExplicitRateAndStateKernels.hpp index d91dff75ef6..350c1fc2c76 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ExplicitRateAndStateKernels.hpp @@ -13,13 +13,10 @@ * ------------------------------------------------------------------------------------------------------------ */ -#ifndef GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ -#define GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ +#ifndef GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_EXPLICITRATEANDSTATEKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_EXPLICITRATEANDSTATEKERNELS_HPP_ -#include "common/DataTypes.hpp" -#include "common/GEOS_RAJA_Interface.hpp" -#include "constitutive/contact/RateAndStateFriction.hpp" -#include "physicsSolvers/inducedSeismicity/rateAndStateFields.hpp" +#include "RateAndStateKernelsBase.hpp" #include "denseLinearAlgebra/denseLASolvers.hpp" namespace geos @@ -28,137 +25,6 @@ namespace geos namespace rateAndStateKernels { -// TBD: Pass the kernel and add getters for relevant fields to make this function general purpose and avoid -// wrappers? -GEOS_HOST_DEVICE -static void projectSlipRateBase( localIndex const k, - real64 const frictionCoefficient, - real64 const shearImpedance, - arrayView2d< real64 const > const traction, - arrayView1d< real64 const > const slipRate, - arrayView2d< real64 > const slipVelocity ) -{ - // Project slip rate onto shear traction to get slip velocity components - real64 const frictionForce = traction[k][0] * frictionCoefficient; - real64 const projectionScaling = 1.0 / ( shearImpedance + frictionForce / slipRate[k] ); - slipVelocity[k][0] = projectionScaling * traction[k][1]; - slipVelocity[k][1] = projectionScaling * traction[k][2]; -} - -/** - * @class ImplicitFixedStressRateAndStateKernel - * - * @brief - * - * @details - */ -class ImplicitFixedStressRateAndStateKernel -{ -public: - - ImplicitFixedStressRateAndStateKernel( SurfaceElementSubRegion & subRegion, - constitutive::RateAndStateFriction const & frictionLaw, - real64 const shearImpedance ): - m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), - m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ), - m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ), - m_traction( subRegion.getField< fields::contact::traction >() ), - m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), - m_shearImpedance( shearImpedance ), - m_frictionLaw( frictionLaw.createKernelUpdates() ) - {} - - /** - * @struct StackVariables - * @brief Kernel variables located on the stack - */ - struct StackVariables - { -public: - - GEOS_HOST_DEVICE - StackVariables( ) - {} - - real64 jacobian[2][2]{}; - - real64 rhs[2]{}; - - }; - - GEOS_HOST_DEVICE - void setup( localIndex const k, - real64 const dt, - StackVariables & stack ) const - { - real64 const normalTraction = m_traction[k][0]; - real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); - // Eq 1: Scalar force balance for slipRate and shear traction magnitude - stack.rhs[0] = shearTractionMagnitude - m_shearImpedance * m_slipRate[k] - - normalTraction * m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); - real64 const dFriction[2] = { -normalTraction * m_frictionLaw.dFrictionCoefficient_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ), - -m_shearImpedance - normalTraction * m_frictionLaw.dFrictionCoefficient_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) }; - - // Eq 2: slip law - stack.rhs[1] = (m_stateVariable[k] - m_stateVariable_n[k]) / dt - m_frictionLaw.stateEvolution( k, m_slipRate[k], m_stateVariable[k] ); - real64 const dStateEvolutionLaw[2] = { 1 / dt - m_frictionLaw.dStateEvolution_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ), - -m_frictionLaw.dStateEvolution_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) }; - - // Assemble Jacobian matrix - stack.jacobian[0][0] = dFriction[0]; // derivative of Eq 1 w.r.t. stateVariable - stack.jacobian[0][1] = dFriction[1]; // derivative of Eq 1 w.r.t. slipRate - stack.jacobian[1][0] = dStateEvolutionLaw[0]; // derivative of Eq 2 w.r.t. stateVariable - stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. slipRate - } - - GEOS_HOST_DEVICE - void solve( localIndex const k, - StackVariables & stack ) const - { - /// Solve 2x2 system - real64 solution[2] = {0.0, 0.0}; - denseLinearAlgebra::solve< 2 >( stack.jacobian, stack.rhs, solution ); - - // Update variables - m_stateVariable[k] -= solution[0]; - m_slipRate[k] -= solution[1]; - } - - GEOS_HOST_DEVICE - void projectSlipRate( localIndex const k ) const - { - real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); - projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity ); - } - - GEOS_HOST_DEVICE - camp::tuple< int, real64 > checkConvergence( StackVariables const & stack, - real64 const tol ) const - { - real64 const residualNorm = LvArray::tensorOps::l2Norm< 2 >( stack.rhs ); - int const converged = residualNorm < tol ? 1 : 0; - camp::tuple< int, real64 > result { converged, residualNorm }; - return result; - } - -private: - - arrayView1d< real64 > const m_slipRate; - - arrayView1d< real64 > const m_stateVariable; - - arrayView1d< real64 const > const m_stateVariable_n; - - arrayView2d< real64 const > const m_traction; - - arrayView2d< real64 > const m_slipVelocity; - - real64 const m_shearImpedance; - - constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw; - -}; - /** * @class ExplicitRateAndStateKernel * @@ -175,7 +41,8 @@ class ExplicitRateAndStateKernel real64 const shearImpedance ): m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ), - m_traction( subRegion.getField< fields::contact::traction >() ), + m_normalTraction( subRegion.getField< fields::rateAndState::normalTraction >() ), + m_shearTraction( subRegion.getField< fields::rateAndState::shearTraction >() ), m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), m_shearImpedance( shearImpedance ), m_frictionLaw( frictionLaw.createKernelUpdates() ) @@ -204,8 +71,8 @@ class ExplicitRateAndStateKernel StackVariables & stack ) const { GEOS_UNUSED_VAR( dt ); - real64 const normalTraction = m_traction[k][0]; - real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); + real64 const normalTraction = m_normalTraction[k]; + real64 const shearTractionMagnitude = LvArray::tensorOps::l2Norm< 2 >( m_shearTraction[k] ); // Slip rate is bracketed between [0, shear traction magnitude / shear impedance] // If slip rate is outside the bracket, re-initialize to the middle value @@ -224,7 +91,7 @@ class ExplicitRateAndStateKernel // Slip rate is bracketed between [0, shear traction magnitude / shear impedance] // Check that the update did not end outside of the bracket. - real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); + real64 const shearTractionMagnitude = LvArray::tensorOps::l2Norm< 2 >( m_shearTraction[k] ); real64 const upperBound = shearTractionMagnitude/m_shearImpedance; if( m_slipRate[k] > upperBound ) m_slipRate[k] = 0.5*upperBound; @@ -245,7 +112,43 @@ class ExplicitRateAndStateKernel void projectSlipRate( localIndex const k ) const { real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); - projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity ); + projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_normalTraction, m_shearTraction, m_slipRate, m_slipVelocity ); + } + + GEOS_HOST_DEVICE + void udpateVariables( localIndex const k ) const + { + projectSlipRate( k ); + } + + GEOS_HOST_DEVICE + void resetState( localIndex const k ) const + { + GEOS_UNUSED_VAR( k ); + } + + /** + * @brief Performs the kernel launch + * @tparam KernelType The Rate-and-state kernel to launch + * @tparam POLICY the policy used in the RAJA kernels + */ + template< typename POLICY > + static real64 + solveRateAndStateEquation( SurfaceElementSubRegion & subRegion, + ExplicitRateAndStateKernel & kernel, + real64 dt, + integer const maxNewtonIter, + real64 const newtonTol ) + { + GEOS_MARK_FUNCTION; + + newtonSolve< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); + + forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + kernel.projectSlipRate( k ); + } ); + return dt; } private: @@ -254,7 +157,9 @@ class ExplicitRateAndStateKernel arrayView1d< real64 > const m_stateVariable; - arrayView2d< real64 const > const m_traction; + arrayView1d< real64 const > const m_normalTraction; + + arrayView2d< real64 const > const m_shearTraction; arrayView2d< real64 > const m_slipVelocity; @@ -264,65 +169,6 @@ class ExplicitRateAndStateKernel }; - -/** - * @brief Performs the kernel launch - * @tparam KernelType The Rate-and-state kernel to launch - * @tparam POLICY the policy used in the RAJA kernels - */ -template< typename KernelType, typename POLICY > -static void -createAndLaunch( SurfaceElementSubRegion & subRegion, - string const & frictionLawNameKey, - real64 const shearImpedance, - integer const maxIterNewton, - real64 const newtonTol, - real64 const time_n, - real64 const dt ) -{ - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( time_n ); - - string const & frictionaLawName = subRegion.getReference< string >( frictionLawNameKey ); - constitutive::RateAndStateFriction const & frictionLaw = subRegion.getConstitutiveModel< constitutive::RateAndStateFriction >( frictionaLawName ); - KernelType kernel( subRegion, frictionLaw, shearImpedance ); - - // Newton loop (outside of the kernel launch) - bool allConverged = false; - for( integer iter = 0; iter < maxIterNewton; iter++ ) - { - RAJA::ReduceMin< parallelDeviceReduce, int > converged( 1 ); - RAJA::ReduceMax< parallelDeviceReduce, real64 > residualNorm( 0.0 ); - forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - typename KernelType::StackVariables stack; - kernel.setup( k, dt, stack ); - kernel.solve( k, stack ); - auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, newtonTol ); - converged.min( elementConverged ); - residualNorm.max( elementResidualNorm ); - } ); - - real64 const maxResidualNorm = MpiWrapper::max( residualNorm.get() ); - GEOS_LOG_RANK_0( GEOS_FMT( "-----iter {} : residual = {:.10e} ", iter, maxResidualNorm ) ); - - if( converged.get() ) - { - allConverged = true; - break; - } - } - if( !allConverged ) - { - GEOS_ERROR( " Failed to converge" ); - } - forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) - { - kernel.projectSlipRate( k ); - } ); -} - /** * @brief Butcher table for embedded RK3(2) method using Kuttas third order * method for the high-order update, and an explicit trapezoidal rule @@ -366,7 +212,8 @@ struct BogackiShampine32Table * * @tparam Butcher table defining the Runge-Kutta method. */ -template< typename TABLE_TYPE > class EmbeddedRungeKuttaKernel +template< typename TABLE_TYPE > +class EmbeddedRungeKuttaKernel { public: @@ -378,8 +225,8 @@ template< typename TABLE_TYPE > class EmbeddedRungeKuttaKernel m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), m_slipVelocity_n( subRegion.getField< fields::rateAndState::slipVelocity_n >() ), - m_deltaSlip( subRegion.getField< fields::rateAndState::deltaSlip >() ), - m_deltaSlip_n( subRegion.getField< fields::rateAndState::deltaSlip_n >() ), + m_deltaSlip( subRegion.getField< fields::contact::deltaSlip >() ), + m_deltaSlip_n( subRegion.getField< fields::contact::deltaSlip_n >() ), m_dispJump( subRegion.getField< fields::contact::dispJump >() ), m_dispJump_n( subRegion.getField< fields::contact::dispJump_n >() ), m_error( subRegion.getField< fields::rateAndState::error >() ), @@ -579,4 +426,4 @@ template< typename TABLE_TYPE > class EmbeddedRungeKuttaKernel } /* namespace geos */ -#endif /* GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ */ +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_EXPLICITRATEANDSTATEKERNELS_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ImplicitRateAndStateKernels.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ImplicitRateAndStateKernels.hpp new file mode 100644 index 00000000000..80d65a59ffd --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/ImplicitRateAndStateKernels.hpp @@ -0,0 +1,226 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_IMPLICITRATEANDSTATEKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_IMPLICITRATEANDSTATEKERNELS_HPP_ + +#include "RateAndStateKernelsBase.hpp" +#include "denseLinearAlgebra/denseLASolvers.hpp" + +namespace geos +{ + +namespace rateAndStateKernels +{ + +/** + * @class ImplicitFixedStressRateAndStateKernel + * + * @brief + * + * @details + */ +class ImplicitFixedStressRateAndStateKernel +{ +public: + + ImplicitFixedStressRateAndStateKernel( SurfaceElementSubRegion & subRegion, + constitutive::RateAndStateFriction const & frictionLaw, + real64 const shearImpedance ): + m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), + m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ), + m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ), + m_slipRate_n( subRegion.getField< fields::rateAndState::slipRate_n >() ), + m_normalTraction( subRegion.getField< fields::rateAndState::normalTraction >() ), + m_shearTraction( subRegion.getField< fields::rateAndState::shearTraction >() ), + m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), + m_shearImpedance( shearImpedance ), + m_frictionLaw( frictionLaw.createKernelUpdates() ) + {} + + /** + * @struct StackVariables + * @brief Kernel variables located on the stack + */ + struct StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables( ) = default; + + real64 jacobian[2][2]{}; + + real64 rhs[2]{}; + + }; + + GEOS_HOST_DEVICE + void setup( localIndex const k, + real64 const dt, + StackVariables & stack ) const + { + real64 const normalTraction = m_normalTraction[k]; + real64 const shearTractionMagnitude = LvArray::tensorOps::l2Norm< 2 >( m_shearTraction[k] ); + + // Eq 1: Scalar force balance for slipRate and shear traction magnitude + stack.rhs[0] = shearTractionMagnitude - m_shearImpedance * m_slipRate[k] + - normalTraction * m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); + real64 const dFriction[2] = { -normalTraction * m_frictionLaw.dFrictionCoefficient_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ), + -m_shearImpedance - normalTraction * m_frictionLaw.dFrictionCoefficient_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) }; + + // Eq 2: slip law + stack.rhs[1] = (m_stateVariable[k] - m_stateVariable_n[k]) / dt - m_frictionLaw.stateEvolution( k, m_slipRate[k], m_stateVariable[k] ); + real64 const dStateEvolutionLaw[2] = { 1.0 / dt - m_frictionLaw.dStateEvolution_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ), + -m_frictionLaw.dStateEvolution_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) }; + + // Assemble Jacobian matrix + stack.jacobian[0][0] = dFriction[0]; // derivative of Eq 1 w.r.t. stateVariable + stack.jacobian[0][1] = dFriction[1]; // derivative of Eq 1 w.r.t. slipRate + stack.jacobian[1][0] = dStateEvolutionLaw[0]; // derivative of Eq 2 w.r.t. stateVariable + stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. slipRate + } + + GEOS_HOST_DEVICE + void solve( localIndex const k, + StackVariables & stack ) const + { + /// Solve 2x2 system + real64 solution[2] = {0.0, 0.0}; + LvArray::tensorOps::scale< 2 >( stack.rhs, -1.0 ); + denseLinearAlgebra::solve< 2 >( stack.jacobian, stack.rhs, solution ); + + // Slip rate is bracketed between [0, shear traction magnitude / shear impedance] + // Check that the update did not end outside of the bracket. + real64 const shearTractionMagnitude = LvArray::tensorOps::l2Norm< 2 >( m_shearTraction[k] ); + real64 const upperBound = shearTractionMagnitude / m_shearImpedance - m_slipRate[k]; + real64 const lowerBound = -m_slipRate[k]; + + real64 scalingFactor = 1.0; + if( solution[1] > upperBound ) + { + scalingFactor = 0.5 * upperBound / solution[1]; + } + else if( solution[1] < lowerBound ) + { + scalingFactor = 0.5 * lowerBound / solution[1]; + } + + LvArray::tensorOps::scale< 2 >( solution, scalingFactor ); + + // Update variables + m_stateVariable[k] += solution[0]; + m_slipRate[k] += solution[1]; + } + + GEOS_HOST_DEVICE + void projectSlipRate( localIndex const k ) const + { + real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); + projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_normalTraction, m_shearTraction, m_slipRate, m_slipVelocity ); + } + + GEOS_HOST_DEVICE + void udpateVariables( localIndex const k ) const + { + projectSlipRate( k ); + m_stateVariable_n[k] = m_stateVariable[k]; + m_slipRate_n[k] = m_slipRate[k]; + } + + GEOS_HOST_DEVICE + camp::tuple< int, real64 > checkConvergence( StackVariables const & stack, + real64 const tol ) const + { + real64 const residualNorm = LvArray::tensorOps::l2Norm< 2 >( stack.rhs ); + int const converged = residualNorm < tol ? 1 : 0; + camp::tuple< int, real64 > result { converged, residualNorm }; + return result; + } + + GEOS_HOST_DEVICE + void resetState( localIndex const k ) const + { + m_stateVariable[k] = m_stateVariable_n[k]; + m_slipRate[k] = m_slipRate_n[k]; + } + + template< typename POLICY > + static real64 solveRateAndStateEquation( SurfaceElementSubRegion & subRegion, + ImplicitFixedStressRateAndStateKernel & kernel, + real64 dt, + integer const maxNewtonIter, + real64 const newtonTol ) + { + bool converged = false; + for( integer attempt = 0; attempt < 5; attempt++ ) + { + if( attempt > 0 ) + { + forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + kernel.resetState( k ); + } ); + } + GEOS_LOG_RANK_0( GEOS_FMT( " Attempt {} ", attempt ) ); + converged = newtonSolve< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); + if( converged ) + { + forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + kernel.udpateVariables( k ); + } ); + return dt; + } + else + { + GEOS_LOG_RANK_0( GEOS_FMT( " Attempt {} failed. Halving dt and retrying.", attempt ) ); + dt *= 0.5; + } + } + if( !converged ) + { + GEOS_ERROR( "Maximum number of attempts reached without convergence." ); + } + return dt; + } + +private: + + arrayView1d< real64 > const m_slipRate; + + arrayView1d< real64 > const m_stateVariable; + + arrayView1d< real64 > const m_stateVariable_n; + + arrayView1d< real64 > const m_slipRate_n; + + arrayView1d< real64 > const m_normalTraction; + + arrayView2d< real64 > const m_shearTraction; + + arrayView2d< real64 > const m_slipVelocity; + + real64 const m_shearImpedance; + + constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw; + +}; + +} /* namespace rateAndStateKernels */ + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_IMPLICITRATEANDSTATEKERNELS_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernelsBase.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernelsBase.hpp new file mode 100644 index 00000000000..c889bb94723 --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernelsBase.hpp @@ -0,0 +1,123 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_RATEANDSTATEKERNELSBASE_HPP_ +#define GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_RATEANDSTATEKERNELSBASE_HPP_ + +#include "common/DataTypes.hpp" +#include "common/GEOS_RAJA_Interface.hpp" +#include "constitutive/contact/RateAndStateFriction.hpp" +#include "physicsSolvers/inducedSeismicity/rateAndStateFields.hpp" + +namespace geos +{ + +namespace rateAndStateKernels +{ + +// TBD: Pass the kernel and add getters for relevant fields to make this function general purpose and avoid +// wrappers? +GEOS_HOST_DEVICE +static void projectSlipRateBase( localIndex const k, + real64 const frictionCoefficient, + real64 const shearImpedance, + arrayView1d< real64 const > const normalTraction, + arrayView2d< real64 const > const shearTraction, + arrayView1d< real64 const > const slipRate, + arrayView2d< real64 > const slipVelocity ) +{ + // Project slip rate onto shear traction to get slip velocity components + real64 const frictionForce = normalTraction[k] * frictionCoefficient; + real64 const projectionScaling = 1.0 / ( shearImpedance + frictionForce / slipRate[k] ); + slipVelocity[k][0] = projectionScaling * shearTraction[k][0]; + slipVelocity[k][1] = projectionScaling * shearTraction[k][1]; +} + +template< typename POLICY, typename KERNEL_TYPE > +static bool newtonSolve( SurfaceElementSubRegion & subRegion, + KERNEL_TYPE & kernel, + real64 const dt, + integer const maxNewtonIter, + real64 const newtonTol ) +{ + bool allConverged = false; + for( integer iter = 0; iter < maxNewtonIter; iter++ ) + { + RAJA::ReduceMin< parallelDeviceReduce, int > converged( 1 ); + RAJA::ReduceMax< parallelDeviceReduce, real64 > residualNorm( 0.0 ); + forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + typename KERNEL_TYPE::StackVariables stack; + kernel.setup( k, dt, stack ); + kernel.solve( k, stack ); + auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, newtonTol ); + converged.min( elementConverged ); + residualNorm.max( elementResidualNorm ); + } ); + + real64 const maxResidualNorm = MpiWrapper::max( residualNorm.get() ); + GEOS_LOG_RANK_0( GEOS_FMT( " Newton iter {} : residual = {:.10e} ", iter, maxResidualNorm ) ); + + if( converged.get() ) + { + allConverged = true; + break; + } + } + return allConverged; +} + +/** + * @brief Performs the kernel launch + * @tparam POLICY the policy used in the RAJA kernels + */ +template< typename KERNEL_TYPE, typename POLICY > +static void +createAndLaunch( SurfaceElementSubRegion & subRegion, + string const & frictionLawNameKey, + real64 const shearImpedance, + integer const maxNewtonIter, + real64 const newtonTol, + real64 const time_n, + real64 const totalDt ) +{ + GEOS_MARK_FUNCTION; + + GEOS_UNUSED_VAR( time_n ); + + string const & frictionaLawName = subRegion.getReference< string >( frictionLawNameKey ); + constitutive::RateAndStateFriction const & frictionLaw = subRegion.getConstitutiveModel< constitutive::RateAndStateFriction >( frictionaLawName ); + KERNEL_TYPE kernel( subRegion, frictionLaw, shearImpedance ); + + real64 dtRemaining = totalDt; + real64 dt = totalDt; + for( integer subStep = 0; subStep < 5 && dtRemaining > 0.0; ++subStep ) + { + real64 dtAccepted = KERNEL_TYPE::template solveRateAndStateEquation< POLICY >( subRegion, kernel, dt, maxNewtonIter, newtonTol ); + dtRemaining -= dtAccepted; + + if( dtRemaining > 0.0 ) + { + dt = dtAccepted; + } + GEOS_LOG_RANK_0( GEOS_FMT( " sub-step = {} completed, dt = {}, remaining dt = {}", subStep, dt, dtRemaining ) ); + } +} + +} /* namespace rateAndStateKernels */ + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCEDSEISMICITY_KERNELS_RATEANDSTATEKERNELSBASE_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp index e060dceb28a..599a6f0e0c4 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp @@ -35,15 +35,23 @@ namespace rateAndState DECLARE_FIELD( slipRate, "slipRate", array1d< real64 >, - 1.0e-6, - NOPLOT, + 0.0, + LEVEL_0, WRITE_AND_READ, "Slip rate" ); +DECLARE_FIELD( slipRate_n, + "slipRate_n", + array1d< real64 >, + 0.0, + NOPLOT, + WRITE_AND_READ, + "Slip rate at timestep n." ); + DECLARE_FIELD( slipVelocity, "slipVelocity", array2d< real64 >, - 0.70710678118e-6, + 0.0, LEVEL_0, WRITE_AND_READ, "Slip velocity" ); @@ -51,7 +59,7 @@ DECLARE_FIELD( slipVelocity, DECLARE_FIELD( slipVelocity_n, "slipVelocity_n", array2d< real64 >, - 0.70710678118e-6, + 0.0, NOPLOT, WRITE_AND_READ, "Slip velocity at previous time step" ); @@ -72,23 +80,53 @@ DECLARE_FIELD( stateVariable_n, WRITE_AND_READ, "Initial rate- and state-dependent friction state variable at this time step" ); +DECLARE_FIELD( normalTraction, + "normalTraction", + array1d< real64 >, + 0.0, + LEVEL_0, + WRITE_AND_READ, + "Normal traction" ); -DECLARE_FIELD( deltaSlip, - "deltaSlip", +DECLARE_FIELD( shearTraction, + "shearTraction", array2d< real64 >, 0.0, LEVEL_0, WRITE_AND_READ, - "Slip increment" ); + "Shear traction" ); + +DECLARE_FIELD( normalTraction_n, + "normalTraction_n", + array1d< real64 >, + 0.0, + LEVEL_0, + WRITE_AND_READ, + "Normal traction at previous timestep n." ); -DECLARE_FIELD( deltaSlip_n, - "deltaSlip_n", +DECLARE_FIELD( shearTraction_n, + "shearTraction_n", array2d< real64 >, 0.0, - NOPLOT, + LEVEL_0, WRITE_AND_READ, - "Initial slip increment at this time step" ); + "Shear traction at previous timestep n." ); +DECLARE_FIELD( backgroundNormalStress, + "backgroundNormalStress", + array1d< real64 >, + 0.0, + LEVEL_0, + WRITE_AND_READ, + "Background Normal Stress" ); + +DECLARE_FIELD( backgroundShearStress, + "backgroundShearStress", + array2d< real64 >, + 0.0, + LEVEL_0, + WRITE_AND_READ, + "Background Shear Stress" ); DECLARE_FIELD( rungeKuttaStageRates, "rungeKuttaStageRates", diff --git a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt index 89a34cbc276..d9faf1fbae4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/multiphysics/CMakeLists.txt @@ -9,6 +9,7 @@ set( physicsSolvers_headers multiphysics/HydrofractureSolver.hpp multiphysics/HydrofractureSolverKernels.hpp multiphysics/MultiphasePoromechanics.hpp + multiphysics/OneWayCoupledFractureFlowContactMechanics.hpp multiphysics/PhaseFieldFractureSolver.hpp multiphysics/PoromechanicsInitialization.hpp multiphysics/PoromechanicsFields.hpp @@ -46,6 +47,7 @@ set( physicsSolvers_sources multiphysics/FlowProppantTransportSolver.cpp multiphysics/HydrofractureSolver.cpp multiphysics/MultiphasePoromechanics.cpp + multiphysics/OneWayCoupledFractureFlowContactMechanics.cpp multiphysics/PhaseFieldFractureSolver.cpp multiphysics/PoromechanicsInitialization.cpp multiphysics/SinglePhasePoromechanics.cpp diff --git a/src/coreComponents/physicsSolvers/multiphysics/OneWayCoupledFractureFlowContactMechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/OneWayCoupledFractureFlowContactMechanics.cpp new file mode 100644 index 00000000000..6f5c291f1ce --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/OneWayCoupledFractureFlowContactMechanics.cpp @@ -0,0 +1,88 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file OneWayCoupledFractureFlowContactMechanics.cpp + */ + +#include "OneWayCoupledFractureFlowContactMechanics.hpp" +#include "physicsSolvers/contact/ContactFields.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "mesh/DomainPartition.hpp" + + +namespace geos +{ +using namespace dataRepository; + +template< typename FLOW_SOLVER > +OneWayCoupledFractureFlowContactMechanics< FLOW_SOLVER >::OneWayCoupledFractureFlowContactMechanics( const string & name, + Group * const parent ) + : Base( name, parent ) +{} + +template< typename FLOW_SOLVER > +void OneWayCoupledFractureFlowContactMechanics< FLOW_SOLVER >::postInputInitialization() +{ + bool const isSequential = this->getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::Sequential; + GEOS_THROW_IF( !isSequential, + "Only sequential coupling is allowed for this solver.", + InputError ); + + Base::postInputInitialization(); +} + + +template< typename FLOW_SOLVER > +real64 OneWayCoupledFractureFlowContactMechanics< FLOW_SOLVER >::sequentiallyCoupledSolverStep( real64 const & time_n, + real64 const & dt, + int const cycleNumber, + DomainPartition & domain ) +{ + forEachArgInTuple( m_solvers, [&]( auto & solver, auto ) + { + solver->solverStep( time_n, dt, cycleNumber, domain ); + } ); + + this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + arrayView2d< real64 > const traction = subRegion.getField< fields::contact::traction >(); + arrayView1d< real64 > const pressure = subRegion.getField< fields::flow::pressure >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + traction( k, 0 ) = traction( k, 0 ) + pressure[k]; + } ); + } ); + } ); + + return dt; +} + +template class OneWayCoupledFractureFlowContactMechanics< SinglePhaseBase >; + +namespace +{ +typedef OneWayCoupledFractureFlowContactMechanics< SinglePhaseBase > OneWayCoupledFractureFlowContactMechanicsSinglePhase; +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, OneWayCoupledFractureFlowContactMechanicsSinglePhase, string const &, dataRepository::Group * const ) +} + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/multiphysics/OneWayCoupledFractureFlowContactMechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/OneWayCoupledFractureFlowContactMechanics.hpp new file mode 100644 index 00000000000..03b6b1a6f4b --- /dev/null +++ b/src/coreComponents/physicsSolvers/multiphysics/OneWayCoupledFractureFlowContactMechanics.hpp @@ -0,0 +1,93 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file OneWayCoupledFractureFlowContactMechanics.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_ONEWAYCOUPLEDFRACTUREFLOWCONTACTMECHANICS_HPP_ +#define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_ONEWAYCOUPLEDFRACTUREFLOWCONTACTMECHANICS_HPP_ + +#include "physicsSolvers/multiphysics/CoupledSolver.hpp" +#include "physicsSolvers/contact/SolidMechanicsLagrangeContactBubbleStab.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" +#include "dataRepository/Group.hpp" + +namespace geos +{ + +template< typename FLOW_SOLVER = SinglePhaseBase > +class OneWayCoupledFractureFlowContactMechanics : public CoupledSolver< FLOW_SOLVER, SolidMechanicsLagrangeContactBubbleStab > +{ +public: + + using Base = CoupledSolver< FLOW_SOLVER, SolidMechanicsLagrangeContactBubbleStab >; + using Base::m_solvers; + using Base::m_dofManager; + using Base::m_localMatrix; + using Base::m_rhs; + using Base::m_solution; + + enum class SolverType : integer + { + Flow = 0, + SolidMechanics = 1 + }; + + /** + * @brief main constructor for OneWayCoupledFractureFlowContactMechanics objects + * @param name the name of this instantiation of OneWayCoupledFractureFlowContactMechanics in the repository + * @param parent the parent group of this instantiation of OneWayCoupledFractureFlowContactMechanics + */ + OneWayCoupledFractureFlowContactMechanics( const string & name, + dataRepository::Group * const parent ); + + /// Destructor for the class + ~OneWayCoupledFractureFlowContactMechanics() override {} + + /** + * @brief name of the node manager in the object catalog + * @return string that contains the catalog name to generate a new OneWayCoupledFractureFlowContactMechanics object through the object + * catalog. + */ + static string catalogName() + { + return "OneWayCoupledFractureFlowContactMechanics"; + } + + /** + * @copydoc PhysicsSolverBase::getCatalogName() + */ + string getCatalogName() const override { return catalogName(); } + + virtual real64 sequentiallyCoupledSolverStep( real64 const & time_n, + real64 const & dt, + int const cycleNumber, + DomainPartition & domain ) override final; + + virtual void postInputInitialization() override final; + + /**@}*/ + +private: + + struct viewKeyStruct : public Base::viewKeyStruct + {}; + +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_MULTIPHYSICS_ONEWAYCOUPLEDFRACTUREFLOWCONTACTMECHANICS_HPP_ */ diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index e130b1e56d4..ff564e050bc 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -363,6 +363,14 @@ + + + + + + + + @@ -371,6 +379,14 @@ + + + + + + + + @@ -383,6 +399,10 @@ + + + + @@ -395,14 +415,6 @@ - - - - - - - - @@ -2282,16 +2294,19 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + + + + - - @@ -3176,6 +3191,74 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -3273,6 +3356,78 @@ Local- Add jump stabilization on interior of macro elements--> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -3395,6 +3550,40 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + + + + + + + + + + + + + + + + + + + + + @@ -3535,78 +3724,6 @@ Level 0 outputs no specific information for this solver. Higher levels require m - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index 3abf5b863f0..12c2814d7e2 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -529,16 +529,19 @@ + + + + + - - @@ -917,6 +920,28 @@ + + + + + + + + + + + + + + + + + + + + + + @@ -945,7 +970,7 @@ - + @@ -956,22 +981,29 @@ - + - - - - - + + + + + + + + + + + + @@ -983,19 +1015,23 @@ + + - + + + - + @@ -1008,7 +1044,7 @@ - + @@ -1019,18 +1055,20 @@ - + + + - +