<jupyter_start><jupyter_text>Benchmark Problem 7: MMS Allen-Cahn<jupyter_code>from IPython.display import HTML

HTML('''{% include jupyter_benchmark_table.html num="[7]" revision=0 %}''')<jupyter_output><empty_output><jupyter_text>* [Overview](Overview)* [Governing equation and manufactured solution](Governing-equation-and-manufactured-solution)* [Domain geometry, boundary conditions, initial conditions, and stopping condition](Domain-geometry,-boundary-conditions,-initial-conditions,-and-stopping-condition)* [Parameter values](Parameter-values)* [Benchmark simulation instructions](Benchmark-simulation-instructions)    * [Part (a)](Part-%28a%29)    * [Part (b)](Part-%28b%29)    * [Part (c)](Part-%28c%29)* [Results](Results)* [Feedback](Feedback)* [Appendix](Appendix)    * [Computer algebra systems](Computer-algebra-systems)    * [Source equation](Source-equation)    * [Code](Code)See the journal publication entitled ["Benchmark problems for numerical implementations of phase field models"][benchmark_paper] for more details about the benchmark problems. Furthermore, read [the extended essay][benchmarks] for a discussion about the need for benchmark problems.[benchmarks]: ../[benchmark_paper]: http://dx.doi.org/10.1016/j.commatsci.2016.09.022  Overview The Method of Manufactured Solutions (MMS) is a powerful technique for verifying the accuracy of a simulation code. In the MMS, one picks a desired solution to the problem at the outset, the "manufactured solution", and then determines the governing equation that will result in that solution. With the exact analytical form of the solution in hand, when the governing equation is solved using a particular simulation code, the deviation from the expected solution can be determined exactly. This deviation can be converted into an error metric to rigously quantify the error for a calculation. This error can be used to determine the order of accuracy of the simulation results to verify simulation codes. It can also be used to compare the computational efficiency of different codes or different approaches for a particular code at a certain level of error. Furthermore, the spatial/temporal distribution can give insight into the conditions resulting in the largest error (high gradients, changes in mesh resolution, etc.).After choosing a manufactured solution, the governing equation must be modified to force the solution to equal the manufactured solution. This is accomplished by taking the nominal equation that is to be solved (e.g.  Allen-Cahn equation, Cahn-Hilliard equation, Fick's second law, Laplace equation) and adding a source term. This source term is determined by plugging the manufactured solution into the nominal governing equation and setting the source term equal to the residual. Thus, the manufactured solution satisfies the MMS governing equation (the nominal governing equation plus the source term). A more detailed discussion of MMS can be found in [the report by <NAME> Knupp][mms_report].In this benchmark problem, the objective is to use the MMS to rigorously verify phase field simulation codes and then provide a basis of comparison for the computational performance between codes and for various settings for a single code, as discussed above. To this end, the benchmark problem was chosen as a balance between two factors: simplicity, to minimize the development effort required to solve the benchmark, and transferability to a real phase field system of physical interest. [mms_report]: http://prod.sandia.gov/techlib/access-control.cgi/2000/001444.pdf  Governing equation and manufactured solutionFor this benchmark problem, we use a simple Allen-Cahn equation as the governing equation$$\begin{equation}\frac{\partial \eta}{\partial t} = - \left[ 4 \eta \left(\eta - 1 \right) \left(\eta-\frac{1}{2} \right) - \kappa \nabla^2 \eta \right] + S(x,y,t) \end{equation}$$where $S(x,y,t)$ is the MMS source term and $\kappa$ is a constant parameter (the gradient energy coefficient). The manufactured solution, $\eta_{sol}$ is a hyperbolic tangent function, shifted to vary between 0 and 1, with the $x$ position of the middle of the interface ($\eta_{sol}=0.5$) given by the function $\alpha(x,t)$:$$\begin{equation}\eta_{sol}(x,y,t) = \frac{1}{2}\left[ 1 - \tanh\left( \frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right) \right] \end{equation}$$$$\begin{equation}\alpha(x,t) = \frac{1}{4} + A_1 t \sin\left(B_1 x \right) + A_2 \sin \left(B_2  x + C_2 t \right)\end{equation}$$where $A_1$, $B_1$, $A_2$, $B_2$, and $C_2$ are constant parameters. This manufactured solution is an equilbrium solution of the governing equation, when $S(x,y,t)=0$ and $\alpha(x,t)$ is constant. The closeness of this manufactured solution to a solution of the nominal governing equation increases the likihood that the behavior of simulation codes when solving this benchmark problem is representive of the solution of the regular Allen-Cahn equation (i.e. without the source term). The form of $\alpha(x,t)$ was chosen to yield complex behavior while still retaining a (somewhat) simple functional form. The two spatial sinusoidal terms introduce two controllable length scales to the interfacial shape. Summing them gives a  "beat" pattern with a period longer than the period of either individual term, permitting a domain size that is larger than the wavelength of the sinusoids without a repeating pattern. The temporal sinusoidal term introduces a controllable time scale to the interfacial shape in addition to the phase transformation time scale, while the linear temporal dependence of the other term ensures that the sinusoidal term can go through multiple periods without $\eta_{sol}$ repeating itself.Inserting the manufactured solution into the governing equation and solving for $S(x,y,t)$ yields:$$\begin{equation}S(x,y,t) = \frac{\text{sech}^2 \left[ \frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right]}{4 \sqrt{\kappa}} \left[-2\sqrt{\kappa} \tanh \left[\frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right] \left(\frac{\partial \alpha(x,t)}{\partial x} \right)^2+\sqrt{2} \left[ \frac{\partial \alpha(x,t)}{\partial t}-\kappa \frac{\partial^2 \alpha(x,t)}{\partial x^2} \right] \right]\end{equation}$$where $\alpha(x,t)$ is given above and where:$$\begin{equation}\frac{\partial \alpha(x,t)}{\partial x} = A_1 B_1 t \cos\left(B_1 x\right) + A_2 B_2 \cos \left(B_2  x + C_2 t \right)\end{equation}$$$$\begin{equation}\frac{\partial^2 \alpha(x,t)}{\partial x^2} = -A_1 B_1^2 t \sin\left(B_1 x\right) - A_2 B_2^2 \sin \left(B_2  x + C_2 t \right)\end{equation}$$$$\begin{equation}\frac{\partial \alpha(x,t)}{\partial t} = A_1 \sin\left(B_1 x\right) + A_2 C_2 \cos \left(B_2  x + C_2 t \right)\end{equation}$$ *N.B.*: Don't transcribe these equations. Please download the appropriate files from the [Appendix](Appendix).  Domain geometry, boundary conditions, initial conditions, and stopping conditionThe domain geometry is a rectangle that spans [0, 1] in $x$ and [0, 0.5] in $y$. This elongated domain was chosen to allow multiple peaks and valleys in $\eta_{sol}$ without stretching the interface too much in the $y$ direction (which causes the thickness of the interface to change) or having large regions where $\eta_{sol}$  never deviates from 0 or 1. Periodic boundary conditions are applied along the $x = 0$ and the $x = 1$ boundaries to accomodate the periodicity of $\alpha(x,t)$. Dirichlet boundary conditions of $\eta$ = 0 and $\eta$ = 1 are applied along the $y = 0$ and the $y = 0.5$ boundaries, respectively. These boundary conditions are chosen to be consistent with $\eta_{sol}(x,y,t)$. The initial condition is the manufactured solution at $t = 0$:$$\begin{equation}\eta_{sol}(x,y,0) = \frac{1}{2}\left[ 1 - \tanh\left( \frac{y-\left(\frac{1}{4}+A_2 \sin(B_2 x) \right)}{\sqrt{2 \kappa}} \right) \right] \end{equation}$$The stopping condition for all calculations is when t = 8 time units, which was chosen to let $\alpha(x,t)$ evolve substantially, while still being slower than the characteristic time for the phase evolution (determined by the CFL condition for a uniform mesh with a reasonable level of resolution of $\eta_{sol}$).  Parameter valuesThe nominal parameter values for the governing equation and manufactured solution are given below. The value of $\kappa$ will change in Part (b) in the following section and the values of $\kappa$ and $C_2$ will change in Part (c).| Parameter | Value ||-----------|-------|| $\kappa$  | 0.0004|| $A_1$     | 0.0075|| $B_1$     | 0.03  || $A_2$     | 8.0   || $B_2$     | 22.0  || $C_2$     | 0.0625|  Benchmark simulation instructionsThis section describes three sets of tests to conduct using the MMS problem specified above. The primary purpose of the first test is provide a computationally inexpensive problem to verify a simulation code. The second and third tests are more computationally demanding and are primarily designed to serve as a basis for performance comparisons.  Part (a)The objective of this test is to verify the accuracy of your simulation code in both time and space. Here, we make use of convergence tests, where either the mesh size (or grid point spacing) or the time step size is systematically changed to determine the response of the error to these quantities. Once a convergence test is completed the order of accuracy can be calculated from the result. The order of accuracy can be compared to the theoretical order of accuracy for the numerical method employed in the simulation. If the two match (to a reasonable degree), then one can be confident that the simulation code is working as expected. The remainder of this subsection will give instructions for convergence tests for this MMS problem.Implement the MMS problem specified above using the simulation code of your choice. Perform a spatial convergence test by running the simulation for a variety of mesh sizes. For each simulation, determine the discrete $L_2$ norm of the error at $t=8$:$$\begin{equation}    L_2 = \sqrt{\sum\limits_{x,y}\left(\eta^{t=8}_{x,y} - \eta_{sol}(x,y,8)\right)^2 \Delta x \Delta y}\end{equation}$$For all of these simulations, verify that the time step is small enough that any temporal error is much smaller that the total error. This can be accomplished by decreasing the time step until it has minimal effect on the error. Ensure that at least three simulation results have $L_2$ errors in the range $[5\times10^{-3}, 1\times10^{-4}]$, attempting to cover as much of that range as possible/practical. This maximum and minimum errors in the range roughly represent a poorly resolved simulation and a very well-resolved simulation.For at least three simulations that have $L_2$ errors in the range $[5\times10^{-3}, 1\times10^{-4}]$, save the effective mesh size and $L_2$ error in a CSV or JSON file. Upload this file to the PFHub website as a 2D data set with the effective mesh size as the x-axis column and the $L_2$ error as the y-axis column. Calculate the effective element size as the square root of the area of the finest part of the mesh for nonuniform meshes. For irregular meshes with continous distributions of element sizes, approximate the effective mesh size as the average of the square root of the area of the smallest 5% of the elements.Next, confirm that the observed order of accuracy is approximately equal to the expected value. Calculate the order of accuracy, $p$, with a least squares fit of the following function:$$\begin{equation}    \log(E)=p \log(R) + b\end{equation}$$where $E$ is the $L_2$ error, $R$ is the effective element size, and b is an intercept. Deviations of ±0.2 or more from the theoretical value are to be expected (depending on the range of errors considered and other factors).Finally, perform a similar convergence test, but for the time step, systematically changing the time step and recording the $L_2$ error. Use a time step that does not vary over the course of any single simulation. Verify that the spatial discretization error is small enough that it does not substantially contribute to the total error. Once again, ensure that at least three simulations have $L_2$ errors in the range $[5\times10^{-3}, 1\times10^{-4}]$, attempting to cover as much of that range as possible/practical. Save the effective mesh size and $L_2$ error for each individual simulation in a CSV or JSON file. [Upload this file to the PFHub website](https://pages.nist.gov/pfhub/simulations/upload_form/) as a 2D data set with the time step size as the x-axis column and the $L_2$ error as the y-axis column. Confirm that the observed order of accuracy is approximately equal to the expected value.  Part (b)Now that your code has been verified in (a), the objective of this part is to determine the computational performance of your code at various levels of error. These results can then be used to objectively compare the performance between codes or settings within the same code. To make the problem more computationally demanding and stress solvers more than in (a), decrease $\kappa$ by a factor of $256$ to $1.5625\times10^{-6}$. This change will reduce the interfacial thickness by a factor of $16$.Run a series of simulations, attempting to optimize solver parameters (mesh, time step, tolerances, etc.) to minimize the required computational resources for at least three levels of $L_2$ error in range  $[5\times10^{-3}, 1\times10^{-5}]$. Use the same CPU and processor type for all simulations. For the best of these simulations, save the wall time, number of computing cores, maximum memory usage, and $L_2$ error for each individual simulation in a CSV or JSON file. [Upload this to the PFHub website](https://pages.nist.gov/pfhub/simulations/upload_form/) as a 3D data set with the wall time as the x-axis column, the number of computing cores as the y-axis column, and the $L_2$ error as the z-axis column. (The PFHub upload system is currently limited to three columns of data. Once this constraint is relaxed, the maximum memory usage data will be incorporated as well.)  Part (c)This final part is designed to stress time integrators even further by increasing the rate of change of $\alpha(x,t)$. Increase $C_2$ to $0.5$. Keep $\kappa= 1.5625\times10^{-6}$ from (b).Repeat the process from (b), uploading the wall time, number of computing cores, processor speed, maximum memory usage, and $L_2$ error at $t=8$ to the PFHub website.  ResultsResults from this benchmark problem are displayed on the [simulation result page]({{ site.baseurl }}/simulations) for different codes.  FeedbackFeedback on this benchmark problem is appreciated. If you have questions, comments, or seek clarification, please contact the [CHiMaD phase field community](https://pages.nist.gov/chimad-phase-field/community/) through the [Gitter chat channel](https://gitter.im/usnistgov/chimad-phase-field) or by [email](https://pages.nist.gov/chimad-phase-field/mailing_list/). If you found an error, please file an [issue on GitHub](https://github.com/usnistgov/chimad-phase-field/issues/new).  Appendix Computer algebra systemsRigorous verification of software frameworks using MMS requires posing the equation and manufacturing the solution with as much complexity as possible. This can be straight-forward, but interesting equations produce complicated source terms. To streamline the MMS workflow, it is strongly recommended that you use a CAS such as SymPy, Maple, or Mathematica to generate source equations and turn it into executable code automatically. For accessibility, we will use [SymPy](http://www.sympy.org/), but so long as vector calculus is supported, and CAS will do.  Source term<jupyter_code>from sympy import Symbol, symbols, simplify
from sympy import Eq, sin, cos, cosh, sinh, tanh, sqrt
from sympy.physics.vector import divergence, gradient, ReferenceFrame, time_derivative
from sympy.printing import pprint
from sympy.abc import kappa, S, t, x, y

