Problem with reading the Trelis model in Pylith

fdiaz · April 9, 2019, 4:35pm

Dear CIG team,
I hope you are well.

I would like to estimate topographic stress and how it affects subduction and the generation of earthquakes and accumulation of energy in this zone, more specific in the Megathrust.

For this estimation, I would like to use Pylith (subduction model in 2d).

I created my mode of subduction in Trelis, with information from continental topographic elevation, oceanic Slab and continental Moho.

However, when I try to run the step 01 in the subduction model of Pylith, it shows me an error with the reading of the model. As seen in the attached picture.

Pylith%20error

Please if you are able to help me in my research, I would be very grateful.

Best Regards.

baagaard · April 9, 2019, 7:17pm

PyLith is trying to read the number of node set groups (num_node_sets) in the Exodus file. PyLith uses the nodesets in the Exodus file to associated pieces of the mesh with boundary conditions. Did you create nodesets in Trelis for your boundary conditions?

fdiaz · April 11, 2019, 2:19am

Thank you for you helping.
I create the Nodesets in my model, but when I try to run my model, anyway its still have problem. This is my new error.

willic3 · April 11, 2019, 4:23am

Dear Felipe,

It seems that there may be a dimension wrong somewhere. Can you provide us with the .cfg files you are using?

Cheers,

Charles

fdiaz · April 11, 2019, 3:12pm

Dear Charles,

[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 = mesh_tri3.exo
filename = modelo1.e
coordsys.space_dim = 2

----------------------------------------------------------------------

problem

----------------------------------------------------------------------

[pylithapp.timedependent]
dimension = 2

[pylithapp.timedependent.formulation.time_step]

Define the total time for the simulation and the time step size.

total_time = 200.0year
dt = 1.0year

----------------------------------------------------------------------

materials

----------------------------------------------------------------------

[pylithapp.timedependent]

Set materials to an array of 4 materials:

‘continent_crust’

‘continent_mantle’

‘ocean_crust’

‘ocean_mantle’

materials = [continent_crust,continent_mantle,ocean_crust,ocean_mantle]

[pylithapp.timedependent.materials]

Set bulk constitutive model for each material.

continent_crust = pylith.materials.ElasticPlaneStrain
ocean_crust = pylith.materials.ElasticPlaneStrain
continent_mantle = pylith.materials.MaxwellPlaneStrain
ocean_mantle = pylith.materials.MaxwellPlaneStrain

Continental crust -----------------------------

[pylithapp.timedependent.materials.continent_crust]
label = Continental crust

The id corresponds to the block number from CUBIT.

id = 1

db_properties.label = Continental crust properties
db_properties.iohandler.filename = mat_concrust.spatialdb

We are doing 2D quadrature for a triangle.

quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

Continental mantle --------------------

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

The id corresponds to the block number from CUBIT.

id = 2

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

We are doing 2D quadrature for a triangle.

quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

Oceanic crust --------------------

[pylithapp.timedependent.materials.ocean_crust]
label = Oceanic crust

The id corresponds to the block number from CUBIT.

id = 3

db_properties.label = Oceanic crust properties
db_properties.iohandler.filename = mat_oceancrust.spatialdb

We are doing 2D quadrature for a triangle.

quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

Oceanic mantle --------------------

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

The id corresponds to the block number from CUBIT.

id = 4

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

We are doing 2D quadrature for a triangle.

quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

----------------------------------------------------------------------

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.formulation]

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
vertex_data_fields = [displacement,velocity]

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_crust.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
cell_data_fields = [stress,total_strain,viscous_strain]

