[pylithapp]

# ----------------------------------------------------------------------
# PROBLEM DESCRIPTION
# ----------------------------------------------------------------------

#
# Run step16.cfg tutorial on a larger mesh 200kmx200kmx100km

[pylithapp.mesh_generator]
reader = pylith.meshio.MeshIOCubit
[pylithapp.mesh_generator.reader]
#filename = mesh/mesh_hex8_10000m_mod.exo
filename = mesh/mesh_hex8_5000m_mod.exo

# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
[pylithapp.timedependent]
# Set bc to an array of 5 boundary conditions: 'x_pos','x_neg', 'y_pos',
# 'y_neg', and 'z_neg'.
bc = [x_pos,x_neg,y_pos,y_neg,z_neg]

# Set gravity field (default is None)
gravity_field = spatialdata.spatialdb.GravityField

[pylithapp.timedependent.implicit]
# 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

# Change time stepping algorithm from uniform time step, to adaptive
# time stepping.
time_step = pylith.problems.TimeStepAdapt

# Change the total simulation time to 200 years, and set the maximum time
# step size to 10 years.
[pylithapp.timedependent.implicit.time_step]
total_time = 200.0*year
max_dt = 10.0*year
stability_factor = 1.0 ; use time step equal to stable value from materials

# ----------------------------------------------------------------------
# materials
# ----------------------------------------------------------------------
# Change material type of lower crust to generalized Maxwell viscoelastic.
[pylithapp.timedependent]
materials.lower_crust = pylith.materials.GenMaxwellIsotropic3D

# Provide a spatial database from which to obtain property values.
# Since there are additional properties and state variables for the
# generalized Maxwell model, we explicitly request that they be output.
# Properties are named in cell_info_fields and state variables are named in
# cell_data_fields.
[pylithapp.timedependent.materials.lower_crust]
db_properties.iohandler.filename = spatialdb/mat_genmaxwell.spatialdb
output.cell_info_fields = [density,mu,lambda,shear_ratio,maxwell_time]
output.cell_data_fields = [total_strain,stress,viscous_strain_1,viscous_strain_2,viscous_strain_3]

# ----------------------------------------------------------------------
# initial stresses
# ----------------------------------------------------------------------
# We must specify initial stresses for each material.
# We provide a filename for the spatial database that gives the stresses,
# and we change the query_type from the default 'nearest' to 'linear'.
[pylithapp.timedependent.materials.upper_crust]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in upper crust
db_initial_stress.iohandler.filename = spatialdb/initial_stress.spatialdb
db_initial_stress.query_type = linear

[pylithapp.timedependent.materials.lower_crust]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in lower crust
db_initial_stress.iohandler.filename = spatialdb/initial_stress.spatialdb
db_initial_stress.query_type = linear

# ----------------------------------------------------------------------
# boundary conditions
# ----------------------------------------------------------------------
# We use the default ZeroDispDB for all boundaries, since we just want
# zero displacements in the direction normal to each face.

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

[pylithapp.timedependent.bc.x_pos]
bc_dof = [0]
label = face_xpos
db_initial.label = Dirichlet BC on +x

[pylithapp.timedependent.bc.x_neg]
bc_dof = [0]
label = face_xneg
db_initial.label = Dirichlet BC on -x

[pylithapp.timedependent.bc.y_pos]
bc_dof = [1]
label = face_ypos
db_initial.label = Dirichlet BC on +y

[pylithapp.timedependent.bc.y_neg]
bc_dof = [1]
label = face_yneg
db_initial.label = Dirichlet BC on -y

[pylithapp.timedependent.bc.z_neg]
bc_dof = [2]
label = face_zneg
db_initial.label = Dirichlet BC on -z

# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
# Give basename for VTK domain output of solution over domain.
[pylithapp.problem.formulation.output.domain]
# We specify that output occurs in terms of a given time frequency, and
# ask for output every 40 years. The time stamps of the output files are
# in years (rather than the default of seconds), and we give a format for
# the time stamp.
output_freq = time_step
time_step = 40.0*year
writer.filename = output1/step16.vtk
writer.time_format = %04.0f
writer.time_constant = 1.0*year

# Give basename for VTK domain output of solution over ground surface.
[pylithapp.problem.formulation.output.subdomain]
# Name of nodeset for ground surface.
label = face_zpos
# We keep the default output frequency behavior (skip every n steps), and
# ask to skip 0 steps between output, so that we get output every time step.
skip = 0
writer.filename = output1/step16-groundsurf.vtk
writer.time_format = %04.0f
writer.time_constant = 1.0*year

# Give basename for VTK output of upper_crust state variables.
[pylithapp.timedependent.materials.upper_crust.output]
# Average values over quadrature points.
cell_filter = pylith.meshio.CellFilterAvgMesh
output_freq = time_step
time_step = 40.0*year
writer.filename = output1/step16-upper_crust.vtk
writer.time_format = %04.0f
writer.time_constant = 1.0*year

# Give basename for VTK output of lower_crust state variables.
[pylithapp.timedependent.materials.lower_crust.output]
# Average values over quadrature points.
cell_filter = pylith.meshio.CellFilterAvgMesh
output_freq = time_step
time_step = 40.0*year
writer.filename = output1/step16-lower_crust.vtk
writer.time_format = %04.0f
writer.time_constant = 1.0*year