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.
My guess is that the scale of the Lagrange multiplier and/or displacement are not correct when you refine. This can lead to the solution converging for one field but the other has significant error.
In PyLith v2.2, the scale for nondimensionalizing the Lagrange multiplier is the nominal shear modulus, and the scale for nondimensionalizing the displacement is the length scale. The shear modulus does not change when you refine, but the length scale should change. For example, if you use uniform refinement by a factor of 2 or 4 then the length scale should be reduced by a factor of 2 or 4.
An additional diagnostic tool to use to verify convergence is to use the PETSc setting ksp_monitor_true_residual. This can be passed via a command line argument --petsc.ksp_monitor_true_residual or in a cfg file:
[pylithapp.petsc]
ksp_monitor_true_residual = True
If the problem is nondimensionalized properly and the mesh is high quality, then the preconditioned residual should be within a factor of about 10 of the true residual. If they differ by more than a factor of 100, then the mesh may have distorted cells, the nondimensionalization is incorrect, or the preconditioner is not working properly.
If these don’t resolve the issue, please provide all of the settings that are different between the coarse simulation and the refined simulation.
Dear Brad,
Thank you for the helpful suggestions. Following your advice, I enabled ksp_monitor_true_residual and found the ratio between the preconditioned and true residuals remains between 2 and 5 throughout the solve, which seems to indicate the nondimensionalization and preconditioner are working properly. I also confirmed that the only difference between my coarse and refined simulations is the .exo mesh file — all configuration files, solver settings, and spatial databases are identical. However, the traction results still differ significantly between the two meshes. I would greatly appreciate any further insights you might have on this.
Best regards,
Xiaopylithapp.cfg (5.4 KB)
finalslip3.txt (5.6 MB)
slurm-3803558.txt (65.1 KB)
Sihai.cfg (7.7 KB)
I also confirmed that the only difference between my coarse and refined simulations is the .exo mesh file — all configuration files, solver settings, and spatial databases are identical.
If the mesh resolution is changing, then you should be changing the length scale accordingly to match the nominal mesh size on the fault as mentioned in my earlier comments.
Some other things to check:
- Do the displacement fields closely agree between the two mesh resolutions? Do they agree only away from the fault or close to the fault as well?
- Do the stress fields in the bulk material closely agree between the two mesh resolutions? Do they agree only away from the fault or close to the fault as well?
If you are not resolving the spatial variations in slip at the coarse resolution, then you could end up with large gradients and traction changes, whereas a fine resolution mesh may.be able to better resolve the spatial variations resulting in smaller gradients and smaller traction changes. As a check, you could adjust the slip distribution to be a pyramid like shape (zero at the edges and increasing linearly into the interior of the fault). I would expect the traction changes to be consistent for the coarse mesh and fine mesh except where the slip is at a maximum (and the gradient in slip changes).