Technical Document : On the Impact of Compute Grids on the 3D GIA Model Seakon

Technical Document : On the Impact of Compute Grids on the 3D GIA Model Seakon

Ryan Love1
1Department of Earth Sciences, University of Ottawa

Correspondence: Ryan Love (rlove@uottawa.ca)

1 abstract

  • For 1D and S40RTS configurations 2 iterations are sufficent for RSL convergence to within %1 for most regions
  • Ring discretization of the compute grid does not significantly impact RSL

2 Introduction & Purpose

The purpose of this document is to quantify the impact of some technical choices upon the scientific results of the ED3D/Seakon glacial isostatis adjustment (GIA) and relative sea level (RSL) model. I focus on the impact of the grid on which the equations from Latychev et al. (2005) are solved, the impact of the distribution of nodes into ‘rings’.

  • Seek to address:
    • Is current basegrid sufficient resolution, too high, or too coarse? Current thinking is that the lithosphere acts as a low-pass filter and so sub-degree resolution is possible overkill. Reducing resolution can also allow for better use of compute resources. Similarly, inputs to the model are generally an order of magnitude (or more) coarse than the numerical grid (excepting surface loading and topography/bathymetry). S40RTS itself specifies a nominal resolution of 1000km as such this should be a lower-limit on the resolution.
    • What impact do choices like ring-discretiziation have on final model outputs? This should have minimal impact (ideally only on numerical noise levels).
    • How many iterations of the Seakon/Ed3D solver are required for ’reasonable’ convergence of solver outputs? Standard practice is 3 as per the 1D spherically symmetric model but is this too many? Does this change with the introduction of horizontally varying viscosity?
  • Future work/unresolved questions:
    • What is the impact of vertical resolution/ Right now the grid-generator is a bit inflexible on layer count but it can be worked with to scale up or scale down. The new grid was generated while getting as close as possible to the original layering within the limits of the grid generator (and user inexperience).
    • Impact of specific layers at interfaces? Certain layers are at specific interfaces in the mantle (e.g. 660km), how important are these interfaces? Note: the reduced resolution grid currently does not respect some of these interfaces. As such, modifications to this and an initial quantification of the impact are high priority.

3 Experimental Design & Methods

Broad-scope model configuration:

  • 3 grids examined so far
    • Basegrid - stock Ed3D/Seakon grid as per K.L.
    • 62L70k - reduced resolution grid. Attemped decrease in nodes by 30%
    • 62L70k100kZL18kZU - reduced resolution grid with modified vertical layering. 100km layer spacing in the lower mantle, 18km spacing in upper mantle.
    • 62L70k70kZL30kZU - reduced resolution grid with modified vertical layering. 70km layer spacing in the lower mantle, 30km spacing in upper mantle.
  • Both grids use the same 3-layer viscosity model
    • UMV = 3×1021Pa s
    • LMV = 90×1021Pa s
    • 120km lithosphere
  • For lateral viscosity structure used S40RTS (Ritsema et al.2010).
  • Loading history is ICE5G (Peltier2004).
  • Unless otherwise noted all runs presented were executed on Archi using compute-0-22/23/24
  • Unless otherwise noted all runs presented were conducted using 48 cores

Table of runs


Table 1: Tableof simulations




n-Rings Grid

Viscosity Structure

n-Iterations




6 basegird

1D-SS

4
6 62L70k

1D-SS

4
6 62L70k100kZL18kZU

1D-SS

2
6 62L70k60kZL30kZU

1D-SS

2
6 62L70k

S40RTS Modulated

3
6 basegird

S40RTS Modulated

2
9 62L70k

1D-SS

3





3.1 Bulk Grid Statistics





Metric 62L70 basegrid



Number of Nodes 5624224 16965334
Number of edges 38644262 112468082
Internodal distance min 9000 9400
Internodal distance max 134530.172 109470.156
Number of elements 31930680 92883260
Size of Work Array 127722720 371533040
Number of surface nodes (top) 529002 2570492
Number of surface nodes (bottom) 31362 49002
Average QF 0.830375612 0.871173441
Maximum QF 0.999297976 0.999945223
Minimum QF 0.44889617 0.615376711
Low QF elements 0 0
Max Element Volume 5.11E+13 4.62E+13
Minimum Element Volume 1.22E+12 2.63E+11
Average number of elements around node 22 21
Maximum number of elements around node 38 40
Minimum number of elements around node 5 5




3.2 Vertical Layering & Horizontal resolution