# Spatial coordinates: x=R[0], y=R[1], z=R[2]
R = ReferenceFrame('R')

# sinusoid amplitudes
A1, A2 = symbols('A1 A2')
B1, B2 = symbols('B1 B2')
C1, C2 = symbols('C1 C2')

# Define interface offset (alpha)
alpha = (1/4 + A1 * t * sin(B1 * R[0]) 
             + A2 * sin(B2 * R[0] + C2 * t)
        ).subs({R[0]: x, R[1]: y})

# Define the solution equation (eta)    
eta = (1/2 * (1 - tanh((R[1] - alpha) /
                        sqrt(2*kappa)))
      ).subs({R[0]: x, R[1]: y})

# Compute the initial condition
eta0 = eta.subs({t: 0, R[0]: x, R[1]: y})

# Compute the source term from the equation of motion
S = simplify(time_derivative(eta, R)
             + 4 * eta * (eta - 1) * (eta - 1/2)
             - divergence(kappa * gradient(eta, R), R)
            ).subs({R[0]: x, R[1]: y})
pprint(Eq(symbols('alpha'), alpha))
pprint(Eq(symbols('eta'), eta))
pprint(Eq(symbols('eta0'), eta0))
pprint(Eq(symbols('S'), S))<jupyter_output>⎛           ⎛√2⋅(A₁⋅t⋅sin(B₁⋅x) + A₂⋅sin(B₂⋅x + C₂⋅t) - y + 0.25)⎞        
    ⎜0.5⋅√κ⋅tanh⎜────────────────────────────────────────────────────⎟ - 0.25⋅
    ⎝           ⎝                        2⋅√κ                        ⎠        
