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}.
Currently MAM uses less accurate Gauss-Hermite quadrature points (defined as double precision but only provided 8 decimal digits). This is expected to introduce a numerical error and thus we will quantify it in this test.
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 uses 2 Gauss-Hermite quadrature points (less accurate), and the roots and associated weights of Hermite polynomials are listed below \cite{zwillinger2002crc}:
\begin{align} x_{_{1}} &= 0.70710678, &\mathbb{W}_{_{1}} &= 0.88622693; \nonumber \\ x_{_{2}} &= -0.70710678, &\mathbb{W}_{_{2}} &= 0.88622693. \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 errors 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}Because MAM uses 2 Gauss-Hermite quadrature points with 8 decimal digits only, an overall PASS is issued if $\phi_{_{MAM}} \lt 10^{-8}$.
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\_rounded\_points}$ and the yaml file is under $\rm \color{blue}{tests/mam/restructed\_tests/yaml\_files\_for\_test\_runner/condensation/gauss\_hermite\_rounded\_points.test.yaml}$.
Note that this test is not expected to run nightly and thus not included in the latest CMDV-Verification GitHub Repo. In order to run this test, please checkout the version with hash 5e0c49ae4b6866742414547a995e085ba92d60d8 first.
If the verification test runs successfully, we could obtain the netcdf file named "condensation_verification_test.nc" under the "postprocess_input" folder with MAM and reference solutions. We could calculate the $\phi$ according to the Sectoin 3.5 and compare with $10^{-8}$ 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'])
#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 < 1.e-8:
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.
[5] Daniel Zwillinger, ``CRC standard mathematical tables and formulae'', 2002.