[pylithapp]

# ----------------------------------------------------------------------
# journal
# ----------------------------------------------------------------------
# Turn on some journals to show progress.
[pylithapp.journal.info]
timedependent = 1
implicit = 1
petsc = 1
meshimporter = 1
meshiocubit = 1
elasticityimplicit = 1
faultcohesivekin = 1

# ----------------------------------------------------------------------
# mesh_generator
# ----------------------------------------------------------------------
[pylithapp.mesh_generator]
# Change the default mesh reader to the CUBIT reader.
reader = pylith.meshio.MeshIOCubit
reader.filename = mesh_tri3.exo
# Reorder cells and vertices using Reverse Cuthill-McKee algorithm
reorder_mesh = True

[pylithapp.mesh_generator.reader]
coordsys.space_dim = 2

# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
[pylithapp.problem]
dimension = 2

normalizer.length_scale = 1.0*km
# Use relaxation time of 1.0*s so timesteps in output are in seconds
normalizer.relaxation_time = 1.0*s

[pylithapp.problem.formulation.time_step]
# Define the total time for the simulation and the time step size.
total_time = 30.0*s
dt = 1.0*s

# ----------------------------------------------------------------------
# materials
# ----------------------------------------------------------------------
# Set array of materials to one material
[pylithapp.problem]
materials = [elastic]
materials.elastic = pylith.materials.ElasticPlaneStrain

[pylithapp.problem.materials.elastic]
label = Elastic material
# The id corresponds to the block id in CUBIT.
id = 1
db_properties.iohandler.filename = mat_elastic.spatialdb
db_properties.label = Elastic properties

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

# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
# Create array of BCs, one for each side of the domain
[pylithapp.problem]
bc = [x_neg,x_pos,y_neg,y_pos]
bc.x_pos = pylith.bc.Neumann ; Change BC type to Neumann
bc.x_neg = pylith.bc.Neumann ; Change BC type to Neumann
bc.y_pos = pylith.bc.Neumann ; Change BC type to Neumann
bc.y_neg = pylith.bc.Neumann ; Change BC type to Neumann

[pylithapp.timedependent.bc.x_pos]
label = bndry_east ; Name of CUBIT nodeset for +x boundary
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC +x edge
db_initial.iohandler.filename = x_tract.spatialdb
db_initial.query_type = nearest
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
quadrature.cell.quad_order = 2

[pylithapp.timedependent.bc.x_neg]
label = bndry_west ; Name of CUBIT nodeset for -x boundary
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC -x edge
db_initial.iohandler.filename = x_tract.spatialdb
db_initial.query_type = nearest
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
quadrature.cell.quad_order = 2

[pylithapp.timedependent.bc.y_pos]
label = face_top ; Name of CUBIT nodeset for +y boundary
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC +y edge
db_initial.iohandler.filename = y_tract.spatialdb
db_initial.query_type = nearest
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
quadrature.cell.quad_order = 2

[pylithapp.timedependent.bc.y_neg]
label = face_bot ; Name of CUBIT nodeset for -y boundary
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC -y edge
db_initial.iohandler.filename = y_tract.spatialdb
db_initial.query_type = nearest
quadrature.cell = pylith.feassemble.FIATLagrange
quadrature.cell.dimension = 1
quadrature.cell.quad_order = 2

#[pylithapp.problem.bc.x_pos]
#label = bndry_east
#bc_dof = [0]
#db_initial = pylith.bc.ZeroDispDB
#db_initial.label = Dirichlet BC on +x

#[pylithapp.problem.bc.x_neg]
#label = bndry_west
#bc_dof = [0]
#db_initial = pylith.bc.ZeroDispDB
#db_initial.label = Dirichlet BC on -x

#[pylithapp.problem.bc.y_pos]
#label = face_top
#bc_dof = [1]
#db_initial = pylith.bc.ZeroDispDB
#db_initial.label = Dirichlet BC on +y
#
#[pylithapp.problem.bc.y_neg]
#label = face_bot
#bc_dof = [1]
#db_initial = pylith.bc.ZeroDispDB
#db_initial.label = Dirichlet BC on -y

# ----------------------------------------------------------------------
# fault
# ----------------------------------------------------------------------
# Create array of faults with one fault.
[pylithapp.problem]
interfaces = [fault]

[pylithapp.problem.interfaces]
fault = pylith.faults.FaultCohesiveKin

#center part of fault
[pylithapp.problem.interfaces.fault]
id = 100
# The label corresponds to the name of the nodeset in CUBIT.
label = fault

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

# The slip time and final slip are defined in spatial databases.
[pylithapp.problem.interfaces.fault.eq_srcs.rupture.slip_function]
slip.iohandler.filename = eqslip.spatialdb
slip.query_type = linear
slip.label = Final slip

# 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 = [1.0*s]


# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
[pylithapp.problem.formulation]
output = [domain,points]
output.points = pylith.meshio.OutputSolnPoints

# Domain
[pylithapp.problem.formulation.output.domain]
writer = pylith.meshio.DataWriterHDF5Mesh
writer.filename = output/eqsim.h5

# Points
[pylithapp.problem.formulation.output.points]
coordsys.space_dim = 2
coordsys.units = km

writer = pylith.meshio.DataWriterHDF5Mesh
writer.filename = output/eqsim-points.h5
reader.filename = output_points.txt

# Fault
[pylithapp.problem.interfaces.fault.output]
writer = pylith.meshio.DataWriterHDF5SubSubMesh
writer.filename = output/eqsim-fault.h5

# Materials
[pylithapp.problem.materials.elastic.output]
cell_info_fields = []
cell_data_fields = []
writer = pylith.meshio.DataWriterHDF5Mesh
writer.filename = output/eqsim-statevars.h5

# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
[pylithapp.petsc]

ksp_type = gmres
ksp_gmres_restart = 100
pc_type = asm
sub_pc_factor_shift_type = nonzero

ksp_rtol = 1.0e-8
ksp_atol = 1.0e-10
ksp_max_it = 500

ksp_monitor = true
ksp_view = true
ksp_converged_reason = true

log_summary = true

# End of file