Hi Luke,
Thanks a lot for the quick response! It’s really good to see the forums are so active.
My pleasure! It’s really exciting from the developer side to see the code being used for so many interesting projects and questions useful to the broader community being asked on the forum.
Likewise, thank you for the detailed response and images. Responses below.
To better clarify than in my original post yes the problem is these anomalously high velocities occur at the front and back boundaries of my model and very little material is inputted through the base of my model:
It’s also worth noting that these high velocities only start to occur after my upper and lower crust have broken up and mantle material is being upwelled to the surface. Before breakup material is largely inputted through the central base of my model:
If I run a model with left and right orthogonal extension and then apply the initial lithostatic pressure condition to the front and back boundaries I get the same problem as to what I get when I run a model with oblique extension.
Based on the above information, I suspect the issues is that the:
-
boundary pressure profile is deviating significantly from the current pressure profile somewhere along the boundary after crustal breakup and/or
-
There is an issue related to the composition not being fixed in areas of inflow along the front/back walls.
I have also tried using “set Fixed composition boundary indicators = bottom, front, back” but this resulted in timesteps becoming really small earlier than my current model and still resulted in these high velocities.
I’m not sure why this is the case, but regardless I think the path forward to having stable simulations is not using pressure boundary conditions on the front/back walls.
I initially tried running a model with tangential boundaries on the front and back, initial lithostatic pressure on the base and oblique extension on the left and right, however, this results in the error message “The iterative Stokes solver in Simulator::solve_stokes did not converge” before timestep 1.
This makes sense, as you then have a situation where on the left/right walls there is a non-zero prescribed y velocity (for oblique extension), but at the front/back walls the y-velocity is fixed to be 0.
Thus, there is effectively a mismatch in prescribed velocities along the lines where the left/right and front/back walls intersect, giving rise to the Stokes solver errors.
You also mentioned using kinematic or periodic boundary conditions on my front and back boundaries. Sorry if I’m sounding naive but what exactly do you mean by this?
No problem at all and not a naive question, I should have said velocity instead of kinematic.
In short, you were close to one of the solutions I (and others) have used when modeling oblique extension/compression.
Option 1 is to prescribe a function for x-y-z velocities on the front/back walls. The function needs to prescribe a change (I usually do linear) in the x and y velocities as one moves from the left to right along the boundaries. This is to make sure the velocities at the front/back and left/right walls match where they intersect. Let me know if you would like an example of how to do this.
Option 2 is to use periodic boundary conditions, which roughly can be thought of balancing things in a manner such that the conditions at the beginning (ex: back boundary) are the same as at the end (ex: front boundary). This allows specifying any velocity (orthogonal, oblique) condition you want on the left/right boundaries, and the velocities/pressure on the front/back walls are then determined internally. The only disadvantage to this is you can’t use periodic BCs with the new GMG Stokes solver (more on this below).
Here is an example of how to obtain periodic BC on the front/back walls, which requires removing the null space in the y-direction. The prm snippet below is for a 500x500x100 km model geometry:
# Model geometry (500x500x100 km, 20 km spacing)
subsection Geometry model
set Model name = box
subsection Box
set Y periodic = true
set X repetitions = 5
set Y repetitions = 5
set Z repetitions = 1
set X extent = 500e3
set Y extent = 500e3
set Z extent = 100e3
end
end
# Remove nullspace in y direction
subsection Nullspace removal
set Remove nullspace = net y translation
end
In your setup (top: free surface, bottom: traction BC), in the Boundary velocity subsection you only need to provide constraints for the left and right walls.
In general, I have found Option 1 to be more stable across a range of complicated multiphase rifting scenarios and it is also compatible with the new GMG Stokes solver developed by Timo Heister and his group. However, Option 2 is easier to implement in the PRM file.
Please let me know if you have questions about any of the options above, and happy to help you implement them in your PRM file if you don’t mind sending it over.
Likewise, happy to share PRM configurations for new faster solver options and active particles (versus) fields that improve performance and accuracy.
Hope this was helpful!
Cheers,
John
,
