Hi,
The output of the Pylith example is almost displacement. I would like to build a 3D fault model to simulate the rate of interseismic deformation at the surface. I have tried several times without success. Can you please help me if you know how to do it?

To compute the rate of interseismic deformation, you need to create a model that has variation in deformation over time. For example, you could impose velocity boundary conditions with a linear elastic rheology. You can compute an average rate by dividing the displacement by time. If you include viscoelastic rheologies, then you may need to run your model for many time steps to find the steady state solution. We might be able to provide more help if you provide a diagram of the boundary value problem you are trying to solve.

It’s unclear exactly what you want. If you want the velocity field at the surface you would need to do a little post processing. PyLith outputs the total displacement field, so to get velocities you would need to compute the displacement differences between one time step and the next, and then divide this by the time step size. This is easy to do with a simple Python script (using h5py and numpy). If you want the strain rate field you would need to do something similar with the total strain output from PyLith.

Here are all the files for my simulation using Pylith, including the finite element model, Pylith control scripts, rheological parameter settings, output files, etc. We want to output the surface deformation rate, and after trying for some time we still can’t succeed.

I think your may be having trouble because you include the velocity field in the solution. The quasistatic formulation does not use the velocity field in the solution; it is used in the dynamic (include inertia). See equation 114 at Quasistatic — PyLith 4.0.0 documentation. As Charles and I have suggested, you should use the displacement field and you can divide by time to get a velocity or rate of deformation.