[pylithapp.timedependent.materials.ocean_crust.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
cell_data_fields = [stress,total_strain,viscous_strain]

----------------------------------------------------------------------

PETSc

----------------------------------------------------------------------

Set the solver options.

[pylithapp.petsc]

Preconditioner settings.

pc_type = asm
sub_pc_factor_shift_type = nonzero

Convergence parameters.

ksp_rtol = 1.0e-8
ksp_atol = 1.0e-12
ksp_max_it = 400
ksp_gmres_restart = 100

Linear solver monitoring options.

ksp_monitor = true
#ksp_view = true
ksp_converged_reason = true
ksp_error_if_not_converged = true

Nonlinear solver monitoring options.

snes_rtol = 1.0e-8
snes_atol = 1.0e-12
snes_max_it = 100
snes_monitor = true
#snes_view = true
snes_converged_reason = true
snes_error_if_not_converged = true

PETSc summary – useful for performance information.

#log_view = true

Uncomment to launch gdb when starting PyLith.

start_in_debugger = true

End of file

And the Step001.

-- Config -- (syntax highlighting)

----------------------------------------------------------------------

Brad T. Aagaard, U.S. Geological Survey

Charles A. Williams, GNS Science

Matthew G. Knepley, University of Chicago

This code was developed as part of the Computational Infrastructure

for Geodynamics (http://geodynamics.org).

Copyright © 2010-2017 University of California, Davis

See COPYING for license information.

----------------------------------------------------------------------

PROBLEM DESCRIPTION

----------------------------------------------------------------------

This simulation involves coseismic slip between the continental

crust and top of the subducting oceanic crust. The slip also

extends down into the top of the mantle below the continental

crust.

----------------------------------------------------------------------

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.

To run the simulation:

pylith step01.cfg

Output will be directed to directory output.

----------------------------------------------------------------------

Output of parameters used and simulation progress.

----------------------------------------------------------------------

[pylithapp]
dump_parameters.filename = output/step01-parameters.json
problem.progress_monitor.filename = output/step01-progress.txt

----------------------------------------------------------------------

problem

----------------------------------------------------------------------

[pylithapp.timedependent.formulation.time_step]
total_time = 0.0year
dt = 5.0year

----------------------------------------------------------------------

boundary conditions

----------------------------------------------------------------------

[pylithapp.timedependent]

Set bc to an array of 4 boundary conditions:

‘boundary_east_crust’

‘boundary_east_mantle’

‘boundary_west’

‘boundary_bottom_mantle’

bc = [boundary_east_crust,boundary_east_mantle,boundary_west,boundary_bottom_mantle]

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.

East boundary (crust)

[pylithapp.timedependent.bc.boundary_east_crust]
bc_dof = [0]
label = bndry_east_crust
db_initial.label = Dirichlet BC on east boundary (crust)

East boundary (mantle)

[pylithapp.timedependent.bc.boundary_east_mantle]
bc_dof = [0]
label = bndry_east_mantle
db_initial.label = Dirichlet BC on east boundary (mantle)

West boundary

[pylithapp.timedependent.bc.boundary_west]
bc_dof = [0]
label = bndry_west
db_initial.label = Dirichlet BC on west boundary

Bottom boundary (mantle)

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

----------------------------------------------------------------------

faults

----------------------------------------------------------------------

[pylithapp.timedependent]
interfaces = [fault]

Set the type of fault interface condition.

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

Set the parameters for the fault interface condition.

[pylithapp.timedependent.interfaces.fault]

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

label = fault_slabtop

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.timedependent.interfaces.fault.eq_srcs.rupture.slip_function]
slip.iohandler.filename = fault_slip_coseismic.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 = [0.0*year]

----------------------------------------------------------------------

output

----------------------------------------------------------------------

Domain

[pylithapp.problem.formulation.output.domain]
writer.filename = output/step01.h5

Ground surface

[pylithapp.problem.formulation.output.subdomain]
writer.filename = output/step01-groundsurf.h5

Fault

[pylithapp.problem.interfaces.fault.output]
writer = pylith.meshio.DataWriterHDF5
writer.filename = output/step01-fault.h5

Materials

[pylithapp.timedependent.materials.continent_crust.output]
writer.filename = output/step01-concrust.h5

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

[pylithapp.timedependent.materials.ocean_crust.output]
writer.filename = output/step01-oceancrust.h5

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

End of file

sorry but I can not upload files because I’m new user

cheers.

baagaard · April 11, 2019, 3:55pm

The error suggests there are likely still problems with your mesh and/or .cfg parameter files. Have your tried generating a mesh without faults and making sure PyLith can read it? Have you run the Trelis mesh quality checks to make sure you do not have distorted cells?

I know Charles suggested sending the .cfg files, but the JSON parameter file generated when PyLith starts is much more informative because it shows all parameters, including the defaults. You will save yourself a lot of time by becoming familiar with it and learning how to diagnose problems using it.

Please use the pylith_parameters tool (it is included in the binary or use the online version at https://geodynamics.github.io/pylith_parameters/) to look at the settings and double check to make sure the dimensions and coordinate system match what you expect. I also suggest comparing the JSON file from your simulation with the one generated for the example. See Section 4.10 of the PyLith manual for more information.

If you still cannot diagnose the problem after doing all of the above, then send a diagram illustrating the boundary value problem you are trying to solve, your Trelis journal files used to create the mesh, and your JSON parameter file along with and how you have setup your mesh.

willic3 · April 11, 2019, 7:06pm

Yes, Brad is correct that the JSON file is much more useful in diagnosing the problem.

Cheers,

Charles

Charles Williams I Geodynamic Modeler
GNS Science I Te
Pῡ Ao

1 Fairway Drive, Avalon 5010, PO Box 30368, Lower Hutt 5040, New Zealand

Ph 0064-4-570-4566 I Mob 0064-22-350-7326 I Fax 0064-4-570-4600 http://www.gns.cri.nz/ I **Email: **C.Williams@gns.cri.nz

Topic		Replies	Views
Topographic stress with Pylith PyLith	1	509	January 28, 2019
My 2D subduction static stress does not converge PyLith generate-stress	4	87	February 27, 2024
Using PyLith to calculate the subduction zone Green's functions PyLith	10	46	September 12, 2024
PyLith v3.0.2 released PyLith	0	227	September 19, 2022
Problem Running Constarined Fault Zone Pylith2 PyLith mesh	5	311	September 30, 2022