From 9689fd8b8e8b46b12640914c26820a430ef6dbd4 Mon Sep 17 00:00:00 2001 From: Trevor Keller Date: Mon, 15 May 2023 18:02:19 -0400 Subject: [PATCH 1/5] Cleaner code generation for MMS source term --- benchmarks/benchmark7.ipynb | 428 ++++++++++----- benchmarks/benchmark7.ipynb.raw.html | 757 ++++++++++++++++----------- 2 files changed, 758 insertions(+), 427 deletions(-) diff --git a/benchmarks/benchmark7.ipynb b/benchmarks/benchmark7.ipynb index 56be1dc2a..44d85c0d4 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" } @@ -396,9 +396,7 @@ }, { "cell_type": "markdown", - "metadata": { - "collapsed": true - }, + "metadata": {}, "source": [ "# Appendix\n", "\n", @@ -417,59 +415,190 @@ "cell_type": "code", "execution_count": 5, "metadata": {}, + "outputs": [], + "source": [ + "# Sympy code to generate expressions for PFHub Problem 7 (MMS)\n", + "\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", + "\n", + "init_printing()\n", + "\n", + "R = ReferenceFrame(\"R\") # spatial coords\n", + "X = R[0]\n", + "Y = R[1]\n", + "\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 α symbolically to simplify the expressions for visual comparison\n", + "\n", + "α = Function(\"alpha\", real=True)(X, t)\n", + "ψ = (Y - α) / sqrt(2 * κ)\n", + "η = (1 - tanh(ψ)) / 2\n", + "\n", + "source = time_derivative(η, R) \\\n", + " + 4 * η * (η - 1) * (η - Rational(1, 2)) \\\n", + " - divergence(κ * gradient(η, R), R)\n", + "\n", + "# 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": "code", + "execution_count": 6, + "metadata": { + "tags": [] + }, "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", + "Compare generated expression for S:\n" + ] + }, + { + "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": [ + "# Print it, for comparison\n", + "print(\"Compare generated expression for S:\")\n", + "Eq(symbols(\"S\"), S.subs({X: x, Y: y}))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Expression for α:\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" + "α = A_1*t*sin(R_x*B_1) + A_2*sin(R_x*B_2 + C_2*t) + 1/4\n" ] } ], "source": [ - "# Sympy code to generate expressions for PFHub Problem 7 (MMS)\n", + "alpha = Rational(1, 4) + A1 * t * sin(B1 * X) + A2 * sin(B2 * X + C2 * t)\n", + "alpha_t = simplify(expand(diff(alpha, t)))\n", + "alpha_x = simplify(expand(diff(alpha, X)))\n", + "alpha_xx = simplify(expand(diff(alpha, X, 2)))\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.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", - "\n", - "# sinusoid amplitudes\n", - "A1, A2 = symbols('A1 A2')\n", - "B1, B2 = symbols('B1 B2')\n", - "C2 = symbols('C2')\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", + "print(\"Expression for α:\\n\")\n", + "print(\"α =\", alpha)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Expression for η:\n", + "\n", + "η = 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": [ + "eta = simplify(η.subs({R[0]: x, R[1]: y, α: alpha.subs({R[0]: x, R[1]: y})}))\n", "\n", - "# Define the solution equation (eta)\n", - "eta = 0.5 * (1 - tanh((R[1] - alpha) / sqrt(2*kappa)))\n", + "print(\"Expression for η:\\n\")\n", + "print(\"η =\", eta)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Expression for η0:\n", + "\n", + "η₀ = tanh(sqrt(2)*(4*A_2*sin(B_2*x) - 4*y + 1)/(8*sqrt(kappa)))/2 + 1/2\n" + ] + } + ], + "source": [ + "eta0 = simplify(eta.subs({t: 0, α: alpha.subs({R[0]: x, R[1]: y, t: 0})}))\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", + "print(\"Expression for η0:\\n\")\n", + "print(\"η₀ =\", eta0)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Expression for S:\n", + "\n", + "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": [ + "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", - "# 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", "\n", - "print(\"alpha =\", alpha, \"\\n\")\n", - "print(\"eta =\", eta, \"\\n\")\n", - "print(\"eta0 =\", eta0, \"\\n\")\n", + "print(\"\\nExpression for S:\\n\")\n", "print(\"S =\", source)" ] }, @@ -492,10 +621,10 @@ "```python\n", "from sympy.utilities.lambdify import lambdify\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", + "fast_alpha = lambdify([x, t], alpha, modules='sympy')\n", + "fast_eta = lambdify([x, y, t], eta, modules='sympy')\n", + "fast_eta0 = lambdify([x, y], eta0, modules='sympy')\n", + "fast_source = lambdify([x, y, t], source, modules='sympy')\n", "```\n", "\n", "> Note: Click \"Code Toggle\" at the top of the page to see the Python expressions." @@ -510,8 +639,11 @@ }, { "cell_type": "code", - "execution_count": 6, - "metadata": {}, + "execution_count": 11, + "metadata": { + "scrolled": true, + "tags": [] + }, "outputs": [ { "name": "stdout", @@ -520,7 +652,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 +663,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 R_x, double t);\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 +675,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 +684,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 R_x, double t) {\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(R_x*B_1) + A_2*sin(R_x*B_2 + 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 +723,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 +743,11 @@ }, { "cell_type": "code", - "execution_count": 7, - "metadata": {}, + "execution_count": 12, + "metadata": { + "scrolled": true, + "tags": [] + }, "outputs": [ { "name": "stdout", @@ -621,66 +756,88 @@ "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, R_x, t)\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) :: R_x\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(R_x*B_1) + A_2*sin(R_x*B_2 + 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_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", - "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_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 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" @@ -692,7 +849,7 @@ "codegen([(\"alpha\", alpha),\n", " (\"eta\", eta),\n", " (\"eta0\", eta),\n", - " (\"S\", S)],\n", + " (\"S\", source)],\n", " language=\"f95\",\n", " prefix=\"manufactured\",\n", " project=\"PFHub\")\n", @@ -710,8 +867,11 @@ }, { "cell_type": "code", - "execution_count": 8, - "metadata": {}, + "execution_count": 13, + "metadata": { + "scrolled": true, + "tags": [] + }, "outputs": [ { "name": "stdout", @@ -719,36 +879,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, R_x, t)\n", "\n", - " out1 = A1.*t.*sin(B1.*x) + A2.*sin(B2.*x + C2.*t) + 0.25\n", + " out1 = A_1 .* t .* sin(R_x .* B_1) + A_2 .* sin(R_x .* B_2 + 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_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 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", @@ -761,7 +921,7 @@ "codegen([(\"alpha\", alpha),\n", " (\"eta\", eta),\n", " (\"eta0\", eta),\n", - " (\"S\", S)],\n", + " (\"S\", source)],\n", " language=\"julia\",\n", " prefix=\"manufactured\",\n", " project=\"PFHub\")\n", @@ -779,20 +939,23 @@ }, { "cell_type": "code", - "execution_count": 9, - "metadata": {}, + "execution_count": 14, + "metadata": { + "scrolled": true, + "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[R_x*B_1] + A_2*Sin[R_x*B_2 + 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 +978,11 @@ }, { "cell_type": "code", - "execution_count": 10, - "metadata": {}, + "execution_count": 15, + "metadata": { + "scrolled": true, + "tags": [] + }, "outputs": [ { "name": "stdout", @@ -825,33 +991,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, R_x, t)\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(R_x.*B_1) + A_2.*sin(R_x.*B_2 + 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_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 = 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" @@ -863,7 +1029,7 @@ "codegen([(\"alpha\", alpha),\n", " (\"eta\", eta),\n", " (\"eta0\", eta),\n", - " (\"S\", S)],\n", + " (\"S\", source)],\n", " language=\"octave\",\n", " project=\"PFHub\")\n", "\n", @@ -876,7 +1042,7 @@ "metadata": { "anaconda-cloud": {}, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -890,7 +1056,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.9.16" }, "toc": { "base_numbering": 1, @@ -912,5 +1078,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..9b3f18d81 100644 --- a/benchmarks/benchmark7.ipynb.raw.html +++ b/benchmarks/benchmark7.ipynb.raw.html @@ -1,16 +1,15 @@ - -
+
In [1]:
%%javascript
-    MathJax.Hub.Config({
-      TeX: { equationNumbers: { autoNumber: "AMS" } }
-    });
+    MathJax.Hub.Config({
+      TeX: { equationNumbers: { autoNumber: "AMS" } }
+    });
-
+
@@ -20,19 +19,19 @@
-
+
-
-
+
@@ -42,12 +41,12 @@
-
+
In [2]:
-
from IPython.display import HTML
+
from IPython.display import HTML
 
 HTML('''<script>
 code_show=true; 
@@ -66,7 +65,7 @@
 <form action="javascript:code_toggle()"><input type="submit" value="Code Toggle"></form>''')
 

 
