Impossibly slow simulation using prescribed velocity. I'm I missing something?

[Francyrad]

Dear users
I wanted to replicate a paper where the authors added an astenosphere fixed velocity of 3 cm/s years. in this way, slabs would have followed this counterclockwise flow changing in angulation:

https://www.nature.com/articles/s41598-017-06551-y

So I installed the shared library “prescribed velocity” and I tried to add a counterclockwise flow of 15 cm/ years in the upper mantle:

subsection Prescribed velocities
  subsection Indicator function
    set Variable names = x,y,t
    set Function constants = r_inner=5870310, r_outer=6300000
    # The following line should be all on one line without the line continuation character
    set Function expression = \
    if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, 1, 0); \
    if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, 1, 0)
  end




# Function to define the velocity values within the annulus
# The velocity is counter-clockwise, hence Vx = -Vy and Vy = Vx
subsection Velocity function
  set Variable names = x,y,t
  set Function constants = r_inner=5870310, r_outer=6300000, velocity_mag=0.15
  set Function expression = if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, -velocity_mag*y, 0); if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, velocity_mag*x, 0)

end


end

In this way everything is inside the upper mantle will have this flow + eventual speeds from simulations and what is outside the upper mantle won’t have it.

The point is that it doesn’t work, and the speeds that I get in general are too big, both inside and outside my annulus:

Obviously this is not what I usually have:

What is not clear to me are the order of magnitude of the speeds. Are they in m/yr? I’m I doing something wrong in my function? Any speed that i insert seems to be too big

Let me know and thank you for the help!

[bobmyhill]

The point is that it doesn’t work

I disagree. It looks like ASPECT is behaving correctly. You’ve just not quite understood the function you’re using. Plug an appropriate x and y value (in meters) into the function:
set Function expression = if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, -velocity_mag*y, 0); if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, velocity_mag*x, 0)

You should obtain the same value that paraview is showing.

Note that “velocity_mag” in the function has units of 1/s. You will need to scale your intended velocity (m/s) by the radius (m) at which you want that velocity.

Bob

[Francyrad]

Thank you for the suggestion! I solved in this way:

subsection Prescribed velocities
  subsection Indicator function
    set Variable names = x,y,t
    set Function constants = r_inner=3481000, r_outer=6300000 #originale 5870310
    # The following line should be all on one line without the line continuation character
    set Function expression = \
    if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, 1, 0); 0
  end


# Function to define the velocity values within the annulus
# The velocity is counter-clockwise, hence Vx = -Vy and Vy = Vx
subsection Velocity function
  set Variable names = x,y,t
  set Function constants = r_inner=5870310, r_outer=6300000, velocity_mag=8.1e-16 #to calculate it, for example, 0.15 m/yr, do (0.15/(3.154e7 * 5870310))
  set Function expression = if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, -velocity_mag*y, 0); if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer, velocity_mag*x, 0)

end

If i want my speed at 0.15 m/yr at 5870 km from the center of my planet, did I calculated it correctly?: (0.15/(3.154e7 * 5870310)) Spees seems reasonable

I only need to hope that my sim won’t crash

[Francyrad]

Anyway there is still something wrong. The red speed should be ok, but as you can see is like the whole mantle of my planet is rotating.

I can give 2 interpretations:

1. My upper mantle is dragging the whole planet, I wouldn’t be surprised
2. My function is still wrong

I ask because my sim can’t converge with a resolution of 7 and 8

[bobmyhill]

Dear @Francyrad,

Your calculation is now correct, and the output looks reasonable. The layer velocity is a bit slow - I suspect this is the solver struggling to find a solution.

Your new observations are not surprising:

• Your upper mantle is indeed dragging its surroundings (that’s what finite viscosity does…), and vice-versa.
• I suspect the solver is failing to converge because the stable solution (if there were no density anomalies and a free slip lower boundary condition, that would be pure rotation of the whole planet) is a long way from the initial guess.

The first issue (unintended-but-physically-necessary consequences of imposed conditions) is normal for simulations involving “Hand of God” velocity manipulations - especially one as extreme as differential shell rotations. That is very much a scientific setup problem, rather than an ASPECT problem.

The second issue (solvers failing) is also not unusual for “Hand of God” velocity manipulations. You are imposing constraints on the system which lead to extreme velocity gradients and stresses, and under such (unnatural) constraints I think any solvers are going to struggle.

I don’t want to impose my scientific views on your research, but I’ll leave you with a question: Can a homogeneous mantle wind exist without the “Hand of God”?

[Francyrad]

Dear @bobmyhill
I remain skeptical about the accuracy of the function and suspect there might be an oversight in my approach.

In utilizing the “Indicator Function,” my intention is to specify the region within the mantle where I aim to apply my flow model. Rather than implementing it throughout the entire mantle, I should focusing solely on the upper mantle. Application across the entire mantle would override the “normal speed” determined by my temperature gradients.

