# Problem with boundary conditions on 3D continental extension model

Hi all,

I’m trying to create a 3D continental extension model that undergoes oblique extension (similar to the 2D continental extension cookbook model), however, I’m having problems with my boundary conditions.

Currently I have:
• Top: free surface
• Left and right: function
set Variable names = x, y, z
set Function expression = if (x<150e3, -0.7071067812cm/year, 0.7071067812cm/year); if
(x<150e3, -0.7071067812cm/year, 0.7071067812cm/year); 0
• Front, back and bottom: initial lithostatic pressure

However, using these boundary conditions I get these ‘triangles’ of material at the front and back of the model and after around 12 Myr of extension I start getting really small timesteps (~700 years). I’m a little unsure as to why this is happening and I feel like I’m perhaps misunderstanding how the initial lithostatic pressure boundary works here, like why is so much material being upwelled from the front of the model compared to the middle of the base of the model? When setting the front and back boundaries of the model to tangential and pulling the box orthogonally the model works as I’d expect it to, so this leads me to believe it’s a problem with the front and back boundaries here, although again I could be wrong on this.

Any help or advice is much appreciated ,
Luke

Density after 8.5 Myr. The middle of my model breaks up how I’d expect it to, but at the front and back this isn’t the case.

X velocity after 8.5 Myr. Again, the centre of my model looks ok but at the front and back it’s a little strange.

Hi Luke,

Welcome and thanks for posting to the forum!

When setting the front and back boundaries of the model to tangential and pulling the box orthogonally the model works as I’d expect it to, so this leads me to believe it’s a problem with the front and back boundaries here, although again I could be wrong on this.

To follow on this finding, does the model also work as you would expect with tangential front and back boundaries, a lithostatic pressure condition at the base, and oblique kinematic boundary conditions on the left and right walls?

I have run a range of oblique 3D extension models, but all with velocity boundary conditions on the left, right, and bottom walls, and either periodic or velocity BCs on the front/back walls.

In most models I normally have a lower resolution and/or strong lithosphere “buffer zone” between the active rift area and the front/back walls.

I would not be surprised at all if the lithostatic pressure BC on the front/back walls is indeed producing issues giving rise to unrealistically high velocities (or some other unrealistic behavior)

Have you tried a very thin 3D model with lithostatic pressure BC on the front/back walls to see if it gives equivalent results to an equivalent model with front/back tangential BCs? That would be my first check for debugging, along rerunning your posted 3D model with kinematic or periodic BCs on the front/back wall.

As the lithostatic pressure boundaries can have inflow/outflow, I wonder if an issue is arising with assigning composition values in areas of inflow? What composition BCs are you assigning on those and other boundaries?

Depending on how much time you have for the project, it may be worth considering a switch to kinematic (or periodic) front/back walls. I would not be surprised if persistent issues arise when 4/6 boundaries are defined as either lithostatic pressure or free surface boundaries.

I hope this helps and lets keep the discussion going.

Cheers,
John

1 Like

Hi John,

Thanks a lot for the quick response! It’s really good to see the forums are so active.

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:

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. After not being able to get this to work I then tried using the initial lithostatic pressure condition on the front and back to create some type of ‘open’ boundary.

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. However, if I run a similar model with orthogonal extension but tangential boundaries on the front and back it works fine.

Currently I use:
subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end
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.

To try and sort out this problem I’ve also tried:
• Removing the free surface but similar results still occur.

• Placing the seed 100 km from the front and back boundaries of the model, which did help in making it take longer for the timesteps to get really small but didn’t fully solve the problem.

• To try and reduce the anomalously high velocities from the z and y axis along the front and back boundaries using both function and initial lithostatic pressure conditions:
subsection Boundary velocity model
set Prescribed velocity boundary indicators = left:function, right:function, front zy:function, back zy:function
set Function expression = if (x<150e3, -0.7071067812cm/year, 0.7071067812cm/year); if (x<=0, -0.7071067812cm/year, \ if (x>=300.e3, 0.7071067812cm/year, 0)); 0

subsection Boundary traction model
set Prescribed traction boundary indicators = bottom:initial lithostatic pressure, front:initial lithostatic pressure, back:initial lithostatic pressure

However, whilst this model did work when extending orthogonally it couldn’t get past timestep 1 when pulled obliquely.

I was considering either starting to build my model from scratch but this time giving my base a set velocity so material is constantly being upwelled from the base and then if my front and back boundaries continue to behave like this I can snip them off when presenting them. Alternatively, I was also thinking (like you mentioned) if placing some stronger material along the front and back boundaries could help, although I’m a little unsure exactly how to implement this.

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?

Big thanks again,
Luke

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:

1. boundary pressure profile is deviating significantly from the current pressure profile somewhere along the boundary after crustal breakup and/or

2. 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.

Cheers,
John

1 Like

Hi John,

Thanks again for the response! I’m glad I didn’t scare you away with my last reply, it was like a two page word document or something crazy that I responded with.

Anyway, I’ve spent the last couple of days playing with the periodic boundary condition and it’s worked great on my small test models as long as there was a gap between the front and back boundaries and the seed (image below). Thanks for suggesting this, I feel a little silly previously trying to use pressure boundaries.