Table 2: Standard resolution grid nominally referred to as basegrid. MRD represents maximal radial deviation. The radius of the Earth utilized in Seakon is 6371000m. The layer at the base of the upper mantle is denoted with a dagger ().










Layer

Nodes

Radius (m)

MRD (m)

Lateral Resolution (m)

Layer

Nodes

Radius (m)

MRD (m)

Lateral Resolution (m)

1

49002

3480000

3.37E-09

59883.4844

34

106092

5125984

1.32E-05

59947.8164

2

50412

3529844.75

9.97E-07

59885.7539

35

108162

5173801

1.13E-05

59925.2461

3

51842

3579689.75

8.24E-07

59887.957

36

110252

5221617.5

9.74E-06

59903.1055

4

53292

3630000

1.30E-06

59897.7812

37

112362

5269434.5

4.89E-06

59881.3828

5

54762

3680600.75

2.69E-06

59912.0664

38

114492

5317251

7.02E-06

59860.0664

6

56252

3731678.5

9.93E-07

59933.6328

39

116642

5365068

4.03E-06

59839.1445

7

57762

3782756.25

8.80E-07

59954.6289

40

118812

5412884.5

9.05E-06

59818.6055

8

59292

3833834

3.07E-06

59975.082

41

121002

5460701

4.39E-06

59798.4375

9

60842

3884911.75

3.25E-06

59995.0078

42

123212

5508518

5.43E-06

59778.6367

10

62412

3935989.5

8.14E-07

60014.4297

43

125442

5553953

4.27E-06

59733.5703

11

64002

3987067.25

2.33E-06

60033.3633

44

127692

5600000

6.86E-07

59695.8281

12

65612

4038145

4.91E-06

60051.832

45

129962

5666000

1.81E-06

59869.582

13

67242

4089222.5

2.54E-06

60069.8477

46

519842

5701000

2.40E-09

30119.9688

14

68892

4140300.25

3.83E-06

60087.4297

47

129962

5735000

3.07E-09

60598.668

15

70562

4191378

7.29E-06

60104.5938

48

132252

5771000

6.55E-06

60448.8203

16

72252

4242456

3.44E-06

60121.3516

49

134562

5820758

1.18E-05

60444.418

17

73962

4293533.5

3.51E-06

60137.7188

50

136892

5871700

1.46E-06

60452.2891

18

75692

4344611.5

8.11E-06

60153.7109

51

139242

5922000

1.10E-05

60453.4688

19

77442

4395689

7.28E-07

60169.3398

52

141612

5971000

4.53E-06

60441.4727

20

79212

4446767

4.19E-06

60184.6172

53

144002

6021872.5

1.02E-05

60448.4688

21

81002

4497844.5

6.17E-07

60199.5547

54

146412

6070243

1.31E-05

60430.4453

22

82812

4548922.5

4.51E-06

60214.1602

55

148842

6118614

1.41E-05

60412.7148

23

84642

4600000

8.67E-07

60228.4531

56

595362

6151000

1.27E-08

30366.4766

24

86492

4647816.5

2.89E-06

60200.1992

57

600252

6171000

2.16E-05

30340.8652

25

88362

4695633.5

5.72E-06

60172.5469

58

605162

6191000

3.73E-05

30315.4629

26

90252

4743450

5.34E-06

60145.4727

59

610092

6211000

3.01E-05

30290.2676

27

92162

4791267

7.56E-06

60118.9648

60

615042

6231000

5.02E-05

30265.2734

28

94092

4839083.5

1.01E-05

60093.0039

61

620012

6251000

5.94E-05

30240.4805

29

96042

4886900.5

4.08E-06

60067.5703

62

625002

6275000

5.56E-05

30235.1602

30

98012

4934717

6.88E-06

60042.6523

63

630012

6298500

5.04E-05

30227.4824

31

100002

4982534

8.04E-06

60018.2344

64

635042

6322000

7.23E-05

30219.8652

32

102012

5030350.5

1.07E-05

59994.2969

65

640092

6346600

5.36E-05

30217.5449

33

104042

5078167.5

1.37E-05

59970.8281

66

2560362

6361600

5.53E-09

15144.5088

67

2570492

6371000

1.59E-04

15136.9717













Table 3: Reduced resolution grid nominally referred to as 62L70 (62-layers, 70km). MRD represents maximal radial deviation. The radius of the Earth utilized in Seakon is 6371000m. The layer at the base of the upper mantle is denoted with a dagger ().










Layer

Nodes

Radius (m)

MRD (m)

Lateral Resolution (m)

Layer

