along interface velocity for subducting slab

[lhy11009]

Add an along surface velocity option for the subducting slab. This will trace the surface of the subducting slab and rotate the velocity to align with the surface. The free parameter is the magnitude of the velocity.

Here is the layout of the test, which connects a horizontal velocity in the plate to the along-surface velocity in the slab.

This is tentative. Eventually, we might want to tweak the velocity in each of the slab segments, but this would be an easy start.

[coveralls]

Coverage Report for CI Build 32893349143

Coverage increased (+0.01%) to 98.314%

Details

• Coverage increased (+0.01%) from the base build.
• Patch coverage: 104 of 104 lines across 3 files are fully covered (100%).
• No coverage regressions found.

Uncovered Changes

No uncovered changes found.

Coverage Regressions

No coverage regressions found.

---

Coverage Stats

Relevant Lines:	12041
Covered Lines:	11838
Line Coverage:	98.31%
Coverage Strength:	815012.44 hits per line

---

💛 - Coveralls

[github-actions]

Benchmark	Main	Feature	Difference (99.9% CI)
Slab interpolation simple none	1.186 ± 0.015 (s=378)	1.192 ± 0.010 (s=381)	+0.3% .. +0.8%
Slab interpolation curved simple none	1.185 ± 0.010 (s=364)	1.192 ± 0.010 (s=396)	+0.3% .. +0.7%
Spherical slab interpolation simple none	1.161 ± 0.006 (s=404)	1.166 ± 0.007 (s=372)	+0.3% .. +0.6%
Slab interpolation simple curved CMS	1.228 ± 0.012 (s=350)	1.238 ± 0.017 (s=382)	+0.6% .. +1.2%
Spherical slab interpolation simple CMS	1.561 ± 0.011 (s=278)	1.580 ± 0.027 (s=297)	+0.8% .. +1.6%
Spherical fault interpolation simple none	1.217 ± 0.012 (s=375)	1.227 ± 0.010 (s=364)	+0.6% .. +1.0%
Cartesian min max surface	2.853 ± 0.020 (s=161)	2.856 ± 0.021 (s=157)	-0.2% .. +0.4%
Spherical min max surface	8.676 ± 0.105 (s=61)	8.589 ± 0.107 (s=46)	-1.8% .. -0.2%

[MFraters]

Cool! Yes, the next step would be have a velocity per segment something like the angles, either [1.5] or [1.5,2.2].

[MFraters]

Looking at it a bit more, I am not sure how to get the result you got, I get the following result.

with the following wb file:

{
  "version":"1.2",
  "coordinate system":{"model":"cartesian"},
  "cross section":[[0,450e3],[10e3,450e3]],
  "features":
  [
   { 
     "model":"subducting plate", "name":"Slab", "dip point":[0,0],
     "coordinates":[[1500e3,1000e3],[1600e3,350e3],[1500e3,0]],
     "segments":
       [
         {"length":300e3, "thickness":[300e3], "top truncation":[-100e3], "angle":[0,60],
          "composition models":[
            {"model":"uniform", "compositions":[3], "max distance slab top":50e3},
            {"model":"uniform", "compositions":[2], "min distance slab top":50e3, "max distance slab top":100e3}]},
         {"length":500e3, "thickness":[300e3], "top truncation":[-100e3], "angle":[60,20]}
       ],
     "sections":
       [
         {"coordinate":0,
          "segments":[{"length":300e3, "thickness":[300e3], "top truncation":[-100e3], "angle":[0,60]},
                      {"length":400e3, "thickness":[300e3], "top truncation":[-100e3], "angle":[60]}],
                       "composition models":[{"model":"uniform", "compositions":[1],"max distance slab top":100e3}]}
       ],
     "temperature models":[{"model":"mass conserving", "density":3300, "spreading velocity":0.02, "subducting velocity":0.02,
                            "ridge coordinates":[[[5e6,-1],[5e6,2000e3]]],"coupling depth":50e3,
                            "min distance slab top":-200e3, "max distance slab top":300e3}],
     "composition models":[{"model":"uniform", "compositions":[2], "max distance slab top":100e3}],
      "velocity models": [{"model":"along surface", "velocity": [1,2,3]}]
   }
  ]
}

and world builder file:

# output variables
grid_type = cartesian
dim = 2
compositions = 6
vtu_output_format = ASCII

# domain of the grid
x_min = -1000e3
x_max = 2000e3 
z_min = 0 
z_max = 600e3 

# low res grid properties
n_cell_x = 600 
n_cell_z = 100

# shown grid properties
# n_cell_x = 2400 
# n_cell_z = 400

I also thought the user would provide a scalar instead a vector, and maybe later on we could add an angle offset.

[lhy11009]

If we add a sign modification to Vx depending on the relative position of the dip point and the trench point, and if dip point have smaller x than the trench point, we change the sign of Vx.
With Menno's help, I also adjusted the angle I use to be the local angle; now the velocity perfectly aligns with the surface.
This angle is added to the distance_from_plane object