By restricting the range to just the upper mantle, only the velocities in that specific area are overwritten and modified according to the flow parameters I have set in velocity function

This approach results in a more constrained flow and maintains the speed of my original flow model. However, I am perplexed by the unexpectedly high velocities within my mantle model, as shown in the attached screenshot.

even if I still don’t understand why the speeds of my mantle are pretty high. I’m also sure that there is the flow in the mapper mantle because i checked x and y magnitudes.

A critical issue arises with the external library: it seems to completely override and eliminate all existing speeds. Despite my adjustments, the simulation remains excessively slow. Implementing a stronger flow results in excessively high speeds within the mantle, a phenomenon that is amplified by a factor of 10,000 or similar, regardless of whether I increase or decrease the velocity. Does another library or another method exist to impose a permanent flow in ASPECT?

Furthermore, I am unsure how to accurately calculate these velocities in cm/year and how to normalize them relative to the radius where I intend to apply my flow. Could a mathematical example be provided for this process?

My expectation was that the library would simply add a velocity or flow, upon which speeds could then vary, similar to adding a perturbation function to an adiabatic temperature profile. I also anticipated that slabs would be dragged in a counterclockwise direction, and that rising plumes in the upper mantle would behave as typically observed in simulations, but in a more permanent manner, as if under a constant net force.

Regarding the question, I just wanted to replicate what the paper did: impose a flow that drags slabs and plumes.
Does the flow exist without the Hand of God? No, because ASPECT doesn’t simulate earth’s rotation, which should produce this flow because. According to the paper, the astenosphere is less viscous than the litosphere and then it should rotating faster producing this flow as a consequence.

[bobmyhill]

In utilizing the “Indicator Function,” my intention is to specify the region within the mantle where I aim to apply my flow model.

By restricting the range to just the upper mantle, only the velocities in that specific area are overwritten and modified according to the flow parameters I have set in velocity function

Yes, that appears to be what is happening.

A critical issue arises with the external library: it seems to completely override and eliminate all existing speeds.

I don’t understand what you mean by “existing speeds”. Velocities in ASPECT are determined by solving the governing equations subject to the imposed boundary, initial, and prescribed conditions. If you change any of these (in your case, by adding prescribed conditions in a part of the domain), velocities throughout the domain can change.

I am unsure how to accurately calculate these velocities in cm/year and how to normalize them relative to the radius where I intend to apply my flow. Could a mathematical example be provided for this process?

As far as I can tell, you did it correctly above.

My expectation was that the library would simply add a velocity or flow, upon which speeds could then vary

The plugin imposes a fixed velocity in a region of the domain. You could impose a velocity that is equal to a vector from another simulation plus another vector (you would need to read the resulting vector in from an ascii file). That would still be a fixed velocity though, and I don’t really see the point.

I also anticipated that slabs would be dragged in a counterclockwise direction

Subduction cannot happen through a layer where you are imposing purely horizontal flow.

I think you are getting confused by this line in the paper you are trying to recreate:

“an imposed throughgoing asthenospheric mantle flow of −3 to +3 cm yr ^{−1} are implemented in the models”

In their simulations, they have fixed inflow and outflow velocities on the sides of the domain. They do not impose a velocity throughout the domain.

In your model, you could do a similar thing by constraining the flow through a narrow circular sector of the shell. You should probably define the flow throughout the sector (from the top to the bottom of the domain), rather than just in a narrow depth domain. You can of course allow the imposed velocity to change direction and speed with radius.

Finally, the origin and magnitude of net westward drift is controversial. I’d recommend you check out the other papers cited in the paragraph that starts

The “westward” drift has been so far considered only from the perspective of an average motion of the lithosphere relative to the mantle…

[Francyrad]

Thank you for clarifying, that makes more sense:

So i should put a free slip on top and bottom (I’m in a stagnant lid model), but a flow of 3cm/yr in left and right boundary (only in the upper mantle in my case), using the “bounday velocity model” subsection without using the external plugin, right?

I tried now to do that, but i cannot use left:function and right: function because i’m in a spherical shell at 360 degrees, Do you know that term should i use? Because i didn’t found them in the User Manual and i can use them only in a shell which has an angulation < of 360°

If i cannot use any terms for my side boundaries, must I use ascii?

Thank you for the preciouses suggestions and let me know

Francesco

[bobmyhill]

In your model, you could do a similar thing by constraining the flow through a narrow circular sector of the shell.

Best wishes,
Bob

[Francyrad]

Thank you for the suggestion. I tried to do that:

subsection Boundary velocity model
  set Prescribed velocity boundary indicators = top:function
  set Tangential velocity boundary indicators = bottom

