New, cutde okada code version seems to fail for near-vertical, strike slip faults.

[alomax]

I am generating a stack of Receiver_Horizontal_Profile's to create a 3D cube of ΔCFS. A vertical slice for a normal fault case looks fine, but the results seem clearly wrong for identical, right-lateral, strike-slip source and receiver faults dipping 89.9 deg or 89.0 deg, but maybe OK for 80.0 deg dip. See figures attached here. The results for 89.9 deg were fine for the previous version of Elastic_stresses_py (bottom figure).

If this is an issue (and I am not doing something wrong), I can give more details of my processing (a complicated bash script running all the layers...).

Thanks, Anthony

[kmaterna]

Well this is interesting... This seems like a real issue that I'd like to fix. But I haven't been able to reproduce this. I ran several examples that have vertical sources, vertical receivers, and/or vertical geometries on the receiver_horizontal_profile. I tried it on the rectangle-okada-wrapper version and the cutde version. Each time I get deformation results that have no diff, and stress results that have a small diff, only in the 5-6th decimal place, or <0.01 Pa.

# Source: strike rake dip magnitude faulting_type fault_corner_lon fault_corner_lat fault_corner_depth
# Receivers: strike rake dip length width updip_corner_lon updip_corner_lat updip_corner_dep
# General: PR1 FRIC lon_min lon_max lon_zero lat_min lat_max lat_zero

General: 0.250 0.40 -125.80 -122.60 -124.50 39.30 41.70 40.30
Source_WC: 228 -2 79.0 6.8 SS -125.134 40.829 16.4 # the catalog case.
Source_WC: 228 -2 89.9 6.8 SS -125.134 40.829 16.4 # vertical case
Receiver: 355 90 89.9 140 120 -125.16 40.3 10.8
Receiver_Horizontal_Profile: 1 15 89 0 -124.8 40.3 100 100 0.5

Would you be able to share more about what you're doing in that script? Looks interesting.

[alomax]

OK. Thanks for checking so thoroughly and quickly.

Attached in test.zip is a directory with README_run_test.txt for running the two versions for a receiver_horizontal_profile and plotting the results. I get clearly different plots for the two versions: horizontal_profile_stresses_???.png.

I am running and plotting with the https://github.com/kmaterna/Elastic_stresses_py/issues/17 fix, and also I modified the plotting of axes and colors in PyCoulomb/output_manager.py -> map_horiz_profile(), so your plots will be different. I ran the new version without the fix and plotting changes, this gives the plot in horizontal_profile_stresses_NEW_NO_LONG_FIX.png, which should match what you get...

Source and receivers faults are the same:

Receiver_Horizontal_Profile: 10.000  140.0 89.9 180.0  -120.45 35.9  80.0 80.0 0.5
Source_FM: 140.0 180.0 89.9  -120.45 35.9  20.0  6.4

I also attach run_elastic_stresses_cube.bash, the complex script that creates a 3D cube of ΔCFS. This should not be relevant to the test, but if you are interesting in what I am doing I would be happy to discuss by e-mail alomax@free.fr

Thanks!

test.zip

run_elastic_stresses_cube.bash.zip

[kmaterna]

Thanks for this! I was able to reproduce the figures.

TLDR: I guess this is why we can't have nice things. Can't simultaneously have triangular dislocation elements on Python 3.11+ and also have point sources with near-vertical dip that don't explode.

Longer version: in migrating to cutde, I also needed to replace the internal concept of a point source (Source_FM) with a very small rectangle (see io_intxt in 93802c24577402d085d54606f6a50895b9c8bfe2); I chose 0.5 m x 0.5 m. I am finding cutde to be unstable when the sources are very small and steeply dipping, i.e., small surface projection. I suspect it's floating point errors or stress concentrations when the plane is small and close to vertical, but I don't know for sure. I tested, and it turns out okada_wrapper is more stable for steeply-dipping-1m-x-1m-rectangles, but still a little finicky compared to the official point source.

The point source was included in dc3d0.f and okada_wrapper but not really carried forward into cutde. In a recent email, Ben suggested that I could translate the point-source function directly into Python, since it's significantly shorter and faster than the rectangular dislocations in dc3d.f. But I'm not able to do that right now, and at 1100 lines of fortran, that file is too long for chatGPT to translate :) Maybe a future student project?

