Numerical breakdown regarding pressure BCs

Hi all,

I have a 3D model Cartesian box for which I want to prescribe traction boundary conditions of the kind of pressure on two of the side walls (left and right) and also at the bottom.

I tried using ‘initial lithostatic pressure’, but it didn’t work, it complained about my initial topography (which is a loaded data). We could discuss that eventually.
After getting tired of not understanding, and also thinking that in my model I prefer to know exactly what pressure profile is being forced on the walls, plus, my model is laterally heterogeneous, thus getting a good pressure profile (rho.g.depth) is not so direct; ..after that, I decided to set up for traction boundaries as functions which I can control. So, I am setting a unique expression that specifies traction vector functions over three different walls.  

I did:

subsection Boundary traction model

set Prescribed traction boundary indicators = bottom:function , left:function , right:function

subsection Function
set Coordinate system = cartesian
set Variable names = x,y,z
set Function constants = Pbott=8.2e9, dens=3400, g=10, H=2.52e5

 set Function expression = if( x<=0, dens*g*(H-z), if(x>=4e5,-dens*g*(H-z),0) ) ; 0 ; if(z<=0,Pbott,0)

set Function expression = if( x<=0, densg(H-z), if(x<4e5,0,-densg(H-z)) ) ; 0 ; if(z<=0,Pbott,0) ###### I also tried this one.


I also tried with that second expression (# set Function expression), which is slightly different, perhaps more correct…

For both of those trial functions:
That vector function with if statements will prescribe signed Px(z) on the lateral walls and Pbott at the bottom.
(I am planning to set more complex expressions for Px(z), regarding the different densities of layers, so this was just a starting point).

For both of them the compiler parses ! Aspects runs !

But, … The iterative Stokes solver did not converge :frowning:

Among other messages it reads: AztecOO::Iterate error code -2: numerical breakdown

The density I used is the maximum of the system, and H is slightly (2 km) larger than the model box height (cause after learning that trying the original value didn’t work, I decided to add the topography height to it); …So there will be an over pressure on the side walls, a pressure profile that is slightly larger than the theoretical lithostatic pressure…

I know that imbalanced forces on the walls may lead to …global acceleration, and a likely breakdown in the Solver. Here I have exactly opposite function forces on the lateral walls, but they are not exactly equal to the lithostatic pressure function, so I would expect some 3D global deformation adjustment in early times…I also tried increasing the bottom pressure…to no avail.

What do you think I should do?

You can find my prm file, the initial topography file, and the logo.txt log file with the run error message, here:



Hi Felipe,

As you stated, the fact that the model did not converge on the very first Stokes solve indicates something is clearly off.

I would not be surprised at all if this is related to having traction boundary conditions for multiple boundaries. In particular, as you noted the fact that you have the left and right boundaries that are not the ‘initial lithostatic pressure’ lends itself to all sorts of issues if inflow/outflow is not balanced correctly.

What do you think I should do?

Simplify your model design and systematically start adding in complexity. This could be done by starting with all free-slip boundaries and then one by one adding in the traction boundaries. Add in one change at a time from a setup that works and make sure you the model results for each change make sense.