Deformation sensitivity to mesh resolution

PyLith
1

I am trying to figure out why the sensitivity of the deformation to the mesh resolution doesn’t seem to be converging.

Trying to run the simplest version of this issue: I’m running the subduction-3d example gmsh mesh, but without slip on faults and instead adding a simple surface load (Neumann BC on the top surface, -8MPa, 60km radius) with a time history: the load ramps up through time, stays steady for a few thousand years, and then is removed through time. Depending on the mesh resolution, the resulting deformation varies pretty significantly. For example, below is an image showing the output after 3000 years, for a 20km, 15km, 10km, and 5km resolution. As the mesh resolution increases, the deformation increases (max subsidence of ~7,8,9,13m). Maybe I’m just not going to a fine enough resolution (?), but it does not appear to be converging at all towards a result that is stable with subsequent mesh refinement. This is using nearest neighbor interpolation for the BC, but very similar results for linear interpolation.

image
image1336×723 102 KB

This sensitivity appears to increase over subsequent timesteps (different runs with different model resolutions start out with about the same deformation, but differ more and more as time goes on). For runs of only a few hundred years, there doesn’t appear to be a difference. Using a fully elastic model, there does not appear to be a sensitivity to mesh resolution (deformation results vary by ~1% but not systematically).

The fact that it diverges more and more through time makes me think that this might be slight differences on a single time step that just get amplified as it accumulates through time? Are there guidelines to help with figuring out an appropriate level of discrepancy that would be acceptable between meshes?

pylithapp.cfg (5.7 KB)

step04_eqcycle_forumQ.cfg (3.8 KB)

2

There are several issues that may be coming into play.

  1. Using nearest interpolation with a variable resolution mesh means the total work done varies with the mesh resolution. To obtain convergence, the load needs to be mesh resolution independent. The load should be spatially continuous, for example a smooth bump rather than a step.
  2. The load should be sufficient far enough from the boundary that the truncation of the domain does not affect the response. If the load has a spatial dimension of 60 km, then the domain boundaries probably need to be 300-600 km from the loadl.
  3. The mesh resolution may be affecting the response as a result of how well the variations in material properties are resolved. For example, for the largest discretization size, the subducting slab may be 1 cell thick.