Nodes

Radius (m)

MRD (m)

Lateral Resolution (m)











1

31362

3480000

1.96E-09

74853.0703

32

75692

4986390

2.78E-06

69039.5234

2

32492

3528593.25

1.21E-06

74566.8594

33

77442

5034983

4.10E-06

68920.1719

3

33642

3577186.5

8.59E-07

74290.5234

34

79212

5083576.5

5.95E-06

68803.5

4

34812

3625779.75

2.27E-06

74023.5469

35

81002

5132169.5

5.45E-06

68689.4219

5

36002

3674373

3.07E-06

73765.4688

36

82812

5180762.5

7.41E-06

68577.8516

6

37212

3722966

1.78E-06

73515.8516

37

84642

5229356

8.79E-06

68468.7031

7

38442

3771559.25

1.32E-06

73274.2812

38

86492

5277949

1.74E-06

68361.9062

8

39692

3820152.5

7.06E-07

73040.375

39

88362

5326542.5

4.01E-06

68257.375

9

40962

3868745.75

9.00E-07

72813.7812

40

90252

5375135.5

5.87E-06

68155.0469

10

42252

3917339

8.48E-09

72594.1562

41

92162

5423729

8.14E-06

68054.8516

11

43562

3965932.25

1.48E-07

72381.1875

42

94092

5472322

8.51E-06

67956.7266

12

44892

4014525.5

9.74E-07

72174.5703

43

96042

5520915.5

1.08E-05

67860.5938

13

46242

4063118.75

3.43E-06

71974.0312

44

98012

5569508.5

1.28E-05

67766.4141

14

47612

4111711.75

2.96E-06

71779.3047

45

100002

5618101.5

1.38E-05

67674.1094

15

49002

4160305

1.40E-06

71590.1328

46

102012

5666695

1.45E-05

67583.6328

16

50412

4208898.5

4.92E-07

71406.2969

47

104042

5701000

9.61E-06

67326.1953

17

51842

4257491.5

2.11E-06

71227.5625

48

106092

5763881.5

1.18E-05

67407.9531

18

53292

4306084.5

2.91E-06

71053.7266

49

108162

5812474.5

1.19E-05

67322.6484

19

54762

4354678

4.56E-06

70884.5859

50

110252

5861068

1.03E-05

67238.9609

20

56252

4403271

2.35E-06

70719.9531

51

112362

5909661

8.96E-06

67156.8594

21

57762

4451864.5

4.26E-06

70559.6484

52

114492

5958254

7.30E-06

67076.2891

22

59292

4500457.5

5.45E-06

70403.5156

53

116642

6006847.5

1.15E-05

66997.2188

23

60842

4549051

7.22E-06

70251.3828

54

118812

6055440.5

9.65E-06

66919.5859

24

62412

4597644

6.26E-06

70103.0938

55

121002

6104034

8.07E-06

66843.375

25

64002

4646237.5

5.65E-06

69958.5156

56

123212

6152627

6.20E-06

66768.5312

26

65612

4694830.5

4.10E-06

69817.5078

57

125442

6201220.5

5.57E-06

66695.0234

27

67242

4743423.5

5.61E-06

69679.9453

58

127692

6249813.5

3.65E-06

66622.8203

28

68892

4792017

3.91E-06

69545.6875

59

129962

6298407

1.95E-06

66551.8828

29

70562

4840610

2.74E-06

69414.625

60

132252

6347000

4.31E-09

66482.1797

30

72252

4889203.5

9.68E-07

69286.6484

61

529002

6362000

2.46E-09

33319.9375

31

73962

4937796.5

9.70E-07

69161.6484

62

529002

6371000

2.45E-09

33367.0742













PIC

Figure 1: Map showing the locations where the USEC is sampled for comparisons between iterations and grid resolutions. The points are equidistantly spaced at 100km intervals between (-65E,47.7N) and (-84.3E,30.4N). This track is roughly parallel with the East Coast of North America between Laméque Island, New Bruinswick, Canada, and Tallahassee, Florida, USA. Points are sampled using GMT5 with the argument ’-E-65/47.7/-84.3/30.4+i100k’ passed to the grdtrack utility.



PIC

Figure 2: Outlined in red is the region of the interior of North America over which meridional averaging is performed in some of the benchmarks. The region was selected as it encompases some of the largest misfits near the southern margin of the North American Ice Complex as well as much of the interior of the ice sheets themselves. The domain encompases (-130E,40N) to (-50E,60N).


4 Results and Discussion


PIC