S = ──────────────────────────────────────────────────────────────────────────
                                                                              

                                          ⎞ ⎛    2⎛√2⋅(A₁⋅t⋅sin(B₁⋅x) + A₂⋅sin
√2⋅(A₁⋅sin(B₁⋅x) + A₂⋅C₂⋅cos(B₂⋅x + C₂⋅t))⎟⋅⎜tanh ⎜───────────────────────────
                                          ⎠ ⎝     ⎝                        2⋅√
──────────────────────────────────────────────────────────────────────────────
                 √κ                                                           

(B₂⋅x + C₂⋅t) - y + 0.25)⎞    ⎞
─────────────────────────⎟ - 1⎟
κ                        ⎠    ⎠
───────────────────────────────<jupyter_text>Code  PythonCopy the first cell under Source Term directly into your program.For a performance boost, convert the expressions into lambda functions:```pythonfrom sympy.utilities.lambdify import lambdifyapy = lambdify([x, y], alpha, modules='sympy')epy = lambdify([x, y], eta,   modules='sympy')ipy = lambdify([x, y], eta0,  modules='sympy')Spy = lambdify([x, y], S,     modules='sympy')``` *N.B.*: You may need to add coefficients to the variables list.  C<jupyter_code>from sympy.utilities.codegen import codegen

