I’ve been playing around with the continental extension model in the manual and was interested if there’s a way of applying a force to the left and right boundaries of the model instead of a set velocity. Maybe I’m completely off the plot here but from looking at the manual it looks like in the Boundary traction model subsection you can specify a function to calculate the traction force for each vector component (x,y,z) for which you wish to have a force in that direction. However, I’m a little unsure how to implement this or if it’ll work like I presume it will.
If anybody has done anything like this before or even knows if you can implement this any help is much appreciated.
Yep, you are on the track. The Boundary traction model can be used to prescribe tractions on boundaries, rather than velocities.
In general, models with both a free surface and boundary tractions on one or more side walls tend to be less stable than equivalent models with kinematic boundaries, but with enough testing various configurations you probably can arrive with a stable suite of simulations.
As lithostatic pressure will dominant for a crustal- or lithospheric-scale model, what you likely need to prescribe on the traction boundaries is the lithostatic pressure plus a constant that defines either inward or outward flow (simple 2D example).
Previously, I believe this combination could only have been done via a carefully designed user-specified function, but now it should be possible to combine multiple tractions models with the following recently merged PR (yay hackathons ) - Problem with update strain rheology - #13 by jbnaliboff
As a starting point, I suggest reading through that PR (particularly the new example test case) and:
Reading through relevant tests cases, including box_lithostatic_pressure_2D.prm. I would search for keywords like lithostatic_pressure and boundary_traction.
Go through the benchmark suite for nonlinear channel flow, which uses boundary tractions to drive deformation.
I hope this helps and it would be great to continue the discussion with any further questions about existing examples or new tests you conduct.
Thanks for the response and the links to all the useful information! It’s good to hear that the traction boundaries work similar to what I thought they would. I’ve also had a read through the links you suggested but there’s still a few things I’m a little unsure about.
Firstly, the closest example I found to what I want to achieve is the attached file and image below which I somehow found on Juliane Dannberg’s Github page. This example uses the box with lithosphere boundary indicators and applies a free surface on top, zero velocity boundary indicators on the bottom boundary, left lithosphere: function, right lithosphere:function, left:initial lithostatic pressure, right:initial lithostatic pressure. It also uses a function that can be summarised as: density * gravity * (height of model – y) + or – stress. I also noticed that the box_lithostatic_pressure_2D.prm has a similar function. My question here is basically why do these models prescribe traction functions in this way? Why not use a constant value across the boundaries? traction_plate.prm (7.4 KB)
I’ve also had a go at applying boundary tractions on a 2D continental extension model, but I’ve been getting some funny results. To try and get my models working I’ve been using zero velocity boundary indicators on the bottom, free surface on the top and x function tractions on the left and right. However, when applying a constant traction on the sides I end up getting really high velocities dipping downwards (image below). I’ve tried playing around with this a bit but haven’t had any luck and I’m not sure what the best thing to try going forward is. Do you have any advice or insight as to what’s going on here?
Haha I guess I asked at a good time then. I also saw you released ASPECT 2.5.0 recently (looks like a lot of really cool changes btw), but I’m still on the past version (2.4.0) and it’ll take me a few weeks at least before I’ll be able to get the newer version going. Would you recommend I keep on trying to get the traction boundaries working on my current version or wait until I can use 2.5.0?
This is a very delayed reply to your last post, but hopefully it may be still be useful for you (or others).
I also noticed that the box_lithostatic_pressure_2D.prm has a similar function. My question here is basically why do these models prescribe traction functions in this way? Why not use a constant value across the boundaries?
Consider a model with no internal variations in density and free slip boundaries on all sides. The Stokes solution to this setup should have effectively zero velocities everywhere and a laterally uniform lithostatic pressure profile.
If one were to replace the side boundaries with traction boundaries, maintaining this state of equilibrium would require having a depth-dependent traction value that matches the lithostatic pressure.
If the traction value at a given depth did not match the lithostatic pressure inside the model, then material would flow in or out as a result of the pressure gradient. As depth increases, this pressure gradient across the boundary would increase and potentially produces extremely large velocities.
Coming back to the models you mentioned in your last post, what they are trying to achieve is adding a tectonic stress (i.e., extension or compression) to a lithostatic pressure profile on the traction boundaries.
This should produce a pressure gradient across the boundary that is equal to the tectonic stress, rather than an additional difference in the lithostatic pressure gradient across the boundary.
Hopefully this makes sense!
Would you recommend I keep on trying to get the traction boundaries working on my current version or wait until I can use 2.5.0?
If you have not already done so, I would switch to using the main branch!
Thanks a lot for the response! You summarised everything I wanted to know really well and after doing some extra reading and playing with the function traction boundaries I’ve realised I really didn’t understand how they worked.
Coming back to my original post, what I’m essentially trying to do is write a function that applies the lithostatic pressure profile as well as an additional force on the traction boundaries. My initial thought on this was to look if there’s any way to apply the initial lithostatic pressure traction condition as well as a function traction condition on the same boundary, but I don’t believe it’s possible to apply multiple boundary conditions like this?
I’m also a little unsure here as to what you mean here by being able “to combine multiple tractions models with the following recently merged PR”. Really sorry if I’m sounding naïve here but I’m not quite sure how something like this would be implemented. I’ve gone through the linked forum post and PR, but I still aren’t entirely sure what you mean. Would you be able to explain how something like this would be implemented or give another example?
I’m a bit late to the party, but I recently got access to the newest version . I’m currently reading through some the new changes that have been made so I’m really excited to play around with it.
I’m also a little unsure here as to what you mean here by being able “to combine multiple tractions models with the following recently merged PR”. Really sorry if I’m sounding naïve here but I’m not quite sure how something like this would be implemented.
Apologies, I linked the wrong PR in the previous post, which should have linked to the following PR:
) - 
. I’m currently reading through some the new changes that have been made so I’m really excited to play around with it.