```{tags} category:benchmark feature:3d feature:spherical feature:analytical-solution ``` (sec:benchmarks:thin-shell-gravity)= # 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}=6371-D-5~\text{ km}$ and the outer radius $R_\text{outer}=6371-D+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}^3-R_\text{inner}^3)\rho_0$. Both the Stokes and energy equations solve are bypassed: ```{literalinclude} thinshell_b.prm ``` 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: ```{literalinclude} thinshell_a.prm ``` 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 longitude-latitude grid spanning the whole Earth, with a $2^\circ\times 2^\circ$ resolution: ```{literalinclude} thinshell_c.prm ``` The Gauss-Legendre 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 ```{math} U({\mathbf r}) = \iiint_\Omega \frac{G \rho({\mathbf r}')}{|{\mathbf r}-{\mathbf r}'|} d{\mathbf r}' = \sum_K \iiint_{\Omega_K} \frac{G \rho({\mathbf r}')}{|{\mathbf r}-{\mathbf r}'|} d{\mathbf r}', ``` 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 ```{math} {\mathbf g}({\mathbf r}) = \iiint_\Omega \frac{G \rho({\mathbf r}')}{|{\mathbf r}-{\mathbf r}'|^2} d{\mathbf r}' = \sum_K \iiint_{\Omega_K} \frac{G \rho({\mathbf r}')}{|{\mathbf r}-{\mathbf r}'|^2} d{\mathbf r}'. ``` 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 $(2n-1)$-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 {cite:t}`turcotte: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. ```{table} Thin shell gravity benchmark: $1\text{ mGal}=10^{-5}\text{ m}/\text{s}^2$. $n_q$ is the number of GLQ points per cell. 'a.v.' stands for analytical value. :name: tab:thin_shell_gravity_benchmark | | | | | | | | | | |:----:|:----:|:-----:|:---------:|:------------:|:---------:|:------------:|:----------------------------------:|:------------:| | 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. air-water, water-crust, 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.