According to the examples in the manual, to avoid large displacements after applying gravity, I set the initial stress and initial conditions of the model based on equilibrium and lithostatic pressure. After simulating for 20,000 years, I found that the model still has a significant velocity field, especially in the high-elevation regions of the Tibetan Plateau, where there is a very large subsidence velocity. Is this correct?
In my understanding, if I have set the correct initial stress, it should be sufficient to support the model, even with high elevation, so large deformation velocities should not occur.
To verify that your are minimizing the deformation due to gravity, I would first run a static simulation with a purely elastic materials to verify the deformation is minimal. If you are computing lithostatic pressure from a layered model, then you might want to test some different options for which layered structure minimizes the deformation for your more complex structure.
This is the type of simulation in which solver tolerances and nondimensionalization are especially important. If the solver tolerances are larger than the deformation, then the deformation will not decay. You may need to adjust ksp_rtol, ksp_atol, snes_rtol, and snes_atol, where ksp refers to the linear solver, snes refers to the nonlinear solver, rtol refers to the relative tolerance (tolerance relative to the initial preconditioned residual for KSP), and atol refers to the absolute tolerance (preconditioned residual for KSP). If the solvers converge due to the relative tolerance and the residuals are greater than your deformation, you should decrease the relative tolerance. If the solvers converge due to the absolute tolerance, and the residuals are greater than the deformation, you should decrease the absolute tolerance. It is usually most efficient to keep the KSP tolerances about 1-2 orders of magnitude smaller than the corresponding SNES tolerances.
I also recommend verifying that your mesh has sufficient resolution in the areas where deformation is occurring.
Thus, I suggest running a suite of simulations to verify convergence in space and time. Visualizing the stress field and ensuring it is smooth (well resolved), and then refine the mesh in regions where it is not smooth.
Thank you for your suggestions. I tested a fully elastic model and found that the displacement field of the model hardly changes, and the velocity is almost zero. Can I conclude that the initial stress I provided is correct? The significant subsidence velocity observed in the previous complex model might be related to material properties: for example, the low viscosity of the lower crust might make it difficult to support high topography?
If the displacement from the fully elastic model is small, then that is a good indication that you have reasonable initial conditions. You are correct that a low viscosity in the lower crust might make it difficult to support high topography.
