diff --git a/benchmarks/benchmark7.ipynb b/benchmarks/benchmark7.ipynb index 56be1dc2a..a88db049b 100644 --- a/benchmarks/benchmark7.ipynb +++ b/benchmarks/benchmark7.ipynb @@ -10,7 +10,7 @@ "application/javascript": [ " MathJax.Hub.Config({\n", " TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n", - " });" + " });\n" ], "text/plain": [ "" @@ -123,7 +123,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 4, "metadata": {}, "outputs": [ { @@ -135,7 +135,7 @@ "" ] }, - "execution_count": 1, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -212,7 +212,7 @@ "\n", "where $A_1$, $B_1$, $A_2$, $B_2$, and $C_2$ are constant parameters. \n", "\n", - "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.\n", + "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_{\\mathrm{sol}}$ repeating itself.\n", "\n", "Inserting the manufactured solution into the governing equation and solving for $S(x,y,t)$ yields:\n", "\n", @@ -234,7 +234,7 @@ "\\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)\n", "\\end{equation}$$\n", "\n", - "** *N.B.*: Don't transcribe these equations. Please download the appropriate files from the [Appendix](#Appendix) **." + "> *N.B.*: Don't transcribe these equations. Please download the appropriate files from the [Appendix](#Appendix)." ] }, { @@ -242,15 +242,15 @@ "metadata": {}, "source": [ "# Domain geometry, boundary conditions, initial conditions, and stopping condition\n", - "The 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$ = 1 and $\\eta$ = 0 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$:\n", + "The 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_{\\mathrm{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_{\\mathrm{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$ = 1 and $\\eta$ = 0 are applied along the $y = 0$ and the $y = 0.5$ boundaries, respectively. These boundary conditions are chosen to be consistent with $\\eta_{\\mathrm{sol}}(x,y,t)$. The initial condition is the manufactured solution at $t = 0$:\n", "\n", "$$\n", "\\begin{equation}\n", - "\\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] \n", + "\\eta_{\\mathrm{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] \n", "\\end{equation}\n", "$$\n", "\n", - "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}$)." + "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_{\\mathrm{sol}}$)." ] }, { @@ -283,15 +283,16 @@ "metadata": {}, "source": [ "## Part (a)\n", + "\n", "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.\n", "\n", "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$:\n", "\n", "$$\\begin{equation}\n", - " L_2 = \\sqrt{\\sum\\limits_{x,y}\\left(\\eta^{t=8}_{x,y} - \\eta_{sol}(x,y,8)\\right)^2 \\Delta x \\Delta y}\n", + " L_2 = \\sqrt{\\sum\\limits_{x,y}\\left(\\eta^{t=8}_{x,y} - \\eta_{\\mathrm{sol}}(x,y,8)\\right)^2 \\Delta x \\Delta y}\n", "\\end{equation}$$\n", "\n", - "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.\n", + "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 $[1 \\times 10^{-4}, 5 \\times 10^{-3}]$, 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.\n", "\n", "Save the effective element size, $h$, and the $L_2$ error for each simulation.\n", "[Archive this data](https://github.com/usnistgov/pfhub/issues/491) in a\n", @@ -300,7 +301,7 @@ "the finest part of the mesh for nonuniform meshes. For irregular meshes\n", "with continuous distributions of element sizes, approximate the effective\n", "element size as the average of the square root of the area of the smallest\n", - "5% of the elements. Then [submit your results on the PFHub website](https://pages.nist.gov/pfhub/simulations/upload_form/) 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.\n", + "5% of the elements. Then [submit your results on the PFHub website](https://pages.nist.gov/pfhub/simulations/upload_form/) 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.\n", "\n", "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:\n", "\n", @@ -308,9 +309,9 @@ " \\log(E)=p \\log(R) + b\n", "\\end{equation}$$\n", "\n", - "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).\n", + "where $E$ is the $L_2$ error, $R$ is the effective element size, and b is an intercept. Deviations of $\\pm 0.2$ or more from the theoretical value are to be expected (depending on the range of errors considered and other factors).\n", "\n", - "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. [Archive the effective mesh size and $L_2$ error](https://github.com/usnistgov/pfhub/issues/491) for each individual simulation in a CSV or JSON file. [Submit your results 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." + "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 $[1\\times10^{-4}, 5\\times10^{-3}]$, attempting to cover as much of that range as possible/practical. [Archive the effective mesh size and $L_2$ error](https://github.com/usnistgov/pfhub/issues/491) for each individual simulation in a CSV or JSON file. [Submit your results 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." ] }, { @@ -318,10 +319,10 @@ "metadata": {}, "source": [ "## Part (b)\n", - "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$.\n", + "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 \\times 10^{-6}$. This change will reduce the interfacial thickness by a factor of $16$.\n", "\n", - "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 (in seconds), number of computing cores, normalized computing cost (wall time in seconds $\\times$ number of cores $\\times$ nominal core speed $/$ 2 GHz), maximum memory usage, and $L_2$ error at $t=8$ for each individual simulation. [Archive this data](https://github.com/usnistgov/pfhub/issues/491) in a\n", - "CSV or JSON file with one column (or key) for each of the quantities mentioned above. [Submit your results to the PFHub website](https://pages.nist.gov/pfhub/simulations/upload_form/) as two 2D data sets. For the first data set use the $L_2$ error as the x-axis column and the normalized computational cost as the y-axis column. For the second data set, use the $L_2$ error as the x-axis column and the wall time as the y-axis column." + "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 $[1 \\times 10^{-5}, 5 \\times 10^{-3}]$. Use the same CPU and processor type for all simulations. For the best of these simulations, save the wall time (in seconds), number of computing cores, normalized computing cost (wall time in seconds $\\times$ number of cores $\\times$ nominal core speed $/$ 2 GHz), maximum memory usage, and $L_2$ error at $t=8$ for each individual simulation. [Archive this data](https://github.com/usnistgov/pfhub/issues/491) in a\n", + "CSV or JSON file with one column (or key) for each of the quantities mentioned above. [Submit your results to the PFHub website](https://pages.nist.gov/pfhub/simulations/upload_form/) as two 2D data sets. For the first data set use the $L_2$ error as the $x$-axis column and the normalized computational cost as the $y$-axis column. For the second data set, use the $L_2$ error as the $x$-axis column and the wall time as the $y$-axis column." ] }, { @@ -342,7 +343,7 @@ "\n", "## Part (a) Guidelines\n", "\n", - "Two data items are required in the \"Data Files\" section of the [upload form]. The data items should be labeled as `spatial` and `temporal` in the `Short name of data` box. The 2D radio button should be checked and the columns corresponding to the x-axis (either $\\Delta t$ or $\\Delta x$) and the y-axis ($e_{L2}$) should be labeled correctly for each CSV file. The CSV file for the spatial data should have the form\n", + "Two data items are required in the \"Data Files\" section of the [upload form]. The data items should be labeled as `spatial` and `temporal` in the `Short name of data` box. The 2D radio button should be checked and the columns corresponding to the $x$-axis (either $\\Delta t$ or $\\Delta x$) and the $y$-axis ($e_{L2}$) should be labeled correctly for each CSV file. The CSV file for the spatial data should have the form\n", "\n", "```\n", "mesh_size,L2_error\n", @@ -364,7 +365,7 @@ "\n", "## Parts (b) and (c) Guidelines\n", "\n", - "Two data items are required in the \"Data Files\" section of the [upload form]. The data items should be labeled as `cost` and `time` in the `Short name of data` box. The 2D radio button should be checked and the columns corresponding to the x-axis ($e_{L2}$) and the y-axis (either $F_{\\text{cost}}$ or $t_{\\text{wall}}$) should be labeled correctly for each CSV file. The CSV file for the cost data should have the form\n", + "Two data items are required in the \"Data Files\" section of the [upload form]. The data items should be labeled as `cost` and `time` in the `Short name of data` box. The 2D radio button should be checked and the columns corresponding to the $x$-axis ($e_{L2}$) and the $y$-axis (either $F_{\\text{cost}}$ or $t_{\\text{wall}}$) should be labeled correctly for each CSV file. The CSV file for the cost data should have the form\n", "\n", "```\n", "cores,wall_time,memory,error,cost\n", @@ -396,14 +397,12 @@ }, { "cell_type": "markdown", - "metadata": { - "collapsed": true - }, + "metadata": {}, "source": [ "# Appendix\n", "\n", "## Computer algebra systems\n", - "Rigorous 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, CAS will do." + "Rigorous 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, any CAS will do." ] }, { @@ -417,88 +416,332 @@ "cell_type": "code", "execution_count": 5, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "alpha = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25 \n", - "\n", - "eta = -0.5*tanh(sqrt(2)*(-A1*t*sin(B1*x) - A2*sin(B2*x + C2*t) + y - 0.25)/(2*sqrt(kappa))) + 0.5 \n", - "\n", - "eta0 = -0.5*tanh(sqrt(2)*(-A2*sin(B2*x) + y - 0.25)/(2*sqrt(kappa))) + 0.5 \n", - "\n", - "S = -(tanh(sqrt(2)*(A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) - y + 0.25)/(2*sqrt(kappa)))**2 - 1)*(0.5*sqrt(kappa)*((A1*B1*t*cos(B1*x) + A2*B2*cos(B2*x + C2*t))**2 + 1)*tanh(sqrt(2)*(A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) - y + 0.25)/(2*sqrt(kappa))) - 0.5*sqrt(kappa)*tanh(sqrt(2)*(A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) - y + 0.25)/(2*sqrt(kappa))) + 0.25*sqrt(2)*kappa*(A1*B1**2*t*sin(B1*x) + A2*B2**2*sin(B2*x + C2*t)) + 0.25*sqrt(2)*(A1*sin(B1*x) + A2*C2*cos(B2*x + C2*t)))/sqrt(kappa)\n" - ] - } - ], + "outputs": [], "source": [ "# Sympy code to generate expressions for PFHub Problem 7 (MMS)\n", "\n", - "from sympy import symbols, simplify\n", - "from sympy import sin, cos, cosh, tanh, sqrt\n", - "from sympy.physics.vector import divergence, gradient, ReferenceFrame, time_derivative\n", + "from sympy import sin, cos, cosh, sech, tanh, sqrt, symbols\n", + "from sympy import Rational, diff, expand, simplify\n", + "from sympy import Eq, Function, init_printing\n", + "from sympy.abc import x, y, t\n", + "from sympy.physics.vector import ReferenceFrame, divergence, gradient, time_derivative\n", "from sympy.utilities.codegen import codegen\n", - "from sympy.abc import kappa, S, x, y, t\n", "\n", - "# Spatial coordinates: x=R[0], y=R[1], z=R[2]\n", - "R = ReferenceFrame('R')\n", + "init_printing()\n", "\n", - "# sinusoid amplitudes\n", - "A1, A2 = symbols('A1 A2')\n", - "B1, B2 = symbols('B1 B2')\n", - "C2 = symbols('C2')\n", + "R = ReferenceFrame(\"R\") # spatial coords\n", + "X = R[0]\n", + "Y = R[1]\n", "\n", - "# Define interface offset (alpha)\n", - "alpha = 0.25 + A1 * t * sin(B1 * R[0]) \\\n", - " + A2 * sin(B2 * R[0] + C2 * t)\n", + "A1, A2 = symbols(\"A_1 A_2\", real=True)\n", + "B1, B2 = symbols(\"B_1 B_2\", real=True)\n", + "C2, κ = symbols(\"C_2 kappa\", real=True)\n", "\n", - "# Define the solution equation (eta)\n", - "eta = 0.5 * (1 - tanh((R[1] - alpha) / sqrt(2*kappa)))\n", + "# Define α symbolically to simplify the expressions for visual comparison\n", "\n", - "# Compute the source term from the equation of motion\n", - "source = simplify(time_derivative(eta, R)\n", - " + 4 * eta * (eta - 1) * (eta - 1/2)\n", - " - kappa * divergence(gradient(eta, R), R))\n", + "α = Function(\"alpha\", real=True)(X, t)\n", + "ψ = (Y - α) / sqrt(2 * κ)\n", + "η = (1 - tanh(ψ)) / 2\n", "\n", - "# Replace R[i] with (x, y)\n", - "alpha = alpha.subs({R[0]: x, R[1]: y})\n", - "eta = eta.subs({R[0]: x, R[1]: y})\n", - "eta0 = eta.subs(t, 0)\n", - "source = source.subs({R[0]: x, R[1]: y})\n", + "source = time_derivative(η, R) \\\n", + " + 4 * η * (η - 1) * (η - Rational(1, 2)) \\\n", + " - divergence(κ * gradient(η, R), R)\n", "\n", - "print(\"alpha =\", alpha, \"\\n\")\n", - "print(\"eta =\", eta, \"\\n\")\n", - "print(\"eta0 =\", eta0, \"\\n\")\n", - "print(\"S =\", source)" + "# Nested `simplify` calls to incorporate the hyperbolic trig\n", + "# identity, which is not implemented in SymPy\n", + " \n", + "S = simplify(\n", + " expand(\n", + " simplify(\n", + " expand(source)\n", + " ).subs({cosh(ψ)**2: 1 / sech(ψ)**2})\n", + " )\n", + ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Code" + "Compare the generated expressions for $S(x, y, \\alpha)$, $\\alpha$, and the derivatives of $\\alpha$:" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 6, + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3oAAABBCAYAAACHBtaLAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dXbLkNLaFDxU1AChGcGEGdDGAGxQP9x2oEVDMoDvqrd8qYAZ0j4CfGdBE9DvVPQPoETScGdRdn47kUjolW3baaWfmUoTTsn62tpa2pL0l2fnOmzdv7uwOEfjrX//6rUJ+0v3Hwxg/bYWA22Qr5OeVq/Z6Vzl/1vWJ/PfzqDiXETACRsAIGIE2BDzvtOGUUlmvSkhczn2OjL9jQ++wgaPgP9H988MYP22FwJJtIlqjKxtK885Wdb2mcoXjR6rPD7r+JL+NvStoXLUjBvzLWJUP4v1Lt+8VNK6rYASuAAHPO22NKJzY0FhE1xUt5vnPxkpWOutWYyA1xE+V8ccNNG8micB7oco+1f1PN1PpnVd0yTaJtL7S/W87r/ZVsCec/63ra1Xm77q8cHIVrXr3tdr0q1QV+VEW/qXrwxTm+3IICN9nooZxjVENxj8o7B+62xkBI1BAQP3D804BlzxIGC2m60ZanIDzHJ+DvKJ/qow/WpGXiyIt4Nh9QCn95KIYv2JmV2iTz0VzkpGn9J/pSjsXu0N7Tf5E+8+1Cg/F5XmUDrzfhc883P6LReCF2hLjIznGzA8UxvhptzwCGHY/6vpGpP+i66flizBFI3BdCKi/eN6pNGkcq5fUdSfrVbAmPqxbFdpIuFT1rjy50jXLuHf03iLHrsPfBJ6PmL3FZGvfYm0SOw+7D3fyNx0/UzpWvX7T/betgSiVvzR/osfODIN2qK/u3+j6Vle3g5P4GIpLabI7K33/UZ5/6HL/yoC5QC+y8PoC+b5UlvPTJSw4uf9cakua73Mj4HmnjPgqehVFaX63blXAXLhsqlvZ0HsQTnYbWJH2bl5BSLcIUsdYuk0+Fc1PY11Gj5/FAYs8uzyOsDR/kd5HuveNWnYUwIvdhL4biuvSKu+9ru8VwASzCZ4qn/7NThSKMpMRBrw/tiQgpjhh1t8Rx/ADy39PoeO0DwiMyaXi8/7IKvwm/cfttS8ExuRmX9xuw40w2nzeWbLmqk+TETVUpmisqVdRtHWrXgPEdttUt7Kh99AoTKDsXqAE2u0DgcXaRO3KVnjYzYtV4/hZ/q4LZRFGZ0wKK2F5nph1N7el+cMISnXvKik82IVj8Ob45UH/GIrrCLz1wO+vysMxv1x5fZui0af8aae16V0lpQ8TpO6dkiw/dbrTZWMv4i4siru3tWZReoxnFId816mWfFfhkXdWWZuc0i/+EQHRbJZLpWUs+ovuR320qQJOtBkCajP6yXNdTUeyYFR5qvKmuGa5gdaNu8XmnR3gOGpENfAIHovoupLDvl5F8datjhthc93q5g09CSuDMEdiXh23j0O2QGDJNomTYr6bR5Vajp89U17SHTiF0WmRmY/lzw0HDMfu+SDThAfRCBOTsjCZv6/nsJOme9/QqvGXlAk+3AAtFPHgF40jA0thyD71JB0roOT5TvdcocSgot4lo2goTlkeHGXrIi1llXYHY8rhm2gEw0L3/s7SUManiuy/Z4nijKJfqtMQrWuOQ5Hgpfq0812tq9KAJ7JyqV9UTbzncl6t70oRTXIprOEVI4/++WddvK9ndzkIpPabPe71qtokN708N/movrLIvLMT8FqMqCqrwmIxXVe00E/6ehVlW7eKLSCMdqNbPYo83fLtpSq/+btDdEJdb6ZcV9xok9qEDqWLXZ6SgxYTbeeUtv8uJoNTd/wMenq+7zJEj8IZ3GgnFC1eJOb5TncGUAyhk5zo/CoCGFnwRxnUKxlu1CM4hdX4I236cMMT+dmhgQ6KATweOcVTbxQQ6vsKv66+8vuL4mrK/1Bcvzw+wVxrp37ao2fxBd4vdae9mp3S078Pdp30jNGL8RzasJnYCQlVFu22iWspO2LCgkWSuSKvkRZywkSP8YGcbla3IpMDgeI19AXd+3I+kGv5KJU/KpdKw2IE7fGH/Pw1zCTZX4JrlbtZ225Z9kLYLS5rwmRUbpbgvUZj6zaZUf5J804Nhw3CW4yoIbZW1asoWG1j3Sq2gLDYjW518zt6apPPdJ198oyykN/2sMKc87Olv7lN1JkwsDBCWOn+XhcGS3Dyo8RjmFVXUhXHREx5uSGAYvO7rr57ogAGMtLTiVNZ8PC6n3jKs2h9q/R8FjpXPn9SGH8HQRj+5Gr8oXSnXTvqnvJ82aOb6OR3cKgdhaSeNWVvKC6njx/6GFcYyXN20piwwWmyU3nJeIFfsEnl06aErepi+bRNap9VyysQxxgD98HdIMUj3xyx7U/YgaTCkQPaAEMPmcExflb7WEixrx8UnoPFn63YE4aDcqn4fFw6O5uRv93L7dmBaS9wFVkbk5t29qal3IE8wHDTWJbV7NR5JyO1nVfY90+xMO52C9QNnJ1Nr4IX8cv8YN3qYZF9U93qpg29KIjIZK0RiAtOaVFwkjLzVH4MgdLuR0g/5SfxoXuu5E8hsWnayD8f2gAjjK3ZhnPCQnRG24RKKz3pWOHE2GJSTW0kb1Dm8mfCOqc88JsM7FzZxxDIn0MepQ9Kuu7PFZAMBeI+1pWMqjvFM8C1KJIYcUnx/0J5+kc/iYNHdk7yetT4yycC8AjtoLyDcqV40mIkHtVZ4ThkveaG4g7yiD6TEnXCMM/xO0hXelA+cHiqe21nsZTtTunB6mdd9NVOpuQP7aV7wr+Yf4lAlcFkxzHcQSNribJqNFQ2fYRFjxe6cjkpZcGQoz/35ZG07DCB6cH7baI5u89D9FxOfCJHKIq5LBCW+teiY3utXip/c7ms8ZbCxeOlyW1ifRf3kqzBWAyfJW/Ku5ncqOzN5SHiN2UsA+/Z885agiSe0BFm60wxP+3RtBAU01OdbtzjoeaUfrZeBU3lZ0y1brUT3eqmDT0JIpM6Qjmo7EWh5Rhcp2TKjxD/izBdTZ2Hsipu8VU/+NOVJpNKscsEqxwMCd7V4ejhkxOpNrVJoYxXCuPYGccPMVqYEFHoikaOwoNyp3to0/h8pzuywEX+mmOQRhlODkOJ8oMTDcrsZCUGj90ory+HGFGEd7QjkUH+VD78dbt7egaHPu1IKtxIH2RY6Uq4gVUt/1BcXkbyU05o4xTQeEeWxwyUEin6KRN9qY+Wwko0ZodFPDlu2jQhzy4oyxjLxBBjbOraTX5ewmfMOtj5zrImLzgzfrD7ep8Cuev5vfz5Av3IUTcuqj7I75pjew2iTeWyxlQKFy6MA5cmt4n9vdwPZA2mFpC3TeRmC3mIeCGHp4xlkMHNnXceci/8Kzxn60xRhpCDKe9Hn0WvAqbIXzjanz0TzlzERZvWnHWrt8gspls9ekvzJn0o40VDoIcGnepgxVpCyyCOEsRxstkudoqDFebZxGLGSHOoM51aRC1/p1TWEjSEt7bJASnVmR0iysdoxtFmnUIXQuJPxAdDDQWPXQ4GF9JiWN3pGZmgkw25lBYjD4W4RY6G6JXKfE4G0U6Ga8hf4k9hyNAfsQDy5W3R4aA08Ao2ueNjLb/EAPDL8xIMFhjxJTcUV0oPHfCe6ni3r9s1nZCZfN8V0jPx5cZ6IckiQWB9jnJyZpFJ2iXIaB4hP7z02/8gieSDcQ15ZJd5Vaeynuma/d7mFOZUDmMiu8K5gQ8Wq4ztI7xtLZcj7AUZuSi5HavQOeMrsgYLp8rbVnKzxTgGXieNZRCIbu68k/Kvde/PtYPlSK4Y14MRJT96AfM+YWNudb0KBiIvjBvWrR5aZBe61aMHXm72l0GkpaORjvdW+sYTCgOKc0tHq4GMEt4p4rVEE8M/m5h+T8lb26TEMzhyNA1D4onuNeOLFULK4Z4u8qHgJscuUK1dKYdjl0y6+GvlJFotd47JQZN3DXmXCn6gzbtSvMcDv7k74E9pkGN2akiLQcTxFmjxDJ3kqBO0cxlhYP5QYUGJ0D3HgXxMErWjlkNx5O270N9URrOxl6V93Sc29Kx8qb8etI/CqTv41eo0RHZq3BcqZ85O5NRy8vS0CfXrt+Nd5AWeEjZ5vtwP1qWjm3maJfzwMcbLEuVAA2W1vzu+5the5DvDfku5LPKWBV6q3GZV2NRbkjUYmi1vG8vNFvIAXkuMZdCZPO+QaU9O7c/cXTWiRnhF7lp03RIZ9Ad0hjG9irzWrXaoWz0uteoNhaFgtAg/Bh07P/cVbGYpKnHgZoX5YEWZMojTVSuvwkbIR2dkkjm3clnlaWJEa5sckRVefHGSuv+sq3pUTmlajp9Bh3bJjaQ75YU/PpoSjATdUYZPXvkWHeTwSLFW+JFsKB3uiL9e2nzX4iGHfpUG5fI93TtDL4YVy1FckG3dj/rJUFxX4LEn0WHSOlB0j5N2IWGSUnmT+gPpdVEGZYVy9Yz/pa5PdK3qVFYwKFctpEwcvIaMWLBgt25ojGAFnDRX4dQWyPEz3ftyPnlsh5auSbKYg0heXZvJZc5LyS/eLlluS1U6a5jwq8kafEyWt8T8VnKzoTxQ9SXGMujMmXfItyeHEYVsce+c2qc/pnVxmYd8CYMseNwr+k16FZSU1rrVDnWrWzf0kE0UmkEn4T1SwGOGsCuh+FaFtV8OyvqrfmB8/ll0+dR++ICD7nRUFFR2XXgP7ejDDgpjgn6uC4dSk46VQqdT/OSHBvRw7z/cAs17/IqnXt2LwvLDZzIMPpb/F6U5Kl/hwSkOZbo5fcyW30bbJE/c82OYPRcPswa1REv52RH7nLr0aIELxxx534lJiJ3DIYU5kVz0rjJr/I2WE/kuGoKFzMmgLER1xmYprhaW2uVJLUEhHJkLslmIGwvCoON9s9BX5UfeP9HzLHrKBx7JIef0xVr/ZzW6irPy0Qfpizjq+KUuaKY+TD/r+q3Cq07pkEUmfGhCg4Up+n+pr8ITvA3JLXXCoDnJqBGNvTjGzrztAl+q35yxffLYXABhUbks0D8IUj3zul+z3B7Ue6OHoqzBy0x5y6uxiNzsVR4iRkuPZZBtnneEDXMFPDBHMJ7i0Cm6xWP56UPoGugqzCmMkwcGl57JS7+DTnLFMT3SG9SZlKbFiErllO6b61UwpXrUdBfrVm9bDbk5kKe3UVVfVcYfV7NceUTsWNTy9zlVVX4Gg9TZJ5NQfgYBjLGjxoxx0GcgudMz/qe6vtWFv6g8Kh1KIasvdGg6U4k2gwlCxIu8QTB0h+7PhOm60zNKXv5xFZTGZHBSPh90oJwkWGRLLtWrNX3Kdyd64Imb1SZkFI2AAf5TnWhxlBIjAWU+DdZghfHHkcj3de8G/1PLm5pfZZf4ayED/0UZyjMrDUbID7oftfNQXE6j71c+djMIRk5aHWlnyQTlKe9RP2gtOKUTHWSTI7FgHrDTHb7+o3voSzENfSUZZ6nPJjL9O7IVeNMdeQoTnfy8g4Gc0U8TrX7eg2elhyf6PP2bcQVDruYYH1BkhlzCm3ozHuzWxTojV0MyzXHoD1sqoXSMcdQ7jL95HsXR5sRPGptzGvhFZxG57NPtP6ucW5PbPgSLPgtP+tdisgZzolmVtz7zSnuS3Cj/ruUh4rH0WAbGTfOO0tG/GZcPxk89d3OI/LRX0JfkTzrUn+VHL0o6FDiz88ZfG4UxXHfG3DCfyp+PqZTJmD1ZZ1K+USe68IJLY/rD04Rf0VhMr6JY0SvpLtatHrCp6l1DTSZMqzL+aCjjlcfRuU5xrJZj7FR3tkaIv1Q8ilzJhSNTop0Ut7BrpGfeu+FLevkgUco/NYy6oKCmASHlZxBjAOqUzazsmqLIUdRup6AhfSqL+6ltktNaxC/+D5Q9PdMmvDfHjt5B3CIFTiQyhwf4byyGd/5qyvNQXAt5VkFb3RMlvG9NvFI6Jm1eMO/wkB+eXutKcoAS2MXLjzwXJ1flZTDP+z+0yJ/ahjrn8XpscigoY+MDPPX7ep94whs+duuEYzKQq1jFNCgRrW5obD/32NzKcy3drcltDYeTw1eSNfgakreT+e4RuBR5gO2lxrIcgrF556kSo8MwducuHz9oL+a/bgFUfvRAdCjGcBxp+u+BJ5ppbA0J9XOKzpRoDN1TuUNpzh4nrNK8GcrWs3Wrh1ZYXLd6fPbW3U+BSYHpd7pRDiWQdHo6cfHYj8Lp7Kes+jHAkZ8dBMpCmVzEiWa3MiM/AwAD28HqVa+gWtm1wWNq+ry4SW0i/t/kmdf0q6w1yZ9Mey3+hujW4hT+TkOFavJTykraosGUEqvMRWUhr4P89MF3dS8t6jDZ04dw7PDeP3jDL/KcP2dRd6+VtlMUFMFKMO9+hvS6F8eWnEDFzwIMfXzIUS6YDrlBvIcyluJUHzAsLQ6FPq/4bsU8yw8eVRwU95HSYhijvLGijqJVMnJZPW7dzRsc21XOpLFZ5S4qlyp/0Km8ru/Jf4tyG/BR3ZFvdl3G5DzHk6P6Jfm5U/jiskbBojsob4o/WX5EI8hELOtSxjHgWWosg1Zyg/IgjDA2GPv+wK87Oli3oK8w5IBFMgzmvkN2Po75SHcwXyi80716GU/RmXqkio9PYuh9MbYXKD5PlrkeycFHlTcYv3XkWvwN0S3FKawb20cwOZLxxyMZrjmazow7AuUhuPwrsFmJZ4etaBwpnBXml7pQplDejlxMwwBfcwxw3+tiVYhBhXdCoMcxgKbOqrRVJxrUmRVwBoDvdP2iaw9uUpuoHq2Cv4e6mYf5CCDzabIqUllZFtjFyXfqch7gjaOw9KdXecSQX+lRCnJHGc3584zJL5r0a8YLlJMhB5bwPeQS3qlPDqUdjRNvJUPuTuEsioHfgVI0SlAJlCdgqDu4pXH3wDBUHOP10FjbFRXTVsf2mHDS2CyaW45Rtyi3oZmEO/JdnH+7Bp/gEb1FZY2iRXNQl4hplpSfi5CHWO8lxzJITnHIDTocYxPjOkc5OcXDGJZOQnBcn/bLHeMMRltK8988ckN/GsPBdNSpXkvK3Gh5TrA+Ards6E1GVx2Ajs9n6DtlQv7QqXVnh+8jxZ+0whzp0SE5y52USwYZ/iONgSUoRHpmZY7JbNQpHQNV2CbXnTrwLhCrl4G+7vBttyME1CYY+bTVqFPaSxyYeU9sFRflGYWc/kFf6h+haS5XtMjPNWQ8Ec9L9vc9wkywxA065YNX0nU7cbHcuwLNIVrQIU8aN/CXxgnKSpN/jV7i+76WYC/hquO9LsbdF7pjNOY7pYx7owq/8gyO7dQV2rqBy6yxWfkZZ0+SS9FIsgIv7FLmvOjxwSkdfHLdmtxGBNa5CdeTZQ3ORGdU3vo1UJ7Z8qO8lyQPVH3JsSyHcnDeEU70q991Z9xIOlNYLFIYhlwaW37Sczde5wUoHKxxTacIHpL691wIqH1uTrd6fC5wd1jOfeTpSQtvEg4GWbblQ+fP8jBgh/dqFLfEqt/RAJeVFbwqh4EHZS7VoZ+kFk46hBzDsVMG9dxhoHDqyWSWBjQ9ns0lvjt+zlbyjgqK7ctE0i0o7Ii9pVhJbd1CD1l8tyWhMCPdyxw7+VnouNNVnJiH6CoPfWEoCXGkKdGG77AQlBNQWnikH8IX/ZCV4n6fow4HYw35dA3hximDMAbpfqe0jE3w0IURLkf/GuvfqQ+OGYSB4A5+WHlnXASzsHsY68+xqyHMwIkxb3BsVzxu9tisMpaSSwy78PW9SJMFwKOFHsXdqtyGhlr5Z7aswZfaplXeumrEth4d15SOI4Us4h70bz1fkjxQ7yXHsg5HeQbHAsXTNozZ3QkDYYe+BD/hXTqwlP+5rqMxX3HPdHH8kzH3qa4jp7jPdB3lPUq4XECqcxrTl6N8YZSEO3PEzelWjy6snZZkNykwTMCDTsJBx0cxQ9HigwzdpTDe/0gdKdCJzxh/dGjy5g5FJBiGeWDmZ0DJjbAsqlPO2FU8GMjzRDF/cZDppUuPDG44sMj5nTowTE0fCs1+mtsky7MLr9oD2UBx50IZ4xqVrQrzTNRDMlLJtv9g1SvJV2rrFqZ/UaJWeSZdKiPRZiWWXey5jragXx441YVJgzYO7azn1I9SOib7j9NDdsdg4GLlmLwHWCiMOOrcuZiO90ZK74akdPS/MC7E9OxuwUPfscNVCs/TgSHK4X0euFe/+KTejJsvYt1hleNXHOusOqWlnq1j+ylj81Jyme9OhjaqVu5hnrk1uR2AY5moubJG6RPlLWd4VH5Em7GE91RrusGljGPUe8mxLOEO3YOxloCCw6AOY3oWx3PSyz6RH4MuLPykNHpmASBhzyIt4y9zROdimnzsPVVn6mgPeFKd+3UayLKfKGFm3aqhOYQT8wEutffDk34fd74b8wiUpMS0bK9z/AUQDzpthCzvtDmKc1f96Iwopn33pQL4LxcUx1J8l15p+Dx7MDYUyMCTKzvh/LnieaeFOOrFBPC+LhQeaHPUDT+KKzjhhwZpw2q57gyGrIJjlJAOpaopvdIVnehMaZMijQ0DMfASNnfygyNKeYt8dWwrH+0y2L5d4sv0pMnmfgL79DEGe67BfIpnMs6V4TuFIc9N+Us8KW/6FDR9Oh39oR6s9HLd6SKOdzLy8eA7PdN3+g4e6XNBUVBe6LPKmGQGA/BgxVfP9A3661A9WERiESqMU7p3q9IKzx3ldrKaR2R+MEyKTRa8ay9tQN0w9sCKD94MyovSTBnbafNS3xwdm8XHInIZ6yU2gqO+1V1/pb1FuY3QrH6bI2swNUXeuko0yg+yn48/XX48FyQPsLvkWAY9+i5ubDwgnrGZMUS34NCNmN8ZU8CRj0QxPqJn0f+YE0jD2JvSMOf8j8LIx5iR5o2QRmGL6EyiO+pUFnMH9Zqki4wSPl8C61ZtWFdl/J03b876gZ02ds+USsKPIo5SdbTquQQLosugzuD7Hh0tljf7z5qX4GnvNCJGq7XJWvUX33Qk/voiKMe6YxQzuPPSdnXy7fOjtCj8q8hjv6wtnlW3z1Quxk/oE608KB9H1FhUGDU+lAZjGcfkxuCH0URbdDvhWRomP5Q2+Ap+xYXJWs8nO9GiXPhuloGhQkUHPjkaRN0mO+VDLpGxwUlf8YyNKCUoPas50ac+sz7GUmIq8k0daUNwX6wtS+VNCRMvi8mlaKE80j6LyFW/HqJ7kXLbr8eaz8KIPnI2WVN5Nfn5X/Hxf7roS4wLjJHfLSkbe5MH1e9OPDWNZTEt2Eyed8h7DU5YrarrromReLdu1QCwcKrK+OOG/Nec5LUq98WKFZy76rciS7snvXabrAUAOyTwPtupozKRH+0YKBw5utOdFUQMF2SWMBYNVlH0RHstx440q51TDZXvlQ8DuGroRWx+Vhpw6tLJz4LLne5ptRWceXcr7VyjMGOkMxmy0rqkcUA7IRtjO2hK0uQ+Fp8YrnMdK+VBnmoERB8ZAwcwX9shB1NlYYinV4pEoaNtl2zHoTIH4yKei8ml6NF+yDiLh/ytRG3ndpCvkciLk9uR+qwRfRZZU/vSH4fk55+K/6fSsaj8SvdTxgeRKLq9yQNMjo5lWU3mzjsZiYv2XqpeBejMndatxsWvKuOPxvNedQoUu3AUao1aasBF2UQRfxmvQQVrDR4ukOaqbbIWHmprjvDlCiuDE8pmkyGmdEzmGBulSZoJ/BfFcf9Cd3ZZFtslEq1zOgyIJkx6TGEAl45O58noX2DeGXlZZB6GgpyMAHAPhqDu/H1JCf+MzDRvbCt2rFh9PsmJBrz+dy6RyAO8jO3SsZAw+hGTuXzk+cRL+JPcPOwUv+jRfsgXSuBe3GJyqfoxPrJQwfuarHQzzizuRBsZuTS5XRyHIYJnlLVW+WFszce5IfYnxe1JHmBc/DCetoxlqZ5z552U/9LvF6lXATqyp8u61bgEVmX81nf0WLFGgUSBXlTBy9rkLKt+WXmX7j1Hm6yKkQYlOhzb6Jzj75zCmZyeMXB1gW89LAbUFgKgh0LHkbuweq/7KhP6W3ZW89HXJiunqi/vRWAUDH2xDEOQdyb67qkCeI8qONHI8e/4oYyUZuE7PLHLdOqR3Bfi8ZTdG8a6Fuwxkk7lVSS2ccLooN9tw8VBqYvJ5ZnrdmlyewD6OR7O1B6j8iM+GMdYwMoV4qUh2Is8UK/WsSxh0I3zKeDG7hevV9Fekm/rVnXBrcr4o3qe64+JgyIK82pKjcrY4wrzbhv3HG2yZuXFP8YcBhvv5nWTrvypE3JMkJ2ZzsVnvpZ2ZLzFfNBB+Wa36lddGJEX52Jd4Puono2VwVjDID5yop0wPTDWIlbgdrSQozAmjW53T8+03eJOdGk/PozBTsxsp/yzjbxYNjykncwiH4pHqURGB9MVMzvwCAHhuFu5PGK2FyDeL0Zue6xfzeME+WEsC+MqeXTxvKjbgzxQIfHBODo6lqXKKz1zL27uvPOQ+4J/Y9tR/9V03bXhUR2sW1VAHpPxmzb0ImastHNUaTWnRkDpt+LUjvDqbdLOSnvKOBDx7gxHMDEguqNP8rMbhbGGIdI3VjAMa0fNwgSuvOxoYayQLuTXc5rAFHQRjkmGeszqC8qH4sk7KEc7nzEObDtjTWH4weoTXcERposPu+Ce68p5qbVBSHzKj8rE2JxtqJ1SNnkpGx6G6CgevDgSvBmfQ/xdYpywRGZ3K5djmCIzW8oDZcPDGJ/XGq+6N8mP6s/HlX6JODDmrYIZbUGbxHI2uVE+fEwo/KR5Z0I5e096kXoVoKq9mZusW9UlbFDGH9Xz3UzM96opK2CXpjRfcwNdXJvEgehbNQq7IezOYaBhOPzeayiO8vIBhXcJj3eMD5TBkqMDp3fIUnyiufiqbSpgpTu7RdR/thNOGLvsapb6KwYdhgr4hlVfPR985VbhKAjfx3hwxQBP6Vcz9FTOJTjqXzr6egm875lHy+WeW2f/vI3Kj6rA3MNXhRlj+eImBqLdAwInzztXAuTF6VXgLlnGyLNuNSyEgzJe/HsFAcvRMBRMXBgwFIY1DTFWdK5qC1z1YYfgme57e7cD/G/SXVqbiF92iYLxljeYwt/Jn/ErjE+X87EL+lQawIqGnuJ5iRrjJRfmjAkAAAfxSURBVKxg6k4Zf9eVjJQpK5vKto0T3xheL3V/bxsOXKoRMAJGwAjcEgKedw5bW3hcnK4rnq1bHTbjwVOLjB8ZesqEAtl99AGKCkO5ZNsXY+hIcSXNpTvVC2G6xM/VXzr0Vf6vtU1ULxZSMNZYqcX4ufpdlNiWi3/Vsio8jjACRsAIGIGbRsDzznHzR0yuUtdV3axbHTf53cHRTYHE7sLRGWyFs6uHoVfcdSjQvcQgji2hfNvtB4GrbBP1J44fcvwy/DfSfuBehxPVl908xpWjD6KsU6KpGgEjYASMwC0j4Hmn2vpXqVdR26hjWLfqNf2Boac4jmZi7JUcR8S+K0VcQ5gEJHxyXXcwsNsBAlfeJgy2vB92EUcv54qD6veB8r7UdfW7lnMxcj4jYASMgBFYDgHPO3Usr1yvouLWrXrN3xl6avz0cYOnvTTpESv5qt7NSxXL7hyl470plFO7fSBwlW0iGeMdvas2flQ/jnxzFJwjm1dt0O6jq5gLI2AEjMBtI+B5p6n9r1KvouZqf+tWPRHoDD2Fv45xfDXwa10HX/TTM59Fv+ajm3eqH0dU+QjND/KjpNptjIDbZOMGOK14jnszlvjI5mk4OrcRMAJGwAi0IeB5ZwQn61UjAO0/epKMH3yMRY1fOrqJksYuV9OKvNJx9DPtDrbChRG5m90N8cKOHn+5cNWGbWvj7CGd22QPrdDOg9qLhZKnul/7KYB2UJzSCBgBI2AEVkPA8840aK1XTcNrD6nnyPiBoUclIhEMPna2ksHGnz/7s+gAZGcEjIARMAJGwAgYASNgBIyAEdg5AkeGXp9fGXjs0AXDT/7NVuejAcpXCqccqeT/x7wr129UPxsBI2AEjIARMAJGwAgYASNw1Qg8pnYyhl7oCl+dLNSWL9hg6E0xsApkTgsSf/eicNIfmovGm9O4cG4jYASMgBEwAkbACBgBI2AEjMD+EXgs4wcDjmOaNUPvSaxG026e6O32HT3xdpV/9r5/MTOHRsAIGAEjYASMgBEwAkbACJwTAXb0eA/v4AubPQbY0ftGRhI7aqNO6b4aTeQERsAIGAEjYASMgBEwAkbACBgBI7AaAvy9Art5dzLQvu6XEsOe6I6xZ2cEjIARMAJGwAgYASNgBIyAETACF4BAekfvQxlzz3Tx3wz8MXo6rvmdjbwLaEWzaASMgBEwAheHgOZXTtM81/XnKcwrn19DmAKY0xoBI2AEbhSB0a9u3igurrYRMAJGwAgYgVURkMH2kwrgf2r9dehVkTZxI2AEjMBtIsDRTTsjYASMgBEwAkZgBAF24HTxwbGTHbQgoruNvJPRNAEjYASMgBEoIWBDr4SKw4yAETACRsAIHCPA6w3p1Ybj2GkhL5X86N34aSSc2ggYASNgBIxAHQEbenVsHGMEjIARMAJGICCgnbdJ79ENwSZaHyj+A92b/rZoiJbjjIARMAJGwAjUELChV0PG4UbACBgBI2AEhIAMMo5Z8hdDTX8z1AAaX7L216wbgHISI2AEjIARmI+ADb352DmnETACRsAI3AYCz2Xs/W2JqorOu6LzVPfibl6MX6Io0zACRsAIGIEbRyD8vcKNY+DqGwEjYASMgBEoIiDDiyObi3yAJRbAe3mvioU9BP6sMvlro2941B3DkPf5Xuh6lcLltzMCRsAIGAEjMIiAd/QG4XGkETACRsAI3CoCMqp4l+5e99+WwCAabfxn7Y8lejGeY6LhS5x6xv+FLgzN17qKu4AKtzMCRsAIGAEjcISAd/SOIHGAETACRsAIGIGAwFcytpZ8l27sS5sYdXcq8x+6PpMXIzMdGf00cOQfI2AEjIARMAKNCHhHrxEoJzMCRsAIGIHbQUAGFkclm49sYpjpejaCEGmS4VZKijGHcccfqePYxbMzAkbACBgBIzALARt6s2BzJiNgBIyAEbhWBGRocWTzXd2bjmwqHe/x/V1X9X/xYpoxwxFD8ftIh//r4329H3Txnp6dETACRsAIGIFJCLzz5s2bSRmc2AgYASNgBIzANSMgwwqD66tCHcNxSoXzrtxvSvcXXbxHh0HIsUwMvj8pLLxjJ3/nFParrg+7gJ5HcRiXv+r6VP7uXTz5/1BY9xEWPWOALvU3Dz0u/GgEjIARMALXhIDf0bum1nRdjIARMAJG4GQEoqHVGVuJYDS6eH/u8ywsfTiFL2li6GHwdfGkU/qWY6Dh2KfSHpULDVykQ7wNvYCIf4yAETACRmAIAR/dHELHcUbACBgBI2AE3iLAEcriMUoZYRhfvH/He3jszuWOD7oMvZtHWt7Pqxl56Qjph6Kd/OSxMwJGwAgYASNQRcA7elVoHGEEjIARMAJGIOyk8W5dMt6eydj6Qc+/6B7+6y7DiHf02L3DsAtHP5WG454/6j62C4cBWXqH70uF84ftH1fiFWxnBIyAETACRuAYAb+jd4yJQ4yAETACRsAIzEJABhlfzOQY5nsYd7r+Jf8n+GcRdCYjYASMgBEwAjMR8NHNmcA5mxEwAkbACBiBAgLpy5svZNyxm/faRl4BJQcZASNgBIzA6gjY0FsdYhdgBIyAETACt4KAjDres+MDLXyUhSsZfvLaGQEjYASMgBE4HwI29M6HtUsyAkbACBiB20CAL3Dyzh1/weCPp9xGm7uWRsAIGIHdIWBDb3dNYoaMgBEwAkbgkhGQcfej+GdXj4+y2BkBI2AEjIAR2ASB/wfbNpWKxDZ7DAAAAABJRU5ErkJggg==", + "text/latex": [ + "$\\displaystyle S = - \\frac{\\left(2 \\sqrt{\\kappa} \\tanh{\\left(\\frac{\\sqrt{2} \\left(y - \\alpha{\\left(x,t \\right)}\\right)}{2 \\sqrt{\\kappa}} \\right)} \\left(\\frac{\\partial}{\\partial x} \\alpha{\\left(x,t \\right)}\\right)^{2} + \\sqrt{2} \\kappa \\frac{\\partial^{2}}{\\partial x^{2}} \\alpha{\\left(x,t \\right)} - \\sqrt{2} \\frac{\\partial}{\\partial t} \\alpha{\\left(x,t \\right)}\\right) \\operatorname{sech}^{2}{\\left(\\frac{\\sqrt{2} \\left(y - \\alpha{\\left(x,t \\right)}\\right)}{2 \\sqrt{\\kappa}} \\right)}}{4 \\sqrt{\\kappa}}$" + ], + "text/plain": [ + " ⎛ 2 2 ⎞ \n", + " ⎜ ⎛√2⋅(y - α(x, t))⎞ ⎛∂ ⎞ ∂ ∂ ⎟ 2⎛√2⋅(y - α(x, t))⎞ \n", + " -⎜2⋅√κ⋅tanh⎜────────────────⎟⋅⎜──(α(x, t))⎟ + √2⋅κ⋅───(α(x, t)) - √2⋅──(α(x, t))⎟⋅sech ⎜────────────────⎟ \n", + " ⎜ ⎝ 2⋅√κ ⎠ ⎝∂x ⎠ 2 ∂t ⎟ ⎝ 2⋅√κ ⎠ \n", + " ⎝ ∂x ⎠ \n", + "S = ───────────────────────────────────────────────────────────────────────────────────────────────────────────\n", + " 4⋅√κ " + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Eq(symbols(\"S\"), S.subs({X: x, Y: y}))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "scrolled": false, + "tags": [] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAArCAYAAAB8dS9hAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAMsElEQVR4Ae2d27XVNhCGN6xTAAkdkA4IqSCkgwQqADoIizfeWKGDJBUk0AGhAi4dhFQAnA7I/+lIxhfJlm1529t7Zi0dWfeZX+MZSfb2ufbly5eDkSFgCBgChsB5IfD06dPbkviFwve6vqxLf1FP2LUhYAgYAobAfhGQA7gh6f5U+KRwR+GWQofMMXQgsQxDwBAwBPaJgN8Z/IJ0uv5VEbuGDl3v5FiGIWAIGAKGwFkjYI7hrKffhDcEDAFDoIuAOYYuJpZjCBgChsBZI2CO4ayn34Q3BAwBQ6CLgDmGLiaWYwgYAobAWSNgjuGsp9+ENwQMAUOgi4A5hi4mlmMIGAKGwFkjYI7hrKffhDcEDAFDoIuAOYYuJpZjCBgChsBZI2CO4cSnX79evLVlEbbO35axK8Gb4V8Cxel9bBz/m16yb9sS7soxMAkK79pC7jUtWZM/ad+QzMwJfG6a9qg7J6Ifm9aLAsxtTv+lFy8Iku2hl8+llRfSh2t7+rqqBMMp3FZ8rcCEzupCPLCSh5/Huv5jVmeRxurzZ2X/oPhxpPigfCaZb6Lc9eUvW/Xgjw9pwd/7VlnRpPrHMVwqLo5DKUbF21F1R+OZfpSavAn9CH8+Jsc98pNC/cuir9BTPz+PFEfvrzFDqo/N639bnt04BoGPocQLHnQ9yzGoPUqDofhJ1x/ocyypHQb5lQJK9mhs+776nr/Xir/vq0eZ6nxW9FYxN0CDlAde4MZnd5d2DuD5o8ap34QNftZKiKdiupMrg8Y0/cgFq3A9YY+hfqLwTIH7s9JJrwvcK8zPb0pnL2ZUN2k3VLZZ/ZecHdrTURKT6VbFmgRWY3MIpQgr6kn9iId/1PAbxUWdgmfmN8W/DzHmcUBZcVAxCn0swWN7PMaC70kkWe4qVFvdSZ2kG5XUnfQotRLJYvpRw2PossT8q48bCtwL6DuLlOcKlVOAB6XDzpr7nzkaQ312Y5b+w4R4W/IeaMh5vZE60YQAw+AQwuoeYziHMBQf1G9DacZ2OLd9z3j31HfOSgZFhVIKjvJD/15Fy/31/ML31Lmh3dS2ScHET2ndSY7VLtDYs/Sr3V8tfXL6UeM9dVli/l+r8zsKQztkdtIcfQZ7kuKpnZ+0G+qL+3WO/jNWCQzaPEfTJ+8YBDjG7aOfxI9eymDwokJnZGJQw8oho/rxqkhOjj1yFRZFPahN6piIlRPG6Viywvc9hU2QcFlCd1aV7cT1YzHshAsLAP73wANdDzlknr39PYGZIbuxKf3vk++iXegBDNncOM+UlzIsod6aMQ9Pw1FIMJid16+GGFQfTCr94JWRm4fYrBzeKH6uuEHKQ8log5LRBrqvfFYjtKct8T9Kh3+MQRv+exL5KB7KiqGHflCIjuVKv/7B2Kd2AF9rXV3BX6eu+IFf+CBmSx1wU3KYfPtwrAPfDxSQ6b5vjRwxZwMv8J+z2/FdLRoV0Z0Yh5Lf9OMKmDH6EYNydp7mAt3kuQKnADG9bI/B/RCOWdtljbT6G2M3tqb/DVnqiYuQ8OC58zddO2OiGMPxn2KMHaACMAYzCa7KABSwxtB7tXPGc0wjtWEceA4UDNx3ISM3Vl/IjBHHUHOW51bbsfYqAxceTDXqKO0clGL4ADMeOFWkNA6WfI5ucF5g6ZyOYgzJO8UvFYIcyuoQW+FBpfX9wedB19wUgbhRmcfflT8ac98JsgdZ6RsngzPA0Iaz1JiOIPdY3fBDlo3EZzHdaXOmvk0/pulHG8pS6bBwHLxvGFDzd6koazGsutl2Q31uRv+Rs48uaoUYMXYHzimQD0AKb3XJK1uAi9HsXe2pfpgEVV2cfmmNx4RCGL6phLEfUgqM8x2NzcOsMCbjtRXvE5kRwvDjfCrDrGucI1UxWH0YY3RS/dI+EP1AvHLXcDRKM49/Mp5C31iug/of1WenwE4nEPLT3wOfgcOrl/tsF8H3nLmp9zX3egndCTyZfgQkrnbUufrxtVXZK8aHYouVq5L5f3Psxpb0v1fiC0p1s2PQMHKdIxPlY1RQdOim6tQN4VXuCn/FBytVfjxSN8YYJQjj2SHVJf+eAsYy9aonBrVXgdSWnQWT/JlrxexaWOnH8FNRlHC4MYryXquIjDlzgKJGH6CJT3hl58CugVdZ34f+dT2EEfXrjgYccWqOJ8WVswt91mLa9cqn9sxncGq1pm6HdVB5bOExasepPkbpjurD8xPPTHBsybNq1d+lfmTiMEc/Dhpjifl3dkF91/W2rluNa+RUCPqcO/eDdkODDOo/jGjsJTBoyDiUcI5BlTCW1U6h1QiAMMCsAp+1ylZJihcmCwAxfg1S3hdlhJu3KlP+bSWCwXHtq0J/oTrk07Z+PNWuFtIYRIwFqxGwce88q4+Y4VLx0WlIUdnWQtRzjiETo8qJuNZXupOrF4NOLYWf8sEZPRzjfD2LXyO1H607al0dndGT+uDGZYfdd2S5O/2QvIM4CJs5+gG20ftH+XPm/5N4j97zym+QHwcDHuTIkTnXbgzqP8wshAH9Yhuz6EKVEYrQZwwp5wz5cqhX1Ul5u76mo1Z86qgxWa2OmVS3Qqjniy8mmnFQsBQ5x6E6lZPUdbV6CI2Uh/P4pJgjNvfLSF1zxMIKnBCUKjQpGQ8qucbHCUJ9cxqcasWr53sIo6ue9Vf1wQvdqHZYyiNNWUxXKIP/NWm07ojZh5LnhULQCxYC5PGMqMIvCKW8XerHBBzG6keAsHTMvDFfLCyGdg18TaDS50yZs+yG+lpV/yVX9g9/cQwcN4jnXqJOHaxkZdWLevxkg5EF6p9JaDzUbXWBQersGFp1UsnGOaHGwonUVw+hHYaXMarVq+ryC0ra31HoGAvllSL4GZLPKarqBUPWGFt8wj+y4QSidRoNfEJ1UWzetsK40o65RjfqN9sTpZ2zVFmbcNj1uu3yRdPia6ruIGfq6C/G8171oxcH4TtXP2JYlshzjlwdofPVPdvuWPzHTkV6ZfZ95NqNVfW/LW9f+rov5AFkWEFW9QUUq2Amm3BQGoVfjTQ+BpGVW98DU7ci9XXH8lpNnNojMyuMlJHHADpcaoOQbhvadh2qM85Ugh+eD/TRfRW2DbarL54xjq8V6OdHl5n/h7YEdkvI5bAOzX3fb0I6ErvnEZH8xbPE22TdUdvGZxPELMaCt/RSuoE8u9OPDBzm6ge4FSfxzWKEZ18saMJr1tU4yrur4JyC4suqQBdK58x9rt1YTf/rMuVcX1BJwvMwFtAAJ5w9c+MDCuHgyz4qr+9mUHF50tjwwkoVxYMXdgyN7+4oj9UAN6yro9ht/ZWfWr2qSoeoy1GQUx7FqdUFyoNzeqg6ihzd1F8wxGDgQHn+wO7hoDS8P1DAOJFPOYab/GcK5MM7hEFhO5t6iPuX6tCuQ74/sKL/g0/X6zEOxpwHp1k7wHpjXeP0kDvMA3rDR8fCmTsOo69f2gU5dbk8iZ+iuqP+wBZd4yZP0e71I4HDXP1I4Tk7H71U+EYdcY9yxMocQdwP2IpBO6E6qbnPtRtH138n4YQ/u/mIXq7smlxuapSj78Fhbner1BPvOG9et1zESS+BkfrEKeFEJuHueZr98HnOhHkZcIJgHwzLnC4XaSveltYP5vKoOKw9/3Pn3refrP8oyjExuFhEM63TpRFgZ8fK+6ir75lCsaqC76mEIV7NGPsbmx/wuSNXnz4oXu2ZSQ+Qi+mHl3sNHFab/0Iyz9V/pnsRDCQfOyF2/5U9sR1Dz9215SJNItthJrO4YVKfRXdV6s+tMBV3nmNtGePAW+Bf6fpxAzcRBnI1ZxX4i8Xiq7h+nCIOMWzG5JWQOfSheJP6L77YYfIiSnV8fTY7Bj853Myc83EkwaqKj++lniOoeNPEJPKsoZiyLYgRxw7VamTTqMaZ45nWDYXG23DCa8syFdcPL/+p4RCf0fzcEnO/Wf2XDvNjzw6d3Y6hg8AJZ3hD/rPizTo3r3g8+Cu+sznhqTsK66egH0cBYsVBtqz/4o0jJF6QYYHDSzPVjsEcw4pKY0MbAoaAIbAWAnIEvITDcSi7ooZjuL4WUzauIWAIGAKGwDoIyBlwhMQRV5TMMURhsUxDwBAwBPaJgJwCL4NEfwAbJDbHEJCw2BAwBAyB80CAtxn7vh5xMMdwHopgUhoChoAhcJBDcB/7HILCHMMQQlZuCBgChsAOEPBHSDcUD74haG8l7WDCTQRDwBAwBIYQkENIfauJH7TyQ02+dcXbSY/NMQyhaeWGgCFgCOwYATmCzxKP/zpY/Y7BjpJ2POEmmiFgCBgCGQjwa3ZCRbZjqKCwC0PAEDAEzgcB7RD4HcMtBY6YID6Z/0b5z/8Hl0I60ddK3dUAAAAASUVORK5CYII=", + "text/latex": [ + "$\\displaystyle \\alpha = A_{1} t \\sin{\\left(B_{1} x \\right)} + A_{2} \\sin{\\left(B_{2} x + C_{2} t \\right)} + \\frac{1}{4}$" + ], + "text/plain": [ + "α = A₁⋅t⋅sin(B₁⋅x) + A₂⋅sin(B₂⋅x + C₂⋅t) + 1/4" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "alpha = Rational(1, 4) + A1 * t * sin(B1 * X) + A2 * sin(B2 * X + C2 * t)\n", + "Eq(symbols(\"alpha\"), alpha.subs({X: x, Y: y}))" + ] + }, + { + "cell_type": "code", + "execution_count": 8, "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAAtCAYAAAAk9oUAAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQt0lEQVR4Ae2d67EdtRKFt10nAAMRwM3AmAgwGRg7ApsMoPyPfy7IABwBNhlgR2DjDC43ApuTge/6dNQqzUMzmvfe0F01RzOSRmotLXVLmtlzbn38+PE0Jj/++OMd5bmv4wsdr3T9buweT3cEHAFHwBFwBPZC4PZYRXJcOLGnOnBgv+q4p7iXCl0cAUfAEXAEHIGzQODW0MosrsieKvwh11bX3+v6K4Xf5vF+7gg4Ao6AI+AIHIHA2MqMFdmT6NRy/VihPVD83TzSzx0BR8ARcAQcgSMQGHNm/5VSL+S0rnPl4jVxbEG6OAKOgCPgCDgChyJwNVS7nBYrsJJ8UMJXpUSPdwQcAUfAEXAE9kJg0JmhhBzaEwVf6rDV2RvF/a5r3mz8S4eLI+AIOAKOgCNwKAKDzkxO6w9p95fC70xLnd/XgYND3Jnd4OB/HQFHwBFwBA5EoPjMLDqyk8LkyNBT168UWNyfB+ruVTsCjoAj4Ag4AgGB3pWZHBYrr/s62F7sE35Ejby4CfyvI+AIOAKOgCNwHAKllRm/K7uWUyt96YPnZb8r3Z6jHdcCr9kRcAQcAUfgX49Ax5nJQbHqwlmxndgRpbNiQ365CfyvI+AIOAKOgCNwLAIdZyZ1Po0qvSmoxqrtnZxar7Mr3OPRjoAj4Ag4Ao7AZgh0nJmcFG8o9m4fKu2B0liZ+WesNusSL9gRcAQcAUdgKgIdZxYLYPVlbyyGKDkyPl31XMeX0eGF+L3+qE62Pg+RI+s+pMFnXOm598W563fGXXtRqnk/H9tdffj3OjNl5MsfPyh8qeMXDl3j3D7XeemlkM1apzr5sPGR34H8IuqwWRvXLlj6ovM/6qcTZ8CDmm66CK78E/lR0zlr5LkQHq7R1HMuozPOBr+afw4tEXHY2uQL/Y0v9++tWyQwb3gOfeJrb7WK9UlPHNldhbeKmVZMUD2snKmTSdDqGKnMQR4onZ+TsP1tLyjxlZpc0I9PsKHfphMylc/k66y5Ih1X54fK5OUx+uEb2q/D5A+l/aqDPvhO4aFj2ZSaE0r3i+HhnPYtvUf47MYB1dUYZ2ftzCIwrxWWfu+2FPvO/bFOBvo3Om984UTXxH+tMB+onTKOjpB+DLjwP+d0vpczw4nwxRiMVmOLeikeKo8BUsUD5f1bed8qxKA2RHFgAjZslW/t0M6WK2r76vxQmRgW/svGMx1wII2RWB/9AUd+0vXqkx2Vu7lI74vj4eagZBUIn004EHEftcm924yZfkef/iQF9v4JAAPOZvHt9qMLOp27YDjCykREoC2bi+rh7dZPFK7qyKLiVTyIbcXg4FT7xLi0hY7t+hZxRW3JPxvXLnvp9Wr8kJ53dIA3mDLR+1lHcmQoqmtbJcPFS34L+hJ5SBcURX2zmGcqY2sOVNnk28VWnkfCQwG19yyOgc73KBsDEjiiLuiEwTxLkW4MOA5bVe6mq+ruYLYSSLU8gPRIyWBiTBH+tdGmsgJX6LfV+24DfryWnvd0jK12WRWz9Wq81OXFycXxsALhNXi2NQeqbPLZOjORnq2QI4iPQbSZZB8X0OlhX8LRccIMY/0+Goz3UR8z4EerN6v+iTyA9CfdU9pCZPWAwx3qX4pYS86KK8JlVX6oPCZNvJj1WOfgOiQ8r7zYz9+pfVPs0bnxcKhfFqXtxIEqm3y1qCUTb44Nt7sYWM8UVzI8EKI0wz7pPmYU9vV+/q/aYx2U+UgHYv+q5uZq4K/KAiwMHWVSBi9OMJOkjJ8V5oJO6Lb3ijHXoXTOyw22hWYTAfsRfOmeqniVi9ECJ4wWOCGPFM+MHMzAi/CVrsPvEOM9/JyDeAwZxg+jgNBnffiGxOzPIA+yfJyiX4cz0gN90YOQrTDDRpfjEu+fw7Vz48pq/BAm9CnPSNjFqJkcgLlt846D3pMj9gMcgoMmjXGe5bHV93+U8aXiEy90XuSyFdoTXjIPe5qzPEo4bsYBlT3ZJl8tb9J4CbHRYV9d54FUCjEs/1MYfrcW8+BEbGCwdTFEfh4kB8OtkEGFsYLYDFh7XmFlKaksyo9OGGEMLXvIELckDBKAHpWoR1XerDC+rhKcQRY3eqp7qAeMTcxgM5gXicqmr8C7gYuuDX/qoh95SJtE10xUiAcznCr9GyYHCjEofyrkG5+mq6I6MsaDcEMsDz1POocPJjhNBh0/MZmMayxkLtequWLKbhWq7WvzwyZNQ2M0NUf144BKE9eUr3Si++lD+MUqMIzr2CYcFQ6VcQOnmFQ1Xt5SPG9TBu4pHOSy7i3JJfOw1Kal8ZtxQP002SZfLW1N5f2QkFVYPjti//yt4u0H2jiSfLUD6dia6IjuY5bMDM2EgcL9j2MEhjNPj9GjAcZ6bMChEwNrVKSndfZo3hUyfNuqD0yQKl1vshb/MpDvqXwe9Fq5ZG4bst7+Uj6cFZOE5Ex0jvFRdJgY5P1OXC5FHuSZdI6xRnj1u+EcdQ03nlOfjqG6QgH5H+VfwrVqruR1bnS+Nj/AFKmaMN5kXfQXJ9VeBcINxDhJnr7JETaGiRP2B9tQw2Vla8gl87DRkBUv9uBAtU3e3JmJQBg8jGB7uw5MMToYSuQz5TFScg3p8mviTHj1OjdYvLqPcQz5FSajaTdUhhjEscFJvRD7bETtZSXCjwhz5wJ+SK+uykv8Qx0Y/8GfPiidVSuG+W/OFbICxGj09amSeoWJS5/06pdlHOJBli1s/fa+YCA90ZUVGqszuPPObtT5GA5LuDbKFdVPn5kjNrUIQ/8pvW9CNGn1rjIm8SNi8jQqY5Oh9nMx0y8fh7n+jXPK1HFNZGX56f6Y/64iGnxTPGM1jFedk46ub3Q0RGk2cWJbnJ2bOVymvUH/RuHdC4zvJB5KHzg4hvcSHp5UxxY824MD1Tb5So382O2P+TEqr/27JgxmWpG1SoYcGGFWUc9aacVL5U/GKGaijur7+wpWmRCKwZBv1fVlrSV1372rx0W9TwoZRA1RHH1rxiilKZ6BbwaUdtcIDo8Bx2yM/gq/F1JZfcZWybvLGOntGQr5An9qcFCeJVwb5UoJP8WDM2OjYcCnoqr7Q/8qrOaH6kjbqtSnezGE7K7kW9Y4hCru6H7agtMzLGvKV/YkxuH3KaZ7Ynmuu0kpBt4jW3J5Mg+lzygewtCwCw3Qn0k2T/f3jtPYN3N5tikHpBv8qrbJOLO28zGwFodRGRQachCk86yrTcIqoHQf5KGMtKKK9Z56ylS2olAO9yTHq/M0m8zuoi50GxXdX5oNDd07adatghoDoVUwBgSD2hDpxcCgHozMqCgfhPqgkO0ajpPO2X4LnzrTOeVtJaM8UP1mpIZ4ZsY86Rr1noLDVK5Vc2Ur8FTuZH7onifCJn9pgskLcTx3MvwYJ8RhCMdWZ3zBJ43PyvKVLYmVnzvTlBhPLA+Yl4Rtyrlc3oyHUrYG79QmtWEqD9O9K59szYFJNvlq5cY1ihPoLLcbcT0X5MmJblkgp822LO6kvJCVvXEGKWAy42gv658qLRhdpQXhPh3Xdt0TYuxsoJ6Utz2btFtwDjZwLK43VBm9s6HezDMiVT6dzYy5JLS3g2Ep80A8zoJy0ipBdfOVBzC7pyPhpvO1pZcHrUoC6RWXJiJ5uvREf/oTx9WbJ89v58o7i2t2v8JqrmT3rHYq/efyA96+HVEkODjlAdfEi/Y90qFv16Wm/FSUymB8wzG41hGlPdDBVjJ8h5MNe6J44wd2Yy6Xt+ThIB7SfykPO5itFLE1BybZ5NsrNWqoGB64o1RD1EHM7OkkjpOuIVkukJfnHG2BmBysFLj3Q55BcaQ19s1jPp73DBn+ZHhifmacfUY6PJ/L6zziXLrhXJg9D73QELCJeZeqyQQh9FVWENdt59DOQ3awnSslHuTlPdJFe0IT0qUzfHitg3K+DpH1fyZzrVX0YVyJfT6LH7q38TkqtQljG94YtPYpD8adZ9NMKhnLDVHcfR3BkSm8zhN1PVp+nj+eUxdjslFXrMPGKf3LD5vbtgQ9+DKJcbWWy7HqEGzGQ+k1hsdSHubtWO1cem/NgUk2+aqvZVISgwSAGEwe/htZ+rIPxuleXjCA8BDKnltQPh3IcYpp7Ifn9fyma2ZSbYGQGHD0417K59Vb29fHybVnZmboWJ1xNAYX5UhYybFtFgaLwtJsk3oZ3IcI+qticLH246Ab34tUHmbL6BjyKAxbRopvrFYVXyvgBeZPVIbd85lO6FeMHMaD52lh5qxr9HusA/4QTzp9QPwzHcQbhhgWtqEwVn1S4sEplgcewXjF67wM6vmgI73OnSdWnE/mWqvM3bkiDFblh8oDW/iEY26I0lgNfaJIeMAWLzxBwBzOjfJNeYrlh5LiH+WDZ5/rkroY62ZLGLMY1ZNCVt7oyUsepgsc4NocGfFFLiutJLvwUHr24bGUh6U2LY6XvltyYJJN7nxoWMoxAFlJQRhIyYNGfrdRMjZK3kZUJ4TlleLcyc2uTOUwKHHORvRJZek+BgaOc2jvflKZR2aOeGAczro90m9VHrQx3wIHlbmIK1GnuQ/m202cdR3bgB1gDM4aM0MVb13+UN1z0qTv1jyEM5vh3ddmtQmbeBjPlnIg3h9s8lXeQCXc0TWOK59RsXpiNcNMa2+HxmqOGbzN4nN155y3H0RPLQNc0MllXwTW5sEe2i/lCs5jdQdS23CNdVvRhEcE8fqkMKyCassp5du6/FK9C+M34+GBeBzGs5XanMZZY2Wmwuksttn4J5xpIOn8juL41xp8zaGxSoppbHXhBFchuspKojLZvuj8EDZlqDyJerJNVto+HCxJ94VZk8LO87/BG884UW1hVnb2KzMglK6r8KCvO9bGQeVdNFdMf2GVT2qZUObbd31QVsVtXX6VEjMzSffVeXjJeMyE8bRGm60MhcEmN1ZmUoxl9AslJkeGslzHOLYgG85M18QxeNmS3EJYDfKsZakTme3IYqNY/q+1QoxFHhOoL+kv2hL6TtdMYvhA8SxHv1Mr1uJBUndDHC6dK0xOmcASJhFea/F/6/KTzhucrM5D6XjJeMyFeI02N8ZZY2U2pJWIjKPjAWtjq1HXFHhf4WbPXVQ2xpfXbw8xtqr3e9XPg87VV54q16USgaN5UKOmc6UGpcvOcwk8vGyEx7XvG2cdZ6ZMbDPyRpCtzt4oDkPO1yR4eaKxQtI1To70fEtCUS6OgCPgCDgCjsA+CDScmRwS+8G8Apu2E3Ru24iswHgZhOdXxJGHrQjOeXUU54fjO2T1pLpdHAFHwBFwBP6lCCRnJieEIzspbKy8Yhz7m/z+AUeWfqSrc14g4C3HzT6JRf0ujoAj4Ag4Ao7AEAJXJMoZsbXICqvzw0jSJazAkBc3QfqL42u/EJIS/cQRcAQcAUfAEdgDgduxEp538cZiyTHxAgbPxew5mulmW4x27aEj4Ag4Ao6AI7A7ArfloFh14azscy8NJZSOw0J4ZpYkuy9sT6YEP3EEHAFHwBFwBHZGgJXZp7HON4W6WbXxSn7b2QUnl8dHB1coxqMdAUfAEXAEHIFtEGBlxm+n2tuHoTal8YIHTqvx27KQePMj5rQtGfOywnNxBBwBR8ARcAR2RYCVGcLqK72OT4ScE28vPtfBJ6z6fizMii7EKz1sVSpMzk1pLo6AI+AIOAKOwC4I5K/mswp7pONDVnPxW2xyXKzCeI7Gp6ZOuk6v7HPt4gg4Ao6AI+AI7IXA/wEn5lxyDlNpLgAAAABJRU5ErkJggg==", + "text/latex": [ + "$\\displaystyle \\frac{\\partial}{\\partial t} \\alpha{\\left(x,t \\right)} = A_{1} \\sin{\\left(B_{1} x \\right)} + A_{2} C_{2} \\cos{\\left(B_{2} x + C_{2} t \\right)}$" + ], + "text/plain": [ + "∂ \n", + "──(α(x, t)) = A₁⋅sin(B₁⋅x) + A₂⋅C₂⋅cos(B₂⋅x + C₂⋅t)\n", + "∂t " + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "### Python\n", + "alpha_t = simplify(expand(diff(alpha, t)))\n", "\n", - "Copy the first cell under Source Term directly into your program.\n", - "For a performance boost, convert the expressions into lambda functions:\n", + "Eq(diff(α, t).subs({X: x, Y: y}), alpha_t.subs({X: x, Y: y}))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdsAAAAtCAYAAADmxsGEAAAACXBIWXMAAA7EAAAOxAGVKw4bAAARA0lEQVR4Ae2d7bXbNhKGZZ9bgDepINkO/FFBnA4cuwI7HWyO//mfT9JB4gripIPYFdhxB5utIM7twPs+EIYGKYAEJFKUrmbOoQDiYzB4McAAIEjd+vTp06aWXrx4cUdpH+r6Wtcb3X+ozevpHAFHwBFwBByBS0Xgdm3FZVgxss91YWB/0XVfYb/JdXIEHAFHwBFwBByBEQRu1axs44r2udwfUl66/4/uH8j9Lg13vyPgCDgCjoAj4Ah8RqB2ZcuK9lk0up9zb1e4jxR+Nw10vyPgCDgCjoAj4Ah8RqDW2P5XWV7LqF5/zrrZxHvC2GJ2cgQcAUfAEXAEHIEMAleZsJ0gGVWe0ZbooyIelCI93BFwBBwBR8ARuHQEqowtIMngPpNzT5etbt8p7HfdczL5L11OjoAj4Ag4Ao6AI5BBoMrYyqj+obx/yf3eeMj/UBcGGHJju8XBfx0BR8ARcAQcgR0EJp/ZRkO7kdsZWrjo/o0cC/tzh7MHOAKOgCPgCDgCjkBAYHRlK4PKyvWhLraPc8RHLqDXW8d/HQFHwBFwBBwBR2CIwNTKlvdqr2V0S1+K4nnt74q357hD/n7vCDgCjoAj4AhcPAJFYysDyqoVY8p28Q4pnhUv9PPW8V9HwBFwBBwBR8ARyCFQNLZK/EXM8C6XUWGsej/I6GaNcSGPBzsCjoAj4Ag4AheHQNHYyohywji7Pay4R4pjZeufabw4lfEKOwKOgCPgCLQiUDS2kRGrVztxHIJkaPk04ytd96JBDuHH/lHZbHGvRmuXv1rFveAiAqeuE6cuXxFYj2hGwNu6GbJZM+TwHzW2ysCXo36Q+5uun7l0j/H9Sv7SoalZhc4xU9n8AcLa32P+OsqRE/EkwyQvMvtrWgu0zono5FTNzkJnXU+nmnE8/kx0cbwS5x+709eq/vXnlOotRWILm38a6v0D0RoyRqXmtPbY5yzXEC1bpuTE0N6VeyuboCFQPHgtjMcIdlCOr4mlxM4Dn/JksrbaxCwVaCm/6jeqk6eElWRhonrSOisZXU/3VNZz0sU9q3hQNuHDwV/Grm91pY9J/1DcL7oYt76Xe7B9EY9eXzsrYyvhAeqt3NJ7v4qen2K5DADfyt/7WpbuCf9Gbtpw8wtxIEfJh0EI/z8s/8HG1sQRr3/kfy8X5e2RwiiPcnnk0GRwlZ62zmLeK2TlmyhnlU4q7SJYtUIgOU5WZyXbWelpK/ZLpj9HXVwSjyFv4YPx4x/sXurCsHZjdtQ7xjAWDz/qvnoBpbTFsUpxXV8b3UZWoadGP0qgNV41ogFspTbEBHmQ69QJRQqrTykAdTmYIh8Ujc955sjaqvfcP5cwEzaGeSb5akFVOrkwVq2VP1hnVZ/0c62t5Y+lPzc9HavLsePOURcnMTpU15T/ji7GKMYhFkY/6eoMLQLo3nbmGBtb37AZG6u6vnabgs6IHguU6hnHjPViAODb0L0Ggn+UB7kwOidJko1OyGWr8rlkRcmgknKaUecvGlupiHkro4XT1+rkklg1VXEmnUWH5tKjIP+Z6mkT9gsnPjtdrMTjUF17q3Lu65raYWMnjkcsNk5Wihe2pCftw9kYWwHA9lIrCLVgTaVjoLSZTy4tcj3ORawdJtwweH9HBfo7ymNG8FDxMIgb8S5tETOTvNY1hh0scjSFeS7PUcNU7xadXBKrfep9UjorLM9VT/fBfvY8Z66Ls+NhDIULiwwO0z6Vn7FojDhjss+nh6fGqtDXrsZKXjouAmHF0NleKqw0cDNYlVZQG+Vj9mP/QsT/6z7VBc8nuiD7S8Dt3cSv+AEgxgK+8OFgETMf+PwkNyXkQr41Vt2pHDk/B5RsG9cmK/bBklz6ljAw2mkTlQdmvB6Gy7aNlavbMildC+aBUSyLDnWdcO61dZLGVtj/VlpO2Heyy0+HpHz4IDf0ROFj5wNGdXLLovudFSvjGuu2j96fms6ejZ4a9kM3toXr4haYg8fgIb6t92oPxm2e07LqrJnwM07Zo6/R4sSvZawKfe1qlONCkRGEsIcufxjw5DLA/U9ueH83psHAGUhsA4wBwUPtYFTkAjCDPYMundj2zY2XosZJeZDrjVxWLzyfYmAtEYM44E9SlKUqbcKML3U1f0BEeSgHnI3M6GFsDiLxxjgFoyQ/eBvRyVByXhVrklnpWzDfKD3lcACBWWtoW7nUGUNKBwM35GSS1DvcpnBOH6JfPL+hHuhPr411b5MURWdpSidDJvGZHatEmn31vlpnk7IW8Qqfs9LTHAiqg+vitg/NNgbncG4Ms/47Zjc6lmrDa92UFntdOjxK2zJWhb521eNwvBsGSFaxwdBSrPzslb+X1z6kgZFLV4oMiCzzd0j5mNkzozQCNPI/jQFfyE3jY3CVwwA81QDIRWebJMlqCjCZdoYE3w3KAxeoStZt0uIvAyTEMXkz4iFA92D/Su5GV9qGIb7ipwZz2GBEh7NW9ASyupKGP8voyagw9OxPhaOD6Md9+TlIYfkUNDq5I56ysjpJZEKLYCVZD9H7ap1N6rGU91z1NMXDdbE/xtKP5hqDU5xb/JQPVS+ytsmbfmvGqtDXjm5sNUAwy2BQG27FUkMGxPt4RF8OBj4GxHQgDIniD6+epIMpW3+sakJ6uU0rrJSx/AyUU41F2Qy8J0OqM6tNXqxOZ3VgCGVlVVrCH+vCgI5tn8IDJcseJlBejBsrXFa3tM0HMkDy15QxiXnkc1cse3qkcNoqtJf8xDOxeKerR4pDPwhjq5jdDzrEP3IxvuwGUIceb4UNaUwn07RLYXWI3lfprDBAf2iPIQVdUnxu8li9E6P8TXqq9OjP8yiMTRrHnsc1Y99aRkx/Y3WxEo9DdDE0p8qZW9dMR1PbMNTj7p566romAL+cGj2bHKvEJ/S1KzH9BPMlSWWk73UymHcr2kG5VBQDwSr05SCueKv03WAeE1FGdf4SY/EFcDp0uhWbS1476Obyzh4W5d7IZaDpkcJobxukujiFM1jYoEq9p2hKyez5KOlC+9SUoTS1mFsd7NBXTl5Lg16ViHpDTC7oXMyG0b/wrp3kyRkTRTfRUlgdovdVOluqv8LBib46NSEpAqW8Qc/kVuupmHXb5jBWXgZodspKj0aasRev1jJMz26kLtbgoXY4RBdVRGjLbF8T73117aPY1oxl6BFlYBStHpM6oDy1Y1Xoaxjb1BCqrOUoCoeAY8aLeJ61Xg8kqQJO+ehc8OhWo7HcTYbnoIidW3iRr5scyN/NfpLUlId8k6T8pdnbWN7qlUJk0lOUAWMUisbvkeRCySgHpRslpTEDNdaONoCa8m4qy6jF3GarpUGWOlga2qdEbEMzWH6Uy9YyF7KyRRs+USp/VwfiEqLNx3jDZ0msOlFUTqveV+tsV8j8nmY9lQjPVNf0cBsTI8J4/t5rJ8KiyE162lJG5G96diN1sRUP4d6qixHG2R3GbXSDSaG1UakQvkrY2YyYb0rPaseq0NeuSiUvEa7KsO04xZo0aaUtPWDZDNLCNkpLRXheQscFXGZH8EjBfa77MIgqLhD5dF3bfcHFYHQdWOmHsx/LhvFKy7PwHVc8srO3nYR7Bog/CsBMv0TUeQfHUuJCeFAyxXWTkDSdZGCQAyuMdzZNmn7gr8JcfGlj2ub+IH+4VdwjXWwFU1949nRK4VYHdAd5waRbpSmeL8yQD/6dDsifUlYn0wTyWzlZHFTGXlgp3156n8hWrbNJntm8kn9fPaX/vK8UZF/sW8rYqC43XRdH8ZhBFyubszlZmIgpF2NR17eHXCR/bid1tM6RR9VYpbShr90eFnyEew7MIGSPVGFWEgwgXBvdMwilxIDHc8Ah0aG4WJmQ92OaQGHE9Z7ZxXQ8nxszSrDpBqSYhxlSbuANz4jTctfwSzYMBrOxsUNJAZ+Ydl8xnyjjcEITeIkveL/VBU7fhMC2n1rM4cqzeNoE3elI93Qeaydk4GX/oT6RhpPIZgSZkAXd6xhtddHik+DOSxk5newSyLMUVs16nwol/2o6K5z31lPl7X1mT/VgUAwnzwf143Yv7BvLsGJvrC5W4HGoLhqGs7qSm8kw7cJCrDdGUJDCHuoKhlbuNWFGuq/Rs9qxKvS1K2OeuiqIQQcA6RS8/mIDV5psL794cfiGylNJe65HeVSOaxPjeP6Rlvur7lmFDInBEOOCvOSFP6912LMcjPBwVWOG4o7iuHpAwycSq2G2EkNDyS3NjiibTr8KUQcVDDaGAZOI3vealYbZHTKGNHLDFonCqWMVKS1lUFYwXPE+zYu+fNTVvYqTRlb6azHfqHwG2a/EF32ivU2faLOw0yCX1TXKziEoa2fk5N4MKeHo0DOFyQn0pX7hO7ZjUdLJjfItjVWz3sd6mXN0nRUms+qp+KGH6DXt29Gc2JfK6AqLHqW7CF0s4HGoLg7hnO1e8rK79S8xpC/zOMHGAMYpxsDJ8U9psnqm/LVjVehrO39EIMZEsPJk8EIgDhvxjiIzhFVJMjCY8ppAaoT3lkl86KhMJqwBmnkpLwM3xn3seU0z37UyRExQzMXqc4wyjoWf6jKrTg7lXgIr8TxYZ6NcBx2QGta15T7WgTGK8WDv/jtW5jHKGCu/NU7yLq2L6M2imOfqrHoxTq+ia4fqQMwf7MNVWjlFMPPEsKbWntUmqz9mAWsbXFbDrM7mWkUOH4qncNT6wQq5nC4Tgbl18hgozqGzGLhFjNwUABqHbGciPI6K9xu5Y7sQU2x78ccoo1fgPDeL6eLKeKyiazPVuetrtwdt/Fz3bKVhdFNii40DJ2H7MI04pl/lIwczHDrbQRTryFb13hTlQB7kcrpABGLbz6KTx4BvLp0VH3aEjq73UX5WVzwq4AQyYxIDGrtws9AxyphF0AETyT3b+JiyXhsPlX90XZujzpFHZx96K1sBzDbEayXqzVi5j2FsMc+yhZs2ZqOf1TXPw3YOWTXyYVJRegZby4pOP9cqu7bMRdIJCyYw9kwXBWGWzB8YHIpRJ+8xyugKO65nLp3spF4Qq3PXWc4jsBjA7Uh4zdkPj1FGJ/vMntl1UfKdMx77wjtHnXt9beeZbUkyKTOGmMMma28lbyQDhoGV9myGoFTvUrjK5ss3PHyfbeuqVJaHnz4Cp6CTUyi5zk4hdDPiz0EXbwbS5Vrk+tqOsVUiTt5yus9Wt+8UhlHhy0Ms5w9dUZYl9BhHwBFwBBwBR+AGItDbRpYh5Wg0R9i7LRn5eRfJ3lHqreIUznaOxT2Q/6kuVp1PdEHBUG+9/usIOAKOgCPgCFwmAret2jKc4ZNmcjtDS5zueYfKwnrPSRTOKyJ8GIDt3He6XunCOHNgAUPsp3QFgpMj4Ag4Ao7AZSMQVrYyjqxOOfzUezk8gQbDCb3eOsEIkyc1pmw78z4Uq1voC11pfAj0H0fAEXAEHAFH4NIQsG1kVqKcOC6dNGZr2L4zaxi9V/p0WxlDzQEqjO5G7uoHqUxQdx0BR8ARcAQcgTURuC2jyKoVY2qfruvJo3hWvBDHmDtS+NAw86WpX7sE7nEEHAFHwBFwBByBgADPbNnuhXjmmiNWvaxYs8aYDNEgY7S7bxArjK9OEebkCDgCjoAj4AhcNAKsbNkKDlu/QyQUxzNYVra9LeFoSPneo616OUDFNnS6rcy/qGT5Dsvxe0fAEXAEHAFH4CYjwMoWYvVqJ45DgAwln0HjdPG9gRElHiPL9VFxrF57n0pTGHGllbKinBwBR8ARcAQcgctBoPuohQwkq9gnulLDmf4tWYdKNLCcNOarUhvd8/oP+fngBa8HYYS7LWXdOzkCjoAj4Ag4AheLwP8B8rYZ1bJfnOwAAAAASUVORK5CYII=", + "text/latex": [ + "$\\displaystyle \\frac{\\partial}{\\partial x} \\alpha{\\left(x,t \\right)} = A_{1} B_{1} t \\cos{\\left(B_{1} x \\right)} + A_{2} B_{2} \\cos{\\left(B_{2} x + C_{2} t \\right)}$" + ], + "text/plain": [ + "∂ \n", + "──(α(x, t)) = A₁⋅B₁⋅t⋅cos(B₁⋅x) + A₂⋅B₂⋅cos(B₂⋅x + C₂⋅t)\n", + "∂x " + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "alpha_x = simplify(expand(diff(alpha, X)))\n", "\n", - "```python\n", - "from sympy.utilities.lambdify import lambdify\n", + "Eq(diff(α, X).subs({X: x, Y: y}), alpha_x.subs({X: x, Y: y}))" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAAuCAYAAADTGpswAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAS6UlEQVR4Ae2d67XdtBLHN2edAnhUcKEDSCog6SCQCnLoIKx8y7cs6ABSAQkdABUk0MHlVkA4HeT+f9qSjmxLtmzL+5WZtbz1Hkl/jTQa+bE/ev/+/a4lPX/+/GPxe+Z5fu7dJ4q/bVmP8TIEDAFDwBAwBAyB3e56AxB+kNL+LvCV/yf5/9T1RYgz1xAwBAwBQ8AQMATaIPDRHIvcW9sPVDWW9u8K/9VvhuIw8R/K/Z00ueT9r66v5B/kJ4+RIWAIGAKGgCFgCCxD4Kq2mJQwCpwjc5Txz7ruKe613D5hjb/tR1rYEDAEDAFDwBAwBNojUGWRS2G7+95yv0+boPBThe/L/SaNT/1K+0HhR3LtaD0FxvyGgCFgCBgChkADBGotcizxGyljFHpKWOYo6S/TyOD38Y8U/irEmWsIGAKGgCFgCBgC7RCoVeTc434lxdx58tyHiePYvUNK49441jj3xjvlOhktYAgYAoaAIWAIGAKLEah6al2KGMu7RO+UcD9N9Er8e7kPifdh3L/TfOY3BAwBQ8AQMAQMgXUIVClyqpASvpHDEXmwrt8o7leFsbyjglYcYV45Q5GHI3cegOvcX1fYyBAwBAwBQ8AQMARWIlClyKWQf1M9f8tN3w9/oDDKHYqKXH7eGf/Yu3L2lJYNceYaAoaAIWAIGAKGwDoEJp9a90p8J9cdk6fVKQ6ljdX9nfxjx+9pMfMbAoaAIWAIGAKGQCMERi1yKWcs7ge6Sk+dY3lDr/aO/RoChoAhYAgYAobAIRGYemqd+9q3UuilL7JxP/xXpYf75odsu9VlCBgChoAhYAh88AgULXIpZ6xtp6hzKCkdSx3iwTYjQ+BkEPCyy7cPIGQYsj/u2eNgvyeAgMnoCQzCBTWhqMjVx099P98U+ou1/pcE0n1TvZDHog2BYyBgf9xzDNStzjkImIzOQcvyjiJQPFqXguZJ9OyRudL4WhsWefHTrKO1WqIhsC0CN5LRcGJETXyY6HPFhdcht63duBsC0wiYjE5jZDkqESgqcl8eqzu+ckacXwxfyssX29LXzkg+OVIbw9Hqwdt2zLoP3tnTqvAi/7jn1OXp1Nt3WiLq1tWz/HMpG+fjSlIO/1FFrgK8UsaHXV7r+olLYRbJ/8hfegDuuL1Malcb+VOXY1phWIG04YMh9Zc+81ri0Uj1/6wrPU1CZvkOwsnLbAk0tf3YslxqWhp/NvIuPI8qp6r/LGX0TOQwlclL9A/m2eR75OeKggSO43/+me2oX5Tzgs+T/x/Ee/bqp/u2gNyP1sqOePD6I7dvwjE5XxJMidMWPhHMZjOrpBXPRo6/293sm/+qg3bQb9rRfJzFc1SWlb4aJ7W9CaktbDhOXt7VziZy2gJ78Wgio+JjcthEipcxEf48IM5c5JsrqSHxm9LYuDE+fHNltU4Sj848u0hFrk4C6B9yS++/K7kt+TpZHB7K37nloDDxX8tNB7dtA06Am/qHwnH/US//akUeuiRe/8r/Vm7uo0TUR70o6o4yV5iJwynSN/JnsVc8spIdN8VXkXiw0eDrh0zWzq2oKgYjmXz7qmRZeRfhNFL9oiS146TlXe1rLqdLsVe5SRmtHQTxMjmsBatxPmGPYuVNmRe6OqctSkPeWLsYHx5yrN7sK29xfVJanGdXYnyJxMNNLOCHJAaJSYmF2CfaQpsunRBWZzVLyMBiNXk+CDOKMkdhnDsK1JfDQmZjhYXIcVSuTWPjlqtvECe+vLnxidxOGwYZl0VUybLv22ycljVpstSpy3tTOV2KvS9XI6OTgJNB/EwOq5DaZxJe6WfGZ5S8yyoeH+tibWLuY6z9qOv2Locbl3CSyPoz9y2vsfUpzrOrtMIL8n8rMKt3PY36zeLAfdjOIMLbt4U2sdBeJKlvKByucBrRqq8IMlSaAEE581e7jtQW4hBynuv4kkt+jrNym6ziuCl/NamOwbhXFx7PWCvLs3Ear3Z5qrBg7p2kvKttW8jpbOzVjjkyWj0YJofVUJGRNWrtOvWHeNzTNTgRVFxKnBxiUIT1MU0b8xfXJ/GK8+ziFLk6xzHGXLDGgKxNYzKHnVeuDG36Npdw7nHCnEXpHy+k//j+ENeCEOSdeHeOzRPG7IRRoin2HDkxHrjhuhGPnLKdGjexOA6pvXNkeQlOW3bs5ORdeG4lp0uwnyOjW47TJO8zl8PJ/i3NIFzYFGIk1HxsCiNiyafMp9YnN8+ul3biUOU8WKE6JuILxZUWdfIxqUrW205l2YHxQALE/6g/0QXfx7qg8Pes+9DIr3gBMooEnvDA+mPnBY8f5aZEm2jboU8K0jZs5ed4MBwrh03Up40qA+PBeKo+MOc1SFyOtEK9O/k/UVyRlD5n3ODHZKXMrS7qgx4rnl04486Y4/6usPu2gi9D+4hnAjPpUcwQcpeTEZeY/IzKcpIP72yceuWzQfWD/i6ZL6co71vJ6WzsheuojOYGw8uUyeHKdTuH7dw4jQXzmvvinMKmRkSJFetTuA1YyuPixY8xnqVXrkc5HjHRA+XuPcjvFnK5LCr/k8sCCoCAifJMgeSYYwyw+EUllWMgWGxZVJnk4Z5Dyk/JeVJ+2sXizQLN/RYW3hJx9MsATZJvR1XehBlf2XNKJInb3Ks6aSfjFCgo1C9CxFJXvFGgTnHKz1gFQhEy9hydz+6zylSPm/JSPzLTGVuF3cZFLv1FHrGwIinMZpN4xp1NDXLqNndy6defcvmfgoCXogY0JcuugOfXHCffmqXzpVreB73eIEIYbSKnG2MfkVA9Jod3c27Vuh1BXecJhsuYrok1aPxuFRgzQNO81euTCrl5dh1Ln56HhRHr2ylxmgcYut7Ky/1OgESB9i1cBJ5jjAGpLJYFllEgwIXHEx/Bgpum++hJh0V+apBoE8pnktTOICSTeU8gA0+Ep+0FUyjbV+VlfL7VxWsYU28VhM0MeTsKT2HG7aXcna6+DCipimrGDWV6T3XwUEvoG8z7Ezgrc8pHu9nkxQ2H/Gy6FO02dmNtB6sSX8oHWoST2jA6FkpfM1+q5T10YmO3Wk49Ls98e4Icl45PF2G/oK8XKYc1WK+UwwVQVxVh/YGqjL591v1vTZ99/pr1yc2zazF9n1ZyLL/aEV9Xkp9FkoWzfzxN81gYEWroM+VJF1fiUMb9OOIhXmFKFQKKhEXV5ZcbF1uXu/6HyTw1oNTLwnk0Uv+o/w9dc9rBAvhXrtGKZ2fM0+CpUgN/aFCH8mGJhoVvkO5KdX8Q5OwDIuKFNYtljlXOuGbb2GU3CE2Om/hy4sJk+Re/XE4fqDsnm4MKfASbzxxNYTAmyym/2Tip/TVjsWa+VMm72gEGzWQyBSX4VccsOVW5eAoBD5VHvjEscqdMs7GH51xSGy5SDoVDDdZr5DCMX1h3UujdWiVsU0MkpE+dcIayqT4JZQeu6kgNgZo+w2NyfVIeN89Q5FGBDmo/XgQWW7TEe824VRjlgeX8opc2GlSZ/mJPPbN49CtggBTHrj09Xu5nI1y7KOfKNolTW8Fuygquqsv3eyeXhaxDimNzGCyZmKZ48GeChN1sTCt4pgSZYyWIfHFsxZ8xYWyLVr/PUzNuYuMww0Kj3cgdE7H5O+Piu5Rm46T2T46Fz5O2ac58qZJ31dFMJtOGBr/4Iws7udVyquw3yv9aV1iDGHPiuD0S5Qy+oiXY0ybkCQrzpGTx73Ptf5m7FyWH6s8k1hnM58jhTuVzipp45jO6ZM6mfD8S+5MyJ1shouT6elC4QXZq+gzvmvXJzbPrUuXHilen6QDXmGIknfvat5l2vvPlM0l3USrLBIRPtKR93bsC37vCXR98KBMmPf509xVyUxdtmySVxwJwfCcz32WY2kHe5Wzj6+wqeywRWrdj7cVXB4XBlz7zmByExTlMELCnXMAOzEtUNW7ix2R6J5fbOVzUwZGz+2Sx/LFu0hrTpCyr/kU4zW2n6gGvOfOlWt7ntmVm/iVyysJfOkWJ1a/AvtMm8WG+lyx+V5/yXKocVmEdQF8gh6Foa5f1/oZx0TVlld9XnqhnVK6mz1Xrk3i5eXbdundr+anDHKVOsSFPCkyaH1DDLjfGKz8d5uliJhGDAJjwSQfhmcJusVaaI8rpug3hjIsyiYu58rLLg2eM82VQbGldPnroiAdtO1lS+xCyzsNdvcaC12AMenmmgk6QlSlukNICagMKDKzZwMQ88oN7jdVfO27UQ1/irl11YI1T/h516dqKsrLcq2wRTj0eg6D6t2i+JIyq5T0p09SrPiySU8a31xDmY+5b/UuxRwHUWvyhKRcph1NYK32tHAb8WrvulEZMWYPi2tCvRO0fnBxP9dnzqF2f3Dy76ldMGPB0PdL1VBcCtJo8T5QoF0LMxSDliIlERzqk/FhClHHlFM61jYWVe6d9YtJxYV1R/l2aQXGkvenFkY97o2NKywFJOc+XHVpucXf341P+5+hX31BqjF1/sUu747D1edP4Of7HytzfaLny4stYcV8VnL92kfN/ascNzmzwnMwl1RCOGwgf389DNPUsJfqXk+WU31Y4zZ4vaaPkP6q8e9lbLafiwxrDYp27JbUU+xqLrAenC160HBawXiuHORxXx6mtbLJ5pgp9hl7qkOIe6HJKXO5tJzEJKK0kX7Xrk5tng2+t0wDVgxLluIcFmfsRfOaSRi8mlecoMlqahMWMzuYeINkpHhCgcB+URdJ9w9aXBRw+QtLZDSkMMEzgDl+FKQ9Px49yupig9BVFjYIfWPmKC/XzKtFgQBSHYqMvWPs7hbMKTvHw4Z5tf/Gn2MmT2g1+9BH5gJwSVXzERH7wZIzTPDyk0z/lIB8ToDNGitspjjqoK/Dojwl4I5fIUz9N0XtSWrEOcii9dtzCfPjHs8b5TBd/hEDfkLdnukJ7Gd8nuuBPPO0AI+Jf6CI+YOTixSM7tzzvgSyr/E5pm+Ik/ozB7PlC2yCVP4q8+3a3ktMgIzzwyVg5kr8J9gk/cMZwGsyHJM+ly2EJ61VyGPDLucKbubn0HrljKR6hfbQ/yAjrE/N2dK1XerbPME7SkDXCo3qlo8iVmUax6+svvE8Vzzl/dsGhoilSWR6AYkPgOieXTjDZUZAohWYkfvAtPm09tyLxY8BZtMNAzWKhcvSVhb84UWcxPPPMHs+sIm/VtUPU0aqtY3zUj6ay3K9rC5zE8+zl3fchfkDGh3dyq26P9XEuhcXPGR5KzxoKpXKHjlc7N5ND8UZeNse6j5nqXa3I+zxrwy367Hk4vXLVqxgL4kYZUOgpsRtgx4jQLSWskLdLC88sxw6X+loRm5hFStw3gI0RbTIyBOYi0FqW59a/JP9Zy7tfIDll49SHJ9VZ9+gTllYz8vUwvietxH2HN5HDQ2FdGDTW9DXreoHteHTDPsd51rfIOetHqAZKUHH8ReILuZ2j7PEml1PFB8EYPU4ql55OEX+edh58SGS6ZDeH+LCpYXOzqN8qx26TBWFwz79b0+WHPBbhWJnFEUwHt0daIKG62G1vavW3aGcND/WliSyndW01Fp7vWcu7+sBax7zvkOKbvarrcTq4Fdrp0MyA2ryFHG6O9cxubp69hXx5+Ynz7DpttRKz5/A+D7vR+2n+pX7VwyLOQpt7gGQp2345bgNwf2GtAl2sxH2D2NkPNkY+7YNyNO4cS7KLNJqHQCtZjrVuOBZnL+/CZvZ30COwFR7xd5t7ZUWRh1NO1ohTnxtbyOGmWFcMx8GzNJKvzjzrWOT0SJU4q1zecOTwRnF8yYp73NwnXqUYVT4IcecBEupuTb4urP5F1vTa9qheni0Au6b31da265LL+zE/iNV/SBx9v44myzV9NXmvQcmtsWdrhZ6DHNaNwvnmys2zjiJXBo5OeF8yWpDy87RkUL6dr1kpzR07e0iw1sPTuo99nNsEeP9O+eFzVsdJoe3mGgKGgCFgCBgCp4jAVWiUlCxKHGUblbgP85R5iOu/T809SF7jwuJ9o+ulrgcKc0SEkuc+uCPFhc3Apg+QhPrMNQQMAUPAEDAEPgQE3D1yKVmO07G8S/esw4MfrwIovkxU1IrnKJ773ljl0Ke60nQ2AfDpbAbEJ2wSlGRkCBgChoAhYAgYAnMQCA+7YUHzFa3S+9xY09zrRVkHmvWPNCr7wT3UEIAy1xAwBAwBQ8AQ2AqBKylYrGQUdfYrNErHUod4Si6S4vtKny/A/RIzmMcQMAQMAUPAEDAENkfgSjVwBA5xjztHWOudP6boZ/LKng1B/GQmGwSufl4LGwKGgCFgCBgChkA7BLDIeTUqPTKP3JXGPW8sct4fjOSVNJ+GC9Y697n7f3DxTOlZvpGReQwBQ8AQMAQMAUNgFQJY5BBWd+ehMylhPlTAU+h86a3/HjQKnKv6n8SU18gQMAQMAUPAEDAEGiMQ3yOXssb6fqwr/Z4w73wPrGrFcWTOE+l8SH+ncNU/iZHXyBAwBAwBQ8AQMATaIfB/u/i5nfT0UREAAAAASUVORK5CYII=", + "text/latex": [ + "$\\displaystyle \\frac{\\partial^{2}}{\\partial x^{2}} \\alpha{\\left(x,t \\right)} = - 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)}$" + ], + "text/plain": [ + " 2 \n", + " ∂ 2 2 \n", + "───(α(x, t)) = - A₁⋅B₁ ⋅t⋅sin(B₁⋅x) - A₂⋅B₂ ⋅sin(B₂⋅x + C₂⋅t)\n", + " 2 \n", + "∂x " + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "alpha_xx = simplify(expand(diff(alpha, X, 2)))\n", "\n", - "apy = lambdify([x, y], alpha, modules='sympy')\n", - "epy = lambdify([x, y], eta, modules='sympy')\n", - "ipy = lambdify([x, y], eta0, modules='sympy')\n", - "Spy = lambdify([x, y], S, modules='sympy')\n", - "```\n", + "Eq(diff(α, X, X).subs({X: x, Y: y}), alpha_xx.subs({X: x, Y: y}))" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "scrolled": false, + "tags": [] + }, + "outputs": [], + "source": [ + "eta = simplify(η.subs({R[0]: x, R[1]: y, α: alpha.subs({R[0]: x, R[1]: y})}))" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "scrolled": false, + "tags": [] + }, + "outputs": [], + "source": [ + "eta0 = simplify(eta.subs({t: 0, α: alpha.subs({R[0]: x, R[1]: y, t: 0})}))" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "scrolled": false, + "tags": [] + }, + "outputs": [], + "source": [ + "source = S.subs(\n", + " {\n", + " α: alpha,\n", + " diff(α, t): alpha_t,\n", + " diff(α, X): alpha_x,\n", + " diff(α, X, 2): alpha_xx\n", + " }\n", + ").subs({X: x, Y: y})\n", + "\n", + "alpha = alpha.subs({X: x, Y: y})\n", + "eta = eta.subs({X: x, Y: y})\n", + "eta0 = eta0.subs({X: x, Y: y})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Code" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Python\n", "\n", - "> Note: Click \"Code Toggle\" at the top of the page to see the Python expressions." + "Substituting the expression for $\\alpha$ into $S(x, y, \\alpha)$, the Python code for the full source term $S(x, y, t)$ is" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "S = -(2*sqrt(kappa)*(A_1*B_1*t*cos(B_1*x) + A_2*B_2*cos(B_2*x + C_2*t))**2*tanh(sqrt(2)*(-A_1*t*sin(B_1*x) - A_2*sin(B_2*x + C_2*t) + y - 1/4)/(2*sqrt(kappa))) + sqrt(2)*kappa*(-A_1*B_1**2*t*sin(B_1*x) - A_2*B_2**2*sin(B_2*x + C_2*t)) - sqrt(2)*(A_1*sin(B_1*x) + A_2*C_2*cos(B_2*x + C_2*t)))*sech(sqrt(2)*(-A_1*t*sin(B_1*x) - A_2*sin(B_2*x + C_2*t) + y - 1/4)/(2*sqrt(kappa)))**2/(4*sqrt(kappa))\n" + ] + } + ], + "source": [ + "print(\"S =\", source)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "α = A_1*t*sin(B_1*x) + A_2*sin(B_2*x + C_2*t) + 1/4\n" + ] + } + ], + "source": [ + "print(\"α =\", alpha)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "η = tanh(sqrt(2)*(4*A_1*t*sin(B_1*x) + 4*A_2*sin(B_2*x + C_2*t) - 4*y + 1)/(8*sqrt(kappa)))/2 + 1/2\n" + ] + } + ], + "source": [ + "print(\"η =\", eta)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "η₀ = tanh(sqrt(2)*(4*A_2*sin(B_2*x) - 4*y + 1)/(8*sqrt(kappa)))/2 + 1/2\n" + ] + } + ], + "source": [ + "print(\"η₀ =\", eta0)" ] }, { @@ -510,8 +753,11 @@ }, { "cell_type": "code", - "execution_count": 6, - "metadata": {}, + "execution_count": 18, + "metadata": { + "scrolled": false, + "tags": [] + }, "outputs": [ { "name": "stdout", @@ -520,7 +766,7 @@ "manufactured.h:\n", "\n", "/******************************************************************************\n", - " * Code generated with sympy 1.2 *\n", + " * Code generated with SymPy 1.11.1 *\n", " * *\n", " * See http://www.sympy.org/ for more information. *\n", " * *\n", @@ -531,10 +777,10 @@ "#ifndef PFHUB__MANUFACTURED__H\n", "#define PFHUB__MANUFACTURED__H\n", "\n", - "double alpha(double A1, double A2, double B1, double B2, double C2, double t, double x);\n", - "double eta(double A1, double A2, double B1, double B2, double C2, double kappa, double t, double x, double y);\n", - "double eta0(double A1, double A2, double B1, double B2, double C2, double kappa, double t, double x, double y);\n", - "double S(double S);\n", + "double alpha(double A_1, double A_2, double B_1, double B_2, double C_2, double t, double x);\n", + "double eta(double A_1, double A_2, double B_1, double B_2, double C_2, double kappa, double t, double x, double y);\n", + "double eta0(double A_2, double B_2, double kappa, double x, double y);\n", + "double S(double A_1, double A_2, double B_1, double B_2, double C_2, double kappa, double t, double x, double y);\n", "\n", "#endif\n", "\n", @@ -543,7 +789,7 @@ "manufactured.c:\n", "\n", "/******************************************************************************\n", - " * Code generated with sympy 1.2 *\n", + " * Code generated with SymPy 1.11.1 *\n", " * *\n", " * See http://www.sympy.org/ for more information. *\n", " * *\n", @@ -552,34 +798,34 @@ "#include \"manufactured.h\"\n", "#include \n", "\n", - "double alpha(double A1, double A2, double B1, double B2, double C2, double t, double x) {\n", + "double alpha(double A_1, double A_2, double B_1, double B_2, double C_2, double t, double x) {\n", "\n", " double alpha_result;\n", - " alpha_result = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25;\n", + " alpha_result = A_1*t*sin(B_1*x) + A_2*sin(B_2*x + C_2*t) + 1.0/4.0;\n", " return alpha_result;\n", "\n", "}\n", "\n", - "double eta(double A1, double A2, double B1, double B2, double C2, double kappa, double t, double x, double y) {\n", + "double eta(double A_1, double A_2, double B_1, double B_2, double C_2, double kappa, double t, double x, double y) {\n", "\n", " double eta_result;\n", - " eta_result = -0.5*tanh((1.0/2.0)*M_SQRT2*(-A1*t*sin(B1*x) - A2*sin(B2*x + C2*t) + y - 0.25)/sqrt(kappa)) + 0.5;\n", + " eta_result = (1.0/2.0)*tanh((1.0/8.0)*M_SQRT2*(4*A_1*t*sin(B_1*x) + 4*A_2*sin(B_2*x + C_2*t) - 4*y + 1)/sqrt(kappa)) + 1.0/2.0;\n", " return eta_result;\n", "\n", "}\n", "\n", - "double eta0(double A1, double A2, double B1, double B2, double C2, double kappa, double t, double x, double y) {\n", + "double eta0(double A_2, double B_2, double kappa, double x, double y) {\n", "\n", " double eta0_result;\n", - " eta0_result = -0.5*tanh((1.0/2.0)*M_SQRT2*(-A1*t*sin(B1*x) - A2*sin(B2*x + C2*t) + y - 0.25)/sqrt(kappa)) + 0.5;\n", + " eta0_result = (1.0/2.0)*tanh((1.0/8.0)*M_SQRT2*(4*A_2*sin(B_2*x) - 4*y + 1)/sqrt(kappa)) + 1.0/2.0;\n", " return eta0_result;\n", "\n", "}\n", "\n", - "double S(double S) {\n", + "double S(double A_1, double A_2, double B_1, double B_2, double C_2, double kappa, double t, double x, double y) {\n", "\n", " double S_result;\n", - " S_result = S;\n", + " S_result = -1.0/4.0*(2*sqrt(kappa)*pow(A_1*B_1*t*cos(B_1*x) + A_2*B_2*cos(B_2*x + C_2*t), 2)*tanh((1.0/2.0)*M_SQRT2*(-A_1*t*sin(B_1*x) - A_2*sin(B_2*x + C_2*t) + y - 1.0/4.0)/sqrt(kappa)) + M_SQRT2*kappa*(-A_1*pow(B_1, 2)*t*sin(B_1*x) - A_2*pow(B_2, 2)*sin(B_2*x + C_2*t)) - M_SQRT2*(A_1*sin(B_1*x) + A_2*C_2*cos(B_2*x + C_2*t)))*pow(1.0/((1.0/2.0)*exp((1.0/2.0)*M_SQRT2*((1.0/2.0)*I*A_1*t*(exp(I*B_1*x) - exp(-I*B_1*x)) + (1.0/2.0)*I*A_2*(-exp(I*(-B_2*x - C_2*t)) + exp(I*(B_2*x + C_2*t))) + y - 1.0/4.0)/sqrt(kappa)) + (1.0/2.0)*exp(-1.0/2.0*M_SQRT2*((1.0/2.0)*I*A_1*t*(exp(I*B_1*x) - exp(-I*B_1*x)) + (1.0/2.0)*I*A_2*(-exp(I*(-B_2*x - C_2*t)) + exp(I*(B_2*x + C_2*t))) + y - 1.0/4.0)/sqrt(kappa))), 2)/sqrt(kappa);\n", " return S_result;\n", "\n", "}\n", @@ -591,8 +837,8 @@ "[(c_name, code), (h_name, header)] = \\\n", "codegen([(\"alpha\", alpha),\n", " (\"eta\", eta),\n", - " (\"eta0\", eta),\n", - " (\"S\", S)],\n", + " (\"eta0\", eta0),\n", + " (\"S\", source)],\n", " language=\"C\",\n", " prefix=\"manufactured\",\n", " project=\"PFHub\")\n", @@ -611,8 +857,11 @@ }, { "cell_type": "code", - "execution_count": 7, - "metadata": {}, + "execution_count": 19, + "metadata": { + "scrolled": false, + "tags": [] + }, "outputs": [ { "name": "stdout", @@ -621,66 +870,83 @@ "manufactured.f:\n", "\n", "!******************************************************************************\n", - "!* Code generated with sympy 1.2 *\n", + "!* Code generated with SymPy 1.11.1 *\n", "!* *\n", "!* See http://www.sympy.org/ for more information. *\n", "!* *\n", "!* This file is part of 'PFHub' *\n", "!******************************************************************************\n", "\n", - "REAL*8 function alpha(A1, A2, B1, B2, C2, t, x)\n", + "REAL*8 function alpha(A_1, A_2, B_1, B_2, C_2, t, x)\n", "implicit none\n", - "REAL*8, intent(in) :: A1\n", - "REAL*8, intent(in) :: A2\n", - "REAL*8, intent(in) :: B1\n", - "REAL*8, intent(in) :: B2\n", - "REAL*8, intent(in) :: C2\n", + "REAL*8, intent(in) :: A_1\n", + "REAL*8, intent(in) :: A_2\n", + "REAL*8, intent(in) :: B_1\n", + "REAL*8, intent(in) :: B_2\n", + "REAL*8, intent(in) :: C_2\n", "REAL*8, intent(in) :: t\n", "REAL*8, intent(in) :: x\n", "\n", - "alpha = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25d0\n", + "alpha = A_1*t*sin(B_1*x) + A_2*sin(B_2*x + C_2*t) + 1.0d0/4.0d0\n", "\n", "end function\n", "\n", - "REAL*8 function eta(A1, A2, B1, B2, C2, kappa, t, x, y)\n", + "REAL*8 function eta(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)\n", "implicit none\n", - "REAL*8, intent(in) :: A1\n", - "REAL*8, intent(in) :: A2\n", - "REAL*8, intent(in) :: B1\n", - "REAL*8, intent(in) :: B2\n", - "REAL*8, intent(in) :: C2\n", + "REAL*8, intent(in) :: A_1\n", + "REAL*8, intent(in) :: A_2\n", + "REAL*8, intent(in) :: B_1\n", + "REAL*8, intent(in) :: B_2\n", + "REAL*8, intent(in) :: C_2\n", "REAL*8, intent(in) :: kappa\n", "REAL*8, intent(in) :: t\n", "REAL*8, intent(in) :: x\n", "REAL*8, intent(in) :: y\n", "\n", - "eta = -0.5d0*tanh(0.70710678118654752d0*kappa**(-0.5d0)*(-A1*t*sin(B1*x &\n", - " ) - A2*sin(B2*x + C2*t) + y - 0.25d0)) + 0.5d0\n", + "eta = (1.0d0/2.0d0)*tanh(0.17677669529663688d0*kappa**(-0.5d0)*(4.0d0* &\n", + " A_1*t*sin(B_1*x) + 4.0d0*A_2*sin(B_2*x + C_2*t) - 4.0d0*y + 1.0d0 &\n", + " )) + 1.0d0/2.0d0\n", "\n", "end function\n", "\n", - "REAL*8 function eta0(A1, A2, B1, B2, C2, kappa, t, x, y)\n", + "REAL*8 function eta0(A_2, B_2, kappa, x, y)\n", "implicit none\n", - "REAL*8, intent(in) :: A1\n", - "REAL*8, intent(in) :: A2\n", - "REAL*8, intent(in) :: B1\n", - "REAL*8, intent(in) :: B2\n", - "REAL*8, intent(in) :: C2\n", + "REAL*8, intent(in) :: A_2\n", + "REAL*8, intent(in) :: B_2\n", "REAL*8, intent(in) :: kappa\n", - "REAL*8, intent(in) :: t\n", "REAL*8, intent(in) :: x\n", "REAL*8, intent(in) :: y\n", "\n", - "eta0 = -0.5d0*tanh(0.70710678118654752d0*kappa**(-0.5d0)*(-A1*t*sin(B1*x &\n", - " ) - A2*sin(B2*x + C2*t) + y - 0.25d0)) + 0.5d0\n", + "eta0 = (1.0d0/2.0d0)*tanh(0.17677669529663688d0*kappa**(-0.5d0)*(4.0d0* &\n", + " A_2*sin(B_2*x) - 4.0d0*y + 1.0d0)) + 1.0d0/2.0d0\n", "\n", "end function\n", "\n", - "REAL*8 function S(S)\n", + "REAL*8 function S(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)\n", "implicit none\n", - "REAL*8, intent(in) :: S\n", + "REAL*8, intent(in) :: A_1\n", + "REAL*8, intent(in) :: A_2\n", + "REAL*8, intent(in) :: B_1\n", + "REAL*8, intent(in) :: B_2\n", + "REAL*8, intent(in) :: C_2\n", + "REAL*8, intent(in) :: kappa\n", + "REAL*8, intent(in) :: t\n", + "REAL*8, intent(in) :: x\n", + "REAL*8, intent(in) :: y\n", "\n", - "S = S\n", + "S = -1.0d0/4.0d0*(2*sqrt(kappa)*(A_1*B_1*t*cos(B_1*x) + A_2*B_2*cos(B_2* &\n", + " x + C_2*t))**2*tanh(0.70710678118654752d0*kappa**(-0.5d0)*(-A_1*t &\n", + " *sin(B_1*x) - A_2*sin(B_2*x + C_2*t) + y - 0.25d0)) + sqrt(2.0d0) &\n", + " *kappa*(-A_1*B_1**2*t*sin(B_1*x) - A_2*B_2**2*sin(B_2*x + C_2*t &\n", + " )) - sqrt(2.0d0)*(A_1*sin(B_1*x) + A_2*C_2*cos(B_2*x + C_2*t)))* &\n", + " 1d0/((1.0d0/2.0d0)*exp(0.70710678118654752d0*kappa**(-0.5d0)*( &\n", + " 0.5d0*cmplx(0,1)*A_1*t*(exp(cmplx(0,1)*B_1*x) - exp(-cmplx(0,1)* &\n", + " B_1*x)) + 0.5d0*cmplx(0,1)*A_2*(-exp(cmplx(0,1)*(-B_2*x - C_2*t &\n", + " )) + exp(cmplx(0,1)*(B_2*x + C_2*t))) + y - 0.25d0)) + (1.0d0/ &\n", + " 2.0d0)*exp(-0.70710678118654752d0*kappa**(-0.5d0)*(0.5d0*cmplx(0, &\n", + " 1)*A_1*t*(exp(cmplx(0,1)*B_1*x) - exp(-cmplx(0,1)*B_1*x)) + 0.5d0 &\n", + " *cmplx(0,1)*A_2*(-exp(cmplx(0,1)*(-B_2*x - C_2*t)) + exp(cmplx(0, &\n", + " 1)*(B_2*x + C_2*t))) + y - 0.25d0)))**2/sqrt(kappa)\n", "\n", "end function\n", "\n" @@ -691,8 +957,8 @@ "[(f_name, code), (f_name, header)] = \\\n", "codegen([(\"alpha\", alpha),\n", " (\"eta\", eta),\n", - " (\"eta0\", eta),\n", - " (\"S\", S)],\n", + " (\"eta0\", eta0),\n", + " (\"S\", source)],\n", " language=\"f95\",\n", " prefix=\"manufactured\",\n", " project=\"PFHub\")\n", @@ -710,8 +976,11 @@ }, { "cell_type": "code", - "execution_count": 8, - "metadata": {}, + "execution_count": 20, + "metadata": { + "scrolled": false, + "tags": [] + }, "outputs": [ { "name": "stdout", @@ -719,36 +988,36 @@ "text": [ "manufactured.jl:\n", "\n", - "# Code generated with sympy 1.2\n", + "# Code generated with SymPy 1.11.1\n", "#\n", "# See http://www.sympy.org/ for more information.\n", "#\n", "# This file is part of 'PFHub'\n", "\n", - "function alpha(A1, A2, B1, B2, C2, t, x)\n", + "function alpha(A_1, A_2, B_1, B_2, C_2, t, x)\n", "\n", - " out1 = A1.*t.*sin(B1.*x) + A2.*sin(B2.*x + C2.*t) + 0.25\n", + " out1 = A_1 .* t .* sin(B_1 .* x) + A_2 .* sin(B_2 .* x + C_2 .* t) + 1 // 4\n", "\n", " return out1\n", "end\n", "\n", - "function eta(A1, A2, B1, B2, C2, kappa, t, x, y)\n", + "function eta(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)\n", "\n", - " out1 = -0.5*tanh(sqrt(2)*(-A1.*t.*sin(B1.*x) - A2.*sin(B2.*x + C2.*t) + y - 0.25)./(2*sqrt(kappa))) + 0.5\n", + " out1 = tanh(sqrt(2) * (4 * A_1 .* t .* sin(B_1 .* x) + 4 * A_2 .* sin(B_2 .* x + C_2 .* t) - 4 * y + 1) ./ (8 * sqrt(kappa))) / 2 + 1 // 2\n", "\n", " return out1\n", "end\n", "\n", - "function eta0(A1, A2, B1, B2, C2, kappa, t, x, y)\n", + "function eta0(A_2, B_2, kappa, x, y)\n", "\n", - " out1 = -0.5*tanh(sqrt(2)*(-A1.*t.*sin(B1.*x) - A2.*sin(B2.*x + C2.*t) + y - 0.25)./(2*sqrt(kappa))) + 0.5\n", + " out1 = tanh(sqrt(2) * (4 * A_2 .* sin(B_2 .* x) - 4 * y + 1) ./ (8 * sqrt(kappa))) / 2 + 1 // 2\n", "\n", " return out1\n", "end\n", "\n", - "function S(S)\n", + "function S(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)\n", "\n", - " out1 = S\n", + " out1 = -(2 * sqrt(kappa) .* (A_1 .* B_1 .* t .* cos(B_1 .* x) + A_2 .* B_2 .* cos(B_2 .* x + C_2 .* t)) .^ 2 .* tanh(sqrt(2) * (-A_1 .* t .* sin(B_1 .* x) - A_2 .* sin(B_2 .* x + C_2 .* t) + y - 1 // 4) ./ (2 * sqrt(kappa))) + sqrt(2) * kappa .* (-A_1 .* B_1 .^ 2 .* t .* sin(B_1 .* x) - A_2 .* B_2 .^ 2 .* sin(B_2 .* x + C_2 .* t)) - sqrt(2) * (A_1 .* sin(B_1 .* x) + A_2 .* C_2 .* cos(B_2 .* x + C_2 .* t))) .* sech(sqrt(2) * (-A_1 .* t .* sin(B_1 .* x) - A_2 .* sin(B_2 .* x + C_2 .* t) + y - 1 // 4) ./ (2 * sqrt(kappa))) .^ 2 ./ (4 * sqrt(kappa))\n", "\n", " return out1\n", "end\n", @@ -760,8 +1029,8 @@ "[(f_name, code)] = \\\n", "codegen([(\"alpha\", alpha),\n", " (\"eta\", eta),\n", - " (\"eta0\", eta),\n", - " (\"S\", S)],\n", + " (\"eta0\", eta0),\n", + " (\"S\", source)],\n", " language=\"julia\",\n", " prefix=\"manufactured\",\n", " project=\"PFHub\")\n", @@ -779,20 +1048,23 @@ }, { "cell_type": "code", - "execution_count": 9, - "metadata": {}, + "execution_count": 21, + "metadata": { + "scrolled": false, + "tags": [] + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "alpha = A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] + 0.25 \n", + "alpha = A_1*t*Sin[B_1*x] + A_2*Sin[B_2*x + C_2*t] + 1/4 \n", "\n", - "eta = -0.5*Tanh[(1/2)*2^(1/2)*(-A1*t*Sin[B1*x] - A2*Sin[B2*x + C2*t] + y - 0.25)/kappa^(1/2)] + 0.5 \n", + "eta = (1/2)*Tanh[(1/8)*2^(1/2)*(4*A_1*t*Sin[B_1*x] + 4*A_2*Sin[B_2*x + C_2*t] - 4*y + 1)/kappa^(1/2)] + 1/2 \n", "\n", - "eta0 = -0.5*Tanh[(1/2)*2^(1/2)*(-A2*Sin[B2*x] + y - 0.25)/kappa^(1/2)] + 0.5 \n", + "eta0 = (1/2)*Tanh[(1/8)*2^(1/2)*(4*A_2*Sin[B_2*x] - 4*y + 1)/kappa^(1/2)] + 1/2 \n", "\n", - "S = -(Tanh[(1/2)*2^(1/2)*(A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] - y + 0.25)/kappa^(1/2)]^2 - 1)*(0.5*kappa^(1/2)*((A1*B1*t*Cos[B1*x] + A2*B2*Cos[B2*x + C2*t])^2 + 1)*Tanh[(1/2)*2^(1/2)*(A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] - y + 0.25)/kappa^(1/2)] - 0.5*kappa^(1/2)*Tanh[(1/2)*2^(1/2)*(A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] - y + 0.25)/kappa^(1/2)] + 0.25*2^(1/2)*kappa*(A1*B1^2*t*Sin[B1*x] + A2*B2^2*Sin[B2*x + C2*t]) + 0.25*2^(1/2)*(A1*Sin[B1*x] + A2*C2*Cos[B2*x + C2*t]))/kappa^(1/2) \n", + "S = -1/4*(2*kappa^(1/2)*(A_1*B_1*t*Cos[B_1*x] + A_2*B_2*Cos[B_2*x + C_2*t])^2*Tanh[(1/2)*2^(1/2)*(-A_1*t*Sin[B_1*x] - A_2*Sin[B_2*x + C_2*t] + y - 1/4)/kappa^(1/2)] + 2^(1/2)*kappa*(-A_1*B_1^2*t*Sin[B_1*x] - A_2*B_2^2*Sin[B_2*x + C_2*t]) - 2^(1/2)*(A_1*Sin[B_1*x] + A_2*C_2*Cos[B_2*x + C_2*t]))*Sech[(1/2)*2^(1/2)*(-A_1*t*Sin[B_1*x] - A_2*Sin[B_2*x + C_2*t] + y - 1/4)/kappa^(1/2)]^2/kappa^(1/2) \n", "\n" ] } @@ -815,8 +1087,11 @@ }, { "cell_type": "code", - "execution_count": 10, - "metadata": {}, + "execution_count": 22, + "metadata": { + "scrolled": false, + "tags": [] + }, "outputs": [ { "name": "stdout", @@ -825,33 +1100,33 @@ "manufactured.nb:\n", "\n", "alpha.m\n", - "function out1 = alpha(A1, A2, B1, B2, C2, t, x)\n", - " %ALPHA Autogenerated by sympy\n", - " % Code generated with sympy 1.2\n", + "function out1 = alpha(A_1, A_2, B_1, B_2, C_2, t, x)\n", + " %ALPHA Autogenerated by SymPy\n", + " % Code generated with SymPy 1.11.1\n", " %\n", " % See http://www.sympy.org/ for more information.\n", " %\n", " % This file is part of 'PFHub'\n", "\n", - " out1 = A1.*t.*sin(B1.*x) + A2.*sin(B2.*x + C2.*t) + 0.25;\n", + " out1 = A_1.*t.*sin(B_1.*x) + A_2.*sin(B_2.*x + C_2.*t) + 1/4;\n", "\n", "end\n", "\n", - "function out1 = eta(A1, A2, B1, B2, C2, kappa, t, x, y)\n", + "function out1 = eta(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)\n", "\n", - " out1 = -0.5*tanh(sqrt(2)*(-A1.*t.*sin(B1.*x) - A2.*sin(B2.*x + C2.*t) + y - 0.25)./(2*sqrt(kappa))) + 0.5;\n", + " out1 = tanh(sqrt(2)*(4*A_1.*t.*sin(B_1.*x) + 4*A_2.*sin(B_2.*x + C_2.*t) - 4*y + 1)./(8*sqrt(kappa)))/2 + 1/2;\n", "\n", "end\n", "\n", - "function out1 = eta0(A1, A2, B1, B2, C2, kappa, t, x, y)\n", + "function out1 = eta0(A_2, B_2, kappa, x, y)\n", "\n", - " out1 = -0.5*tanh(sqrt(2)*(-A1.*t.*sin(B1.*x) - A2.*sin(B2.*x + C2.*t) + y - 0.25)./(2*sqrt(kappa))) + 0.5;\n", + " out1 = tanh(sqrt(2)*(4*A_2.*sin(B_2.*x) - 4*y + 1)./(8*sqrt(kappa)))/2 + 1/2;\n", "\n", "end\n", "\n", - "function out1 = S(S)\n", + "function out1 = S(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)\n", "\n", - " out1 = S;\n", + " out1 = -(2*sqrt(kappa).*(A_1.*B_1.*t.*cos(B_1.*x) + A_2.*B_2.*cos(B_2.*x + C_2.*t)).^2.*tanh(sqrt(2)*(-A_1.*t.*sin(B_1.*x) - A_2.*sin(B_2.*x + C_2.*t) + y - 1/4)./(2*sqrt(kappa))) + sqrt(2)*kappa.*(-A_1.*B_1.^2.*t.*sin(B_1.*x) - A_2.*B_2.^2.*sin(B_2.*x + C_2.*t)) - sqrt(2)*(A_1.*sin(B_1.*x) + A_2.*C_2.*cos(B_2.*x + C_2.*t))).*sech(sqrt(2)*(-A_1.*t.*sin(B_1.*x) - A_2.*sin(B_2.*x + C_2.*t) + y - 1/4)./(2*sqrt(kappa))).^2./(4*sqrt(kappa));\n", "\n", "end\n", "\n" @@ -862,8 +1137,8 @@ "code = \\\n", "codegen([(\"alpha\", alpha),\n", " (\"eta\", eta),\n", - " (\"eta0\", eta),\n", - " (\"S\", S)],\n", + " (\"eta0\", eta0),\n", + " (\"S\", source)],\n", " language=\"octave\",\n", " project=\"PFHub\")\n", "\n", @@ -876,7 +1151,7 @@ "metadata": { "anaconda-cloud": {}, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -912,5 +1187,5 @@ } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/benchmarks/benchmark7.ipynb.raw.html b/benchmarks/benchmark7.ipynb.raw.html index 5e638ea69..7271011f2 100644 --- a/benchmarks/benchmark7.ipynb.raw.html +++ b/benchmarks/benchmark7.ipynb.raw.html @@ -26,13 +26,14 @@ -
+
@@ -158,7 +159,7 @@

Benchmark Problem 7: MMS Allen-Cahn
-
In [1]:
+
In [4]:
from IPython.display import HTML
@@ -176,7 +177,7 @@ 
Benchmark Problem 7: MMS Allen-Cahn
 
 

 
-
Out[1]:

+
Out[4]:

 
 
 
@@ -243,7 +244,7 @@ 
Governing equation and man
 \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.

+
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_{\mathrm{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]
@@ -258,7 +259,8 @@ 
Governing equation and man
 
$$\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 .

+
N.B.: Don't transcribe these equations. Please download the appropriate files from the Appendix.

+
@@ -267,13 +269,13 @@

Governing equation and man

-

Domain geometry, boundary conditions, initial conditions, and stopping condition

The 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$ = 1 and $\eta$ = 0 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$:

+

Domain geometry, boundary conditions, initial conditions, and stopping condition

The 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_{\mathrm{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_{\mathrm{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$ = 1 and $\eta$ = 0 are applied along the $y = 0$ and the $y = 0.5$ boundaries, respectively. These boundary conditions are chosen to be consistent with $\eta_{\mathrm{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] +\eta_{\mathrm{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}$).

+

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_{\mathrm{sol}}$).

@@ -336,9 +338,9 @@

Benchmark simulation instructionsPart (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} + L_2 = \sqrt{\sum\limits_{x,y}\left(\eta^{t=8}_{x,y} - \eta_{\mathrm{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 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 $[1 \times 10^{-4}, 5 \times 10^{-3}]$, 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.

Save the effective element size, $h$, and the $L_2$ error for each simulation. Archive this data in a CSV or JSON file, using one column (or key) each for $h$ and $L_2$. @@ -346,13 +348,13 @@

Part (a)

submit your results on 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.

+5% of the elements. Then submit your results on 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.

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. Archive the effective mesh size and $L_2$ error for each individual simulation in a CSV or JSON file. Submit your results to the PFHub website 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.

+

where $E$ is the $L_2$ error, $R$ is the effective element size, and b is an intercept. Deviations of $\pm 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 $[1\times10^{-4}, 5\times10^{-3}]$, attempting to cover as much of that range as possible/practical. Archive the effective mesh size and $L_2$ error for each individual simulation in a CSV or JSON file. Submit your results to the PFHub website 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.

@@ -361,9 +363,9 @@

Part (a)

-

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 (in seconds), number of computing cores, normalized computing cost (wall time in seconds $\times$ number of cores $\times$ nominal core speed $/$ 2 GHz), maximum memory usage, and $L_2$ error at $t=8$ for each individual simulation. Archive this data in a -CSV or JSON file with one column (or key) for each of the quantities mentioned above. Submit your results to the PFHub website as two 2D data sets. For the first data set use the $L_2$ error as the x-axis column and the normalized computational cost as the y-axis column. For the second data set, use the $L_2$ error as the x-axis column and the wall time as the y-axis column.

+

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 \times 10^{-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 $[1 \times 10^{-5}, 5 \times 10^{-3}]$. Use the same CPU and processor type for all simulations. For the best of these simulations, save the wall time (in seconds), number of computing cores, normalized computing cost (wall time in seconds $\times$ number of cores $\times$ nominal core speed $/$ 2 GHz), maximum memory usage, and $L_2$ error at $t=8$ for each individual simulation. Archive this data in a +CSV or JSON file with one column (or key) for each of the quantities mentioned above. Submit your results to the PFHub website as two 2D data sets. For the first data set use the $L_2$ error as the $x$-axis column and the normalized computational cost as the $y$-axis column. For the second data set, use the $L_2$ error as the $x$-axis column and the wall time as the $y$-axis column.

@@ -382,7 +384,7 @@

Part (c)

-

Submission Guidelines

Part (a) Guidelines

Two data items are required in the "Data Files" section of the upload form. The data items should be labeled as spatial and temporal in the Short name of data box. The 2D radio button should be checked and the columns corresponding to the x-axis (either $\Delta t$ or $\Delta x$) and the y-axis ($e_{L2}$) should be labeled correctly for each CSV file. The CSV file for the spatial data should have the form

+

Submission Guidelines

Part (a) Guidelines

Two data items are required in the "Data Files" section of the upload form. The data items should be labeled as spatial and temporal in the Short name of data box. The 2D radio button should be checked and the columns corresponding to the $x$-axis (either $\Delta t$ or $\Delta x$) and the $y$-axis ($e_{L2}$) should be labeled correctly for each CSV file. The CSV file for the spatial data should have the form

mesh_size,L2_error
 0.002604167,2.55E-06
@@ -394,7 +396,7 @@ 
Submission Guidelines
-

Parts (b) and (c) Guidelines

Two data items are required in the "Data Files" section of the upload form. The data items should be labeled as cost and time in the Short name of data box. The 2D radio button should be checked and the columns corresponding to the x-axis ($e_{L2}$) and the y-axis (either $F_{\text{cost}}$ or $t_{\text{wall}}$) should be labeled correctly for each CSV file. The CSV file for the cost data should have the form

+

Parts (b) and (c) Guidelines

Two data items are required in the "Data Files" section of the upload form. The data items should be labeled as cost and time in the Short name of data box. The 2D radio button should be checked and the columns corresponding to the $x$-axis ($e_{L2}$) and the $y$-axis (either $F_{\text{cost}}$ or $t_{\text{wall}}$) should be labeled correctly for each CSV file. The CSV file for the cost data should have the form

cores,wall_time,memory,error,cost
 1,1.35,25800,0.024275131,1.755
@@ -427,7 +429,7 @@ 
Feedback

 

 

-
Appendix
Computer algebra systems
Rigorous 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, but so long as vector calculus is supported, CAS will do.

+
Appendix
Computer algebra systems
Rigorous 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, but so long as vector calculus is supported, any CAS will do.

 
 

 

@@ -447,42 +449,65 @@ 
Source term

 
# Sympy code to generate expressions for PFHub Problem 7 (MMS)
 
-from sympy import symbols, simplify
-from sympy import sin, cos, cosh, tanh, sqrt
-from sympy.physics.vector import divergence, gradient, ReferenceFrame, time_derivative
+from sympy import sin, cos, cosh, sech, tanh, sqrt, symbols
+from sympy import Rational, diff, expand, simplify
+from sympy import Eq, Function, init_printing
+from sympy.abc import x, y, t
+from sympy.physics.vector import ReferenceFrame, divergence, gradient, time_derivative
 from sympy.utilities.codegen import codegen
-from sympy.abc import kappa, S, x, y, t
-
-# 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')
-C2 = symbols('C2')
-
-# Define interface offset (alpha)
-alpha = 0.25 + A1 * t * sin(B1 * R[0]) \
-             + A2     * sin(B2 * R[0] + C2 * t)
-
-# Define the solution equation (eta)
-eta = 0.5 * (1 - tanh((R[1] - alpha) / sqrt(2*kappa)))
-
-# Compute the source term from the equation of motion
-source = simplify(time_derivative(eta, R)
-             + 4 * eta * (eta - 1) * (eta - 1/2)
-             - kappa * divergence(gradient(eta, R), R))
-
-# Replace R[i] with (x, y)
-alpha = alpha.subs({R[0]: x, R[1]: y})
-eta = eta.subs({R[0]: x, R[1]: y})
-eta0 = eta.subs(t, 0)
-source = source.subs({R[0]: x, R[1]: y})
-
-print("alpha =", alpha, "\n")
-print("eta =", eta,     "\n")
-print("eta0 =", eta0,  "\n")
-print("S =", source)
+
+init_printing()
+
+R = ReferenceFrame("R")  # spatial coords
+X = R[0]
+Y = R[1]
+
+A1, A2 = symbols("A_1 A_2",   real=True)
+B1, B2 = symbols("B_1 B_2",   real=True)
+C2, κ  = symbols("C_2 kappa", real=True)
+
+# Define α symbolically to simplify the expressions for visual comparison
+
+α = Function("alpha", real=True)(X, t)
+ψ = (Y - α) / sqrt(2 * κ)
+η = (1 - tanh(ψ)) / 2
+
+source = time_derivative(η, R) \
+       + 4 * η * (η - 1) * (η - Rational(1, 2)) \
+       - divergence(κ * gradient(η, R), R)
+
+# Nested `simplify` calls to incorporate the hyperbolic trig
+# identity, which is not implemented in SymPy
+        
+S = simplify(
+        expand(
+            simplify(
+                expand(source)
+            ).subs({cosh(ψ)**2: 1 / sech(ψ)**2})
+        )
+)
+

+
+

+
+
+ +
+
+
+
+
+

Compare the generated expressions for $S(x, y, \alpha)$, $\alpha$, and the derivatives of $\alpha$:

+ +
+
+
+
+
+
In [6]:
+
+
+
Eq(symbols("S"), S.subs({X: x, Y: y}))
@@ -495,24 +520,209 @@

Source term
-
+
Out[6]:
-
-
alpha = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25 
 
-eta = -0.5*tanh(sqrt(2)*(-A1*t*sin(B1*x) - A2*sin(B2*x + C2*t) + y - 0.25)/(2*sqrt(kappa))) + 0.5 
 
-eta0 = -0.5*tanh(sqrt(2)*(-A2*sin(B2*x) + y - 0.25)/(2*sqrt(kappa))) + 0.5 
+

+$\displaystyle S = - \frac{\left(2 \sqrt{\kappa} \tanh{\left(\frac{\sqrt{2} \left(y - \alpha{\left(x,t \right)}\right)}{2 \sqrt{\kappa}} \right)} \left(\frac{\partial}{\partial x} \alpha{\left(x,t \right)}\right)^{2} + \sqrt{2} \kappa \frac{\partial^{2}}{\partial x^{2}} \alpha{\left(x,t \right)} - \sqrt{2} \frac{\partial}{\partial t} \alpha{\left(x,t \right)}\right) \operatorname{sech}^{2}{\left(\frac{\sqrt{2} \left(y - \alpha{\left(x,t \right)}\right)}{2 \sqrt{\kappa}} \right)}}{4 \sqrt{\kappa}}$
+

+
+
+ +
+

+ +
+
+
+
In [7]:
+
+
+
alpha = Rational(1, 4) + A1 * t * sin(B1 * X) + A2 * sin(B2 * X + C2 * t)
+Eq(symbols("alpha"), alpha.subs({X: x, Y: y}))
+
+ +
+
+
+ +
+
+ + +
+ +
Out[7]:
+ + + + +
+$\displaystyle \alpha = A_{1} t \sin{\left(B_{1} x \right)} + A_{2} \sin{\left(B_{2} x + C_{2} t \right)} + \frac{1}{4}$ +
+ +
+ +
+
+ +
+
+
+
In [8]:
+
+
+
alpha_t  = simplify(expand(diff(alpha, t)))
+
+Eq(diff(α, t).subs({X: x, Y: y}), alpha_t.subs({X: x, Y: y}))
+
-S = -(tanh(sqrt(2)*(A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) - y + 0.25)/(2*sqrt(kappa)))**2 - 1)*(0.5*sqrt(kappa)*((A1*B1*t*cos(B1*x) + A2*B2*cos(B2*x + C2*t))**2 + 1)*tanh(sqrt(2)*(A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) - y + 0.25)/(2*sqrt(kappa))) - 0.5*sqrt(kappa)*tanh(sqrt(2)*(A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) - y + 0.25)/(2*sqrt(kappa))) + 0.25*sqrt(2)*kappa*(A1*B1**2*t*sin(B1*x) + A2*B2**2*sin(B2*x + C2*t)) + 0.25*sqrt(2)*(A1*sin(B1*x) + A2*C2*cos(B2*x + C2*t)))/sqrt(kappa) -
+
+
+ +
+
+ + +
+ +
Out[8]:
+ + + + +
+$\displaystyle \frac{\partial}{\partial t} \alpha{\left(x,t \right)} = A_{1} \sin{\left(B_{1} x \right)} + A_{2} C_{2} \cos{\left(B_{2} x + C_{2} t \right)}$ +
+
+
+
+
+
In [9]:
+
+
+
alpha_x  = simplify(expand(diff(alpha, X)))
+
+Eq(diff(α, X).subs({X: x, Y: y}), alpha_x.subs({X: x, Y: y}))
+
+ +
+
+
+ +
+
+ + +
+ +
Out[9]:
+ + + + +
+$\displaystyle \frac{\partial}{\partial x} \alpha{\left(x,t \right)} = A_{1} B_{1} t \cos{\left(B_{1} x \right)} + A_{2} B_{2} \cos{\left(B_{2} x + C_{2} t \right)}$ +
+ +
+ +
+
+ +
+
+
+
In [10]:
+
+
+
alpha_xx = simplify(expand(diff(alpha, X, 2)))
+
+Eq(diff(α, X, X).subs({X: x, Y: y}), alpha_xx.subs({X: x, Y: y}))
+
+ +
+
+
+ +
+
+ + +
+ +
Out[10]:
+ + + + +
+$\displaystyle \frac{\partial^{2}}{\partial x^{2}} \alpha{\left(x,t \right)} = - 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)}$ +
+ +
+ +
+
+ +
+
+
+
In [11]:
+
+
+
eta = simplify(η.subs({R[0]: x, R[1]: y, α: alpha.subs({R[0]: x, R[1]: y})}))
+
+ +
+
+
+ +
+
+
+
In [12]:
+
+
+
eta0 = simplify(eta.subs({t: 0, α: alpha.subs({R[0]: x, R[1]: y, t: 0})}))
+
+ +
+
+
+ +
+
+
+
In [13]:
+
+
+
source = S.subs(
+    {
+        α: alpha,
+        diff(α, t):    alpha_t,
+        diff(α, X):    alpha_x,
+        diff(α, X, 2): alpha_xx
+    }
+).subs({X: x, Y: y})
+
+alpha = alpha.subs({X: x, Y: y})
+eta = eta.subs({X: x, Y: y})
+eta0 = eta0.subs({X: x, Y: y})
+
+ +
+
+
+
@@ -526,20 +736,134 @@

Code

-

Python

Copy the first cell under Source Term directly into your program. -For a performance boost, convert the expressions into lambda functions:

-
from sympy.utilities.lambdify import lambdify
-
-apy = 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')
+
Python
Substituting the expression for $\alpha$ into $S(x, y, \alpha)$, the Python code for the full source term $S(x, y, t)$ is

+
+
+
+
+
+
+
In [14]:
+
+
+
print("S =", source)
-

Note: Click "Code Toggle" at the top of the page to see the Python expressions.

-
+
+ +
+
+ + +
+ +
+ + +
+
S = -(2*sqrt(kappa)*(A_1*B_1*t*cos(B_1*x) + A_2*B_2*cos(B_2*x + C_2*t))**2*tanh(sqrt(2)*(-A_1*t*sin(B_1*x) - A_2*sin(B_2*x + C_2*t) + y - 1/4)/(2*sqrt(kappa))) + sqrt(2)*kappa*(-A_1*B_1**2*t*sin(B_1*x) - A_2*B_2**2*sin(B_2*x + C_2*t)) - sqrt(2)*(A_1*sin(B_1*x) + A_2*C_2*cos(B_2*x + C_2*t)))*sech(sqrt(2)*(-A_1*t*sin(B_1*x) - A_2*sin(B_2*x + C_2*t) + y - 1/4)/(2*sqrt(kappa)))**2/(4*sqrt(kappa))
+
+
+
+ +
+
+ +
+
+
+
In [15]:
+
+
+
print("α =", alpha)
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
α = A_1*t*sin(B_1*x) + A_2*sin(B_2*x + C_2*t) + 1/4
+
+
+
+ +
+
+ +
+
+
+
In [16]:
+
+
+
print("η =", eta)
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
η = tanh(sqrt(2)*(4*A_1*t*sin(B_1*x) + 4*A_2*sin(B_2*x + C_2*t) - 4*y + 1)/(8*sqrt(kappa)))/2 + 1/2
+
+
+
+ +
+
+ +
+
+
+
In [17]:
+
+
+
print("η₀ =", eta0)
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
η₀ = tanh(sqrt(2)*(4*A_2*sin(B_2*x) - 4*y + 1)/(8*sqrt(kappa)))/2 + 1/2
+
+
+
+ +
+
+
@@ -551,14 +875,14 @@

C

-
In [6]:
+
In [18]:
[(c_name, code), (h_name, header)] = \
 codegen([("alpha", alpha),
          ("eta",   eta),
-         ("eta0",  eta),
-         ("S",     S)],
+         ("eta0",  eta0),
+         ("S",     source)],
          language="C",
          prefix="manufactured",
          project="PFHub")
@@ -585,7 +909,7 @@ 
C

 
manufactured.h:
 
 /******************************************************************************
- *                       Code generated with sympy 1.2                        *
+ *                      Code generated with SymPy 1.11.1                      *
  *                                                                            *
  *              See http://www.sympy.org/ for more information.               *
  *                                                                            *
@@ -596,10 +920,10 @@ 
C

 #ifndef PFHUB__MANUFACTURED__H
 #define PFHUB__MANUFACTURED__H
 
-double alpha(double A1, double A2, double B1, double B2, double C2, double t, double x);
-double eta(double A1, double A2, double B1, double B2, double C2, double kappa, double t, double x, double y);
-double eta0(double A1, double A2, double B1, double B2, double C2, double kappa, double t, double x, double y);
-double S(double S);
+double alpha(double A_1, double A_2, double B_1, double B_2, double C_2, double t, double x);
+double eta(double A_1, double A_2, double B_1, double B_2, double C_2, double kappa, double t, double x, double y);
+double eta0(double A_2, double B_2, double kappa, double x, double y);
+double S(double A_1, double A_2, double B_1, double B_2, double C_2, double kappa, double t, double x, double y);
 
 #endif
 
@@ -608,7 +932,7 @@ 
C

 manufactured.c:
 
 /******************************************************************************
- *                       Code generated with sympy 1.2                        *
+ *                      Code generated with SymPy 1.11.1                      *
  *                                                                            *
  *              See http://www.sympy.org/ for more information.               *
  *                                                                            *
@@ -617,34 +941,34 @@ 
C

 #include "manufactured.h"
 #include <math.h>
 
-double alpha(double A1, double A2, double B1, double B2, double C2, double t, double x) {
+double alpha(double A_1, double A_2, double B_1, double B_2, double C_2, double t, double x) {
 
    double alpha_result;
-   alpha_result = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25;
+   alpha_result = A_1*t*sin(B_1*x) + A_2*sin(B_2*x + C_2*t) + 1.0/4.0;
    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(double A_1, double A_2, double B_1, double B_2, double C_2, double kappa, double t, double x, double y) {
 
    double eta_result;
-   eta_result = -0.5*tanh((1.0/2.0)*M_SQRT2*(-A1*t*sin(B1*x) - A2*sin(B2*x + C2*t) + y - 0.25)/sqrt(kappa)) + 0.5;
+   eta_result = (1.0/2.0)*tanh((1.0/8.0)*M_SQRT2*(4*A_1*t*sin(B_1*x) + 4*A_2*sin(B_2*x + C_2*t) - 4*y + 1)/sqrt(kappa)) + 1.0/2.0;
    return eta_result;
 
 }
 
-double eta0(double A1, double A2, double B1, double B2, double C2, double kappa, double t, double x, double y) {
+double eta0(double A_2, double B_2, double kappa, double x, double y) {
 
    double eta0_result;
-   eta0_result = -0.5*tanh((1.0/2.0)*M_SQRT2*(-A1*t*sin(B1*x) - A2*sin(B2*x + C2*t) + y - 0.25)/sqrt(kappa)) + 0.5;
+   eta0_result = (1.0/2.0)*tanh((1.0/8.0)*M_SQRT2*(4*A_2*sin(B_2*x) - 4*y + 1)/sqrt(kappa)) + 1.0/2.0;
    return eta0_result;
 
 }
 
-double S(double S) {
+double S(double A_1, double A_2, double B_1, double B_2, double C_2, double kappa, double t, double x, double y) {
 
    double S_result;
-   S_result = S;
+   S_result = -1.0/4.0*(2*sqrt(kappa)*pow(A_1*B_1*t*cos(B_1*x) + A_2*B_2*cos(B_2*x + C_2*t), 2)*tanh((1.0/2.0)*M_SQRT2*(-A_1*t*sin(B_1*x) - A_2*sin(B_2*x + C_2*t) + y - 1.0/4.0)/sqrt(kappa)) + M_SQRT2*kappa*(-A_1*pow(B_1, 2)*t*sin(B_1*x) - A_2*pow(B_2, 2)*sin(B_2*x + C_2*t)) - M_SQRT2*(A_1*sin(B_1*x) + A_2*C_2*cos(B_2*x + C_2*t)))*pow(1.0/((1.0/2.0)*exp((1.0/2.0)*M_SQRT2*((1.0/2.0)*I*A_1*t*(exp(I*B_1*x) - exp(-I*B_1*x)) + (1.0/2.0)*I*A_2*(-exp(I*(-B_2*x - C_2*t)) + exp(I*(B_2*x + C_2*t))) + y - 1.0/4.0)/sqrt(kappa)) + (1.0/2.0)*exp(-1.0/2.0*M_SQRT2*((1.0/2.0)*I*A_1*t*(exp(I*B_1*x) - exp(-I*B_1*x)) + (1.0/2.0)*I*A_2*(-exp(I*(-B_2*x - C_2*t)) + exp(I*(B_2*x + C_2*t))) + y - 1.0/4.0)/sqrt(kappa))), 2)/sqrt(kappa);
    return S_result;
 
 }
@@ -667,14 +991,14 @@ 
Fortran
-
In [7]:
+
In [19]:
[(f_name, code), (f_name, header)] = \
 codegen([("alpha", alpha),
          ("eta",   eta),
-         ("eta0",  eta),
-         ("S",     S)],
+         ("eta0",  eta0),
+         ("S",     source)],
          language="f95",
          prefix="manufactured",
          project="PFHub")
@@ -700,66 +1024,83 @@ 
Fortran

 
manufactured.f:
 
 !******************************************************************************
-!*                       Code generated with sympy 1.2                        *
+!*                      Code generated with SymPy 1.11.1                      *
 !*                                                                            *
 !*              See http://www.sympy.org/ for more information.               *
 !*                                                                            *
 !*                        This file is part of 'PFHub'                        *
 !******************************************************************************
 
-REAL*8 function alpha(A1, A2, B1, B2, C2, t, x)
+REAL*8 function alpha(A_1, A_2, B_1, B_2, C_2, t, x)
 implicit none
-REAL*8, intent(in) :: A1
-REAL*8, intent(in) :: A2
-REAL*8, intent(in) :: B1
-REAL*8, intent(in) :: B2
-REAL*8, intent(in) :: C2
+REAL*8, intent(in) :: A_1
+REAL*8, intent(in) :: A_2
+REAL*8, intent(in) :: B_1
+REAL*8, intent(in) :: B_2
+REAL*8, intent(in) :: C_2
 REAL*8, intent(in) :: t
 REAL*8, intent(in) :: x
 
-alpha = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25d0
+alpha = A_1*t*sin(B_1*x) + A_2*sin(B_2*x + C_2*t) + 1.0d0/4.0d0
 
 end function
 
-REAL*8 function eta(A1, A2, B1, B2, C2, kappa, t, x, y)
+REAL*8 function eta(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)
 implicit none
-REAL*8, intent(in) :: A1
-REAL*8, intent(in) :: A2
-REAL*8, intent(in) :: B1
-REAL*8, intent(in) :: B2
-REAL*8, intent(in) :: C2
+REAL*8, intent(in) :: A_1
+REAL*8, intent(in) :: A_2
+REAL*8, intent(in) :: B_1
+REAL*8, intent(in) :: B_2
+REAL*8, intent(in) :: C_2
 REAL*8, intent(in) :: kappa
 REAL*8, intent(in) :: t
 REAL*8, intent(in) :: x
 REAL*8, intent(in) :: y
 
-eta = -0.5d0*tanh(0.70710678118654752d0*kappa**(-0.5d0)*(-A1*t*sin(B1*x &
-      ) - A2*sin(B2*x + C2*t) + y - 0.25d0)) + 0.5d0
+eta = (1.0d0/2.0d0)*tanh(0.17677669529663688d0*kappa**(-0.5d0)*(4.0d0* &
+      A_1*t*sin(B_1*x) + 4.0d0*A_2*sin(B_2*x + C_2*t) - 4.0d0*y + 1.0d0 &
+      )) + 1.0d0/2.0d0
 
 end function
 
-REAL*8 function eta0(A1, A2, B1, B2, C2, kappa, t, x, y)
+REAL*8 function eta0(A_2, B_2, kappa, x, y)
 implicit none
-REAL*8, intent(in) :: A1
-REAL*8, intent(in) :: A2
-REAL*8, intent(in) :: B1
-REAL*8, intent(in) :: B2
-REAL*8, intent(in) :: C2
+REAL*8, intent(in) :: A_2
+REAL*8, intent(in) :: B_2
 REAL*8, intent(in) :: kappa
-REAL*8, intent(in) :: t
 REAL*8, intent(in) :: x
 REAL*8, intent(in) :: y
 
-eta0 = -0.5d0*tanh(0.70710678118654752d0*kappa**(-0.5d0)*(-A1*t*sin(B1*x &
-      ) - A2*sin(B2*x + C2*t) + y - 0.25d0)) + 0.5d0
+eta0 = (1.0d0/2.0d0)*tanh(0.17677669529663688d0*kappa**(-0.5d0)*(4.0d0* &
+      A_2*sin(B_2*x) - 4.0d0*y + 1.0d0)) + 1.0d0/2.0d0
 
 end function
 
-REAL*8 function S(S)
+REAL*8 function S(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)
 implicit none
-REAL*8, intent(in) :: S
+REAL*8, intent(in) :: A_1
+REAL*8, intent(in) :: A_2
+REAL*8, intent(in) :: B_1
+REAL*8, intent(in) :: B_2
+REAL*8, intent(in) :: C_2
+REAL*8, intent(in) :: kappa
+REAL*8, intent(in) :: t
+REAL*8, intent(in) :: x
+REAL*8, intent(in) :: y
 
-S = S
+S = -1.0d0/4.0d0*(2*sqrt(kappa)*(A_1*B_1*t*cos(B_1*x) + A_2*B_2*cos(B_2* &
+      x + C_2*t))**2*tanh(0.70710678118654752d0*kappa**(-0.5d0)*(-A_1*t &
+      *sin(B_1*x) - A_2*sin(B_2*x + C_2*t) + y - 0.25d0)) + sqrt(2.0d0) &
+      *kappa*(-A_1*B_1**2*t*sin(B_1*x) - A_2*B_2**2*sin(B_2*x + C_2*t &
+      )) - sqrt(2.0d0)*(A_1*sin(B_1*x) + A_2*C_2*cos(B_2*x + C_2*t)))* &
+      1d0/((1.0d0/2.0d0)*exp(0.70710678118654752d0*kappa**(-0.5d0)*( &
+      0.5d0*cmplx(0,1)*A_1*t*(exp(cmplx(0,1)*B_1*x) - exp(-cmplx(0,1)* &
+      B_1*x)) + 0.5d0*cmplx(0,1)*A_2*(-exp(cmplx(0,1)*(-B_2*x - C_2*t &
+      )) + exp(cmplx(0,1)*(B_2*x + C_2*t))) + y - 0.25d0)) + (1.0d0/ &
+      2.0d0)*exp(-0.70710678118654752d0*kappa**(-0.5d0)*(0.5d0*cmplx(0, &
+      1)*A_1*t*(exp(cmplx(0,1)*B_1*x) - exp(-cmplx(0,1)*B_1*x)) + 0.5d0 &
+      *cmplx(0,1)*A_2*(-exp(cmplx(0,1)*(-B_2*x - C_2*t)) + exp(cmplx(0, &
+      1)*(B_2*x + C_2*t))) + y - 0.25d0)))**2/sqrt(kappa)
 
 end function
 
@@ -781,14 +1122,14 @@ 
Julia
-
In [8]:
+
In [20]:
[(f_name, code)] = \
 codegen([("alpha", alpha),
          ("eta",   eta),
-         ("eta0",  eta),
-         ("S",     S)],
+         ("eta0",  eta0),
+         ("S",     source)],
          language="julia",
          prefix="manufactured",
          project="PFHub")
@@ -813,36 +1154,36 @@ 
Julia

 

 
manufactured.jl:
 
-#   Code generated with sympy 1.2
+#   Code generated with SymPy 1.11.1
 #
 #   See http://www.sympy.org/ for more information.
 #
 #   This file is part of 'PFHub'
 
-function alpha(A1, A2, B1, B2, C2, t, x)
+function alpha(A_1, A_2, B_1, B_2, C_2, t, x)
 
-    out1 = A1.*t.*sin(B1.*x) + A2.*sin(B2.*x + C2.*t) + 0.25
+    out1 = A_1 .* t .* sin(B_1 .* x) + A_2 .* sin(B_2 .* x + C_2 .* t) + 1 // 4
 
     return out1
 end
 
-function eta(A1, A2, B1, B2, C2, kappa, t, x, y)
+function eta(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)
 
-    out1 = -0.5*tanh(sqrt(2)*(-A1.*t.*sin(B1.*x) - A2.*sin(B2.*x + C2.*t) + y - 0.25)./(2*sqrt(kappa))) + 0.5
+    out1 = tanh(sqrt(2) * (4 * A_1 .* t .* sin(B_1 .* x) + 4 * A_2 .* sin(B_2 .* x + C_2 .* t) - 4 * y + 1) ./ (8 * sqrt(kappa))) / 2 + 1 // 2
 
     return out1
 end
 
-function eta0(A1, A2, B1, B2, C2, kappa, t, x, y)
+function eta0(A_2, B_2, kappa, x, y)
 
-    out1 = -0.5*tanh(sqrt(2)*(-A1.*t.*sin(B1.*x) - A2.*sin(B2.*x + C2.*t) + y - 0.25)./(2*sqrt(kappa))) + 0.5
+    out1 = tanh(sqrt(2) * (4 * A_2 .* sin(B_2 .* x) - 4 * y + 1) ./ (8 * sqrt(kappa))) / 2 + 1 // 2
 
     return out1
 end
 
-function S(S)
+function S(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)
 
-    out1 = S
+    out1 = -(2 * sqrt(kappa) .* (A_1 .* B_1 .* t .* cos(B_1 .* x) + A_2 .* B_2 .* cos(B_2 .* x + C_2 .* t)) .^ 2 .* tanh(sqrt(2) * (-A_1 .* t .* sin(B_1 .* x) - A_2 .* sin(B_2 .* x + C_2 .* t) + y - 1 // 4) ./ (2 * sqrt(kappa))) + sqrt(2) * kappa .* (-A_1 .* B_1 .^ 2 .* t .* sin(B_1 .* x) - A_2 .* B_2 .^ 2 .* sin(B_2 .* x + C_2 .* t)) - sqrt(2) * (A_1 .* sin(B_1 .* x) + A_2 .* C_2 .* cos(B_2 .* x + C_2 .* t))) .* sech(sqrt(2) * (-A_1 .* t .* sin(B_1 .* x) - A_2 .* sin(B_2 .* x + C_2 .* t) + y - 1 // 4) ./ (2 * sqrt(kappa))) .^ 2 ./ (4 * sqrt(kappa))
 
     return out1
 end
@@ -865,7 +1206,7 @@ 
Mathematica

 

 

-
In [9]:

+
In [21]:

 

     

 
from sympy.printing import mathematica_code
@@ -890,13 +1231,13 @@ 
Mathematica

-
alpha = A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] + 0.25 
+
alpha = A_1*t*Sin[B_1*x] + A_2*Sin[B_2*x + C_2*t] + 1/4 
 
-eta = -0.5*Tanh[(1/2)*2^(1/2)*(-A1*t*Sin[B1*x] - A2*Sin[B2*x + C2*t] + y - 0.25)/kappa^(1/2)] + 0.5 
+eta = (1/2)*Tanh[(1/8)*2^(1/2)*(4*A_1*t*Sin[B_1*x] + 4*A_2*Sin[B_2*x + C_2*t] - 4*y + 1)/kappa^(1/2)] + 1/2 
 
-eta0 = -0.5*Tanh[(1/2)*2^(1/2)*(-A2*Sin[B2*x] + y - 0.25)/kappa^(1/2)] + 0.5 
+eta0 = (1/2)*Tanh[(1/8)*2^(1/2)*(4*A_2*Sin[B_2*x] - 4*y + 1)/kappa^(1/2)] + 1/2 
 
-S = -(Tanh[(1/2)*2^(1/2)*(A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] - y + 0.25)/kappa^(1/2)]^2 - 1)*(0.5*kappa^(1/2)*((A1*B1*t*Cos[B1*x] + A2*B2*Cos[B2*x + C2*t])^2 + 1)*Tanh[(1/2)*2^(1/2)*(A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] - y + 0.25)/kappa^(1/2)] - 0.5*kappa^(1/2)*Tanh[(1/2)*2^(1/2)*(A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] - y + 0.25)/kappa^(1/2)] + 0.25*2^(1/2)*kappa*(A1*B1^2*t*Sin[B1*x] + A2*B2^2*Sin[B2*x + C2*t]) + 0.25*2^(1/2)*(A1*Sin[B1*x] + A2*C2*Cos[B2*x + C2*t]))/kappa^(1/2) 
+S = -1/4*(2*kappa^(1/2)*(A_1*B_1*t*Cos[B_1*x] + A_2*B_2*Cos[B_2*x + C_2*t])^2*Tanh[(1/2)*2^(1/2)*(-A_1*t*Sin[B_1*x] - A_2*Sin[B_2*x + C_2*t] + y - 1/4)/kappa^(1/2)] + 2^(1/2)*kappa*(-A_1*B_1^2*t*Sin[B_1*x] - A_2*B_2^2*Sin[B_2*x + C_2*t]) - 2^(1/2)*(A_1*Sin[B_1*x] + A_2*C_2*Cos[B_2*x + C_2*t]))*Sech[(1/2)*2^(1/2)*(-A_1*t*Sin[B_1*x] - A_2*Sin[B_2*x + C_2*t] + y - 1/4)/kappa^(1/2)]^2/kappa^(1/2) 
 
 

 

@@ -916,14 +1257,14 @@ 
Matlab

 

 

 

-
In [10]:

+
In [22]:

 

     

 
code = \
 codegen([("alpha", alpha),
          ("eta",   eta),
-         ("eta0",  eta),
-         ("S",     S)],
+         ("eta0",  eta0),
+         ("S",     source)],
          language="octave",
          project="PFHub")
 
@@ -949,33 +1290,33 @@ 
Matlab

 
manufactured.nb:
 
 alpha.m
-function out1 = alpha(A1, A2, B1, B2, C2, t, x)
-  %ALPHA  Autogenerated by sympy
-  %   Code generated with sympy 1.2
+function out1 = alpha(A_1, A_2, B_1, B_2, C_2, t, x)
+  %ALPHA  Autogenerated by SymPy
+  %   Code generated with SymPy 1.11.1
   %
   %   See http://www.sympy.org/ for more information.
   %
   %   This file is part of 'PFHub'
 
-  out1 = A1.*t.*sin(B1.*x) + A2.*sin(B2.*x + C2.*t) + 0.25;
+  out1 = A_1.*t.*sin(B_1.*x) + A_2.*sin(B_2.*x + C_2.*t) + 1/4;
 
 end
 
-function out1 = eta(A1, A2, B1, B2, C2, kappa, t, x, y)
+function out1 = eta(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)
 
-  out1 = -0.5*tanh(sqrt(2)*(-A1.*t.*sin(B1.*x) - A2.*sin(B2.*x + C2.*t) + y - 0.25)./(2*sqrt(kappa))) + 0.5;
+  out1 = tanh(sqrt(2)*(4*A_1.*t.*sin(B_1.*x) + 4*A_2.*sin(B_2.*x + C_2.*t) - 4*y + 1)./(8*sqrt(kappa)))/2 + 1/2;
 
 end
 
-function out1 = eta0(A1, A2, B1, B2, C2, kappa, t, x, y)
+function out1 = eta0(A_2, B_2, kappa, x, y)
 
-  out1 = -0.5*tanh(sqrt(2)*(-A1.*t.*sin(B1.*x) - A2.*sin(B2.*x + C2.*t) + y - 0.25)./(2*sqrt(kappa))) + 0.5;
+  out1 = tanh(sqrt(2)*(4*A_2.*sin(B_2.*x) - 4*y + 1)./(8*sqrt(kappa)))/2 + 1/2;
 
 end
 
-function out1 = S(S)
+function out1 = S(A_1, A_2, B_1, B_2, C_2, kappa, t, x, y)
 
-  out1 = S;
+  out1 = -(2*sqrt(kappa).*(A_1.*B_1.*t.*cos(B_1.*x) + A_2.*B_2.*cos(B_2.*x + C_2.*t)).^2.*tanh(sqrt(2)*(-A_1.*t.*sin(B_1.*x) - A_2.*sin(B_2.*x + C_2.*t) + y - 1/4)./(2*sqrt(kappa))) + sqrt(2)*kappa.*(-A_1.*B_1.^2.*t.*sin(B_1.*x) - A_2.*B_2.^2.*sin(B_2.*x + C_2.*t)) - sqrt(2)*(A_1.*sin(B_1.*x) + A_2.*C_2.*cos(B_2.*x + C_2.*t))).*sech(sqrt(2)*(-A_1.*t.*sin(B_1.*x) - A_2.*sin(B_2.*x + C_2.*t) + y - 1/4)./(2*sqrt(kappa))).^2./(4*sqrt(kappa));
 
 end