Figure 3: USEC RSL through time for the basegrid run at the 4th iteration. Due to the scale and the way the RSL field converges this figure is indistinguishable from the 2nd and 3rd iterations (shown in Fig. 4).


4.1 Impact of Seakon Run Iterations


PIC PIC PIC

Figure 4: RSL differences shown are between successive iterations of the Seakon code along the USEC. Differences are for basegrid grid configuration in a 1D spherically symmetric Earth structure configuration. Each iteration results in an order of magnitude reduction in differences between the previous iteration.



PIC PIC PIC

Figure 5: Top two panels show RSL differences between successive iterations of the Seakon code for the interior of North America. Bottom panel shows the difference between the 3rd and 2nd iteration as a percentage of the 4th iteration’s RSL value. The RSL field for the 4th iteration are contoured in black. RSL differences are calculated using the following CDO commands: cdo -L -sqrt -mermean -sqr -sellonlatbox,230,310,40,60 Differences are for basegrid grid configuration in a 1D spherically symmetric Earth structure configuration. As for USEC RSL, each iteration results in an order of magnitude reduction in differences between the previous iteration.


  • Run times between iterations are very similar, for a basegrid run with 6-rings on a super node run times are between 35.5 and 37 hours.
  • Run times may decrease with iteration count, for basegrid run with 6-rings on a super node run times (in order of shortest->longest) were as follows: Iteration 3,2,4,1. However, without many repeated trials this is within the noise-floor of runtimes on a super-node.
  • Point of diminishing returns on iteration for USEC RSL is 2-iterations as difference between iteration 2 and iteration 3 is O (0.1m) (see Fig. 4) and an order of magnitude smaller between iteration 3 and iteration 4. By comparison the RSL signal itself is O(100m) (see Fig. 3).
  • Fig. 5 demonstrates that interior of ice sheet RSL field shows similar results to USEC RSL results.

4.2 Impact of Ring number and Ring Layout

  • Ring layout is decided in sftdd code with a suggested layout provided in the interactive version of the software
  • Suggested ring layout does not always make best use of resources and generally requires some degree of hand-tuning
  • Suggested ring number has thusfar been reflective of fastest configuration
  • SFTDD may not be able to create valid configurations for every ring number (too high and rings may not get assigned nodes and will fail)
  • Figs. 6 and 7 show RSL differences which arise due to different ring configurations are O(1cm) or less.
  • Ring configuration can be considered as a run-time impact, not a numerical one given scale of variability between different ring configurations


PIC
Figure 6: USEC transect RSL anomaly using the 62L70k grid for 6-Ring vs. 9-Ring decomposition of the grid.



PIC

Figure 7: RSL anomaly at LGM using the 62L70k grid for 6-Ring vs. 9-Ring decomposition of the grid.


4.3 Impact of Grid Resolution on RSL field


PIC

Figure 8: 62L70k - Basegrid RSL values for the 3rd iteration of each run along the USEC transect.



PIC PIC

Figure 9: Top panel shows the difference in RSL at LGM for 62L70 - Basegrid, RSL values are contoured at 50m increments. Bottom panel shows RSL at LGM for the Basegrid configuration.



PIC

Figure 10: RSL for 62L70k-basegrid results expressed as a percent of the basegrid RSL value. Contours reflect the RSL values from the basegrid simulation. RSL values sampled along the USEC transect. Runs shown use 1D spherically symmetric configuration.



PIC

Figure 11: Difference in the RSL field between the 62L70-Basegrid configurations, as in Figure 9, zoomed in on the North American Ice Complex and the Eurasian icesheet. RSL values from the Basegrid configuration are contoured in black at 50m intervals. The average rate of change in ice thickness between 21kaBP and 25kaBP, expressed in m/ka, are contoured in red at 10m/ka intervals.


4.4 Impact of Grid Resolution on Walltime and Memory

  • Default basegrid configuration utilizes 250GB
  • On super nodes the best configurations obtained result in a total wall time of 36h per iteration

4.5 Impact of Vertical Resolution on RSL field

Concepts:

  • In order to determine the ideal layering need to find the vertical length scales of changes in the mantle
  • Length scales set by fastest varying (in the vertical dimension) component
  • PREM varies slowly outside discontinuities (already baked into model and tomography)
  • Current input to Seakon uses 55 layers for S40RTS input, S40RTS has 24 layers (potentially only 21 being used)
  • Need some sort of mathematical means of determining the ’length scale’ of changes or the minimum resolution to represent a specific amount of variance
  • Possible avenues of exploration: Fourier analysis (also Aliasing and sampling theory) and EOFs of the S40RTS(z)
  • UM is 25% the thickness of the LM (670km ignoring lithosphere; 500km assuming 120km lithosphere; LM is 2200km). Therefore it seems reasonable that decreasing resolution in LM can greatly enhance UM resolution at little cost.