subsection Function
  set Variable names = x,y
  set Function constants = r_inner=5870310, r_outer=6300000, velocity_mag=8e-16 
  set Function expression = if(sqrt(x*x + y*y) > r_inner && sqrt(x*x + y*y) <= r_outer, -velocity_mag*y, 0); if(sqrt(x*x + y*y) > r_inner && sqrt(x*x + y*y) <= r_outer, velocity_mag*x, 0)

But honestly I don’t know if i did it correctly, this is the result:

This is for a planet in an adiabatic profile. If I use an ASCII function to impose my temperatures, this won’t even happen (why?). It won’t happen also if I decide to harmonic average my material…
Did I do correctly? I have some issues. For example, it is not important what velocity magnitude i use. It can be 1e-16 or 1e16, I will get only these velocity. Also, the velocity boundary is only along a line, with very weird circumferences which are upward and downward. The flow is present only on that line and not in every part of my upper mantle.

I got better results in a box, but here i can choose “left and right” boundaries, and i cannot choose periodic boundaries:

subsection Boundary velocity model
  set Tangential velocity boundary indicators = bottom, top
  set Prescribed velocity boundary indicators = left: function, right: function
  
  
  subsection Function
    set Variable names      = x,y
    set Function constants  = cm=0.01, year=1
    set Function expression =  if (y>2.2e6 && y<2.8e6, 3*cm/year, 0) ; 0
  end
end

I’m I doing something wrong? Thank you for help and suggestions
Francesco

PS: I’ll try to do that in python and check if the ASCII method will work

[bobmyhill]

My suggestion was to constrain the velocities in a narrow circular sector: Circular sector - Wikipedia.
If that is what you want to do, you would set the Indicator Function to be 1 inside the sector and 0 outside.

The Boundary velocity model cannot be used to prescribe velocities inside a domain, but you could run a model with a chunk that was almost closed (359 degrees) and apply the Boundary velocity model to the left and right sides of that model. That would be an alternative to the “circular sector” approach, and might run faster (IIRC the Prescribed Velocities function doesn’t yet work well with MPI).

The Velocity function can have any value where velocities are not prescribed (you shouldn’t need the if/else blocks in that function as you have done so far).

The artefacts in your box simulation (last image) might go away if you had a smoother function - at the moment there are velocity gradient singularities at y = 2.2e6 and y= 2.8e6.

[Francyrad]

Hi @bobmyhill , i solved using this function:

set Additional shared libraries = ./libprescribed_velocity.so

## Turn prescribed velocities on
set Prescribe internal velocities = true


subsection Prescribed velocities
  subsection Indicator function
    set Variable names = x,y,t
    set Function constants = r_inner=3481000, r_outer=6300000, pi=3.141592653589793, tol=1e-5
    # This function expression identifies the edges of a 90 degree sector, both top and bottom with a tolerance tol
    # and apply the velocity only along those edges.
    set Function expression = \
    if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer && \
       (((atan2(y,x) >= -tol && atan2(y,x) <= tol) || (atan2(y,x) >= pi/2 - tol && atan2(y,x) <= pi/2 + tol)) || \
        ((atan2(y,x) >= -pi-tol && atan2(y,x) <= -pi+tol) || (atan2(y,x) >= -pi/2 - tol && atan2(y,x) <= -pi/2 + tol))), 1, 0); 0
  end

  subsection Velocity function
    set Variable names = x,y,t
    set Function constants = r_inner=5870310, r_outer=6300000, velocity_mag=1e-9, pi=3.141592653589793, tol=1e-5
    # The velocity function applies the prescribed velocity along the edges of the sector identified by the indicator function, within a tolerance tol.
    # It will apply a counter-clockwise velocity along both the top and bottom edges of the sector.
    set Function expression = \
    if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer && \
       (((atan2(y,x) >= -tol && atan2(y,x) <= tol) || (atan2(y,x) >= pi/2 - tol && atan2(y,x) <= pi/2 + tol)) || \
        ((atan2(y,x) >= -pi-tol && atan2(y,x) <= -pi+tol) || (atan2(y,x) >= -pi/2 - tol && atan2(y,x) <= -pi/2 + tol))), -velocity_mag*y/sqrt(x*x + y*y), 0); \
    if(sqrt(x*x + y*y) >= r_inner && sqrt(x*x + y*y) <= r_outer && \
       (((atan2(y,x) >= -tol && atan2(y,x) <= tol) || (atan2(y,x) >= pi/2 - tol && atan2(y,x) <= pi/2 + tol)) || \
        ((atan2(y,x) >= -pi-tol && atan2(y,x) <= -pi+tol) || (atan2(y,x) >= -pi/2 - tol && atan2(y,x) <= -pi/2 + tol))), velocity_mag*x/sqrt(x*x + y*y), 0)
  end
end

It didn’t have the effects that i wanted, but the idea it was.

So thank you for the help and suggestions!