[CIG-SHORT] Speed up nonlinear solver for faults with friction

[Jiangzhi Chen]

Hi Brad,

    Is there some way to speed up the nonlinear solver when calculating 
the interseismic fault deformation with friction on the fault? In my 
current simulation, when I use kinematic constraints on faults, the 
whole calculation takes about 10 hours, but if I use friction instead of 
fault slip, the calculation seems to run forever (it has been running on 
server during the whole AGU week but not even the first time step is 
completed). I attached the configuration file.

cheers,
Jiangzhi
-------------- next part --------------
[pylithapp]
# initialize_only = True
# ----------------------------------------------------------------------
# PROBLEM DESCRIPTION
# ----------------------------------------------------------------------
#
# This simulation involves aseismic creep along the interfaces between
# the subducting oceanic lithosphere and the mantle. The subduction is
# modeled by a velocity condition at the west boundary. Upper interface
# slips follow a distribution with overall NE slip rate 6 cm/year. The
# slip time is 40 years for quick test. A locked zone exists until 50 km,
# and the continental plate moves 1 cm/year NE (prescribed as surface BC).
# The top interface has uniform friction
#
# ----------------------------------------------------------------------
# RUNNING THE SIMULATION
# ----------------------------------------------------------------------
#
# 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 specified 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_litho,boundary_west_mantle,boundary_east,boundary_bottom_mantle,boundary_side,ocean_surf,con_surf]

# The label corresponds to the name of the nodeset in CUBIT.

# west boundary (lithosphere)
[pylithapp.timedependent.bc.boundary_west_litho]
bc_dof = [0,1,2]
label = bndry_west_litho
db_rate = spatialdata.spatialdb.UniformDB
db_rate.label = Dirichlet rate BC on west lithosphere boundary
db_rate.values = [displacement-rate-x,displacement-rate-y,displacement-rate-z,rate-start-time]
db_rate.data = [4.24*cm/year,4.24*cm/year,0*cm/year,0*year]

# 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)

# ocean plate surface
[pylithapp.timedependent.bc.ocean_surf]
bc_dof = [0,1,2]
label = oceansurf
db_rate = spatialdata.spatialdb.UniformDB
db_rate.label = Dirichlet rate BC on ocean plate surface
db_rate.values = [displacement-rate-x,displacement-rate-y,displacement-rate-z,rate-start-time]
db_rate.data = [4.24*cm/year,4.24*cm/year,0*cm/year,0*year]

# continental plate surface
[pylithapp.timedependent.bc.con_surf]
bc_dof = [0,1,2]
label = consurf
db_rate = spatialdata.spatialdb.SimpleDB
db_rate.label = Dirichlet rate BC on continental plate surface
db_rate.iohandler.filename = continent_bc.spatialdb

# front and back sides, simplified treatment
[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]

[pylithapp.timedependent.implicit]
solver = pylith.problems.SolverNonlinear

# 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-7
# We must define the quadrature information for fault cells.
# The fault cells are 2D (surface).
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

# # Kinematic fault, removed when using friction
# # 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_slabtop.eq_srcs.rupture]
# slip_function = pylith.faults.ConstRateSlipFn

# # The slip time and final slip are defined in spatial databases.
# [pylithapp.timedependent.interfaces.fault_slabtop.eq_srcs.rupture.slip_function]
# slip_rate.iohandler.filename = sliprate_upper_locked.spatialdb
# # slip_rate.query_type = linear
# slip_rate.label = Slip of upper interface

# # 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]

# Set static friction model parameters using a uniform DB. Set the
# static coefficient of friction to 0.6 and cohesion to 0.0 Pa.
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]

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

# 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.iohandler.filename = sliprate_lower.spatialdb
# slip_rate.query_type = linear
slip_rate.label = Slip of lower interface

# 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]

# ----------------------------------------------------------------------
# PETSc settings
# ----------------------------------------------------------------------
# NOTE: There are additional settings specific to fault friction.
[pylithapp.petsc]

# Friction sensitivity solve 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

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

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

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

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

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

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

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

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

# End of file