Results:

  • Changing grid spacing to not be as equidistant as default leads to decreased quality factor
  • 62L70k grid modified with 100km spacing in lower mantle and 18km in upper mantle reduced quality factor average from 0.83 to 0.61
  • Fig. 12 show that modifying the vertical distribution of nodes such that the upper mantle has higher vertical resolution than the lower mantle results in greater misfits to Basegrid relative to 62L70k. As such changing from relatively equidistant spacing in the mantle to higher resolution UM at the expense of LM is not worth it for at least spherically symmetric or equidistant tomography data. If dealing with ensemble of incorporating high-resolution UM tomography may be worth revisiting in the future.


PIC PIC PIC
Figure 12: RMS of RSL anomalies relative to basegrid for 62L70k, 62L70k60kZL30kZU, and 62L70k60kZL30kZU respectively. Values are calculated using meriodonal means for formerly glaciated regions of North America (specifically, 230 E ,40 N and 310E,60N).


5 Conclusions

  • Seakon model RSL field converges to O(10cm) after only two solver iterations. This result holds for both basegrid and 62L70k grid configurations. See Fig. 4.
  • Choice of ring decomposition does not significantly affect RSL results, RSL anomalies between the optimal ring decomposition (6-ring) and an inefficent ring (9-ring) are O(1cm) or less. See Figs. 6 and 7.
  • A reduction in lateral resolution of 33% (i.e. the difference between basegrid and 62L70k) does not significantly alter the resulting RSL field for most regions. Differences are generally 2% except in regions where RSL is very small (close to zero). See Fig. 10.
  • The largest differences between the grids correlates (visual inspection) more strongly with gradients of ice sheet height. See Fig. 11. Examination of the gradient field and the dz
dt field are to be conducted next.

References

   Latychev, K., Mitrovica, J. X., Tromp, J., Tamisiea, M. E., Komatitsch, D., and Christara, C. C.: Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation, Geophysical Journal International, 161, 421–444, https://doi.org/10.1111/j.1365-246x.2005.02536.x, 2005.

   Peltier, W.: GLOBAL GLACIAL ISOSTASY AND THE SURFACE OF THE ICE-AGE EARTH: The ICE-5G (VM2) Model and GRACE, Annual Review of Earth and Planetary Sciences, 32, 111–149, https://doi.org/10.1146/annurev.earth.32.082503.144359, 2004.

   Ritsema, J., Deuss, A., van Heijst, H. J., and Woodhouse, J. H.: S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, Geophysical Journal International, 184, 1223–1236, https://doi.org/10.1111/j.1365-246x.2010.04884.x, 2010.

Appendix A: M-Ring Program This program was derived from the sftdd program and modified to run as a standalone utility for examining the impact of requested processor count on the construction of the mring array. Note that this is not designed as an example of GOOD Fortran programming (it is not) but mearly to provide the maths behind the default ring layout.

PROGRAM mRingTesting 
REAL*8 :: pi, dt12, dt, dm 
INTEGER, DIMENSION(1000) :: mring 
INTEGER :: np, nring, j, k 
REAL*8, DIMENSION(1001) :: tring 
REAL*8, DIMENSION(100 ) :: sring 
WRITE(6,*) time() 
READ(5,*) np 
pi=2.0d0*dasin(1.0d0) 
dt12=dacos(1.0d0-2.0d0/dble(np)) 
dt=dt12+dt12 
nring=nint((pi-dt)/dt)+2 
 
mring(1)=1 
mring(nring)=1 
 
tring(1)=0.0d0 
tring(2)=dt12 
tring(nring)=pi-dt12 
tring(nring+1)=pi 
 
do j=3,nring-1 
        tring(j)=dt12+dble(j-2)*dt 
end do 
do j=1,nring 
        sring(j)=dcos(tring(j+1))-dcos(tring(j)) 
end do 
dm=sring(1) 
do j=2,nring-1 
        k=nring+1-j 
        if(k < j) exit 
        mring(j)=nint(sring(j)/dm) 
        if(k > j) mring(k)=mring(j) 
end do 
DO j=1,nring,1 
  WRITE(6,*) j, mring(j) 
ENDDO 
WRITE(6,*) time() 
END