Dear Baagaard,
I hope this message finds you well. I encountered an issue when using time_history. The error message I received is as follows:
mpinemesis: …/…/…/pylith-4.1.3/libsrc/pylith/bc/TimeDependentAuxiliaryFactory.cc:320: static void pylith::bc::TimeDependentAuxiliaryFactory::updateAuxiliaryField(pylith::topology::Field*, PylithReal, PylithReal, spatialdata::spatialdb::TimeHistory*): Assertion `auxiliaryFieldArray’ failed.
file.zip (1.6 MB)
I have tried various approaches, but I’m still uncertain about the root cause of the problem. Could you kindly take a look and assist me?
Thank you very much!
I think the problem is that this assert() is incorrect when a process does not have any points for the boundary condition. I will fix the code and get this into v4.2.0 which I hope to get released this week.
I noticed you are running on 50 processes with a mesh that is not very big. You can avoid this bug if you run on a single process.
Also, PyLith is setup to run in parallel using ``-nodes=NPROCSwhereNPROCSis the number of processes. Usingmpirun -np NPROCS` is not as portable.
1 Like
Thank you very much for your help. Everything is now running smoothly. I look forward to your new achievements being published. Thank you for your outstanding contributions. I wish you good health and all the best!
Sorry to bother you again. The results from this simulation show that the fault near the reservoir exhibits negative pore pressure, while regions farther away display higher pore pressure. This seems to be an unreasonable situation. Could it be that there is an issue with my simulation conditions? Could it be that the conditions applied via time_history cannot effectively simulate the diffusion of water?
Please explain exactly what the plot is showing. Also please provide a simple diagram showing the boundary value problem you are trying to solve. It could be you are misinterpreting the output, the simulation parameters to not match what you have in mind, or the results are more complex than your intuition.
The white trajectory line at the top of the figure represents the two sections of the reservoir area I applied, where I imposed pressure conditions equivalent to the water level on the surface. The total simulation time is over 1000 days. I simulated the change in pore pressure of surrounding faults under the continuous diffusion of water. However, I found that the pore pressure near the reservoir area decreased in the figure, which theoretically should not be the case. So, I am wondering if there is a problem with the pressure conditions I applied. I believe the pressure should be in the negative direction along the z-axis, so I set the pressure conditions in the impulse.timedb file to negative values. From my personal understanding, the pore pressure caused by water diffusion should be positive in any area. Or is there something wrong with my understanding? However, it shouldn’t be the case that the closer the area is to the reservoir, the more negative the pore pressure becomes, right?
al.zip (1.6 MB)
It looks like you are applying Neumann boundary conditions on the surface as a proxy for surface loading of a reservoir. However, you do not have any boundary conditions for the fluid pressure. I think you want at least an initial condition for the fluid pressure that is consistent with your initial value of the Neumann BC.
In your Neumann BC, you have an amplitude for the normal component of the time history of -6.78 MPa and negative values for the time history, so the resulting value will be positive, meaning an outward normal traction. For simplicity, just consider the normal component of the traction. The applied traction (outward normal is positive) is A(x) * f(t) where A(x) is the time_history_normal_amplitude given in a spatial database and f(t) are the values given in the time history file. We usually use f(t) to specify the “shape” of the time history, so it has a normalized amplitude, such as a peak value of 1.0. We use A(x) to specify the amplitude, so it encompasses the direction and amplitude. It looks to me like you need to normalize your time history values and make them positive, so that the applied normal traction is negative.
Thank you for your patient explanation; I now understand the issue. By the way, I’d like to ask why, in PyLith 4, reducing the time step sometimes leads to iteration failure. This seems different from PyLith 2, where a smaller time step typically ensured the calculation’s stability and accuracy.
At the same time, I would like to ask you how to set a fault as locked in Pylith4. I remember that in Pylith2, this was done through parameters like friction, but I did not see this setting in the manual for the new version
Effect of time step size on solve
Please provide detailed information about the boundary value problem (for example, the governing equation and bulk rheology) in which reducing the time step results in failure of the linear solve. For linear elasticity, the time step size should have little affect on the solve. There can be an effect if the time scale used for nondimensionalization becomes too large or too small. We have changed most of the preconditioners between PyLith v2 and PyLith v4; in most cases the ones used in PyLith v4 are better and should be more robust.
Locked faults
Fault friction is not implemented (yet) in PyLith v3 and later. For a locked fault, you can use prescribed slip with zero slip. This is equivalent to imposing a high coefficient of friction that prevents the fault from slipping.
Thank you very much! Actually, the convergence failure regarding the dt setting occurred during the fluid diffusion simulation I shared with you earlier. Is it because of an issue with my normalizer.relaxation_time setting that prevents dt from being set too small? Or is there another reason?
The mechanics of the solver for different governing equations are quite different, so it isn’t appropriate to compare the robustness of changing the time step in a quasistatic linear elastic simulation to quasistatic poroelasticity.
The time scale used for nondimensionalization (normalizer.relaxation_time) is important for poroelasticity. As mentioned in PETSC ERROR when use simulations in Pylith 4.1.3 - #6 by baagaard, for a well-conditioned system for poroelasticity, the ratio of the nondimensional permeability to viscosity should be about 1.
Another potential issue is that the default time stepping algorithm for quasistatic simulations is backward Euler. A higher order time stepping method may work better. This is something we have not yet played around with.
Examining the preconditioned residual and true residual in the linear solve (ksp_monitor_true_residual = true) is often helpful in understanding the convergence of the linear solve. If the true residual is decreasing much slower than the true residual, it means the preconditioner is not very effective or the scales for nondimensionalization need adjustment. This is mentioned in a recently added section of the PyLith manual on diagnosing solver issues (Troubleshooting — PyLith 4.2.0dev documentation).
I feel a bit confused, sir. You mentioned that the ratio of 1 refers to the dimensionless permeability and viscosity ratio. So, are these permeability and viscosity values the ones I specify in the material file of the model, or do I need to define them in the pylith.app file? Also, you mentioned that it involves the time scale, so I would like to ask what the calculation formula for this content is, or which section of the manual discusses it? Please forgive my ignorance, as I thought I understood it before, but now I am still a bit confused.
Given the poroelastic material properties that you specify in the spatial databases, you may need to adjust the time scale used for nondimensionalization (normalizer.relaxation_time) so that the ratio of the nondimensional permeability to the nondimensional viscosity should be about 1.
From pylith/libsrc/pylith/materials/AuxiliaryFactoryPoroelastic.cc at main · geodynamics/pylith · GitHub and pylith/libsrc/pylith/materials/AuxiliaryFactoryPoroelasticity.cc at main · geodynamics/pylith · GitHub
permeabilityScale = lengthScale*lengthScale
viscosityScale = pressureScale*timeScale
1 Like
Thank you so much for your continued help! I sincerely appreciate it! You’ve truly enlightened me, and I thank you once again!