[lhy11009]

reference_on_side_of_line doesn't seem to give me the structure I want. I return it in the "distance_from_curved_plane" function and multiply it to the x component:

As the result, I can get one test right:

But not the other one, which has a slab dipping in the opposite direction:

[MFraters]

hmm, can you take the closest_trench_point from PointDistanceFromCurvedPlanes to compute a direction vector?

[lhy11009]

Yes, I think that works. Now the two tests are merged into one, where two slabs dip in different directions, and the velocity is along the surface. of the features. I think this fulfills our purpose here.

[lhy11009]

@MFraters can you take another look at this? I think all tests have passed.

[MFraters]

Very cool! Thanks for working on this.

The code looks good to me. Have you also tested this in 3D and in spherical, to make sure that also works as expected (and maybe add tests for that as well)?

[lhy11009]

@MFraters , I follow the cookbook of 3d_spherical_subduction, and add a velocity of 1.0 with the "along surface" model.
In the figure below, I plot the contour of temperature as well as the velocity vectors.
The slab dip decreases to 0 at the end, but notice that the vector doesn't track this change:

[lhy11009]

It seems this angle I use is not tracking the angle of the slab curvature in this case, which worked fine in 2-D test.

[MFraters]

Can you test 3D cartesian, 2D spherical and 3D spherical seperately? That might make it easier to track down whether the issue is 3D or spherical coordinates.

[lhy11009]

I seem to come across another bug testing the Cartesian model. When I take the 3d_cartesian_curved_subduction cookbook and run the model, initially it ran through and produced a vtu outputs.

But when I add the "along surface" velocity to the subducting plate feature, I received an Arithmetic exception.
I output the variables in gdb and the "angle" I computed is "nan" value, and the point this happens (650000, 1000000, 1000000).

Here is the layout of the point (white point) with a 1473 K iso-volume plot of the slab.

, and from another angle:

[lhy11009]

The "angle" ultimately comes from this "interpolated_angle_top)

[lhy11009]

I tried to bypass this error if the error is related to a few points with this added line.

[lhy11009]

@MFraters , But then I got 0 velocity everywhere, meaning the angle is nan everywhere for this case. Is it this wb have multiple sections in the slab, and that causes the problem? I can try a different test file later, but maybe you could first comment on this issue.

[MFraters]

I read through the comments, and it looks like we need to take a look at when and why this angle is actually computed. Since I am at a conference right now I don't have the time to really look deep into the code right now, but I will try to do this after the conference.

[lhy11009]

@MFraters, would you have a time next week to look at this? Perhaps we could get this and my other PR merged before the Hackathon.

[MFraters]

Sorry it took me a bit to get around to this. The issue with the segfault has to do with an early return case around line 615 in utilities.cc. You need to add return_values.angle = total_average_angle; there.

Your current approach to get the directly only works in 2D though. I think you should export the x_axis variable (around line 684). I think that might get you the right direction both in Cartesian and spherical.

[lhy11009]

@MFraters Thanks for your reply. I would need to get back to this after the ASPECT Hack.

[lhy11009]

Good, I seem to have fixed the problem, at least for Cartesian geometry. I also added a new angle to the PointDistanceFromCurvedPlanes struct that tracks the direction of the x_axis. This allows us to predict the velocity direction in 3D as well. Below are test results in 2-d and in 3-d.

I will then test the implementation in Spherical geometry.

@MFraters I have to get to the bottom of the distance_point_from_curved_planes function to reach this point. This function is indeed a very thoughtful implementation and really helps me understand how GWB works in general. I do still have a few small questions to ask, will you have a 30 mins meeting time next week?

[lhy11009]

Here are the results from the 3-D spherical test:

To achieve this, I added two additional functions to the utilities. I found that the utilities already contain coordinate transformations between Cartesian and spherical geometries, but seem to be missing the corresponding vector transformations.

These new functions transform a vector between its local (r, north, east) components and its global Cartesian (x, y, z) components at a given position in spherical geometry. The implementation is essentially a matrix multiplication using the local spherical basis vectors.

I think these functions could also be generally useful beyond this particular application.

Why are these functions needed here?

The first step is to determine the dipping angle and the azimuth angle, where the azimuth angle is derived from the a_axis. The natural way to define an azimuth angle is as an angle on a surface of constant radius. In Cartesian geometry, this is simply the x-y plane, while in spherical geometry it corresponds to the local north-east (n-e) plane. This is where the transformation from a Cartesian vector to its local spherical components is needed.

When computing the velocity, it is also natural to first express it in the local coordinate system. In Cartesian geometry, these are simply (Vx, Vy, Vz), while in spherical geometry they are (Vr, Vn, Ve). Therefore, after computing the velocity in the local spherical coordinate system, we need the inverse transformation from the local spherical vector back to the global Cartesian (x, y, z) vector to obtain the final velocity.

[MFraters]

very cool pictures! Yes, lets discuss.