[CIG-SHORT] Fault opening with nonzero traction

Tue Oct 7 10:08:16 PDT 2014

Hi everyone,

     I am working on a simulation of the interseismic deformation of a 
subduction zone. The upper subduction slab interface is a dynamic fault 
with friction 0.6 and zero cohesion, and the lower interface is a 
kinematic fault. The simulation stops with error

    RuntimeError: WARNING! Fault opening with nonzero traction.,
    v_fault: 16935, opening: 1.37307e-06, normal traction: -2.32876e-11

The faults are supposed to rupture without opening. I tried the 
open_free_interface but it does not change anything. Is there a way to 
force no opening for dynamic rupture? The .cfg files are attached. Thanks.

Jiangzhi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20141007/ac2c9be7/attachment-0001.html>
-------------- next part --------------
[pylithapp]

# ----------------------------------------------------------------------
# PROBLEM DESCRIPTION
# ----------------------------------------------------------------------
#
# This simulation involves aseismic creep along the interfaces between
# the subducting oceanic litho and the mantle. Both interfaces have constant slip rate 8 cm/yr. The slip time is 40 years for quick test. The slab top interface has friction.
#
# ----------------------------------------------------------------------
# RUNNING THE SIMULATON
# ----------------------------------------------------------------------
#
# This is not a self-contained simulation configuration file. This
# file specifies only the boundary conditions and earthquake
# parameters for the simulation. The general quasi-static and mesh
# parameters are specificed in the pylithapp.cfg file which PyLith
# reads by default.
#
# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
[pylithapp.timedependent.formulation.time_step]
total_time = 40.0*year
dt = 10.0*year

# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
[pylithapp.timedependent]
# Set bc to an array of 3 boundary conditions:
# 'boundary_west_mantle'
# 'boundary_east'
# 'boundary_bottom_mantle'
# 'boundary_side'
bc = [boundary_west_mantle,boundary_east,boundary_bottom_mantle,boundary_side]

# For all boundaries, we fix the displacement normal to the boundary
# (roller boundary condition) by retaining the default ZeroDispDB,
# which specifies a zero value.
#
# The label corresponds to the name of the nodeset in CUBIT.

# west boundary (mantle)
[pylithapp.timedependent.bc.boundary_west_mantle]
bc_dof = [0]
label = bndry_west_mantle
db_initial.label = Dirichlet BC on west boundary (mantle)

# east boundary
[pylithapp.timedependent.bc.boundary_east]
bc_dof = [0]
label = bndry_east
db_initial.label = Dirichlet BC on east boundary

# Bottom boundary (mantle)
[pylithapp.timedependent.bc.boundary_bottom_mantle]
bc_dof = [2]
label = bndry_bot_mantle
db_initial.label = Dirichlet BC on bottom boundary (mantle)

# front and back sides
[pylithapp.timedependent.bc.boundary_side]
bc_dof = [1]
label = bndry_side
db_initial.label = Dirichlet BC on side boundary (all)

# ----------------------------------------------------------------------
# faults
# ----------------------------------------------------------------------
[pylithapp.timedependent]
interfaces = [fault_slabtop,fault_slabbot]

# Set the type of fault interface condition.
[pylithapp.timedependent.interfaces]
fault_slabtop = pylith.faults.FaultCohesiveDyn
fault_slabbot = pylith.faults.FaultCohesiveKin

# Slab top --------------------
[pylithapp.timedependent.interfaces.fault_slabtop]
# The label corresponds to the name of the nodeset in CUBIT.
label = fault_slabtop
id = 100
up_dir = [0,0,1]
zero_tolerance = 1.0e-11
# We must define the quadrature information for fault cells.
# The fault cells are 2D (surface).
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

# Output direction information for debugging fault problems.
output.vertex_info_fields = [normal_dir,strike_dir,dip_dir]

# Set static friction model parameters using a uniform DB. Set the
# static coefficient of friction to 0.6 and cohesion to 0.0 Pa.
open_free_surface = True
friction = pylith.friction.StaticFriction
friction.label = Static friction

friction.db_properties = spatialdata.spatialdb.UniformDB
friction.db_properties.label = Static friction
friction.db_properties.values = [friction-coefficient,cohesion]
friction.db_properties.data = [0.6,0.0*Pa]

# Slab bottom --------------------
[pylithapp.timedependent.interfaces.fault_slabbot]
# The label corresponds to the name of the nodeset in CUBIT.
label = fault_slabbot
id = 101
up_dir = [-0.04,0,0.9992]

# We must define the quadrature information for fault cells.
# The fault cells are 2D (surface).
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

