Feature #133

surface fails on prime grid dimensions

Added by Florian almost 6 years ago. Updated about 5 years ago.

Status:FeedbackStart date:2012-08-22
Priority:NormalDue date:
Assignee:-% Done:

0%

Category:-
Target version:-
Platform:

Description

Consider this example:

# create a test grid:
grdmath -R-21/21/-15/15 -I0.3 X Y HYPOT DUP 2 MUL PI MUL 8 DIV COS EXCH NEG 10 DIV EXP MUL = sombrero.nc
grd2cpt -E50 sombrero.nc > sombrero.cpt
grdimage -E300 -Csombrero.cpt -JX18c/0 sombrero.nc > sombrero.ps

# create subset of test grid:
grdmath -R-21/21/-15/15 -I0.3/1.5 1 = mask_1.nc
grdmath -R-21/21/-15/15 -I1.5/0.3 1 = mask_2.nc
grd2xyz mask_1.nc mask_2.nc | xyz2grd -R-21/21/-15/15 -I0.3 -Gmask.nc
grdmath sombrero.nc mask.nc OR = sombrero_sub.nc
grdimage -E300 -Csombrero.cpt -JX18c/0 sombrero_sub.nc > sombrero_sub.ps

# This works:
grd2xyz -s sombrero_sub.nc | surface -Gsombrero_sfc.nc -R-21/21/-15/15 -I0.01 -C0.01 -T0.2 -V
surface (surface.c:1881): Grid domain: W: -21 E: 21 S: -15 N: 15 nx: 4200 ny: 3000 [gridline registration]
surface (gmtapi_util.c:401): Reading Data Table from Standard Input stream
surface (surface.c:659): Minimum value of your dataset x,y,z at: surface (surface.c:660): -2.40000009537 3 -0.675758719444
surface (surface.c:661): Maximum value of your dataset x,y,z at: surface (surface.c:662): 0 0 1
surface (surface.c:1937): Grid    Mode    Iteration    Max Change    Conv Limit    Total Iterations
surface (surface.c:1021):  600    D           1    -0.229151425536    1.66666666667e-05             1
surface (surface.c:1021):  120    I          69    8.26775888864e-05    8.33333333333e-05            70
surface (surface.c:1021):  120    D          20    3.93378150874e-05    8.33333333333e-05            90
surface (surface.c:1021):   24    I          42    0.000401469439945    0.000416666666667           132
surface (surface.c:1021):   24    D          41    0.000375375806668    0.000416666666667           173
surface (surface.c:1021):    8    I           5    0.000979155245687    0.00125           178
surface (surface.c:1021):    8    D           7    0.00118705111161    0.00125           185
surface (surface.c:1021):    4    I           1    0.00228155047846    0.0025           186
surface (surface.c:1021):    4    D           4    0.00197947568642    0.0025           190
surface (surface.c:1021):    2    I           1    0.00135280234289    0.005           191
surface (surface.c:1021):    2    D           1    0.00159548455998    0.005           192
surface (surface.c:1021):    1    I           1    0.000659909000875    0.01           193
surface (surface.c:1021):    1    D           1    0.000333866291778    0.01           194
surface (surface.c:1165): Fit info: N data points  N nodes    mean error    rms error    curvature
surface (surface.c:1167):         5281    12607201    0    0    6898.058036
surface (gmtapi_util.c:1839): Writing grid to file sombrero_sfc.nc

# Visual check:
grdimage -R-21/21/-15/15 -E300 -Csombrero.cpt -JX18c/0 sombrero_sfc.nc > sombrero_sfc.ps

# This fails:
grd2xyz -s sombrero_sub.nc | surface -Gsombrero_sfc.nc -R-21/21/-15/15 -I0.009 -C0.01 -T0.2 -V
surface (surface.c:1881): Grid domain: W: -21 E: 21 S: -15 N: 15 nx: 4667 ny: 3333 [gridline registration]
surface (surface.c:1883): Warning: Your grid dimensions are mutually prime.
surface (surface.c:1459): Hint: Choosing nx = 4800, ny = 3456 might cut run time by a factor of 1640.6031
surface (surface.c:1459): Hint: Choosing nx = 4860, ny = 3456 might cut run time by a factor of 1592.2457
surface (surface.c:1459): Hint: Choosing nx = 4800, ny = 3600 might cut run time by a factor of 1573.3261
surface (surface.c:1459): Hint: Choosing nx = 5120, ny = 3456 might cut run time by a factor of 1537.19
surface (surface.c:1459): Hint: Choosing nx = 4860, ny = 3600 might cut run time by a factor of 1529.3905
surface (surface.c:1459): Hint: Choosing nx = 5184, ny = 3456 might cut run time by a factor of 1519.4014
surface (surface.c:1459): Hint: Choosing nx = 5000, ny = 3600 might cut run time by a factor of 1489.9449
surface (surface.c:1459): Hint: Choosing nx = 4800, ny = 3840 might cut run time by a factor of 1476.8402
surface (surface.c:1459): Hint: Choosing nx = 5120, ny = 3600 might cut run time by a factor of 1466.4519
surface (surface.c:1459): Hint: Choosing nx = 5184, ny = 3600 might cut run time by a factor of 1456.3589
surface (gmtapi_util.c:401): Reading Data Table from Standard Input stream
surface (surface.c:659): Minimum value of your dataset x,y,z at: surface (surface.c:660): -2.40000009537 3 -0.675758719444
surface (surface.c:661): Maximum value of your dataset x,y,z at: surface (surface.c:662): 0 0 1
surface (surface.c:1937): Grid    Mode    Iteration    Max Change    Conv Limit    Total Iterations
surface (surface.c:1021):    1    D          23    0.00990246274785    0.01            23
surface (surface.c:1165): Fit info: N data points  N nodes    mean error    rms error    curvature
surface (surface.c:1167):         5281    15563112    0.00139725731454    0.122792414764    85740.6146404
surface (gmtapi_util.c:1839): Writing grid to file sombrero_sfc.nc

