Dear Aagaard,
I’m developing a simulation to model the evolution of pore pressure following an earthquake using pylith 4.1.3. From my understanding, the pore pressure example in PyLith only allows the input of initial displacement conditions. As a result, the actual earthquake process is not included in the simulation.
To incorporate co-seismic deformation into the post-seismic diffusion process, I modified the example magma-2d as a base. I first ensured that step01_inflation runs correctly under a typical model.
However, when I changed the solution from:
solution = pylith.problems.SolnDispPresTracStrain
to:
solution = pylith.problems.SolnDispPresTracStrainLagrange
the simulation begin fails. If I only change the solution (without adding a fault to the model) and run pressure diffusion, the simulation fails to converge and eventually stops at the maximum number of iterations ( Figure 1).
If I then follow the coseismic example and add a fault, the simulation crashes immediately with a zero pivot error during LU factorization (Figure 2). I’m unable to locate the cause of the singular matrix or identify which part of the setup is responsible.
Am I on the right track by trying to use SolnDispPresTracStrainLagrange and add a fault in order to combine coseismic slip with postseismic pore pressure simulation? Is this the proper method?
If the information I provided is not sufficient to diagnose the problem, I would be happy to provide more details or files. I would greatly appreciate any advice or suggestions.
Best regards,
If you want to combine poroelasticity with a fault, then you are correct in using the SolnDispPresTracStrainLagrange solution field container. Some things to note:
- The solver will fail if you have the Lagrange multiplier field in the solution but do not have a fault. I will create a GitHub issue to note that this causes the solver to fail, so that we can generate an informative error message.
- I don’t think we have an example or test that uses poroelasticity with a fault. Robert Walker and Matt Knepley created an example for a paper in GJI based on strikeslip-2d. However, we have changed some solver defaults since then, so the settings they use may not be appropriate.
I think the following combination of the elasticity+fault solver settings and poroelasticity solver settings should work (but we haven’t tested these yet):
[pylithapp.problem.petsc_defaults]
# Turn off default solver settings
solver = False
[pylithapp.petsc]
# These two settings reorder the degrees of freedom so that we can use the variable point block Jacobi preconditioner
dm_reorder_section = True
dm_reorder_section_type = cohesive
# These are the poroelasticity solver settings
pc_type = fieldsplit
pc_fieldsplit_type = multiplicative
pc_fieldsplit_0_fields = 2
pc_fieldsplit_1_fields = 1
pc_fieldsplit_2_fields = 0
fieldsplit_trace_strain_pc_type = bjacobi
fieldsplit_pressure_pc_type = bjacobi
# These are the elasticity+fault solver settings
fieldsplit_displacement_pc_type = gamg
fieldsplit_displacement_ksp_type = gmres
mg_fine_pc_type = vpbjacobi
We will try to find time to set up appropriate automated tests and an example for poroelasticity with a fault in the next few weeks.
Dear aagaard,
Thank you for your prompt reply. I have tried the PETSc solver configuration you suggested, but it doesn’t seem to resolve the issue. I am currently working in the field for a while, but I will make time to further adjust the solver configuration. I will post an update in the forum once I make progress.
Thanks again for your assistance!
Best regards,
Dear Aagaard，
I’ve returned from fieldwork and completed the coupled modeling some time ago, but I’ve been held up in responding because one issue has been bothering me.
First, regarding the error in the program: as I posted in my earlier message with screenshots, the issue seems to stem from the fieldsplit configuration. I was able to run the simulation successfully with the following minimal settings, but I’m not sure whether using the default PETSc solver and preconditioner affects performance. When I slightly refined the mesh of my 3D model, the computation stopped due to failure to converge within 10,000 iterations.
[pylithapp.problem.petsc_defaults]
# Turn off default solver settings
solver = False
[pylithapp.petsc]
# These two settings reorder the degrees of freedom so that we can use the variable point block Jacobi preconditioner
dm_reorder_section = True
dm_reorder_section_type = cohesive
What’s been puzzling me more, though, is the behavior of the pressure field after fault rupture. As shown in the attached figure, negative pressure values appear near the fault tips. I initially assumed that, like mass, pressure should be greater than or equal to zero, and so I thought there might be something wrong with the parameters passed to the pressure field when handling the fault.
However, I’ve started to revise my understanding. Since the fault rupture itself is not a physical source of pressure, the pressure field should remain a conservative field. This would mean that the pressure values do not represent absolute values but rather changes or perturbations, which would explain why negative values can occur.
Am I correct in this line of reasoning?
Best regards,
