[WIP] Allow for Oblique Slabs/Faults Relative to Surface

[danieldouglas92]

Currently, faults and slabs dip from the surface at an angle orthogonal to the surface coordinates. In this PR, I want to add the ability for a user to specify an "obliquity vector" which tells the world builder what direction to dip the feature relative to the surface coordinates. An example of why this is useful can be seen in the images below:

The convergence of two plates on Earth is rarely perfectly orthogonal to the plate boundary, and in subduction zones there is varying degrees of oblique convergence of a plate relative to the plate margin. Currently, if a user wants to initialize a subducting plate that has undergone oblique convergence, the subducting plate will look like this:

The trench is rotated 45 degrees relative to the convergence direction, and so WorldBuilder by default also rotates the slab so that it remains orthogonal to the trench. If I specify an "obliquity vector" that is in the positive x-direction, then the slab looks like this:

which maintains the same convergence direction even with an oblique trench.

This is implemented by modifying the distance_point_from_curved_planes function in utilities.cc, and consists of two main chunks of code:

1. Most of the time the current implementation will use a bezier_curve to find a point on the trench/fault that is closest to the check point. In this case, I use the check point and the obliquity_vector to create a line, which I then incrementally shift towards the bezier curve. If it intersects the bezier_curve, then this intersection point becomes the new "closest point" on the bezier_curve. If it doesn't converge, then the "closest point" is set to be a NAN.

2. The bezier_curve misses points that are near the terminus of the surface coordinates. The way this is currently handled is definitely not generalized. It won't work for spherical and I don't think it will even work well for more complicated Cartesian geometries. I think this will require modifying the closest_point_on_curve_segment function in bezier_curve.cc. Right now, my idea is to check what the min/max values of the surface coordinates are, draw a line between the two extrema, and see if the check point (with the obliquity vector) will intersect this line. If it does, I check to see if the intersection point lies within the min/max values.

This is still very much a work in progress, I still need to add the "obliquity vector" as an input parameter that gets passed to distance_point_from_curved_planes, but I just wanted to open this PR so that I can get feedback on it!

@MFraters

[danieldouglas92]

I've added an argument to distance_point_from_curved_planes that defines the direction of obliquity. I made the argument have a default value that was (DSNAN, DSNAN) so that all of the unit_tests should pass without modification. This works for simple linear trenches, but does not work for curved trenches yet.

[MFraters]

Cool, nice minimal code change and well documented! A few initial comments and questions.

