GenericMappingTools / GenericMappingTools/gmt

gravfft problems computing moho geopotential

Open
#7,255 4 comments 0 reactions 0 assignees View on GitHub
bug
Dominant language
C
Stars
979
Forks
414
Avg merge
17h 26m
Merged PRs (30d)
54

Description

**Description of the problem**

The gravfft module with the -T+m option produces an output moho geopotential grid with artifacts. I've also noticed the module behavior with this option activated is not deterministic, and may result in a segfault.

The module seems to work fine with the -D option or with the -Q and -T (no moho) options.

**Full script that generated the error**

To generate the moho geopotential grid:
```
gmt grdcut @earth_relief_01m_p -R-120/-110/-40/-34 -Gbat.nc
gmt gravfft bat.nc -T1000/2700/3300/1035+m -Z10000 -E1 -Gmoho_g.nc
```
The geographic range, flexural parameters, number of Parker expansion terms, and output potential field (-F option) are not important, and similar artifacts are produced regardless.

**Actual outcome**

Here is a visualization of the moho_g.nc grid:
![grav_comp_bad](https://user-images.githubusercontent.com/38174025/220214887-39e0c520-5772-40a0-9a8a-505867399b3a.png)

Here is an example of the occasional segfault (running the above gravfft command with -V):
```
gravfft [INFORMATION]: Allocates memory and read data file
gravfft [INFORMATION]: netCDF grid bat.nc has a default CPT geo.
gravfft [INFORMATION]: Round-off patrol changed geographic grid increment for longitude from 0.0166666666666666456 to 0.0166666666666666664
gravfft [INFORMATION]: Round-off patrol changed geographic grid increment for latitude from 0.0166666666666666699 to 0.0166666666666666664
gravfft [INFORMATION]: Reading grid from file bat.nc
gravfft [INFORMATION]: Data dimension 600 360 time factor 986.691 rms error 1.18196082e-05 bytes 1728080
gravfft [INFORMATION]: Highest speed 640 384 time factor 709.75788 rms error 1.00090156e-05 bytes 1966160
gravfft [INFORMATION]: Most accurate 768 384 time factor 836.96326 rms error 9.44472165e-06 bytes 2359344
gravfft [INFORMATION]: Least storage 600 360 time factor 986.691 rms error 1.18196082e-05 bytes 1728080
gravfft [INFORMATION]: Selected FFT dimensions for the overall best solution (fastest).
gravfft [INFORMATION]: Grid dimensions (n_rows by n_columns): 360 x 600 FFT dimensions: 384 x 640
gravfft [INFORMATION]: Extend grid via copy onto larger memory-aligned grid
gravfft [INFORMATION]: Mid value removed from real component: -2653.25 Variance reduction: 92.52
gravfft [INFORMATION]: Grid (real component) extended via edge-point symmetry at all edges, then tapered to zero over 100 % of extended area
gravfft [INFORMATION]: Level used for upward continuation: 2653.25
gravfft [INFORMATION]: Forward FFT...
gravfft [INFORMATION]: 2-D FFT using FFTW
gravfft [INFORMATION]: Picking a (probably sub-optimal) FFTW plan quickly.
gravfft [INFORMATION]: Inverse FFT...
gravfft [INFORMATION]: 2-D FFT using FFTW
gravfft [INFORMATION]: Picking a (probably sub-optimal) FFTW plan quickly.
gravfft [INFORMATION]: Demultiplexing complex grid before writing can take place.
gravfft [INFORMATION]: Evaluating Parker for term = 1
gravfft [INFORMATION]: 2-D FFT using FFTW
gravfft [INFORMATION]: Picking a (probably sub-optimal) FFTW plan quickly.
gravfft [INFORMATION]: 2-D FFT using FFTW
gravfft [INFORMATION]: Picking a (probably sub-optimal) FFTW plan quickly.
gravfft [INFORMATION]: Write Output...
gravfft [INFORMATION]: Writing grid to file moho_g.nc
gravfft [INFORMATION]: No valid values in grid [moho_g.nc]
ERROR: Caught signal number 11 (Segmentation fault: 11) at
0 libhdf5.200.dylib 0x000000010142cd1b H5FL__blk_gc_list + 71
1 ??? 0xffffffff7fffffff 0x0 + 18446744071562067967
Stack backtrace:
0 libgmt.6.4.0.dylib 0x00000001007205a5 sig_handler_unix + 210
1 libsystem_platform.dylib 0x00007ff81f4bbdfd _sigtramp + 29
2 libgmt.6.4.0.dylib 0x00000001009ade1d gmt_token_check.var_token + 156979
3 libhdf5.200.dylib 0x000000010142cdec H5FL__blk_gc + 53
4 libhdf5.200.dylib 0x000000010142c486 H5FL_garbage_coll + 42
5 libhdf5.200.dylib 0x000000010142c2ac H5FL_term_package + 32
6 libhdf5.200.dylib 0x0000000101361424 H5_term_library + 3119
7 libsystem_c.dylib 0x00007ff81f39edeb __cxa_finalize_ranges + 416
8 libsystem_c.dylib 0x00007ff81f39ebfe exit + 35
9 libdyld.dylib 0x00007ff81f4b2375 _ZNK5dyld416LibSystemHelpers4exitEi + 11
10 dyld 0x0000000101c75558 start + 504
```

**Expected outcome**

Current workaround is to compute the flexural topography of the moho (-T and -Q options) then compute the geopotential with the -D option:
```
gmt grdcut @earth_relief_01m_p -R-120/-110/-40/-34 -Gbat.nc
gmt gravfft bat.nc -T1000/2700/3300/1035 -Z10000 -Q -Gmoho_d.nc
gmt gravfft moho_d.nc -D600 -E1 -Gmoho_g.nc
```
Here is a visualization of that grid:
![grav_comp_good](https://user-images.githubusercontent.com/38174025/220214995-7d9c26d8-4ff2-4db9-b6c5-266ad4c3ed2f.png)
note that both grids are shown with the same colorscale (-40 - 40 mgal).

**System information**

* Operating system: macOS 12.6
* GMT version (`gmt --version`): 6.4.0

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.