-
+
@@ -76,12 +75,11 @@
-
Out[2]:
+
Out[2]:
-
- -
-
+
@@ -104,12 +101,12 @@
-
+
In [3]:
-
from IPython.display import HTML
+
from IPython.display import HTML
 
 HTML('''
 <a href="https://github.com/usnistgov/pfhub/raw/master/benchmarks/benchmark7.ipynb"
@@ -119,7 +116,7 @@
 ''')
 

 
-
+
@@ -129,17 +126,15 @@
-
Out[3]:
+
Out[3]:
- -
@@ -148,25 +143,24 @@
-
-
-
+
+

Benchmark Problem 7: MMS Allen-Cahn

-
+
-
In [1]:
+
In [4]:
-
from IPython.display import HTML
+
from IPython.display import HTML
 
 HTML('''{% include jupyter_benchmark_table.html num="[7]" revision=0 %}''')
 

 
-
+
@@ -176,13 +170,11 @@

Benchmark Problem 7: MMS Allen-Cahn
-
Out[1]:
+
Out[4]:
-
-{% include jupyter_benchmark_table.html num="[7]" revision=0 %} -
+
{% include jupyter_benchmark_table.html num="[7]" revision=0 %}
@@ -190,129 +182,117 @@

Benchmark Problem 7: MMS Allen-Cahn

-
-
-
+
+
-

Table of Contents

-
+

Table of Contents

+
+
-
-
+
+
-

See the journal publication entitled "Benchmark problems for numerical implementations of phase field models" for more details about the benchmark problems. Furthermore, read the extended essay for a discussion about the need for benchmark problems.

+

See the journal publication entitled "Benchmark problems for numerical implementations of phase field models" for more details about the benchmark problems. Furthermore, read the extended essay for a discussion about the need for benchmark problems.

-
-
-
+
+

Overview

-
-
-
+
+
-

The Method of Manufactured Solutions (MMS) is a powerful technique for verifying the accuracy of a simulation code. In the MMS, one picks a desired solution to the problem at the outset, the "manufactured solution", and then determines the governing equation that will result in that solution. With the exact analytical form of the solution in hand, when the governing equation is solved using a particular simulation code, the deviation from the expected solution can be determined exactly. This deviation can be converted into an error metric to rigously quantify the error for a calculation. This error can be used to determine the order of accuracy of the simulation results to verify simulation codes. It can also be used to compare the computational efficiency of different codes or different approaches for a particular code at a certain level of error. Furthermore, the spatial/temporal distribution can give insight into the conditions resulting in the largest error (high gradients, changes in mesh resolution, etc.).

+

The Method of Manufactured Solutions (MMS) is a powerful technique for verifying the accuracy of a simulation code. In the MMS, one picks a desired solution to the problem at the outset, the "manufactured solution", and then determines the governing equation that will result in that solution. With the exact analytical form of the solution in hand, when the governing equation is solved using a particular simulation code, the deviation from the expected solution can be determined exactly. This deviation can be converted into an error metric to rigously quantify the error for a calculation. This error can be used to determine the order of accuracy of the simulation results to verify simulation codes. It can also be used to compare the computational efficiency of different codes or different approaches for a particular code at a certain level of error. Furthermore, the spatial/temporal distribution can give insight into the conditions resulting in the largest error (high gradients, changes in mesh resolution, etc.).

