In the aerosol condensation process, MAM needs to perform an evaluation of integral with the following formula:
\begin{equation} \label{goveq} \int_{-\infty}^{\infty} f(x) e^{-x^2} dx \,, \end{equation}where f(x) is a function of x. The value of Eq. \eqref{goveq} can be approximated by the Gauss-Hermite quadrature rule \cite{GH-WIKI}, which takes the following formula:
\begin{equation} \label{ghqeq} \int_{-\infty}^{\infty} f(x) e^{-x^2} dx \approx \sum_{k=1}^{K} \mathbb{W}_{_{k}} f(\mathbb{R}_{_{k}}) \,, \end{equation}where $\mathbb{R}_{_{k}}$ is the k-th root of Hermite polynomials (k = 1,2,...,K) and the $\mathbb{W}_{_{k}}$ is the associated weight. Using K quadrature points, Gauss-Hermite quadrature rule can evaluate the Eq. \eqref{goveq} exactly if f(x) is a polynomial of degree up to 2K-1 \cite{steen-1969-mc}.
In this test, we first aim at examining the correctness of Gauss-Hermite quadrature rule implemented in MAM (e.g., no bug in the index of array or variable name).
This verification test uses a box model that inherits the data structure and condensation subroutines from MAM in E3SM. The box model serves as a driver program to specify the necessary parameters in order to calculate the MAM and reference solutions.
According to Eq. \eqref{goveq}, it is trival to find out that:
\begin{itemize} \item If f(x) is an odd function (i.e., f(x) = -f(-x)), \begin{align} \int_{-\infty}^{\infty} f(x) e^{-x^2} dx &= 0 \,. \\ \end{align} \item If f(x) is an even function (i.e., f(x) = f(-x)), \begin{equation} \int_{-\infty}^{\infty} f(x) e^{-x^2} dx = 2 \int_{0}^{\infty} f(x) e^{-x^2} dx \,. \end{equation} \end{itemize}Moreover, the integral of $\int_{0}^{\infty} f(x) e^{-x^2} dx$ can be evaluated analytically if f(x) has a polynomial formula such as \cite{GH-Ingr}:
\begin{align} \int_{0}^{\infty} x^{2N} \exp^{-x^2} dx &= \sqrt{\pi} \frac{(2N-1)!!}{2^{N+1}} \,, \\ \int_{0}^{\infty} x^{2N+1} \exp^{-x^2} dx &= \frac{N!}{2} \,, \end{align}where N! is the factorial of a positive integer N (i.e., $N! = N \times (N-1) \times ... \times 2 \times 1$); (2N-1)!! is the double factorial that calculates the product of all the odd positive integers up to 2N-1 (i.e., $(2N-1)!! = (2N-1) \times (2N-3) \times ... \times 2 \times 1$). By definition, 0! = 1, 0!! = (-1)!! = 1 \cite{arfken-2005}.
Therefore, in this verification test, a polynomial $f(x) = 4x^3 + 3x^2 + 2x + 1$ is used to evaluate the integral in Eq. \eqref{goveq}.
We use the original MAM subroutines to evaluate the Eq. \eqref{goveq} by setting f(x) to the polynomial formula in Section 3.1. MAM originally uses 2 Gauss-Hermite quadrature points with 8 decimal digits but in this test, we would use more accurate Gauss-Hermite points, which are listed below:
\begin{align} x_{_{1}} &= 0.7071067811865475, &\mathbb{W}_{_{1}} &= 0.8862269254527578; \nonumber \\ x_{_{2}} &= -0.7071067811865475, &\mathbb{W}_{_{2}} &= 0.8862269254527578. \nonumber \end{align}This solution is denoted as $C_{_{MAM}}$.
Given the formula of f(x) in Section 3.1, the exact value of Eq. \eqref{goveq} is calculated by:
\begin{align} \int_{-\infty}^{\infty} f(x) e^{-x^2} dx = 2 \left( 3 \int_0^{\infty} x^2 e^{-x^2} dx + \int_0^{\infty} e^{-x^2} dx \right) = \frac{3}{2}\sqrt{\pi} + \sqrt{\pi} = \frac{5}{2}\sqrt{\pi} \label{refsol} \,. \end{align}This exact value then serves as the reference solution and is denoted as $C_{_{REF}}$.
In this verification test, the following quantities are specified:
\begin{enumerate} \item The formula of f(x), which is shown in Section 3.1. \item Number of quadrature points to be used for the test solution, which is 2 throughout this verification test. \item Gauss-Hermite quadrature roots and their associated weights, which are already listed in Section 3.2. \end{enumerate}The relative error between $C_{_{MAM}}$ and $C_{_{REF}}$ (denoted as $\phi_{_{MAM}}$) is first calculated by:
\begin{equation} \phi_{_{MAM}} = \frac{ \left| C_{_{MAM}} - C_{_{REF}} \right| }{C_{_{REF}}} \,. \end{equation}A tolerance level $\epsilon$ is further defined as 10 times the machine epsilon returned by the Fortran intrinsic function $\sf{epsilon}$. This machine epsilon depends on the working precision of the code, double or single (for this verification test, we choose the double precision that is used in MAM).
Because 2 Gauss-Hermite quadrature points are used in this verification test, Eq. \eqref{goveq} can be evaluated exactly for a polynomail f(x) of degree up to 3. Thus an overall PASS is issued if $\phi_{_{MAM}} \lt \epsilon$.
The verification test requires a test driver to perform the verification test and the configuration files for compilation and execution of cmdv-test-runner (e.g., the yaml file). If we go the root directory of CMDV-Verification folder that is cloned from GitHub Repo, the test driver code (driver.F90) is under $\rm \color{blue}{tests/mam/restructed\_tests/F90\_SRC/test\_drivers/condensation/gauss\_hermite\_accuracy}$ and the yaml file is under $\rm \color{blue}{tests/mam/restructed\_tests/yaml\_files\_for\_test\_runner/condensation/gauss\_hermite\_accuracy.test.yaml}$.
If the verification test runs successfully, we could obtain the netcdf file named "condensation_verification_test.nc" under the "postprocess_input" folder with the MAM and reference solutions. We could calculate the $\phi$ according to the Sectoin 3.5 and compare with $\epsilon$ (also stored in the condensation_verification_test.nc) to see whether an overall PASS/FAIL should be issued.
#load library and modules
import numpy as np
import xarray as xr
#this cell is used to read in the input_data
#define the file path and name
nc_file = './postprocess_input/condensation_verification_test.nc'
#read in the MAM and reference solutions, tolerance, number of polynomials and polynomial coefficients
ds = xr.open_dataset(nc_file)
mam = np.array(ds['mam_sol'])
ref = np.array(ds['ref_sol'])
tol = ds['tol'].values
#check the relative error
mam_err = np.abs(mam - ref) / ref
print("relative error of mam solution: %.6e" % mam_err)
#run this cell to issue the PASS/FAIL
if mam_err < tol:
print("PASS")
else:
print("FAIL")
[1] , ``Gauss–Hermite quadrature'', . online
[2] Steen NM, Byrne GD and Gelbard EM, ``Gaussian Quadratures for the Integrals∫∞ 0 exp (-x 2) f (x) dx and∫ b 0 exp (-x 2) f (x) dx'', Mathematics of Computation, vol. , number , pp. 661--671, 1969.
[3] , ``Gauss Integral'', . online
[4] George B Arfken and Hans J Weber, ``Mathematical methods for physicists, 6th ed.'', 2005.