[(c_name, c_code), (h_name, c_header)] = codegen([('alpha', alpha),
                                                  ('eta', eta),
                                                  ('eta0', eta0),
                                                  ('S', S)],
                                                 language='C', prefix='MMS', project='PFHub')
print(c_code)<jupyter_output>/******************************************************************************
 *                      Code generated with sympy 1.1.1                       *
 *                                                                            *
 *              See http://www.sympy.org/ for more information.               *
 *                                                                            *
 *                        This file is part of 'PFHub'                        *
 ******************************************************************************/
#include "MMS.h"
#include <math.h>

double alpha(double A1, double A2, double B1, double B2, double C2, double t, double x) {

   double alpha_result;
   alpha_result = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25;
   return alpha_result;

}

double eta(double A1, double A2, double B1, double B2, double C2, double kappa, double t, double x, double y) {

   double eta_result;
   eta_result = -0.5*tanh((1.0L/2.0L)*sqrt(2)*(-A1*t*sin(B1*x)[...]<jupyter_text>C++<jupyter_code>from sympy.printing.cxxcode import cxxcode
print("α:")
cxxcode(alpha)
print("η:")
cxxcode(eta)
print("η₀:")
cxxcode(eta0)
print("S:")
cxxcode(S)<jupyter_output>S:<jupyter_text>Fortran<jupyter_code>from sympy.printing import fcode
print("α:")
fcode(alpha)
print("η:")
fcode(eta)
print("η₀:")
fcode(eta0)
print("S:")
fcode(S)<jupyter_output>S:<jupyter_text>Julia<jupyter_code>from sympy.printing import julia_code
print("α:")
julia_code(alpha)
print("η:")
julia_code(eta)
print("η₀:")
julia_code(eta0)
print("S:")
julia_code(S)<jupyter_output>S:<jupyter_text>Mathematica<jupyter_code>from sympy.printing import mathematica_code
print("α:")
mathematica_code(alpha)
print("η:")
mathematica_code(eta)
print("η₀:")
mathematica_code(eta0)
print("S:")
mathematica_code(S)<jupyter_output>S:
{'fraction_non_alphanumeric': 0.1227452705675319, 'fraction_numerical': 0.016375812680256147, 'mean_word_length': 3.9237063778580024, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '7602', 'license': 'MIT', 'n_tokens_mistral': 6009, 'n_tokens_neox': 5624, 'n_words': 2458, 'path': 'benchmarks/benchmark7.ipynb', 'repo_name': 'stvdwtt/chimad-phase-field', 'size': '20461'}