After choosing a manufactured solution, the governing equation must be modified to force the solution to equal the manufactured solution. This is accomplished by taking the nominal equation that is to be solved (e.g. Allen-Cahn equation, Cahn-Hilliard equation, Fick's second law, Laplace equation) and adding a source term. This source term is determined by plugging the manufactured solution into the nominal governing equation and setting the source term equal to the residual. Thus, the manufactured solution satisfies the MMS governing equation (the nominal governing equation plus the source term). A more detailed discussion of MMS can be found in the report by Salari and Knupp.

In this benchmark problem, the objective is to use the MMS to rigorously verify phase field simulation codes and then provide a basis of comparison for the computational performance between codes and for various settings for a single code, as discussed above. To this end, the benchmark problem was chosen as a balance between two factors: simplicity, to minimize the development effort required to solve the benchmark, and transferability to a real phase field system of physical interest.

-
-
-
+
+

Governing equation and manufactured solution

For this benchmark problem, we use a simple Allen-Cahn equation as the governing equation

-

$$\begin{equation} +$$\begin{equation} \frac{\partial \eta}{\partial t} = - \left[ 4 \eta \left(\eta - 1 \right) \left(\eta-\frac{1}{2} \right) - \kappa \nabla^2 \eta \right] + S(x,y,t) -\end{equation}$$

-

where $S(x,y,t)$ is the MMS source term and $\kappa$ is a constant parameter (the gradient energy coefficient).

+\end{equation}$$

where $S(x,y,t)$ is the MMS source term and $\kappa$ is a constant parameter (the gradient energy coefficient).

The manufactured solution, $\eta_{sol}$ is a hyperbolic tangent function, shifted to vary between 0 and 1, with the $x$ position of the middle of the interface ($\eta_{sol}=0.5$) given by the function $\alpha(x,t)$:

-

$$\begin{equation} +$$\begin{equation} \eta_{sol}(x,y,t) = \frac{1}{2}\left[ 1 - \tanh\left( \frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right) \right] -\end{equation}$$

-

$$\begin{equation} +\end{equation}$$$$\begin{equation} \alpha(x,t) = \frac{1}{4} + A_1 t \sin\left(B_1 x \right) + A_2 \sin \left(B_2 x + C_2 t \right) -\end{equation}$$

-

where $A_1$, $B_1$, $A_2$, $B_2$, and $C_2$ are constant parameters.

-

This manufactured solution is an equilbrium solution of the governing equation, when $S(x,y,t)=0$ and $\alpha(x,t)$ is constant. The closeness of this manufactured solution to a solution of the nominal governing equation increases the likihood that the behavior of simulation codes when solving this benchmark problem is representive of the solution of the regular Allen-Cahn equation (i.e. without the source term). The form of $\alpha(x,t)$ was chosen to yield complex behavior while still retaining a (somewhat) simple functional form. The two spatial sinusoidal terms introduce two controllable length scales to the interfacial shape. Summing them gives a "beat" pattern with a period longer than the period of either individual term, permitting a domain size that is larger than the wavelength of the sinusoids without a repeating pattern. The temporal sinusoidal term introduces a controllable time scale to the interfacial shape in addition to the phase transformation time scale, while the linear temporal dependence of the other term ensures that the sinusoidal term can go through multiple periods without $\eta_{sol}$ repeating itself.

+\end{equation}$$

where $A_1$, $B_1$, $A_2$, $B_2$, and $C_2$ are constant parameters.

+

This manufactured solution is an equilbrium solution of the governing equation, when $S(x,y,t)=0$ and $\alpha(x,t)$ is constant. The closeness of this manufactured solution to a solution of the nominal governing equation increases the likihood that the behavior of simulation codes when solving this benchmark problem is representive of the solution of the regular Allen-Cahn equation (i.e. without the source term). The form of $\alpha(x,t)$ was chosen to yield complex behavior while still retaining a (somewhat) simple functional form. The two spatial sinusoidal terms introduce two controllable length scales to the interfacial shape. Summing them gives a "beat" pattern with a period longer than the period of either individual term, permitting a domain size that is larger than the wavelength of the sinusoids without a repeating pattern. The temporal sinusoidal term introduces a controllable time scale to the interfacial shape in addition to the phase transformation time scale, while the linear temporal dependence of the other term ensures that the sinusoidal term can go through multiple periods without $\eta_{sol}$ repeating itself.

Inserting the manufactured solution into the governing equation and solving for $S(x,y,t)$ yields:

-

$$\begin{equation} +$$\begin{equation} S(x,y,t) = \frac{\text{sech}^2 \left[ \frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right]}{4 \sqrt{\kappa}} \left[-2\sqrt{\kappa} \tanh \left[\frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right] \left(\frac{\partial \alpha(x,t)}{\partial x} \right)^2+\sqrt{2} \left[ \frac{\partial \alpha(x,t)}{\partial t}-\kappa \frac{\partial^2 \alpha(x,t)}{\partial x^2} \right] \right] -\end{equation}$$

-

where $\alpha(x,t)$ is given above and where:

-

$$\begin{equation} +\end{equation}$$

where $\alpha(x,t)$ is given above and where:

+$$\begin{equation} \frac{\partial \alpha(x,t)}{\partial x} = A_1 B_1 t \cos\left(B_1 x\right) + A_2 B_2 \cos \left(B_2 x + C_2 t \right) -\end{equation}$$

-

$$\begin{equation} +\end{equation}$$$$\begin{equation} \frac{\partial^2 \alpha(x,t)}{\partial x^2} = -A_1 B_1^2 t \sin\left(B_1 x\right) - A_2 B_2^2 \sin \left(B_2 x + C_2 t \right) -\end{equation}$$

-

$$\begin{equation} +\end{equation}$$$$\begin{equation} \frac{\partial \alpha(x,t)}{\partial t} = A_1 \sin\left(B_1 x\right) + A_2 C_2 \cos \left(B_2 x + C_2 t \right) -\end{equation}$$

-

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

+\end{equation}$$

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

-
-
-
+
+

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$:

-

$$ +$$ \begin{equation} \eta_{sol}(x,y,0) = \frac{1}{2}\left[ 1 - \tanh\left( \frac{y-\left(\frac{1}{4}+A_2 \sin(B_2 x) \right)}{\sqrt{2 \kappa}} \right) \right] \end{equation} -$$

-

The stopping condition for all calculations is when t = 8 time units, which was chosen to let $\alpha(x,t)$ evolve substantially, while still being slower than the characteristic time for the phase evolution (determined by the CFL condition for a uniform mesh with a reasonable level of resolution of $\eta_{sol}$).

+$$

The stopping condition for all calculations is when t = 8 time units, which was chosen to let $\alpha(x,t)$ evolve substantially, while still being slower than the characteristic time for the phase evolution (determined by the CFL condition for a uniform mesh with a reasonable level of resolution of $\eta_{sol}$).

-
-
-
+
+

Parameter values

The nominal parameter values for the governing equation and manufactured solution are given below. The value of $\kappa$ will change in Part (b) in the following section and the values of $\kappa$ and $C_2$ will change in Part (c).

- - - + + + + - - + + - - + + - - + + - - + + - - + + - - + +
ParameterValue
ParameterValue
$\kappa$0.0004$\kappa$0.0004
$A_1$0.0075$A_1$0.0075
$B_1$$8.0 \pi$$B_1$$8.0 \pi$
$A_2$0.03$A_2$0.03
$B_2$$22.0 \pi$$B_2$$22.0 \pi$
$C_2$$0.0625 \pi$$C_2$$0.0625 \pi$
@@ -320,25 +300,22 @@

Parameter values
-
-
+
+

Benchmark simulation instructions

This section describes three sets of tests to conduct using the MMS problem specified above. The primary purpose of the first test is provide a computationally inexpensive problem to verify a simulation code. The second and third tests are more computationally demanding and are primarily designed to serve as a basis for performance comparisons.

-
-
-
+
+

Part (a)

The objective of this test is to verify the accuracy of your simulation code in both time and space. Here, we make use of convergence tests, where either the mesh size (or grid point spacing) or the time step size is systematically changed to determine the response of the error to these quantities. Once a convergence test is completed the order of accuracy can be calculated from the result. The order of accuracy can be compared to the theoretical order of accuracy for the numerical method employed in the simulation. If the two match (to a reasonable degree), then one can be confident that the simulation code is working as expected. The remainder of this subsection will give instructions for convergence tests for this MMS problem.

Implement the MMS problem specified above using the simulation code of your choice. Perform a spatial convergence test by running the simulation for a variety of mesh sizes. For each simulation, determine the discrete $L_2$ norm of the error at $t=8$:

-

$$\begin{equation} +$$\begin{equation} L_2 = \sqrt{\sum\limits_{x,y}\left(\eta^{t=8}_{x,y} - \eta_{sol}(x,y,8)\right)^2 \Delta x \Delta y} -\end{equation}$$

-

For all of these simulations, verify that the time step is small enough that any temporal error is much smaller that the total error. This can be accomplished by decreasing the time step until it has minimal effect on the error. Ensure that at least three simulation results have $L_2$ errors in the range $[5\times10^{-3}, 1\times10^{-4}]$, attempting to cover as much of that range as possible/practical. This maximum and minimum errors in the range roughly represent a poorly resolved simulation and a very well-resolved simulation.

+\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.

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$. @@ -348,18 +325,16 @@

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.

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} +$$\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).