# Visual check:
grdimage -R-21/21/-15/15 -E300 -Csombrero.cpt -JX18c/0 sombrero_sfc.nc > sombrero_sfc_fail.ps

History

#1 Updated by Paul almost 6 years ago

  • Status changed from New to Resolved

No, this is user error. We are pretty clear that one has to be thoughtful in selecting -R -I so that it results in many intermediate grid spacings. In this case you get all that information:

surface (surface.c:1883): Warning: Your grid dimensions are mutually prime.
surface (surface.c:1459): Hint: Choosing nx = 4800, ny = 3456 might cut run time by a factor of 1640.6031
surface (surface.c:1459): Hint: Choosing nx = 4860, ny = 3456 might cut run time by a factor of 1592.2457
etc etc

There is no way surface can reach a stable convergence without those intermediate steps in a reasonable time. You might be able to reach a solution by making the convergence limit very small and wait a very very long time for the changes to wash across the large grid at dx=dy=1. Much better to take the hints.

#2 Updated by Florian almost 6 years ago

  • Status changed from Resolved to Feedback

Making the convergence limit very small and waiting very long does not change the result in any way, I'm afraid.

Here are some things that should be considered:
  1. I don't find any comments or warnings on how to choose the correct -R and -I in the surface manpage
  2. Without -V there is no warning that anything might go wrong:
    grd2xyz -s sombrero_sub.nc | surface -Gsombrero_sfc.nc -R-21/21/-15/15 -I0.009 -C0.01 -T0.2
    
    This means that if you are not looking for the problem you might not even be aware that your model did not converge yet and that you have lots of artifacts in the output grid.
  3. The warning just states that there is room for increasing the runtime. It does not say that the computation might fail.
  4. Let me say that I created this example on purpose to show severe artifacts. However, the artifacts might go unnoticed if they are very subtle. Sometimes you see them after calculating the gradient. I have seen these typical checkerboard patterns in published datasets before I even knew how they are caused. This means that there are other people oblivious to the the problem.
  5. What harm does it do if surface temporarily increases the top and right boundary such that an optimal nx and ny guarantees a stable solution? Before writing the grid, modify the padding such that the output file is written with the user defined region. I believe this has also the benefit, that changing -R and -I does not change the solution anymore.

#3 Updated by Paul over 5 years ago

1. We cannot automatically change -R since then the boundary conditions are applied at the adjusted boundary. Thus it has to be an option, e.g., -E "extend region to optimize performance". I think that is a good compromise.
2. The horrible case of prime and that message should be elevated to GMT_MSG_NORMAL from GMT_MSG_VERBOSE so that it is issued without -V.
3. The -E option will allow us to put more language into man page and program usage to inform the user about the issue.
4. I dont think we want to change -I, it would mean a 2nd resampling in the end.
Seems reasonable?

#4 Updated by Florian over 5 years ago

1. We cannot automatically change -R since then the boundary conditions are applied at the adjusted boundary. Thus it has to be an option, e.g., -E "extend region to optimize performance". I think that is a good compromise.
2. The horrible case of prime and that message should be elevated to GMT_MSG_NORMAL from GMT_MSG_VERBOSE so that it is issued without -V.
3. The -E option will allow us to put more language into man page and program usage to inform the user about the issue.
4. I dont think we want to change -I, it would mean a 2nd resampling in the end.
Seems reasonable?

I'm not sure what the boundary condition is. Can you set the boundary to a specific value?

Can you explain, why the region auto-extension is/needs an extra case? I mean if I want a stable result I need to:
  1. choose -I and extend -R
  2. grdcut to original -R

So why not doing this automatically in every case? If you don't extend the region the result is bogus anyway. So why the extra option -E? My intention was not changing the spacing (-I) only the region should be extended so that for a given -I the model is guaranteed to converge. After that ofc we need to auto-crop to the region given with -R on the command line.

I understand that auto-extending a periodic region is a different case though.

#5 Updated by Paul about 5 years ago

Picking up this thread again... Since then I have added periodic boundary condition in longitude if grid is 360. Otherwise the boundary conditions are always the natural conditions for splines in tension. I guess I resist auto-changing the region because some people want to set things a particular way. It also would not be backwards compatible. We have -Q that gives the various hints to the user on better choices for nx/ny. Perhaps we could add a modifier, like -Qa, to automatically pick the top choice and extend region, etc, except in the case of 360 grids with periodic condition?

#6 Updated by Paul about 5 years ago

  • Tracker changed from Bug to Feature

Based on the comments so far I am changing this from a bug to a feature request.

Also available in: Atom PDF