gravfft problems computing moho geopotential

[hughaharper]

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:

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:

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

[welcome]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.

[joa-quim]

Yep, and I'm surprised how you even managed to have your first figure. gcc is really permissive in accessing invalid memory.

What happens is that at gravfft.c#L693

		GMT_Report (API, GMT_MSG_INFORMATION, "Inverse FFT...\n");
		if (GMT_FFT (API, Grid[0], GMT_FFT_INV, GMT_FFT_COMPLEX, K))
			Return (GMT_RUNTIME_ERROR);

we do the inverse transform and the grid Grid[0] looses the imaginary part and hence shrinks to half size. And when it later is fed to the fft function at line 765 it crashes because Grid[0] is too small (lost the imaginary part)

if (GMT_FFT (API, Grid[0], GMT_FFT_INV, GMT_FFT_COMPLEX, K)) {

gravfft went into several cycles of rewriting/improving but this case seems to have escaped our tests and now it's not obvious on what to do to fix it.

[PaulWessel]

I will have to find some time on the weekend to look at this and fix the issue. Would be nice to have a simple test script to aim for.

[joa-quim]

I have actually found what seemed to be a fix, but the two solutions differ by ~13 mGal and the difference is not constant.