+\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.

-
-
-
+
+

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 @@ -368,9 +343,8 @@

Part (b)

-
-
-
+
+

Part (c)

This final part is designed to stress time integrators even further by increasing the rate of change of $\alpha(x,t)$. Increase $C_2$ to $0.5$. Keep $\kappa= 1.5625\times10^{-6}$ from (b).

Repeat the process from (b), uploading the wall time, number of computing cores, processor speed, normalized computing cost, maximum memory usage, and $L_2$ error at $t=8$ to the PFHub website.

@@ -378,115 +352,127 @@

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
 0.00390625,6.26E-06
-...
+... +

and the CSV file for the temporal data should have the form

- 
time_step,L2_error
 5.00E-04,5.80162E-06
 4.00E-04,4.69709E-06
-...
-

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
 1,4.57,39400,0.010521502,5.941
-...
+... +

Only one CSV file is required with the same link in both data sections.

-
-
-
+
+
-

Results

Results from this benchmark problem are displayed on the simulation result page for different codes.

+

Results

Results from this benchmark problem are displayed on the [simulation result page]({{ site.baseurl }}/simulations) for different codes.

-
-
-
+
+
-

Feedback

Feedback on this benchmark problem is appreciated. If you have questions, comments, or seek clarification, please contact the CHiMaD phase field community through the Gitter chat channel or by email. If you found an error, please file an issue on GitHub.

+

Feedback

Feedback on this benchmark problem is appreciated. If you have questions, comments, or seek clarification, please contact the [CHiMaD phase field community]({{ site.baseurl }}/community/) through the [Gitter chat channel]({{ site.links.chat }}) or by [email]({{ site.baseurl }}/mailing_list/). If you found an error, please file an [issue on GitHub]({{ site.links.github }}/issues/new).

-
-
-
+
+

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.

-
-
-
+
+

Source term

-
+
In [5]:
# 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.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)
+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
+
+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})
+        )
+)
+
+ +
+
+
+
In [6]:
+
+
+
# Print it, for comparison
+print("Compare generated expression for S:")
+Eq(symbols("S"), S.subs({X: x, Y: y}))
+
+ +
+
@@ -495,17 +481,63 @@

Source term
-
+
-
alpha = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25 
+
Compare generated expression for S:
+

+
+
+ +
+ +
Out[6]:
+ + + + +
+$\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)
+alpha_t  = simplify(expand(diff(alpha, t)))
+alpha_x  = simplify(expand(diff(alpha, X)))
+alpha_xx = simplify(expand(diff(alpha, X, 2)))
+
+print("Expression for α:\n")
+print("α =", alpha)
+
+ +
+
+
+ +
+
+ -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 -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) +
+
Expression for α:
+
+α = A_1*t*sin(R_x*B_1) + A_2*sin(R_x*B_2 + C_2*t) + 1/4
@@ -514,51 +546,166 @@

Source term

-
+
+
+
In [8]:
+
+
+
eta = simplify(η.subs({R[0]: x, R[1]: y, α: alpha.subs({R[0]: x, R[1]: y})}))
+
+print("Expression for η:\n")
+print("η =", eta)
+
+ +
+
+ +
+
+ + +
+ +
+ + +
+
Expression for η:
+
+η = 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 [9]:
-
-

Code