# Output direction information for debugging fault problems.
output.vertex_info_fields = [normal_dir,strike_dir,dip_dir,final_slip_rupture,slip_time_rupture]

# Switch to constant slip rate time function.
[pylithapp.timedependent.interfaces.fault_slabbot.eq_srcs.rupture]
slip_function = pylith.faults.ConstRateSlipFn

# The creep rate and slip time are uniform, so use UniformDB for simplicity.
[pylithapp.timedependent.interfaces.fault_slabbot.eq_srcs.rupture.slip_function]
slip_rate = spatialdata.spatialdb.UniformDB
slip_rate.label = Slip rate
slip_rate.values = [left-lateral-slip,reverse-slip,fault-opening]
slip_rate.data = [0.0*cm/year, -8.0*cm/year, 0.0*cm/year]

# Slip time is uniform, so use UniformDB for convenience
slip_time = spatialdata.spatialdb.UniformDB
slip_time.label = Slip time
slip_time.values = [slip-time]
slip_time.data = [0.0*year]

# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
# Domain
[pylithapp.problem.formulation.output.domain]
writer.filename = output/interseismic1/interseismic.h5

# Ground surface
[pylithapp.problem.formulation.output.subdomain]
writer.filename = output/interseismic1/interseismic-groundsurf.h5

# Faults
[pylithapp.problem.interfaces.fault_slabtop.output]
writer = pylith.meshio.DataWriterHDF5
writer.filename = output/interseismic1/interseismic-fault-slabtop.h5

[pylithapp.problem.interfaces.fault_slabbot.output]
writer = pylith.meshio.DataWriterHDF5
writer.filename = output/interseismic1/interseismic-fault-slabbot.h5

# Materials
[pylithapp.timedependent.materials.continent_litho.output]
writer.filename = output/interseismic1/interseismic-conlitho.h5

[pylithapp.timedependent.materials.continent_mantle.output]
writer.filename = output/interseismic1/interseismic-conmantle.h5

[pylithapp.timedependent.materials.ocean_litho.output]
writer.filename = output/interseismic1/interseismic-oceanlitho.h5

[pylithapp.timedependent.materials.ocean_mantle.output]
writer.filename = output/interseismic1/interseismic-oceanmantle.h5

# End of file
-------------- next part --------------
[pylithapp]

# This is not a self-contained simulation configuration file. This
# file only specifies the general parameters common to the simulations
# in this directory.

# ----------------------------------------------------------------------
# journal
# ----------------------------------------------------------------------
# Turn on some journals to show progress.
[pylithapp.journal.info]
timedependent = 1
implicit = 1
petsc = 1
solverlinear = 1
meshiocubit = 1
implicitelasticity = 1
faultcohesivekin = 1
fiatsimplex = 1
pylithapp = 1
materials = 1

# ----------------------------------------------------------------------
# mesh_generator
# ----------------------------------------------------------------------
[pylithapp.mesh_generator]
# Change the default mesh reader to the CUBIT reader.
reader = pylith.meshio.MeshIOCubit

[pylithapp.mesh_generator.reader]
filename = mesh3d_block.exo
coordsys.space_dim = 3

# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
[pylithapp.timedependent]
dimension = 3

[pylithapp.timedependent.formulation.time_step]
# Define the total time for the simulation and the time step size.
total_time = 200.0*year
dt = 10.0*year

# ----------------------------------------------------------------------
# materials
# ----------------------------------------------------------------------

dt = 10.0*year
# materials
# ----------------------------------------------------------------------
[pylithapp.timedependent]

# Set materials to an array of 4 materials:
#   'ocean_mantle'
#   'continent_litho'
#   'continent_mantle'
#   'ocean_litho'
materials = [ocean_mantle,continent_litho,continent_mantle,ocean_litho]

[pylithapp.timedependent.materials]
# Set bulk constitutive model for each material.
continent_litho = pylith.materials.ElasticIsotropic3D
ocean_litho = pylith.materials.ElasticIsotropic3D
continent_mantle = pylith.materials.MaxwellIsotropic3D
ocean_mantle = pylith.materials.MaxwellIsotropic3D

# Oceanic mantle
[pylithapp.timedependent.materials.ocean_mantle]
label = Oceanic mantle

# The id corresponds to the block number from CUBIT.
id = 1

db_properties.label = Oceanic mantle properties
db_properties.iohandler.filename = mat_oceanmantle.spatialdb

# We are doing 3D quadrature for a tetrahedron.
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

# Continental litho
[pylithapp.timedependent.materials.continent_litho]
label = Continental litho

# The id corresponds to the block number from CUBIT.
id = 2

