Thin shell gravity benchmark#
This section was contributed by Cedric Thieulot, Bart Root and Paul Bremner.
The idea behind this benchmark is to test the accuracy of the gravity postprocessor: the domain is a thin shell of constant density somewhere in the Earth mantle and we wish to compute the resulting gravity field and potential at satellite height.
The domain is a spherical shell of 10 km radius centered at a depth \(D\), i.e. the inner radius is \(R_{inner}=6371D5~\text{ km}\) and the outer radius \(R_\text{outer}=6371D+5~ \text{ km}\). It is filled with a fluid of constant density \(\rho_0=3300~\text{ kg}/\text{m}^3\) with total mass \(M=\frac43 \pi (R_\text{outer}^3R_\text{inner}^3)\rho_0\).
Both the Stokes and energy equations solve are bypassed:
set Nonlinear solver scheme = no Advection, no Stokes
We make use of the Custom mesh subdivision option to generate a mesh where a single cell is used in the radial direction (parameterized with the number of ‘slices’) while 6 blocks (default value) of \(2^5\times 2^5\) cells are used in the lateral direction. This gives a total of 6,144 mesh cells. For \(D=100~\text{ km}\) the parameterization of the mesh is then as follows:
subsection Geometry model
set Model name = spherical shell
subsection Spherical shell
set Inner radius = 6266e3
set Outer radius = 6276e3
set Custom mesh subdivision = number of slices
set Number of slices = 1
set Initial lateral refinement = 5
end
end
We wish to compute the gravity potential \(U\) and the gravity vector \({\mathbf g}\) at a radius of \(R_{s}=6371+250=6621~\text{ km}\) (i.e. the gravity experienced by a satellite flying at a height of \(250~\text{ km}\) above the surface of the Earth) on a regular longitudelatitude grid spanning the whole Earth, with a \(2^\circ\times 2^\circ\) resolution:
subsection Postprocess
set List of postprocessors = gravity calculation,visualization
subsection Gravity calculation
set Sampling scheme = map
set Minimum radius = 6621e3
set Minimum longitude = 180
set Maximum longitude = 179
set Number points longitude = 180
set Minimum latitude = 90
set Maximum latitude = +90
set Number points latitude = 90
set Quadrature degree increase = 1
end
end
The GaussLegendre Quadrature (GLQ) algorithm is central to the code as it is used to compute the elemental integrals stemming from the discretization of the weak form of the PDEs. Logically, the gravity postprocessor is also based on the GLQ since the potential \(U\) at any location in space is given by the integral
where the sum runs over all cells of the mesh and \({G}=6.67430\times 10^{11}~\text{ m}^3/\text{kg}/\text{s}^2\) is the gravitational constant. The gravity acceleration vector is obtained via \({\mathbf g}={\mathbf \nabla} U\), or
The default number of quadrature points for each cell in the postprocessor is \(2^\text{ndim}\), where ndim is the number of dimensions. The \(n\)point GLQ allows exact integration of \((2n1)\)order polynomials. However, the integrand of the Newton integral is not a polynomial (it contains a term \(\sim r^{m}\)), so there is not an optimal number of quadrature points to use. Therefore, the postprocessor allows the user to choose how many additional points per dimension are used with an expectation that an increase in the number of quadrature points inside the cells leads to a more accurate calculation. This increase number \(I\) is set to 1 in the example above, although it can also be chosen to be 0, 1, 2 or even 1. The case \(I=1\) approximates the cell as a point mass since there is a single quadrature point in the middle of the cell.
Given the symmetry of the problem, the values of the potential and
acceleration depend solely on the distance \(r\) from the origin, see Turcotte and Schubert [2014]. Their analytical values are given
by \(U(r) =  GM/r\) and \({\mathbf g}(r)= g_r(r)=GM/r^2\). Minimum, maximum,
and average values of both the potential and the acceleration are printed in
the statistics file while measurements at all the desired points are written
in the output_gravity
folder inside the output folder. We ran the input file
for \(I\in \{1,0,1,2\}\) for \(D=0,100,500,1500,3000~\text{ km}\) and the results
are presented in Table [1] alongside the analytical values.
D 
\(I\) 
\(n_q\) 
\(g_r\) (mGal) 
\(U\) (J kg^{−1}) 