+
+
eta0 = simplify(eta.subs({t: 0, α: alpha.subs({R[0]: x, R[1]: y, t: 0})}))
+
+print("Expression for η0:\n")
+print("η₀ =", eta0)
+
+ +
+ +
+
+ + +
+ +
+ + +
+
Expression for η0:
+
+η₀ = tanh(sqrt(2)*(4*A_2*sin(B_2*x) - 4*y + 1)/(8*sqrt(kappa)))/2 + 1/2
+
+
-
+ +
+
+
+
+
+
In [10]:
+
+
source = S.subs(
+    {
+        α: alpha,
+        diff(α, t):    alpha_t,
+        diff(α, X):    alpha_x,
+        diff(α, X, 2): alpha_xx
+    }
+).subs({X: x, Y: y})
+
+
+print("\nExpression for S:\n")
+print("S =", source)
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
+Expression for S:
+
+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))
+
+
+
+ +
+
+ +
+
+
+
+

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')
+fast_alpha  = lambdify([x, t],    alpha,  modules='sympy')
+fast_eta    = lambdify([x, y, t], eta,    modules='sympy')
+fast_eta0   = lambdify([x, y],    eta0,   modules='sympy')
+fast_source = lambdify([x, y, t], source, modules='sympy')
-

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

+
+

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

-
-
-
+
+

C

-
+
-
In [6]:
+
In [11]:
[(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")
@@ -568,7 +715,7 @@ 
C

 print(code)
-
+
@@ -578,14 +725,14 @@

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 +743,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 R_x, double t);
+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 +755,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 +764,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 R_x, double t) {
 
    double alpha_result;
-   alpha_result = A1*t*sin(B1*x) + A2*sin(B2*x + C2*t) + 0.25;
+   alpha_result = A_1*t*sin(R_x*B_1) + A_2*sin(R_x*B_2 + 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;
 
 }
@@ -657,24 +804,23 @@ 
C
-
-
-
+
+

Fortran

-
+
-
In [7]:
+
In [12]:
[(f_name, code), (f_name, header)] = \
 codegen([("alpha", alpha),
          ("eta",   eta),
          ("eta0",  eta),
-         ("S",     S)],
+         ("S",     source)],
          language="f95",
          prefix="manufactured",
          project="PFHub")
@@ -683,7 +829,7 @@ 
Fortran

 print(code)
-
+
@@ -693,73 +839,95 @@

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, R_x, t)
 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) :: R_x
 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(R_x*B_1) + A_2*sin(R_x*B_2 + 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_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
 
-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_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 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
 
@@ -771,24 +939,23 @@ 
Fortran
-
-
-
+
+

Julia

-
+
-
In [8]:
+
In [13]:
[(f_name, code)] = \
 codegen([("alpha", alpha),
          ("eta",   eta),
          ("eta0",  eta),
-         ("S",     S)],
+         ("S",     source)],
          language="julia",
          prefix="manufactured",
          project="PFHub")
@@ -797,7 +964,7 @@ 
Julia

 print(code)
-
+
@@ -807,42 +974,42 @@

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, R_x, t)
 
-    out1 = A1.*t.*sin(B1.*x) + A2.*sin(B2.*x + C2.*t) + 0.25
+    out1 = A_1 .* t .* sin(R_x .* B_1) + A_2 .* sin(R_x .* B_2 + 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_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 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
@@ -855,20 +1022,19 @@ 
Julia
-
-
-
+
+

Mathematica

-
+
-
In [9]:
+
In [14]:
-
from sympy.printing import mathematica_code
+
from sympy.printing import mathematica_code
 
 print("alpha =", mathematica_code(alpha), "\n")
 print("eta =", mathematica_code(eta), "\n")
@@ -876,7 +1042,7 @@ 
Mathematicaprint("S =", mathematica_code(source), "\n")
 

 
-
+
@@ -886,17 +1052,17 @@

Mathematica
-
+
-
alpha = A1*t*Sin[B1*x] + A2*Sin[B2*x + C2*t] + 0.25 
+
alpha = A_1*t*Sin[R_x*B_1] + A_2*Sin[R_x*B_2 + 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)
@@ -906,24 +1072,23 @@

Mathematica

-
-
-
+
+

Matlab

-
+
-
In [10]:
+
In [15]:
code = \
 codegen([("alpha", alpha),
          ("eta",   eta),
          ("eta0",  eta),
-         ("S",     S)],
+         ("S",     source)],
          language="octave",
          project="PFHub")
 
@@ -932,7 +1097,7 @@ 
Matlab

     print(f)