I still have quite a good amount of time to work on my project so I believe it’s worth considering both options particularly if using these function boundaries is more stable. Currently, I’m struggling a little bit to imagine how I would implement function boundaries on the front and back wall, could you perhaps show me an example of how this has been implemented before?

I also noticed that you mentioned using the GMG Stokes solver a couple of times, so I did a bit of reading on it and had a bit of a play with it in some basic 2D models continental extension models. Is using the GMG stokes solver one of the fastest and/or most accurate way of modelling a 3D continental extension model?

Wow thanks a lot for the offer! I’ll attach the .prm file I used with periodic boundary condition down below and if you have any tips or tricks for improving my model’s performance or accuracy let me know.

Thanks again you’ve been a big help,

Luke
frontback_test_03.prm (12.2 KB)

Hi Luke,

Thanks again for the response! I’m glad I didn’t scare you away with my last reply, it was like a two page word document or something crazy that I responded with.

More information is always good, you certainly didn’t scare me away. Thanks for the helpful response to all the points!

Great to hear the periodic boundary conditions are making the simulations more stable and you now have at least one viable options to move forward with.

Rather than update your PRM file directly, I’ve attached an updated version of the continental extension cookbook that uses both the GMG solver and active particles with particle interpolation methods. You already use a significant number of initial global refinements (5), which will help the solvers quite a bit.

My suggestion would be to first test the solver settings I use (GMG, iterated defection correction Stokes), and then add in the particles to see how it affects the run times and accuracy (less diffusion compared to fields). In some 3D extension models we have that found the GMG solver decreases runtimes by ~ a factor of 2 (2x faster) for models where the Stokes solver dominates the total wall clock time. It should also decrease total memory usage. However, the performance changes is somewhat on a case by case basis.

The GMG solver does not currently work with periodic boundary conditions, so you will need to do the testing with either free-slip front/back walls or prescribed velocities on the front/back walls.

For the latter method (prescribed velocities on front/back), here is a PRM snippet from a 700x700x175 km model with prescribed velocities on all walls (constant `theta` determines the obliquity).

``````subsection Boundary velocity model
set Tangential velocity boundary indicators =
set Prescribed velocity boundary indicators = left xy: function, right xy:function, front xy:function, back xy:function, bottom z: function
subsection Function
set Variable names      = x,y,z
set Function constants  = v=0.004, d=175.e3, w=700.e3, theta = 0., xa = 350.e3
set Function expression = if( x<w/2., -v*cos(theta*pi/180)*(1 - x/xa), v*cos(theta*pi/180)*(x/xa  - 1)); \
if( x<w/2., v*sin(theta*pi/180)*(1 - x/xa), -v*sin(theta*pi/180)*(x/xa  - 1)); \
v*cos(theta*pi/180)*2*d/w
end
end
``````

Let us know how the solver and particle testing goes, and likewise if you have questions about any of the above!

Cheers,
John

1 Like

Hi again John,

Bit of a late response since I’ve had some time off and have given myself a good amount of time to test and learn how everything you’ve suggested works.

I first started testing 2D continental extension models and implemented your solver settings (GMG, iterated defection correction Stokes) and then tried adding particles. The GMG solver you’ve suggested using has worked great in both my 2D and 3D models, making models run a lot quicker and letting me really up the resolution in the higher resolution models. The particles have in some cases also made things slightly faster but normally it takes roughly the same length of time to run models with and without particles, nevertheless, they help accuracy so have generally been a good addition to my current models.

Next, I started playing with the function boundary conditions you sent me and have got them working. Currently I have a periodic front and back boundary model as well as a seemingly more stable function boundary condition model which uses your solver settings and particles. Here I’ve decided that for my models I’m going to be moving forward with the function boundary conditions.

I only have one concern with my models going forward. I’m implementing an offset seed so in order to try and mitigate edge effects from the front and back boundary I’m only going to be visualizing a small area, furthermore, to better runtime I’ve also reduced the resolution above and below the regions with a seed (image below). I figured I’d ask you if this sounds like a sensible solution? Or do you have a better idea on how to reduce the edge effects? I’m just a little concerned with the velocities along the front and back since they’re near zero around the centre which doesn’t quite link up with the velocities in the middle of my model.

Big thanks again! I really couldn’t have got all this working without your guidance.

Luke

Hi Luke,

This is great news, thank you for doing those tests and reporting back!

The GMG solver you’ve suggested using has worked great in both my 2D and 3D models, making models run a lot quicker and letting me really up the resolution in the higher resolution models.

If you have time, could you check and let us know how much the speedup was when switching from AMG to GMG for both the 2D and 3D models? This will help us document guidelines for performance improvements users can expect. Also, in those tests was AMG-GMG the only thing that changed, or did the changes also include things like the switch to defect Correction Stokes nonlinear solver, different material averaging, etc?

Boundary conditions -

I figured I’d ask you if this sounds like a sensible solution?

Yes, definitely! This is the approach I also take in similar studies.

For example, in this recent paper we refined the grid twice toward the model center in order to avoid edge effects.

In ongoing work that looks at rifting through breakup with oblique deformation we also use this approach, but the manner in which we refine changed a bit.

This all to say that the manner in which you refine the mesh to maximize resolution in regions of interest and dampen boundary effects may change quite a bit between studies, but I’ve generally found is a useful approach.

Thank you again for performing the solver tests!

Cheers,
John