db_properties.label = Continental litho properties
db_properties.iohandler.filename = mat_conlitho.spatialdb

# We are doing 3D quadrature for a tetrahedron.
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

# Continental mantle
[pylithapp.timedependent.materials.continent_mantle]
label = Continental mantle

# The id corresponds to the block number from CUBIT.
id = 3

db_properties.label = Continental mantle properties
db_properties.iohandler.filename = mat_conmantle.spatialdb

# We are doing 3D quadrature for a tetrahedron.
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

# Oceanic litho
[pylithapp.timedependent.materials.ocean_litho]
label = Oceanic litho

# The id corresponds to the block number from CUBIT.
id = 4

db_properties.label = Oceanic litho properties
db_properties.iohandler.filename = mat_oceanlitho.spatialdb

# We are doing 3D quadrature for a tetrahedron.
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
# Names of output files are set in stepXX.cfg. We consolidate all of the
# output settings that are common to all of the simulations here.

[pylithapp.timedependent.implicit]
# Fault friction is a nonlinear problem so we need to use the nonlinear
# solver.
solver = pylith.problems.SolverNonlinear

# Set the output to an array of 2 output managers.
# We will output the solution over the domain and the ground surface.
output = [domain,subdomain]

# Set subdomain component to OutputSolnSubset (subset of domain).
output.subdomain = pylith.meshio.OutputSolnSubset

# Domain
[pylithapp.problem.formulation.output.domain]
output_freq = time_step
time_step = 9.99999*year
writer = pylith.meshio.DataWriterHDF5

# Ground surface
[pylithapp.problem.formulation.output.subdomain]
label = groundsurf ; # Name of CUBIT nodeset for ground surface.
writer = pylith.meshio.DataWriterHDF5

# Materials
[pylithapp.timedependent.materials.continent_litho.output]
cell_filter = pylith.meshio.CellFilterAvg
output_freq = time_step
time_step = 9.99999*year
writer = pylith.meshio.DataWriterHDF5

[pylithapp.timedependent.materials.continent_mantle.output]
cell_filter = pylith.meshio.CellFilterAvg
output_freq = time_step
time_step = 9.99999*year
writer = pylith.meshio.DataWriterHDF5

[pylithapp.timedependent.materials.ocean_litho.output]
cell_filter = pylith.meshio.CellFilterAvg
output_freq = time_step
time_step = 9.99999*year
writer = pylith.meshio.DataWriterHDF5

[pylithapp.timedependent.materials.ocean_mantle.output]
cell_filter = pylith.meshio.CellFilterAvg
output_freq = time_step
time_step = 9.99999*year
writer = pylith.meshio.DataWriterHDF5

# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
# Set the solver options.

[pylithapp.timedependent.formulation]
split_fields = True
matrix_type = aij
use_custom_constraint_pc = True

[pylithapp.petsc]
fs_pc_type = fieldsplit
fs_pc_use_amat = true
fs_pc_fieldsplit_type = schur
fs_pc_fieldsplit_schur_factorization_type = full
fs_pc_fieldsplit_schur_precondition = user
fs_fieldsplit_displacement_pc_type = gamg
fs_fieldsplit_displacement_ksp_type = preonly
fs_fieldsplit_displacement_ksp_rtol = 5.0e-10
fs_fieldsplit_lagrange_multiplier_pc_type = none
fs_fieldsplit_lagrange_multiplier_ksp_type = minres
fs_fieldsplit_lagrange_multiplier_ksp_rtol = 5.0e-10

# Convergence parameters.
ksp_rtol = 1.0e-20
ksp_atol = 1.0e-12
ksp_max_it = 300
ksp_gmres_restart = 50

# Linear solver monitoring options.
ksp_monitor = true
ksp_view = true
ksp_converged_reason = true

# Nonlinear solver monitoring options.
snes_rtol = 1.0e-20
snes_atol = 1.0e-10
snes_max_it = 300
snes_monitor = true
snes_view = true
snes_converged_reason = true

# Friction sensitivity solver used to compute the increment in slip
# associated with changes in the Lagrange multiplier imposed by the
# fault constitutive model.
friction_pc_type = asm
friction_sub_pc_factor_shift_type = nonzero
friction_ksp_max_it = 25
friction_ksp_gmres_restart = 30
# Uncomment to view details of friction sensitivity solve.
friction_ksp_monitor = true
friction_ksp_view = true
friction_ksp_converged_reason = true

# PETSc summary -- useful for performance information.
# log_summary = true

# Uncomment to launch gdb when starting PyLith.
# start_in_debugger = true

# End of file