-
+
@@ -942,40 +1107,40 @@

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, R_x, t)
+  %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(R_x.*B_1) + A_2.*sin(R_x.*B_2 + 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_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 = 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
 
@@ -987,5 +1152,5 @@ 
Matlab
- + From e6ea7e7618f574cfb23ab4e7c22b55c27c1bf013 Mon Sep 17 00:00:00 2001 From: Trevor Keller Date: Tue, 16 May 2023 17:29:23 -0400 Subject: [PATCH 2/5] Improve formatting & typesetting, address notes from @wd15 --- benchmarks/benchmark7.ipynb | 111 ++++++++++++++++---------- benchmarks/benchmark7.ipynb.raw.html | 115 +++++++++++++++------------ 2 files changed, 132 insertions(+), 94 deletions(-) diff --git a/benchmarks/benchmark7.ipynb b/benchmarks/benchmark7.ipynb index 44d85c0d4..9544a1161 100644 --- a/benchmarks/benchmark7.ipynb +++ b/benchmarks/benchmark7.ipynb @@ -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", @@ -401,7 +402,7 @@ "# 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." ] }, { @@ -458,20 +459,21 @@ ")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compare generated expression for $S(x, y, \\alpha)$:" + ] + }, { "cell_type": "code", "execution_count": 6, "metadata": { + "scrolled": true, "tags": [] }, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Compare generated expression for S:\n" - ] - }, { "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==", @@ -494,22 +496,28 @@ } ], "source": [ - "# Print it, for comparison\n", - "print(\"Compare generated expression for S:\")\n", "Eq(symbols(\"S\"), S.subs({X: x, Y: y}))" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Expression for $\\alpha$:" + ] + }, { "cell_type": "code", "execution_count": 7, - "metadata": {}, + "metadata": { + "scrolled": false, + "tags": [] + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Expression for α:\n", - "\n", "α = A_1*t*sin(R_x*B_1) + A_2*sin(R_x*B_2 + C_2*t) + 1/4\n" ] } @@ -520,21 +528,28 @@ "alpha_x = simplify(expand(diff(alpha, X)))\n", "alpha_xx = simplify(expand(diff(alpha, X, 2)))\n", "\n", - "print(\"Expression for α:\\n\")\n", "print(\"α =\", alpha)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Expression for $\\eta$:" + ] + }, { "cell_type": "code", "execution_count": 8, - "metadata": {}, + "metadata": { + "scrolled": false, + "tags": [] + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Expression for η:\n", - "\n", "η = 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" ] } @@ -542,21 +557,28 @@ "source": [ "eta = simplify(η.subs({R[0]: x, R[1]: y, α: alpha.subs({R[0]: x, R[1]: y})}))\n", "\n", - "print(\"Expression for η:\\n\")\n", "print(\"η =\", eta)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Expression for $\\eta_0$:" + ] + }, { "cell_type": "code", "execution_count": 9, - "metadata": {}, + "metadata": { + "scrolled": false, + "tags": [] + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Expression for η0:\n", - "\n", "η₀ = tanh(sqrt(2)*(4*A_2*sin(B_2*x) - 4*y + 1)/(8*sqrt(kappa)))/2 + 1/2\n" ] } @@ -564,15 +586,21 @@ "source": [ "eta0 = simplify(eta.subs({t: 0, α: alpha.subs({R[0]: x, R[1]: y, t: 0})}))\n", "\n", - "print(\"Expression for η0:\\n\")\n", "print(\"η₀ =\", eta0)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Full expression for $S(x, y, t)$:" + ] + }, { "cell_type": "code", "execution_count": 10, "metadata": { - "scrolled": true, + "scrolled": false, "tags": [] }, "outputs": [ @@ -580,9 +608,6 @@ "name": "stdout", "output_type": "stream", "text": [ - "\n", - "Expression for S:\n", - "\n", "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" ] } @@ -597,8 +622,6 @@ " }\n", ").subs({X: x, Y: y})\n", "\n", - "\n", - "print(\"\\nExpression for S:\\n\")\n", "print(\"S =\", source)" ] }, @@ -941,7 +964,7 @@ "cell_type": "code", "execution_count": 14, "metadata": { - "scrolled": true, + "scrolled": false, "tags": [] }, "outputs": [ diff --git a/benchmarks/benchmark7.ipynb.raw.html b/benchmarks/benchmark7.ipynb.raw.html index 9b3f18d81..b1fb1312e 100644 --- a/benchmarks/benchmark7.ipynb.raw.html +++ b/benchmarks/benchmark7.ipynb.raw.html @@ -25,9 +25,9 @@ -
+
-
+
+
@@ -101,12 +105,12 @@
-
+
In [3]:
-
from IPython.display import HTML
+
from IPython.display import HTML
 
 HTML('''
 <a href="https://github.com/usnistgov/pfhub/raw/master/benchmarks/benchmark7.ipynb"
@@ -116,7 +120,7 @@
 ''')
 

 
-
+
@@ -126,15 +130,17 @@
-
Out[3]:
+
Out[3]:
+ +
@@ -143,24 +149,25 @@
-
-
+
+
+

Benchmark Problem 7: MMS Allen-Cahn

-
+
In [4]:
-
from IPython.display import HTML
+
from IPython.display import HTML
 
 HTML('''{% include jupyter_benchmark_table.html num="[7]" revision=0 %}''')
 

 
-
+
@@ -170,11 +177,13 @@

Benchmark Problem 7: MMS Allen-Cahn
-
Out[4]:
+
Out[4]:
-
{% include jupyter_benchmark_table.html num="[7]" revision=0 %}
+
+{% include jupyter_benchmark_table.html num="[7]" revision=0 %} +
@@ -182,119 +191,130 @@

Benchmark Problem 7: MMS Allen-Cahn

-
-
+
+
+
-

Table of Contents

-
- +

Table of Contents

-
-
+
+
+
-

See the journal publication entitled "Benchmark problems for numerical implementations of phase field models" for more details about the benchmark problems. Furthermore, read the extended essay for a discussion about the need for benchmark problems.

+

See the journal publication entitled "Benchmark problems for numerical implementations of phase field models" for more details about the benchmark problems. Furthermore, read the extended essay for a discussion about the need for benchmark problems.

-
-
+
+
+

Overview

-
-
+
+
+
-

The Method of Manufactured Solutions (MMS) is a powerful technique for verifying the accuracy of a simulation code. In the MMS, one picks a desired solution to the problem at the outset, the "manufactured solution", and then determines the governing equation that will result in that solution. With the exact analytical form of the solution in hand, when the governing equation is solved using a particular simulation code, the deviation from the expected solution can be determined exactly. This deviation can be converted into an error metric to rigously quantify the error for a calculation. This error can be used to determine the order of accuracy of the simulation results to verify simulation codes. It can also be used to compare the computational efficiency of different codes or different approaches for a particular code at a certain level of error. Furthermore, the spatial/temporal distribution can give insight into the conditions resulting in the largest error (high gradients, changes in mesh resolution, etc.).

+

The Method of Manufactured Solutions (MMS) is a powerful technique for verifying the accuracy of a simulation code. In the MMS, one picks a desired solution to the problem at the outset, the "manufactured solution", and then determines the governing equation that will result in that solution. With the exact analytical form of the solution in hand, when the governing equation is solved using a particular simulation code, the deviation from the expected solution can be determined exactly. This deviation can be converted into an error metric to rigously quantify the error for a calculation. This error can be used to determine the order of accuracy of the simulation results to verify simulation codes. It can also be used to compare the computational efficiency of different codes or different approaches for a particular code at a certain level of error. Furthermore, the spatial/temporal distribution can give insight into the conditions resulting in the largest error (high gradients, changes in mesh resolution, etc.).

After choosing a manufactured solution, the governing equation must be modified to force the solution to equal the manufactured solution. This is accomplished by taking the nominal equation that is to be solved (e.g. Allen-Cahn equation, Cahn-Hilliard equation, Fick's second law, Laplace equation) and adding a source term. This source term is determined by plugging the manufactured solution into the nominal governing equation and setting the source term equal to the residual. Thus, the manufactured solution satisfies the MMS governing equation (the nominal governing equation plus the source term). A more detailed discussion of MMS can be found in the report by Salari and Knupp.

In this benchmark problem, the objective is to use the MMS to rigorously verify phase field simulation codes and then provide a basis of comparison for the computational performance between codes and for various settings for a single code, as discussed above. To this end, the benchmark problem was chosen as a balance between two factors: simplicity, to minimize the development effort required to solve the benchmark, and transferability to a real phase field system of physical interest.

-
-
+
+
+

Governing equation and manufactured solution

For this benchmark problem, we use a simple Allen-Cahn equation as the governing equation

-$$\begin{equation} +

$$\begin{equation} \frac{\partial \eta}{\partial t} = - \left[ 4 \eta \left(\eta - 1 \right) \left(\eta-\frac{1}{2} \right) - \kappa \nabla^2 \eta \right] + S(x,y,t) -\end{equation}$$

where $S(x,y,t)$ is the MMS source term and $\kappa$ is a constant parameter (the gradient energy coefficient).

+\end{equation}$$

+

where $S(x,y,t)$ is the MMS source term and $\kappa$ is a constant parameter (the gradient energy coefficient).

The manufactured solution, $\eta_{sol}$ is a hyperbolic tangent function, shifted to vary between 0 and 1, with the $x$ position of the middle of the interface ($\eta_{sol}=0.5$) given by the function $\alpha(x,t)$:

-$$\begin{equation} +

$$\begin{equation} \eta_{sol}(x,y,t) = \frac{1}{2}\left[ 1 - \tanh\left( \frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right) \right] -\end{equation}$$$$\begin{equation} +\end{equation}$$

+

$$\begin{equation} \alpha(x,t) = \frac{1}{4} + A_1 t \sin\left(B_1 x \right) + A_2 \sin \left(B_2 x + C_2 t \right) -\end{equation}$$

where $A_1$, $B_1$, $A_2$, $B_2$, and $C_2$ are constant parameters.

-

This manufactured solution is an equilbrium solution of the governing equation, when $S(x,y,t)=0$ and $\alpha(x,t)$ is constant. The closeness of this manufactured solution to a solution of the nominal governing equation increases the likihood that the behavior of simulation codes when solving this benchmark problem is representive of the solution of the regular Allen-Cahn equation (i.e. without the source term). The form of $\alpha(x,t)$ was chosen to yield complex behavior while still retaining a (somewhat) simple functional form. The two spatial sinusoidal terms introduce two controllable length scales to the interfacial shape. Summing them gives a "beat" pattern with a period longer than the period of either individual term, permitting a domain size that is larger than the wavelength of the sinusoids without a repeating pattern. The temporal sinusoidal term introduces a controllable time scale to the interfacial shape in addition to the phase transformation time scale, while the linear temporal dependence of the other term ensures that the sinusoidal term can go through multiple periods without $\eta_{\mathrm{sol}}$ repeating itself.

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

Inserting the manufactured solution into the governing equation and solving for $S(x,y,t)$ yields:

-$$\begin{equation} +

$$\begin{equation} S(x,y,t) = \frac{\text{sech}^2 \left[ \frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right]}{4 \sqrt{\kappa}} \left[-2\sqrt{\kappa} \tanh \left[\frac{y-\alpha(x,t)}{\sqrt{2 \kappa}} \right] \left(\frac{\partial \alpha(x,t)}{\partial x} \right)^2+\sqrt{2} \left[ \frac{\partial \alpha(x,t)}{\partial t}-\kappa \frac{\partial^2 \alpha(x,t)}{\partial x^2} \right] \right] -\end{equation}$$

where $\alpha(x,t)$ is given above and where:

-$$\begin{equation} +\end{equation}$$

+

where $\alpha(x,t)$ is given above and where:

+

$$\begin{equation} \frac{\partial \alpha(x,t)}{\partial x} = A_1 B_1 t \cos\left(B_1 x\right) + A_2 B_2 \cos \left(B_2 x + C_2 t \right) -\end{equation}$$$$\begin{equation} +\end{equation}$$

+

$$\begin{equation} \frac{\partial^2 \alpha(x,t)}{\partial x^2} = -A_1 B_1^2 t \sin\left(B_1 x\right) - A_2 B_2^2 \sin \left(B_2 x + C_2 t \right) -\end{equation}$$$$\begin{equation} +\end{equation}$$

+

$$\begin{equation} \frac{\partial \alpha(x,t)}{\partial t} = A_1 \sin\left(B_1 x\right) + A_2 C_2 \cos \left(B_2 x + C_2 t \right) -\end{equation}$$

-

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

+\end{equation}$$

+

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

-
-
+
+
+

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

-
-
+
+
+

Parameter values

The nominal parameter values for the governing equation and manufactured solution are given below. The value of $\kappa$ will change in Part (b) in the following section and the values of $\kappa$ and $C_2$ will change in Part (c).

- - - - + + + - - + + - - + + - - + + - - + + - - + + - - + +
ParameterValue
ParameterValue
$\kappa$0.0004$\kappa$0.0004
$A_1$0.0075$A_1$0.0075
$B_1$$8.0 \pi$$B_1$$8.0 \pi$
$A_2$0.03$A_2$0.03
$B_2$$22.0 \pi$$B_2$$22.0 \pi$
$C_2$$0.0625 \pi$$C_2$$0.0625 \pi$
@@ -302,22 +322,25 @@

Parameter values
-
+
+
+

Benchmark simulation instructions

This section describes three sets of tests to conduct using the MMS problem specified above. The primary purpose of the first test is provide a computationally inexpensive problem to verify a simulation code. The second and third tests are more computationally demanding and are primarily designed to serve as a basis for performance comparisons.

-
-
+
+
+

Part (a)

The objective of this test is to verify the accuracy of your simulation code in both time and space. Here, we make use of convergence tests, where either the mesh size (or grid point spacing) or the time step size is systematically changed to determine the response of the error to these quantities. Once a convergence test is completed the order of accuracy can be calculated from the result. The order of accuracy can be compared to the theoretical order of accuracy for the numerical method employed in the simulation. If the two match (to a reasonable degree), then one can be confident that the simulation code is working as expected. The remainder of this subsection will give instructions for convergence tests for this MMS problem.

Implement the MMS problem specified above using the simulation code of your choice. Perform a spatial convergence test by running the simulation for a variety of mesh sizes. For each simulation, determine the discrete $L_2$ norm of the error at $t=8$:

-$$\begin{equation} +

$$\begin{equation} 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 $[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.

+\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 $[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$. @@ -327,16 +350,18 @@

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.

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} +

$$\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 $\pm 0.2$ or more from the theoretical value are to be expected (depending on the range of errors considered and other factors).

+\end{equation}$$

+

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.

-
-
+
+
+

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 @@ -345,8 +370,9 @@

Part (b)

-
-
+
+
+

Part (c)

This final part is designed to stress time integrators even further by increasing the rate of change of $\alpha(x,t)$. Increase $C_2$ to $0.5$. Keep $\kappa= 1.5625\times10^{-6}$ from (b).

Repeat the process from (b), uploading the wall time, number of computing cores, processor speed, normalized computing cost, maximum memory usage, and $L_2$ error at $t=8$ to the PFHub website.

@@ -354,77 +380,81 @@

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
 0.00390625,6.26E-06
-...
-
+...

and the CSV file for the temporal data should have the form

+ 
time_step,L2_error
 5.00E-04,5.80162E-06
 4.00E-04,4.69709E-06
-...
+...
+

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
 1,4.57,39400,0.010521502,5.941
-...
-
+...

Only one CSV file is required with the same link in both data sections.

-
-
+
+
+
-

Results

Results from this benchmark problem are displayed on the [simulation result page]({{ site.baseurl }}/simulations) for different codes.

+

Results

Results from this benchmark problem are displayed on the simulation result page for different codes.

-
-
+
+
+
-

Feedback

Feedback on this benchmark problem is appreciated. If you have questions, comments, or seek clarification, please contact the [CHiMaD phase field community]({{ site.baseurl }}/community/) through the [Gitter chat channel]({{ site.links.chat }}) or by [email]({{ site.baseurl }}/mailing_list/). If you found an error, please file an [issue on GitHub]({{ site.links.github }}/issues/new).

+

Feedback

Feedback on this benchmark problem is appreciated. If you have questions, comments, or seek clarification, please contact the CHiMaD phase field community through the Gitter chat channel or by email. If you found an error, please file an issue on GitHub.

-
-
+
+
+

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.

-
-
+
+
+

Source term

-
+
In [5]:
# Sympy code to generate expressions for PFHub Problem 7 (MMS)
 
-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 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
 
 init_printing()
 
@@ -458,20 +488,21 @@ 
Source term)
-
+
-
-
+
+
+

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

-
+
In [6]:
@@ -479,7 +510,7 @@

Source term
Eq(symbols("S"), S.subs({X: x, Y: y}))
-

+
@@ -489,7 +520,7 @@

Source term
-
Out[6]:
+
Out[6]:
@@ -504,7 +535,7 @@

Source term

-
+
In [7]:
@@ -513,7 +544,7 @@

Source termEq(symbols("alpha"), alpha.subs({X: x, Y: y}))

-
+
@@ -523,7 +554,7 @@

Source term
-
Out[7]:
+
Out[7]:
@@ -538,7 +569,7 @@

Source term

-
+
In [8]:
@@ -548,7 +579,7 @@

Source termEq(diff(α, t).subs({X: x, Y: y}), alpha_t.subs({X: x, Y: y}))

-
+
@@ -558,7 +589,7 @@

Source term
-
Out[8]:
+
Out[8]:
@@ -573,7 +604,7 @@

Source term

-
+
In [9]:
@@ -583,7 +614,7 @@

Source termEq(diff(α, X).subs({X: x, Y: y}), alpha_x.subs({X: x, Y: y}))

-
+
@@ -593,7 +624,7 @@

Source term
-
Out[9]:
+
Out[9]:
@@ -608,7 +639,7 @@

Source term

-
+
In [10]:
@@ -618,7 +649,7 @@

Source termEq(diff(α, X, X).subs({X: x, Y: y}), alpha_xx.subs({X: x, Y: y}))

-
+
@@ -628,7 +659,7 @@

Source term
-
Out[10]:
+
Out[10]:
@@ -643,7 +674,7 @@

Source term

-
+
In [11]:
@@ -651,12 +682,12 @@

Source term
eta = simplify(η.subs({R[0]: x, R[1]: y, α: alpha.subs({R[0]: x, R[1]: y})}))
-

+
-
+
In [12]:
@@ -664,12 +695,12 @@

Source term
eta0 = simplify(eta.subs({t: 0, α: alpha.subs({R[0]: x, R[1]: y, t: 0})}))
-

+
-
+
In [13]:
@@ -688,27 +719,29 @@

Source termeta0 = eta0.subs({X: x, Y: y})

-
+
-
-
+
+
+

Code

-
-
+
+
+

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]:
@@ -716,7 +749,7 @@

Python

Su 

print("S =", source)
-
+
@@ -726,7 +759,7 @@

Python

Su

-
+
@@ -739,7 +772,7 @@

Python

Su

-
+
In [15]:
@@ -747,7 +780,7 @@

Python

Su 

print("α =", alpha)
-
+
@@ -757,7 +790,7 @@

Python

Su

-
+
@@ -770,7 +803,7 @@

Python

Su

-
+
In [16]:
@@ -778,7 +811,7 @@

Python

Su 

print("η =", eta)
-
+
@@ -788,7 +821,7 @@

Python

Su

-
+
@@ -801,7 +834,7 @@

Python

Su

-
+
In [17]:
@@ -809,7 +842,7 @@

Python

Su 

print("η₀ =", eta0)
-
+
@@ -819,7 +852,7 @@

Python

Su

-
+
@@ -832,14 +865,15 @@

Python

Su

-
-
+
+
+

C

-
+
In [18]:
@@ -858,7 +892,7 @@

C

print(code)
-
+
@@ -868,7 +902,7 @@

C

-
+
@@ -947,14 +981,15 @@

C

-
-
+
+
+

Fortran

-
+
In [19]:
@@ -972,7 +1007,7 @@

Fortran

print(code)
-
+
@@ -982,7 +1017,7 @@

Fortran

-
+
@@ -1077,14 +1112,15 @@

Fortran

-
-
+
+
+

Julia

-
+
In [20]:
@@ -1102,7 +1138,7 @@

Julia

print(code)
-
+
@@ -1112,7 +1148,7 @@

Julia

-
+
@@ -1160,19 +1196,20 @@

Julia

-
-
+
+
+

Mathematica

-
+
In [21]:
-
from sympy.printing import mathematica_code
+
from sympy.printing import mathematica_code
 
 print("alpha =", mathematica_code(alpha), "\n")
 print("eta =", mathematica_code(eta), "\n")
@@ -1180,7 +1217,7 @@ 
Mathematicaprint("S =", mathematica_code(source), "\n")
 

 
-
+
@@ -1190,7 +1227,7 @@

Mathematica
-
+
@@ -1210,14 +1247,15 @@

Mathematica

-
-
+
+
+

Matlab

-
+
In [22]:
@@ -1235,7 +1273,7 @@

Matlab

print(f)
-
+
@@ -1245,7 +1283,7 @@

Matlab

-
+
@@ -1290,5 +1328,5 @@

Matlab

- +