Rotational component in the 3D

[HanyMKhalil]

Dear Romain,

I am trying to implement varying velocity in a 3D layer to account for rotational extension, is it possible in UWGeodynamics or not yet implement so that I will need to write a function for it?following the paper Money et al., 2017?

Thanks a lot

[HanyMKhalil]

I am sorry the paper is Mondy et al., 2017

[LukeMondy]

The model in Mondy et al., uses velocity conditions that are perpendicular to wall boundary (as opposed to a 2D rotational one), and increases linearly from 0.5 cm/yr (close to the Euler pole) to 5 cm/yr (far from the pole). So the function should just be a linear one.

[HanyMKhalil]

Thanks Mondy so much So I just write a linear function to increase the velocity on oneside, for example left side in 3D, and increase the velocity along y axis, correct?

[LukeMondy]

We applied the velocity to both walls, but essentially yes.

For example,

import numpy

y_in_km = numpy.linspace(0, 1000)

velocity_along_y_axis = y/1000*4.5+0.5

print(velocity_along_y_axis)
array([0.5       , 0.59183673, 0.68367347, 0.7755102 , 0.86734694,
       0.95918367, 1.05102041, 1.14285714, 1.23469388, 1.32653061,
       1.41836735, 1.51020408, 1.60204082, 1.69387755, 1.78571429,
       1.87755102, 1.96938776, 2.06122449, 2.15306122, 2.24489796,
       2.33673469, 2.42857143, 2.52040816, 2.6122449 , 2.70408163,
       2.79591837, 2.8877551 , 2.97959184, 3.07142857, 3.16326531,
       3.25510204, 3.34693878, 3.43877551, 3.53061224, 3.62244898,
       3.71428571, 3.80612245, 3.89795918, 3.98979592, 4.08163265,
       4.17346939, 4.26530612, 4.35714286, 4.44897959, 4.54081633,
       4.63265306, 4.7244898 , 4.81632653, 4.90816327, 5.        ])

I'm not entirely sure how to translate it into UW2 code, but hopefully you get the idea.

[HanyMKhalil]

Dear Mondy,

Thank you so much for that really 💥💥

[HanyMKhalil]

Hello Romain,

I tried the linear function in this way (I calculated the slope and intercept)

conditions = [(True, Model.y GEO.nd(0.00125 u.centimeter / u.year) + GEO.nd(0.5 * u.centimeter / u.year))]

function = GEO.uw.function.branching.conditional(conditions)

Model.set_velocityBCs(left=[-1.0 * function, None , None], right=[function , None , None], front=[None, 0.0, None], back=[None, 0.0, None], bottom=GEO.LecodeIsostasy(reference_mat=Crust))

but in the output it always give me a constant velocity? although the velocity BC line excited well? what is the problem you think? why it couldn't change the velocity along the wall?

Thanks a lot

[lmoresi]

Hany, thanks for getting in touch - I will reply here rather than email as it is a better way to track the issues and other people can read the answers. I suggest that you send through a stripped-back notebook that Romain and others can quickly look at. The resolution can be very coarse as this is just to debug the syntax of the boundary condition parsing. The first impression I have is that this is already rather complicated for what appears to be a constant boundary value.

I am now trying to implement the rotational extension following your paper in Mondy et al, 2017, I contacted Mondy and Romain in underworld and they told me to make it as a linear equation however in the code it always gave the same velocity?  I calculated the slope and the intercept but I do not know what is the problem

conditions = [(True, Model.y * GEO.nd(0.00125 * u.centimeter / u.year) + GEO.nd(0.5 * u.centimeter / u.year))]

function = GEO.uw.function.branching.conditional(conditions)

Model.set_velocityBCs(left=[-1.0 * function, None , None],
                        right=[function , None , None],
                       front=[None, 0.0, None], back=[None, 0.0, None],
                        bottom=GEO.LecodeIsostasy(reference_mat=Crust))

I do not understand what is the problem? I am starting by a simple viscoplastic layer in the 3D I can send you my code if you have time?

[HanyMKhalil]

Dear Dr Moresi, Thanks a lot I attached the file, I downloaded it as notebook.ipynb but I added .txt to be able to upload it.

Rotation.ipynb.txt

[rbeucher]

Hi @HanyMKhalil

Thank you for posting your model.

I did run that set up a while ago. Here is what I did:

velFunc = -Model.y / GEO.nd(Model.maxCoord[1]) * GEO.nd(2.0 * u.centimeter / u.year) + GEO.nd(2.5 * u.centimetre / u.year)

Model.set_velocityBCs(left=[-1.0 * velFunc, None, None],
                      right=[velFunc, None, None],
                      back=[None, 0., None],
                      front=[None, 0., None],
                      bottom=GEO.LecodeIsostasy(reference_mat=mantle))

Note that the values may be slightly different.

[HanyMKhalil]

Dear Romain,

Yes it works perfect for me, I changed the values and plot it. just to make sure so in your line the 2.0 represent the (maximum velocity-minimum velocity), while the 2.5 represents the minimum velocity. so that it goes from 2.5 to 4.0, right?

[HanyMKhalil]

sorry it goes from 2.5 to 4.5 right?

[rbeucher]

I have added a Jupyter notebook to the @underworld-community/contributions repository. It is the setup I used to reproduced Mondy et al, 2018, I hope this help:

https://github.com/underworld-community/contributions

Romain

[HanyMKhalil]

Perfect, Thanks Romain so much 💥💥

[rbeucher]

It might need some work but that will get you started.

[HanyMKhalil]

Thanks a lot , this will surely help as I am implementing some scenarios to test💥💥