Dear Baagaard ,
I am running quasi-static simulations of coseismic slip-induced stress changes (traction) on a planar fault using PyLith 2.2.2. I have encountered a puzzling issue: the computed traction distribution changes substantially when I refine the mesh, and the coarser mesh actually appears to produce more physically reasonable results.
My guess is that the differences are related to nondimensionalization or solver tolerances.
Are you using fault friction? Do the simulations have similar results at the first time step? Do the solutions diverge slowly from each other or suddenly?
We did not account for fault friction because we focused on the coseismic slip effect. Therefore, we set the total simulation time to 0 s, and the specific parameter settings are shown below.
Sihai.cfg (7.7 KB)
pylithapp.cfg (5.3 KB)
You seem to be mixing solver settings for v4.x and later and those with v2.2. I suspect this is causing the solver to give results with substantial errors. If you are using prescribed slip, I strongly recommend using PyLith v5.0.1. The latest release gives good parallel scalability with good default solver parameters, so you should only need to set the length scale and possibly the displacement scale without adjusting the tolerances or preconditioner settings.
If you want to continue using v2.2.2, I recommend starting with the solver settings in share/settings/solver_fault_fieldsplit.cfg, the recommended tolerances (ksp_atol, ksp_rtol, snes_rtol, and snes_atol), and a length scale equal to the discretization size on the fault.
Dear Brad,
Thank you for the clarification. I am indeed still using PyLith v2.2.2. I have already incorporated the length scale setting as you suggested:
pylithapp.timedependent.normalizer
length_scale = 0.2*km
Additionally, I have tightened the tolerances to 1e-18. However, the issue persists despite these adjustments. I would greatly appreciate any further guidance you might have on resolving this.
Tolerances of 1.0e-18 are too small in a case like this. Double precision will likely only provide accuracy to around 1.0e-15. The tolerances are relative to the nondimensionlized displacement and Lagrange multipliers in the solution field, so 1.0e-8 for the SNES absolute tolerance is reasonable (0.1 mm relative to the length scale 0f 0.2 km).
Looking at your configuration files, the interpolation scheme for the final slip spatial database is nearest, so the resolution of your spatial database is controlling the spatial variation of the slip, not the resolution of your mesh. You should consider how best to represent the slip distribution- nearest interpolation implies uniform values between points in the spatial database whereas linear interpolation implies linear variations. Then, select a mesh resolution that can represent the slip distribution. With a linear interpolation of the slip distribution, refining the mesh should lead to reduced discretization error.
Thank you very much for your detailed and helpful suggestions.
I have followed your advice and made the following changes: (1) adjusted the tolerances to more reasonable values, (2) switched the interpolation scheme in the final slip spatial database from nearest to linear, and (3) refined the mesh resolution to better represent the slip distribution.
However, the results remain the same as before. Interestingly, the coarse mesh now gives improved results, but the refined mesh still produces the same outcome as previously.