Boundary velocity (flow in at bottom) leads to a failure of convergence

ASPECT
1

Hi all,

Hope you’re doing well. I’m pretty new to ASPECT tbh and I’m trying to reproduce plume cases in an early paper of Dr. Dannberg (Compressible magma/mantle dynamics: 3-D, adaptive simulations in ASPECT, GJI 2016). The repo could be accessed at GJI2016 of Dr. Dannberg.

The problem is that for subcases plume_compressible and plume_hydrous_compressible, the following boundary condition would lead to a failure of convergence:

subsection Boundary velocity model
  set Tangential velocity boundary indicators = left, right, top
  set Prescribed velocity boundary indicators = bottom:function
  set Zero velocity boundary indicators       = 

  subsection Function
    set Function constants  = b=100000, c=20000
    set Variable names      = x, y
    set Function expression = 0.0; -0.024995 + 0.1 * exp(-((x-b)*(x-b)+y*y)/(2*c*c))
  end
end

And I believe it is the vertical velocity at the bottom (the math expression) that caused the problem. A full prm file is attached
plume_compressible.prm (6.6 KB)
. This prm is modified accordingly since I’m using 2.6. You can see the log
log.txt (1.4 KB)
and solver history
solver_history.txt (15.6 KB).

What I’m thinking is that this y component velocity is a source of mass loss or gain of the system, since there is no other boundary condition for mass balance that I can think of for this case. Surprisingly, in Table 3 of Dr. Dannberg’s paper, u_s (solid velocity) is omitted, I interpreted it as (0, 0). But in section 4.6, a description of velocity boundary condition is added:
the velocity boundary conditions are free slip everywhere except for the bottom boundary layer, where the hydrostatic pressure is applied, but material is allowed to flow in and out.
So the last piece says y component of us is actually free. So the Table, the boundary condition in the context and boundary condition in prm are not consistent? Which is right?

When I set u_s (0, 0), i.e., in the prm set Function expression = 0.0; 0.0, it becomes better but again fails during the first step. The log file
log.txt (2.3 KB)
and solver history file
solver_history.txt (11.6 KB)
are attached.

Another interesting thing is that velocity field for step 0 is not as expected:

Screen Shot 2024-09-02 at 12.10.08
Screen Shot 2024-09-02 at 12.10.081920×1070 130 KB

u_s,x, u_s,y, u_s,z at bottom boundary are consistent with what is given in the prm. But the velocity over the domain should be (0, 0), isn’t it? Or is this the result of some initial guess/adjustment?

Any thoughts? Really appreciate your help!

Mingming

2

PS: smaller initial perturbation like set Amplitude = 50 plus u_s = (0, 0) would help.

3

@optimux - Welcome to the community and thank you for posting your question to the forum! Second, apologies it took awhile for a reply.

The issue is that the Boundary velocity model you prescribed in the updated PRM is not actually what was used in the PRM files from GJI2016 of Dr. Dannberg.

In the “plume” series of PRM files, the prescribed velocities constraints are:

  set Prescribed velocity boundary indicators =

  set Tangential velocity boundary indicators = 0,1,3
  set Zero velocity boundary indicators       =

There is a prescribed velocity section elsewhere in the PRM file, but the function from that section is not actually applied to the bottom boundary.

Per the code snippet above and what you noted above from the paper, the velocity boundary conditions are free-slip on the left (0), right (1), and top (3) boundaries. A condition is not applied for the bottom (2) boundary, which means the default boundary setting (pressure) is applied.

Can you try re-running the model after adding the following section to produce a lithostatic pressure condition on the bottom boundary?:

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

  subsection Initial lithostatic pressure
    set Number of integration points = 1000
    set Representative point         = 2000000,0
  end
end

Note - the above snippet was taken and modified from this cookbook.

Cheers,
John

4

Hi John,

Thank you so much for the reply, I was so wrong! The syntax is quite misleading.

The original PRM from Juliane2016 has the following:

subsection Boundary velocity model
  subsection Function
    set Function constants  = b=100000, c=20000
    set Variable names      = x,y
    set Function expression = 0.0; -0.024995 + 0.1 * exp(-((x-b)*(x-b)+y*y)/(2*c*c))
  end
end

and

subsection Model settings
  ...
  set Prescribed velocity boundary indicators =

  set Tangential velocity boundary indicators = 0,1,3
  set Zero velocity boundary indicators       = 
  ...
end

So the first question is the long math expression, what is it describing? y component of solid velocity?

The second question is the role that subsection Boundary velocity model plays. To which boundary shall this Boundary velocity model be applied? It’s like a body without a head (not a good metaphor though). I was thinking this Boundary velocity model would be applied to the bottom boundary because the Prescribed velocity boundary indicators is left blank/default, i.e., no 2 up there.

The third question is how do you know Prescribed velocity boundary indicators is DEFAULT to pressure. Is it in the source code or cookbook? In manual 2.1.0, I don’t see any clues.

The snippet you provided is helpful, it’s working. Thank you !

Cheers,
Mingming

5

Hi Mingming,

So the first question is the long math expression, what is it describing? y component of solid velocity?

Yes, it was describing the y-component of velocity. However, in this case that function was not actually used (i.e., it was just extra code from previous testing).

The second question is the role that subsection Boundary velocity model plays. …

I don’t recall the exact syntax for this earlier version of ASPECT, but if one changed the model setting to have set Prescribed velocity boundary indicators = 2 and indicated that the velocity would be prescribed through a function, then the aforementioned function would have utilized.

Quite some time ago many of the parameters in subsection Model settings were moved to more relevant subsection to make the PRM files easier to ready. For example, all of the parameters related to velocity boundary indicators under subsection Model settings are now specified under subsection Boundary velocity model.

The third question is how do you know Prescribed velocity boundary indicators is DEFAULT to pressure. Is it in the source code or cookbook? In manual 2.1.0, I don’t see any clues.

That is a great question. I know this anecdotally, but a quick search in the manual did not reveal anything. This could certainly be sorted out by looking at the source code, but it should be in the manual. I made an issue about this on the github page.

The snippet you provided is helpful, it’s working. Thank you !

Great, glad that it is working!

Cheers,
John

6

Hi John,

Thank you so much! BTW, in Section 4.6 of Juliane2016 paper, she explained: the velocity boundary conditions are free slip everywhere except for the bottom boundary layer, where the hydrostatic pressure is applied, but material is allowed to flow in and out. But unfortunately in her original repo Juliane2016_repo, there is no subsection of Initial lithostatic pressure. May I ask her to update the repo? @jdannberg

Cheers,
Mingming