(km) 
avrg. 
min 
max 
avrg. 
min 
max 

0 
1 
\(1^3\) 
2563.6541 
2530.6859 
2607.5210 
169764.4978 
169832.3647 
169744.3149 
0 
\(2^3\) 
2562.8680 
2553.5683 
2571.9157 
169676.4060 
169681.1065 
169671.6199 

1 
\(3^3\) 
2562.6878 
2561.9847 
2563.7170 
169676.3142 
169676.8060 
169675.9148 

2 
\(4^3\) 
2562.6993 
2562.6215 
2562.7395 
169676.3211 
169676.3361 
169676.2903 

a.v. 
2562.6993 
169676.3210 

100 
1 
\(1^3\) 
2484.0821 
2479.3619 
2499.2301 
164477.2574 
164520.0786 
164472.4899 
0 
\(2^3\) 
2482.9027 
2481.7033 
2484.0834 
164391.6151 
164392.2246 
164390.9963 

1 
\(3^3\) 
2482.8799 
2482.7731 
2482.9959 
164391.6031 
164391.6623 
164391.5473 

2 
\(4^3\) 
2482.8819 
2482.8759 
2482.8865 
164391.6041 
164391.6067 
164391.6012 

a.v. 
2482.8818 
164391.6041 

500 
1 
\(1^3\) 
2177.3562 
2177.2415 
2179.8737 
144163.9859 
144178.2736 
144162.2361 
0 
\(2^3\) 
2176.2392 
2176.2253 
2176.2404 
144088.7939 
144088.7971 
144088.7538 

1 
\(3^3\) 
2176.2391 
2176.2391 
2176.2392 
144088.7937 
144088.7938 
144088.7936 

2 
\(4^3\) 
2176.2391 
2176.2391 
2176.2391 
144088.7937 
144088.7937 
144088.7937 

a.v. 
2176.2391 
144088.7937 

1500 
1 
\(1^3\) 
1498.7997 
1498.7569 
1499.0522 
99236.0403 
99238.4122 
99235.3167 
0 
\(2^3\) 
1498.0239 
1498.0236 
1498.0240 
99184.1673 
99184.1678 
99184.1647 

1 
\(3^3\) 
1498.0240 
1498.0240 
1498.0240 
99184.1672 
99184.1672 
99184.1672 

2 
\(4^3\) 
1498.0240 
1498.0240 
1498.0240 
99184.1672 
99184.1672 
99184.1672 

a.v. 
1498.0240 
99184.1672 

3000 
1 
\(1^3\) 
717.8387 
717.8315 
717.8533 
47528.1777 
47528.3488 
47528.0725 
0 
\(2^3\) 
717.4641 
717.4640 
717.4641 
47503.2981 
47503.2982 
47503.2980 

1 
\(3^3\) 
717.4641 
717.4641 
717.4641 
47503.2981 
47503.2981 
47503.2981 

2 
\(4^3\) 
717.4641 
717.4641 
717.4641 
47503.2981 
47503.2981 
47503.2981 

a.v. 
717.4641 
47503.2981 
The accuracy of the calculations increases with \(I\) but so does the time spent in the postprocessor: for \(I\in\{1,0,1,2\}\) this time was about 18, 132, 440 and 1040 s respectively. This is easily explained when one realizes that increasing \(I\) from 0 to 1 means that the number of GLQ points per cell goes from \(2^3=8\) to \(3^3=27\), i.e. a \(3.375\) increase in operations for the same number of cells and measurement points. The time spent in the postprocessor increases by a similar factor (\(440/132\sim 3.4\)). The accuracy obtained with lower \(I\) values increases with increasing anomaly depth (or increasing distance from the observation point). Note that this benchmark has uniform density so that the projection of the density from the nodes onto the quadrature points is exact and it does not introduce any smoothing of data which might occur in the presence of density discontinuities (e.g. airwater, watercrust, moho, etc.) inside a cell. \(I=1\) seems to be the best compromise between accuracy (gravity acceleration errors are less than 0.01 mGal for all shells) and compute time for this experiment.