1. I think the algorithm can maybe be simplified, which may help with the edges. I might be wrong, but I think you should use the closest point on the curve just as a guide. If it is not on the curve, use the endpoint which is closest to compute the distance. That would bet most of the code outside of the if (!std::isnan(closest_point_on_line_2d[0]) statement and should fix the corner cases in complex Cartesian and spherical cases. I haven't worked it out, so might be wrong ;)
2. Can you show what happens in spherical?
3. I think in the final version this parameter should become part of the sections. Do you agree? Just to note, that doesn't need to happen in this pull request, but we need to keep that in mind when making this.
4. The slab seems to extend a bit beyond the plate, or is that just a visualization thing?

[MFraters]

why 100? What if we don't converge before?

[danieldouglas92]

Yea I don't really have a good reason for picking 100, I was torn between making the number of iterations an input parameter or not. I just chose 100 as a starting number, but if you think it should be higher/a user input parameter I would be ok with changing this.

[MFraters]

I don't think this should be a user defined value, since this is deep in the internals of the code. We need to make sure we always get a good solution within an acceptable tolerance.

I don't know what the right number is, but I think that is we do not converge, it should fail and not silently continue. This lets us (and the user) know that something went wrong, and that they get the wrong solution, so we can investigate it.

[danieldouglas92]

I tested adding an assert but the iteration actually fails quite often because there are points that shouldn't converge just by the nature of migrating check points along the obliquity vector

[danieldouglas92]

@MFraters Thanks for taking a look at this so quickly! I addressed the comments on the code you left so far, and for the 4 points you outlined I'll work on addressing in the next couple days. Hopefully you are right that using the end of the bezier curve will work and simplify the code significantly.

For your points 2-4:

I actually haven't tried doing this in spherical, the current logic (in the loop outside of if (!std::isnan(closest_point_on_line_2d[0])) implicitly assumes a cartesian geometry since it just makes a straight line, whereas in spherical it would have to make an arc. If what you suggest works though, then what I have should work in spherical.

I do agree that this should become part of the sections, this would allow a lot more flexibility in the geometry of the slab along-strike.

I'm not sure what you mean by the slab seems to extend beyond the plate?

[MFraters]

Great, thanks!

I'm not sure what you mean by the slab seems to extend beyond the plate?

[danieldouglas92]

I changed the default value for the obliquity vector to be [inf, inf] (I was having issues with NANs). I also updated the flow of the code to first check to see if the point is near the edges of the Bezier curve. I'll work on testing this for curved trenches in Cartesian, and then start testing it in spherical geometries.

The slab extending beyond the plate is a visualization artifact!

[danieldouglas92]

@MFraters

This seems to be working with a curved trench as well!

I'll post what happens for spherical slabs soon

[danieldouglas92]

Currently, the spherical case is not working. Because I only personally needed the Cartesian case I've been less motivated to fix this, but I'm hoping to have time soon to try to figure out what's going on. This is what the slab looks like for the spherical case when trying to specify that the obliquity vector is to the right:

The slab at the very bottom of the trench looks correct (below the tear), but the tear and everything above it is incorrect

[github-actions]

Benchmark	Main	Feature	Difference (99.9% CI)
Slab interpolation simple none	1.294 ± 0.011 (s=344)	1.308 ± 0.007 (s=350)	+0.9% .. +1.2%
Slab interpolation curved simple none	1.290 ± 0.010 (s=341)	1.303 ± 0.009 (s=356)	+0.8% .. +1.2%
Spherical slab interpolation simple none	1.259 ± 0.009 (s=345)	1.268 ± 0.007 (s=369)	+0.5% .. +0.9%
Slab interpolation simple curved CMS	1.326 ± 0.007 (s=347)	1.336 ± 0.006 (s=332)	+0.6% .. +0.9%
Spherical slab interpolation simple CMS	1.622 ± 0.025 (s=271)	1.629 ± 0.008 (s=285)	+0.1% .. +0.8%
Spherical fault interpolation simple none	1.273 ± 0.008 (s=365)	1.282 ± 0.009 (s=342)	+0.6% .. +0.9%
Cartesian min max surface	3.049 ± 0.011 (s=149)	3.071 ± 0.011 (s=148)	+0.6% .. +0.9%
Spherical min max surface	7.856 ± 0.055 (s=62)	7.811 ± 0.059 (s=55)	-1.0% .. -0.1%

[danieldouglas92]

@MFraters So the code that I had from before the hackathon worked in Cartesian but not robustly, I had only tested rightward dipping slabs, and after doing some more tests the old code failed if the obliquity vector and the reference point were not "compatible" enough. The point of failure seemed to be where I tried to determine which direction along the obliquity vector (positive or negative) the check point should be migrated (see commented lines 574-576 in utilities.cc. Honestly, I feel like what I had on lines 574-576 should have worked, and I'm not sure why it isn't robust, but iterating over both positive/negative directions seems to work for all the cases I've thrown at it so far.

I still can't get spherical to work, and for now I'm not sure what else to try.

[MFraters]

Could you describe what you mean with not working robustly? I have some ideas what could potentially cause issues, but it would be good to have a chat I think.

[danieldouglas92]

@MFraters The robustness was that if the reference point changed sides of the trench, then the obliquity vector would have to change signs in order to still produce a slab.

For example if you have a vertical trench at x = 0, then if the obliquity vector was [1, 0] and the reference point was [10e3, 0], this would work, but if the obliquity vector was [1, 0] and the reference point was [-10e3, 0], this would result in no slab being present.

I had tried adding additional checks to modify the direction that the check point is migrated along the obliquity vector based on which side of the feature the reference point was on, but this didn't seem to change anything. And I don't think that it should change anything because fundamentally we just want to move towards the trench regardless of where the reference point is. There must be something wrong with the logic that I implemented to determine whether the obliquity vector is pointing towards the trench.

What I currently do is create a vector from the closest point on the curve to the check point, and then check the sign of the dot product between this vector and the obliquity vector. Positive dot product means parallel, so move negatively along the obliquity vector, and vice versa.

[danieldouglas92]

@MFraters Also I messaged you on element, but it would be great to chat about this and also the spherical issue

[danieldouglas92]

@MFraters I think I've got something working that is much cleaner. The issue that I discussed with you on zoom for why I couldn't always determine the orientation of the vector is because the distance that was being returned was signed. After I used the absolute value of the distance, the method that I had originally worked.

There is just one thing that is pretty hacky right now that it would be great to get your feedback on. The distance that I'm using to determine convergence is also tricky, if I make it too strict then in Cartesian the iteration won't reliably converge creating gaps in the slab. For spherical though, since the distances returned are not in meters, the tolerance has to be very small. So to get around this, I set the tolerance differently for either cartesian or spherical geometries. Is there a better way to set this tolerance for both geometries?

I'll write up proper tests and demonstrate that this works with a complicated trench geometry after I get some lunch, and then also add the functionality for faults. But I think this is basically good now.

[MFraters]

There is just one thing that is pretty hacky right now that it would be great to get your feedback on. The distance that I'm using to determine convergence is also tricky, if I make it too strict then in Cartesian the iteration won't reliably converge creating gaps in the slab. For spherical though, since the distances returned are not in meters, the tolerance has to be very small. So to get around this, I set the tolerance differently for either cartesian or spherical geometries. Is there a better way to set this tolerance for both geometries?

I think it is fine to have different tolerances for spherical and cartesian. They are in different units, so it doesn't surprise me that they need different tolerances. I think something similar is done in several places in the library.