If PyLith gets “stuck” it usually happens when running in parallel, one process aborts, and the others hang. Are you running in parallel? If so, can you run in serial for 1 time step. This might allow you to see an error. Have you tried the same geometry with a coarser mesh?
I have resolved the previous configuration issue, but now when writing outputs I encounter the following error:
Fatal error. Calling MPI_Abort() to abort PyLith application.
Traceback (most recent call last):
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/apps/PetscApplication.py", line 74, in onComputeNodes
self.main(*args, **kwds)
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/apps/PyLithApp.py", line 143, in main
self.problem.finalize()
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/TimeDependent.py", line 224, in finalize
self.formulation.finalize()
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/Implicit.py", line 278, in finalize
Formulation.finalize(self)
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/Formulation.py", line 280, in finalize
output.close()
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/meshio/OutputManager.py", line 202, in close
self._close()
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/meshio/OutputManager.py", line 464, in _close
self.writer.close()
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/meshio/DataWriterHDF5.py", line 83, in close
xdmf.write(ModuleDataWriterHDF5.hdf5Filename(self), verbose=False)
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/meshio/Xdmf.py", line 75, in write
self.h5 = h5py.File(filenameH5, "r")
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/h5py-2.9.0-py2.7-linux-x86_64.egg/h5py/_hl/files.py", line 394, in __init__
swmr=swmr)
File "/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/h5py-2.9.0-py2.7-linux-x86_64.egg/h5py/_hl/files.py", line 170, in make_fid
fid = h5f.open(name, flags, fapl=fapl)
File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
File "h5py/h5f.pyx", line 85, in h5py.h5f.open
IOError: Unable to open file (truncated file: eof = 3756584, sblock->base_addr = 0, stored_eof = 4240840)
application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0
/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/bin/nemesis: mpirun: exit 255
/home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/bin/pylith: /home/zhengyong/users/xiaoyang/xiaoy/xiaoy/Pylith/pylith-2.2.2-linux-x86_64/bin/nemesis: exit 1
I have searched the community for solutions but could not find an answer. What is causing this?
The backtrace shows that it is trying to create the Xdmf file from the HDF5 file. It looks like h5py is having trouble accessing the HDF5 file and thinks it is corrupted. My guess is that the file is corrupt. It can get corrupt if multiple processes write to it and it doesn’t get closed or one process aborts before the write is complete or you run out of disk space. If you are running in parallel, then I recommend using the DataWriterHDF5External which uses external raw datasets written using MPI I/O and just writes an HDF5 with metadata. It produces a lot of files, but it is more robust.
Thank you very much for your help，I have resolved the issue! Wishing you all the best.
Hi Dr. Aagaard,
I have a related question regarding using viscoelastic material to compute ground deformation due to a back slip rate in Cascadia. The 3D model (360 km in z direction) I used and a simulation result (by applying unit slip impulse on one triangular fault patch, red dots in figure 1) are attached. The rigidity of the mantle and crust is set as 64GPa and 48GPa, and continental mantle has viscosity of 10^19 Pa.S and oceanic mantle has viscosity of 10^20 Pa.S, with elastic oceanic and continental crust. Theoretically, the Maxwell relation time of the continental mantle should be around 5 yrs (viscosity/rigidity)? But after 300 yrs of simulation, the ground surface deformation still did not get stable results, especially for the z direction, which keeps decreasing near 300 yr (subplot 3, figure 2). The ground surface horizontal and vertical deformation near 300 yr look not very correct to me (subplots 1,2, figure 2). I wonder what could cause this problem. Your experience and judgement could be very helpful for me to find out the problem fast. Thanks in advance.
Sincerely,
Qingjun
3D_model.pdf (391.0 KB)
You applied a static fault slip impulse at time 0 and there is no other time dependence to the forcing? If the deformation does not decrease as expected, then it could be your materials are not in equilibrium with the boundaries, so something else besides the slip is driving the deformation. What are the boundary conditions on the lateral and bottom sides of the domain? If you don’t apply any fault slip, do you get deformation?
Dear Dr. Aagaard,
Thank you for your response. I applied a static reverse fault slip impulse at time 9.9999 yr, and before this time the deformation is totally 0, which could be found in figure 2 (subplot 3). Thus, the deformation should be zero without any fault slip? I used Dirichlet roller boundary conditions on four side boundaries and the bottom boundary, keeping the top surface free:
-
xneg/xpos：u_x=0
-
yneg/ypos：u_y=0
-
zneg：u_z=0
When I tried to change two side boundaries, for example the northern and southern side boundaries, into free boundaries, Pylith gives error message during simulation.
By the way, I combine all input parameters I used into the file combined_input_files.txt for your reference, as attached, which includes information from files: pylithapp.cfg, step02_visco_reverse_318.cfg and mat_viscoelastic.cfg.combined_input_files.txt (14.4 KB)
Hi Dr. Aagaard,
Now I start to consider the model mesh size could be the potential problem that causes long-term deformation instability. Though my current model uses basis_order=2 during calculation and the result has been tested to be consistent with the analytical results for the elastic GF simulation. But, it might be a problem for viscoelastic simulation, as the strain and stress deformation is accumulating through time, which may amplify the error? I found the Pylith will crush when total element number reaches about 0.8 million using basis_order=2, I wonder is there a way to fix this problem to make it run successfully even with many elements? Also the coseismic simulation seems different with GF calculation, which could select to do element based calculation by setting slip_basis_order =0, but coseismic simulation seems do not have this option and can be only node based? By the way, Pylith 4.0 is not able to do time adaptive calculation but 2.0 version can, right? Thanks in advance and I really look forward to your reply.
Best,
Qingjun
Now I start to consider the model mesh size could be the potential problem that causes long-term deformation instability. Though my current model uses basis_order=2 during calculation and the result has been tested to be consistent with the analytical results for the elastic GF simulation. But, it might be a problem for viscoelastic simulation, as the strain and stress deformation is accumulating through time, which may amplify the error?
A basis order of 2 should reduce numerical errors, not increase them. Is the slip in or very close to a viscoelastic bulk rheology? If so, some shear stress may persist, and the viscoelastic bulk rheology will continue to relax, but the strain rate will not decrease to zero.
I found the Pylith will crush when total element number reaches about 0.8 million using basis_order=2, I wonder is there a way to fix this problem to make it run successfully even with many elements?
Do you get any sort of error message when PyLith crashes? You may be running out of memory.
Also the coseismic simulation seems different with GF calculation, which could select to do element based calculation by setting slip_basis_order =0, but coseismic simulation seems do not have this option and can be only node based?
I do not understand what you are asking. Please provide more details, including small sections of the parameter file if necessary. You can control the discretization of the slip fields. We strongly recommend using a basis order of 1 for Green’s functions and forward coseismic slip simulations to maintain continuity in the slip.
By the way, Pylith 4.0 is not able to do time adaptive calculation but 2.0 version can, right?
In PyLith 3 and later, we use PETSc time-stepping algorithms. These include adaptive schemes. We have added example use cases and parameter files in the PyLith development version. Using adaptive time stepping does not require any changes to the PyLith code, so you can use these features now. Refer to Adaptive Time-Stepping Options section in the PyLith manual for more information.