So what about your example? I got stable results when I changed the 89.9°-dip-point-source into a similar rectangle with 200+ m dimensions. I would do this with Source_Patch. But, the finite patch size required to avoid a numerical explosion depends on dip, as mentioned. At 89.0° dip, the source patch needs to be about ~4 m to achieve stability. Below 80° dip, the tiny Source_FM still works as expected. This is really annoying. Perhaps you could live with an assumption of a moderate fault patch instead of a true point. It also means that you'll need to calculate your own Source_Patch slip using the shear modulus. We could have a new source format with length, width, and magnitude for convenience. Is there a plausible finite-fault assumption for your Parkfield study?

Long term, I think it's unreasonable for me to simply not support vertical strike-slip beach balls, which are common in science. This needs to be a documented caveat or eventually get fixed with future improvements.

Perhaps I can write an automatic check for very small, pathological, vertically-dipping sources and warn the user? Or we stay on okada_wrapper forever? (not my goal)

[kmaterna]

tbenthompson thread might be of interest to you. the point source strikes again (although maybe our fault for using tiny triangles at 89.9° dip). In general, I do appreciate moving to cutde btw -- it made the code faster.

[alomax]

Hi @kmaterna

Thanks for all the further investigation and details. I think I understand the problem(s), and all pretty much agrees with what I see.

My only additional comment would be that I recall singularities with rectangles, perhaps at their edges and corners - is this not why a tapered slip is used in many papers? (With a point source, there is also a singularity, at the source point, I have been masking this.) For this reason, for calculation efficiency, and perhaps for theoretical elegance, for my application* I think a point source is best. So I will probably stay with the old version for the time being.

I looked at the DC3D.f point source code under DC3D0() - it is not alot, but tedious enough FORTRAN that one would have to be systematic and careful translating by hand. Maybe the code already exists in C somewhere? Otherwise, is it better to create a C version or might pure Python be fast enough? In the FORTRAN I see no nested loops, only many loops from 1 to only 12 - these might be fine in Python or could perhaps be converted to vector operations with numpy.

• I will send you an e-mail describing my application, for your interest and because your comments would be very helpful!

Thanks, Anthony

[kmaterna]

@alomax if it's not too much trouble, could you switch to the current commit and test it out on your workflow?
I just used chatGPT to translate DC3D0.f into Python and added the code into this repo, bringing back the internal concept of a point source. It passes my tests on point-source displacements and strains versus the former Okada_wrapper function. I'm curious to hear how this version works for your applications.
I also streamlined the locations where geographic to cartesian coordinate transformations get done, so we can then easily take your other changes (assuming this works).

[alomax]

Hi @kmaterna It seems to work better, though there is distortion and missing negative lobes, e.g. in the horizontal plot at 2km above the point source:

Previous version: Current commit:

Vertical sections (as in first comment above for this issue)

Previous version: Current commit:

Note that I re-added my lat/lon fixes for the plots, but a conflict here with you new code might add distortion, but not loose lobes...

What cases do your tests on point-sources cover? Maybe some issue with strike not being along an axis, or dip not exactly 90.0 (I am using 89.9)?

I also had to install pip install okada_wrapper since https://github.com/kmaterna/Elastic_stresses_py/blob/master/PyCoulomb/init.py imports run_okada_wrapper which imports okada_wrapper - is this necessary?

Thanks, Anthony

[kmaterna]

thank you for pointing this out! I am essentially reproducing it -- a couple of the strain tensor components have some sign errors etc., and I didn't catch them in the first version of the tests. Now the tests are more sensitive and should tell me when I've caught the bugs, hopefully pretty soon. thanks!

[alomax]

Good to hear! Thanks for the update. Anthony

[kmaterna]

I found it -- and now the tests pass, for displacement AND strain. I had forgotten to convert the displacement gradient tensor to the strain tensor!!! Does it work?

[alomax]

It seems to work perfectly!

I generated a CFS cube with the new code and plotted the vertical section as above, it looks identical (pixel to pixel) to the FORTRAN version. I also calculated further results using the stress cube and they are identical.

Thanks,

Anthony