{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Differential testing\n",
"\n",
"Test multiple different methods for implementation of finite difference differentiation.\n",
"This includes using `findiff`"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import findiff as fd\n",
"import numpy as np\n",
"import sympy as sp\n",
"from termcolor import colored\n",
"\n",
"import discretisedfield as df\n",
"import discretisedfield.util as dfu"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create testing functions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def test_function(fun, verbose=1):\n",
" print(colored(f\"Testing {fun.__name__}...\", \"blue\"))\n",
" p1 = (0.0, 0.0, 0.0)\n",
" p2 = (10.0, 10.0, 10.0)\n",
"\n",
" def f(pos):\n",
" x, y, z = pos\n",
" return (x**4, y**4, z**4)\n",
"\n",
" total = 0\n",
" success = 0\n",
"\n",
" for i in range(1, 7):\n",
" for order in range(1, 5):\n",
" total += 1\n",
" print(f\"Test {i} cell(s) in z - order {order}\", end=\" \")\n",
" mesh = df.Mesh(p1=p1, p2=p2, n=(10, 10, i))\n",
" field = df.Field(mesh, dim=3, value=f)\n",
" try:\n",
" out = fun(field, \"z\", order=order)\n",
" if np.allclose(out.array, 0):\n",
" print(colored(\"Success with assumption\", \"yellow\"))\n",
" else:\n",
" print(colored(\"Success\", \"green\"))\n",
" success += 1\n",
" except Exception as e:\n",
" if verbose:\n",
" print(colored(\"Failed\", \"red\"), end=\" \")\n",
" print(e)\n",
" else:\n",
" print(colored(\"Failed\", \"red\"))\n",
"\n",
" print(colored(f\"Success rate: {success}/{total}\", \"blue\"))\n",
"\n",
"\n",
"def test_value(fun, verbose=1):\n",
" # Test results against sympy\n",
" print(colored(f\"Testing {fun.__name__}...\", \"blue\"))\n",
" p1 = (0, 0, 0)\n",
" p2 = (10, 10, 10)\n",
"\n",
" x = sp.symbols(\"x\")\n",
" fx = sp.sqrt(x)\n",
"\n",
" lam_f = sp.lambdify(x, fx, \"numpy\")\n",
"\n",
" def value_fun(point):\n",
" x, y, z = point\n",
" return lam_f(x)\n",
"\n",
" total = 0\n",
" success = 0\n",
"\n",
" for i in range(1, 5):\n",
" for order in range(1, 5):\n",
" value = float(fx.diff(x, order).subs(x, 5))\n",
" total += 1\n",
" print(f\"Test {i} cell(s) in x - order {order}\", end=\" \")\n",
" mesh = df.Mesh(p1=p1, p2=p2, n=(i, 1, 1))\n",
" field = df.Field(mesh, dim=1, value=value_fun)\n",
" try:\n",
" fun_value = fun(field, \"x\", order=order)((5, 5, 5))\n",
" if np.allclose(fun_value, value, rtol=0.05):\n",
" if verbose:\n",
" print(colored(\"Success\", \"green\"), end=\" \")\n",
" print(f\"Expected {value}, got {fun_value}\")\n",
" else:\n",
" print(colored(\"Success\", \"green\"))\n",
" success += 1\n",
" else:\n",
" if verbose:\n",
" print(colored(\"Failed\", \"red\"), end=\" \")\n",
" print(f\"Expected {value}, got {fun_value}\")\n",
" else:\n",
" print(colored(\"Failed\", \"red\"))\n",
" except:\n",
" print(\"\")\n",
"\n",
" print(colored(f\"Success rate: {success}/{total}\", \"blue\"))\n",
"\n",
"\n",
"def test_accuracy(fun, verbose=1):\n",
" print(colored(f\"Testing {fun.__name__}...\", \"blue\"))\n",
" # Test results against sympy\n",
" p1 = (0, 0, 0)\n",
" p2 = (10, 10, 10)\n",
" n = (30, 1, 1)\n",
" mesh = df.Mesh(p1=p1, p2=p2, n=n)\n",
"\n",
" x = sp.symbols(\"x\")\n",
" fx = x**6\n",
"\n",
" lam_f = sp.lambdify(x, fx, \"numpy\")\n",
"\n",
" def value_fun(point):\n",
" x, y, z = point\n",
" return lam_f(x)\n",
"\n",
" field = df.Field(mesh, nvdim=1, value=value_fun)\n",
"\n",
" total = 0\n",
" success = 0\n",
"\n",
" for order in [1, 2, 3, 4, 5]:\n",
" value = float(fx.diff(x, order).integrate((x, 0, 10)) / 10)\n",
" for acc in [2, 4, 6, 8, 10]:\n",
" print(f\"Test {order=}, {acc=}\", end=\" \")\n",
" total += 1\n",
" try:\n",
" fun_value = fun(field, \"x\", order=order).mean()\n",
" if np.allclose(fun_value, value, rtol=0.05):\n",
" if verbose:\n",
" print(colored(\"Success\", \"green\"), end=\" \")\n",
" print(f\"Expected {value}, got {fun_value}\")\n",
" else:\n",
" print(colored(\"Success\", \"green\"))\n",
" success += 1\n",
" else:\n",
" if verbose:\n",
" print(colored(\"Failed\", \"red\"), end=\" \")\n",
" print(f\"Expected {value}, got {fun_value}\")\n",
" else:\n",
" print(colored(\"Failed\", \"red\"))\n",
" except:\n",
" print(\"\")\n",
" print(colored(f\"Success rate: {success}/{total}\", \"blue\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Original df.Field differential"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def diff_original(self, direction, order=1):\n",
" direction = dfu.axesdict[direction]\n",
"\n",
" # If there are no neighbouring cells in the specified direction, zero\n",
" # field is returned.\n",
" if self.mesh.n[direction] <= 1:\n",
" return self.zero\n",
"\n",
" padded_array = self.array\n",
"\n",
" if order not in (1, 2):\n",
" msg = f\"Derivative of the n={order} order is not implemented.\"\n",
" raise NotImplementedError(msg)\n",
"\n",
" elif order == 1:\n",
" if self.nvdim == 1:\n",
" derivative_array = np.gradient(\n",
" padded_array[..., 0], self.mesh.cell[direction], axis=direction\n",
" )[..., np.newaxis]\n",
" else:\n",
" derivative_array = np.gradient(\n",
" padded_array, self.mesh.cell[direction], axis=direction\n",
" )\n",
"\n",
" elif order == 2:\n",
" derivative_array = np.zeros_like(padded_array)\n",
" for i in range(padded_array.shape[direction]):\n",
" if i == 0:\n",
" i1, i2, i3 = i + 2, i + 1, i\n",
" elif i == padded_array.shape[direction] - 1:\n",
" i1, i2, i3 = i, i - 1, i - 2\n",
" else:\n",
" i1, i2, i3 = i + 1, i, i - 1\n",
" index1 = dfu.assemble_index(slice(None), 4, {direction: i1})\n",
" index2 = dfu.assemble_index(slice(None), 4, {direction: i2})\n",
" index3 = dfu.assemble_index(slice(None), 4, {direction: i3})\n",
" index = dfu.assemble_index(slice(None), 4, {direction: i})\n",
" derivative_array[index] = (\n",
" padded_array[index1] - 2 * padded_array[index2] + padded_array[index3]\n",
" ) / self.mesh.cell[direction] ** 2\n",
"\n",
" return self.__class__(\n",
" self.mesh, dim=self.nvdim, value=derivative_array, vdims=self.vdims\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[34mTesting diff_original...\u001b[0m\n",
"Test 1 cell(s) in z - order 1 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 2 cell(s) in z - order 2 \u001b[31mFailed\u001b[0m index 2 is out of bounds for axis 2 with size 2\n",
"Test 2 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 2 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"Test 3 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 3 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"Test 4 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 4 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"Test 5 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 5 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"Test 6 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 6 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"\u001b[34mSuccess rate: 13/24\u001b[0m\n"
]
}
],
"source": [
"test_function(diff_original)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Edit Original df.Field differential\n",
"Accounts for `order=2` and number of cells `2`."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def diff_original_edit(self, direction, order=1):\n",
" direction_idx = dfu.axesdict[direction]\n",
"\n",
" # If there are no neighbouring cells in the specified direction, zero\n",
" # field is returned.\n",
" if self.mesh.n[direction_idx] <= order:\n",
" return self.zero\n",
"\n",
" padded_array = self.array\n",
"\n",
" if order not in (1, 2):\n",
" msg = f\"Derivative of the n={order} order is not implemented.\"\n",
" raise NotImplementedError(msg)\n",
"\n",
" elif order == 1:\n",
" if self.nvdim == 1:\n",
" derivative_array = np.gradient(\n",
" padded_array[..., 0], self.mesh.cell[direction_idx], axis=direction_idx\n",
" )[..., np.newaxis]\n",
" else:\n",
" derivative_array = np.gradient(\n",
" padded_array, self.mesh.cell[direction_idx], axis=direction_idx\n",
" )\n",
"\n",
" elif order == 2:\n",
" derivative_array = np.zeros_like(padded_array)\n",
" for i in range(padded_array.shape[direction_idx]):\n",
" if i == 0:\n",
" i1, i2, i3 = i + 2, i + 1, i\n",
" elif i == padded_array.shape[direction_idx] - 1:\n",
" i1, i2, i3 = i, i - 1, i - 2\n",
" else:\n",
" i1, i2, i3 = i + 1, i, i - 1\n",
" index1 = dfu.assemble_index(slice(None), 4, {direction_idx: i1})\n",
" index2 = dfu.assemble_index(slice(None), 4, {direction_idx: i2})\n",
" index3 = dfu.assemble_index(slice(None), 4, {direction_idx: i3})\n",
" index = dfu.assemble_index(slice(None), 4, {direction_idx: i})\n",
" derivative_array[index] = (\n",
" padded_array[index1] - 2 * padded_array[index2] + padded_array[index3]\n",
" ) / self.mesh.cell[direction_idx] ** 2\n",
"\n",
" return self.__class__(\n",
" self.mesh, dim=self.nvdim, value=derivative_array, vdims=self.vdims\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[34mTesting diff_original_edit...\u001b[0m\n",
"Test 1 cell(s) in z - order 1 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 2 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 4 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 4 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 5 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 5 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"Test 6 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 6 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"\u001b[34mSuccess rate: 19/24\u001b[0m\n"
]
}
],
"source": [
"test_function(diff_original_edit)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Edit Fast Original df.Field differential\n",
"Accounts for `order=2` and number of cells `2` and faster implementation of `order=2`.\n",
"\n",
"Pad array so the same stencil can be used everywhere\n",
"\\begin{align}\n",
" f(1) + f(-1) - 2f(0) &= f(2) + f(0) - 2 f(1) \\\\ \n",
" f(-1) &= f(2) - 3 f(1) + 3f(0)\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def diff_original_fast(self, direction, order=1):\n",
" direction_idx = dfu.axesdict[direction]\n",
"\n",
" # If there are no neighbouring cells in the specified direction, zero\n",
" # field is returned.\n",
" if self.mesh.n[direction_idx] <= order:\n",
" return self.zero\n",
"\n",
" padded_array = self.array\n",
"\n",
" if order not in (1, 2):\n",
" msg = f\"Derivative of the n={order} order is not implemented.\"\n",
" raise NotImplementedError(msg)\n",
"\n",
" elif order == 1:\n",
" if self.nvdim == 1:\n",
" derivative_array = np.gradient(\n",
" padded_array[..., 0], self.mesh.cell[direction_idx], axis=direction_idx\n",
" )[..., np.newaxis]\n",
" else:\n",
" derivative_array = np.gradient(\n",
" padded_array, self.mesh.cell[direction_idx], axis=direction_idx\n",
" )\n",
"\n",
" elif order == 2:\n",
"\n",
" def pad_fun(vector, pad_width, iaxis, kwargs):\n",
" if iaxis == direction_idx:\n",
" vector[0] = vector[3] - 3 * vector[2] + 3 * vector[1]\n",
" vector[-1] = vector[-4] - 3 * vector[-3] + 3 * vector[-2]\n",
"\n",
" pad_width = [(0, 0)] * 4\n",
" pad_width[direction_idx] = (1, 1)\n",
" padded_array = np.pad(padded_array, pad_width, pad_fun)\n",
" # padded_array = self.pad({'x': (1, 1)}, mode='constant').array\n",
" # padded_array[0] = padded_array[3] - 3*padded_array[2] + 3*padded_array[1]\n",
" # padded_array[-1] = padded_array[-4] - 3*padded_array[-3] + 3*padded_array[-2]\n",
" # derivative_array = np.empty_like(padded_array)\n",
" index_p1 = slice(2, None)\n",
" index_0 = slice(1, -1)\n",
" index_m1 = slice(None, -2)\n",
" derivative_array = (\n",
" padded_array[index_p1] - 2 * padded_array[index_0] + padded_array[index_m1]\n",
" ) / self.mesh.cell[direction_idx] ** 2\n",
"\n",
" # derivative_array[0] = derivative_array[1]\n",
" # derivative_array[-1] = derivative_array[-2]\n",
"\n",
" return self.__class__(\n",
" self.mesh, dim=self.nvdim, value=derivative_array, vdims=self.vdims\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'field' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m/home/sam/repos/ubermag-devtools/repos/discretisedfield/dev/diff.ipynb Cell 13\u001b[0m in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> <a href='vscode-notebook-cell://ssh-remote%2B10.15.40.246/home/sam/repos/ubermag-devtools/repos/discretisedfield/dev/diff.ipynb#X50sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a>\u001b[0m get_ipython()\u001b[39m.\u001b[39;49mrun_cell_magic(\u001b[39m'\u001b[39;49m\u001b[39mtimeit\u001b[39;49m\u001b[39m'\u001b[39;49m, \u001b[39m'\u001b[39;49m\u001b[39m'\u001b[39;49m, \u001b[39m\"\u001b[39;49m\u001b[39mdiff_original_fast(field, \u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39mx\u001b[39;49m\u001b[39m'\u001b[39;49m\u001b[39m, order=2)\u001b[39;49m\u001b[39m\\n\u001b[39;49;00m\u001b[39m\"\u001b[39;49m)\n",
"File \u001b[0;32m~/miniconda3/envs/ubermagdev/lib/python3.8/site-packages/IPython/core/interactiveshell.py:2358\u001b[0m, in \u001b[0;36mInteractiveShell.run_cell_magic\u001b[0;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[1;32m 2356\u001b[0m \u001b[39mwith\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mbuiltin_trap:\n\u001b[1;32m 2357\u001b[0m args \u001b[39m=\u001b[39m (magic_arg_s, cell)\n\u001b[0;32m-> 2358\u001b[0m result \u001b[39m=\u001b[39m fn(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m 2359\u001b[0m \u001b[39mreturn\u001b[39;00m result\n",
"File \u001b[0;32m~/miniconda3/envs/ubermagdev/lib/python3.8/site-packages/IPython/core/magics/execution.py:1162\u001b[0m, in \u001b[0;36mExecutionMagics.timeit\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n\u001b[1;32m 1160\u001b[0m \u001b[39mfor\u001b[39;00m index \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(\u001b[39m0\u001b[39m, \u001b[39m10\u001b[39m):\n\u001b[1;32m 1161\u001b[0m number \u001b[39m=\u001b[39m \u001b[39m10\u001b[39m \u001b[39m*\u001b[39m\u001b[39m*\u001b[39m index\n\u001b[0;32m-> 1162\u001b[0m time_number \u001b[39m=\u001b[39m timer\u001b[39m.\u001b[39;49mtimeit(number)\n\u001b[1;32m 1163\u001b[0m \u001b[39mif\u001b[39;00m time_number \u001b[39m>\u001b[39m\u001b[39m=\u001b[39m \u001b[39m0.2\u001b[39m:\n\u001b[1;32m 1164\u001b[0m \u001b[39mbreak\u001b[39;00m\n",
"File \u001b[0;32m~/miniconda3/envs/ubermagdev/lib/python3.8/site-packages/IPython/core/magics/execution.py:156\u001b[0m, in \u001b[0;36mTimer.timeit\u001b[0;34m(self, number)\u001b[0m\n\u001b[1;32m 154\u001b[0m gc\u001b[39m.\u001b[39mdisable()\n\u001b[1;32m 155\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m--> 156\u001b[0m timing \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49minner(it, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mtimer)\n\u001b[1;32m 157\u001b[0m \u001b[39mfinally\u001b[39;00m:\n\u001b[1;32m 158\u001b[0m \u001b[39mif\u001b[39;00m gcold:\n",
"File \u001b[0;32m<magic-timeit>:1\u001b[0m, in \u001b[0;36minner\u001b[0;34m(_it, _timer)\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'field' is not defined"
]
}
],
"source": [
"%%timeit\n",
"diff_original_fast(field, \"x\", order=2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"5.02 ms ± 19.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"%%timeit\n",
"diff_original_edit(field, \"x\", order=2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[34mTesting diff_original_fast...\u001b[0m\n",
"Test 1 cell(s) in z - order 1 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 2 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 4 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 4 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 5 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 5 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"Test 6 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Derivative of the n=3 order is not implemented.\n",
"Test 6 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Derivative of the n=4 order is not implemented.\n",
"\u001b[34mSuccess rate: 19/24\u001b[0m\n"
]
}
],
"source": [
"test_function(diff_original_fast)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Findiff implementation"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def diff_a(self, direction, order=1, acc=None):\n",
" direction_idx = dfu.axesdict[direction]\n",
"\n",
" # If there are no neighbouring cells in the specified direction, zero\n",
" # field is returned.\n",
"\n",
" # ASSUME that needs more cells than curvature!\n",
" if self.mesh.n[direction_idx] <= order:\n",
" return self.zero\n",
"\n",
" # Create FinDiff object.\n",
" diff_fd = fd.FinDiff(direction_idx, self.mesh.cell[direction_idx], order, acc=acc)\n",
"\n",
" derivative_array = diff_fd(self.array)\n",
"\n",
" return self.__class__(\n",
" self.mesh,\n",
" nvdim=self.nvdim,\n",
" value=derivative_array,\n",
" vdims=self.vdims,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[34mTesting diff_a...\u001b[0m\n",
"Test 1 cell(s) in z - order 1 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 1 \u001b[31mFailed\u001b[0m Shift slice out of bounds\n",
"Test 2 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 2 \u001b[31mFailed\u001b[0m Shift slice out of bounds\n",
"Test 3 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 4 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Shift slice out of bounds\n",
"Test 4 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 5 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Shift slice out of bounds\n",
"Test 5 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Shift slice out of bounds\n",
"Test 6 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 3 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Shift slice out of bounds\n",
"\u001b[34mSuccess rate: 18/24\u001b[0m\n"
]
}
],
"source": [
"test_function(diff_a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Findiff implementation\n",
"Allowing for mixed accuracy default behaviour"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"def fd_diff_edit_b(array, h, order, dim, acc=None):\n",
" \"\"\"Edit of findiff.diff.diff to allow differnt accuracies\"\"\"\n",
" if acc is None:\n",
" coefs = {}\n",
" coefs[\"center\"] = fd.coefficients(order, acc=2)[\"center\"]\n",
" coefs[\"forward\"] = fd.coefficients(order, offsets=range(0, order + 1))\n",
" coefs[\"backward\"] = fd.coefficients(order, offsets=range(-order, 1))\n",
" else:\n",
" coefs = fd.coefficients(order, acc=acc)\n",
"\n",
" try:\n",
" npts = array.shape[dim]\n",
" except AttributeError as err:\n",
" raise ValueError(\n",
" \"FinDiff objects can only be applied to arrays or evaluated(!) functions\"\n",
" \" returning arrays\"\n",
" ) from err\n",
" yd = np.zeros_like(array)\n",
" scheme = \"center\"\n",
" weights = coefs[scheme][\"coefficients\"]\n",
" offsets = coefs[scheme][\"offsets\"]\n",
"\n",
" num_bndry_points = len(weights) // 2\n",
" ref_slice = slice(num_bndry_points, npts - num_bndry_points, 1)\n",
" off_slices = [\n",
" fd.diff.Diff._shift_slice(..., ref_slice, offsets[k], npts)\n",
" for k in range(len(offsets))\n",
" ]\n",
"\n",
" fd.diff.Diff._apply_to_array(..., yd, array, weights, off_slices, ref_slice, dim)\n",
"\n",
" scheme = \"forward\"\n",
" weights = coefs[scheme][\"coefficients\"]\n",
" offsets = coefs[scheme][\"offsets\"]\n",
"\n",
" ref_slice = slice(0, num_bndry_points, 1)\n",
" off_slices = [\n",
" fd.diff.Diff._shift_slice(..., ref_slice, offsets[k], npts)\n",
" for k in range(len(offsets))\n",
" ]\n",
"\n",
" fd.diff.Diff._apply_to_array(..., yd, array, weights, off_slices, ref_slice, dim)\n",
"\n",
" scheme = \"backward\"\n",
" weights = coefs[scheme][\"coefficients\"]\n",
" offsets = coefs[scheme][\"offsets\"]\n",
"\n",
" ref_slice = slice(npts - num_bndry_points, npts, 1)\n",
" off_slices = [\n",
" fd.diff.Diff._shift_slice(..., ref_slice, offsets[k], npts)\n",
" for k in range(len(offsets))\n",
" ]\n",
"\n",
" fd.diff.Diff._apply_to_array(..., yd, array, weights, off_slices, ref_slice, dim)\n",
"\n",
" h_inv = 1.0 / h**order\n",
"\n",
" return yd * h_inv"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def diff_b(self, direction, order=1, acc=None):\n",
" direction_idx = dfu.axesdict[direction]\n",
"\n",
" # If there are no neighbouring cells in the specified direction, zero\n",
" # field is returned.\n",
"\n",
" # ASSUME that needs more cells than curvature!\n",
" if self.mesh.n[direction_idx] <= order:\n",
" return self.zero\n",
"\n",
" # Create FinDiff object.\n",
" derivative_array = fd_diff_edit_b(\n",
" self.array, self.mesh.cell[direction_idx], order, direction_idx, acc=acc\n",
" )\n",
"\n",
" return self.__class__(\n",
" self.mesh,\n",
" nvdim=self.nvdim,\n",
" value=derivative_array,\n",
" vdims=self.vdims,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[34mTesting diff_b...\u001b[0m\n",
"Test 1 cell(s) in z - order 1 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 2 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 4 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m Shift slice out of bounds\n",
"Test 4 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 5 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 3 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m Shift slice out of bounds\n",
"Test 6 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 3 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 4 \u001b[32mSuccess\u001b[0m\n",
"\u001b[34mSuccess rate: 22/24\u001b[0m\n"
]
}
],
"source": [
"test_function(diff_b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Findiff implementation\n",
"Generalise defaults to work for all orders"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def fd_diff_edit_c(array, h, order, dim, acc=None):\n",
" \"\"\"Edit of findiff.diff.diff to allow differnt accuracies\"\"\"\n",
" if acc is None:\n",
" coefs = {}\n",
" coefs[\"center\"] = fd.coefficients(order, acc=2)[\"center\"]\n",
" coefs[\"forward\"] = fd.coefficients(order, offsets=range(0, order + 1))\n",
" coefs[\"backward\"] = fd.coefficients(order, offsets=range(-order, 1))\n",
" else:\n",
" coefs = fd.coefficients(order, acc=acc)\n",
"\n",
" try:\n",
" npts = array.shape[dim]\n",
" except AttributeError as err:\n",
" raise ValueError(\n",
" \"FinDiff objects can only be applied to arrays or evaluated(!) functions\"\n",
" \" returning arrays\"\n",
" ) from err\n",
"\n",
" yd = np.zeros_like(array)\n",
"\n",
" len_needed = len(coefs[\"center\"][\"coefficients\"]) // 2 - 1 + order\n",
" if npts <= len_needed:\n",
" for i in range(npts):\n",
" coefs = fd.coefficients(order, offsets=(np.arange(npts) - i).tolist())\n",
" weights = coefs[\"coefficients\"]\n",
" offsets = coefs[\"offsets\"]\n",
"\n",
" ref_slice = slice(i, i + 1, 1)\n",
" off_slices = [\n",
" fd.diff.Diff._shift_slice(..., ref_slice, offsets[k], npts)\n",
" for k in range(len(offsets))\n",
" ]\n",
" fd.diff.Diff._apply_to_array(\n",
" ..., yd, array, weights, off_slices, ref_slice, dim\n",
" )\n",
"\n",
" else:\n",
" scheme = \"center\"\n",
" weights = coefs[scheme][\"coefficients\"]\n",
" offsets = coefs[scheme][\"offsets\"]\n",
"\n",
" num_bndry_points = len(weights) // 2\n",
"\n",
" ref_slice = slice(num_bndry_points, npts - num_bndry_points, 1)\n",
" off_slices = [\n",
" fd.diff.Diff._shift_slice(..., ref_slice, offsets[k], npts)\n",
" for k in range(len(offsets))\n",
" ]\n",
"\n",
" fd.diff.Diff._apply_to_array(\n",
" ..., yd, array, weights, off_slices, ref_slice, dim\n",
" )\n",
"\n",
" scheme = \"forward\"\n",
" weights = coefs[scheme][\"coefficients\"]\n",
" offsets = coefs[scheme][\"offsets\"]\n",
"\n",
" ref_slice = slice(0, num_bndry_points, 1)\n",
" off_slices = [\n",
" fd.diff.Diff._shift_slice(..., ref_slice, offsets[k], npts)\n",
" for k in range(len(offsets))\n",
" ]\n",
"\n",
" fd.diff.Diff._apply_to_array(\n",
" ..., yd, array, weights, off_slices, ref_slice, dim\n",
" )\n",
"\n",
" scheme = \"backward\"\n",
" weights = coefs[scheme][\"coefficients\"]\n",
" offsets = coefs[scheme][\"offsets\"]\n",
"\n",
" ref_slice = slice(npts - num_bndry_points, npts, 1)\n",
" off_slices = [\n",
" fd.diff.Diff._shift_slice(..., ref_slice, offsets[k], npts)\n",
" for k in range(len(offsets))\n",
" ]\n",
"\n",
" fd.diff.Diff._apply_to_array(\n",
" ..., yd, array, weights, off_slices, ref_slice, dim\n",
" )\n",
"\n",
" h_inv = 1.0 / h**order\n",
"\n",
" return yd * h_inv"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def diff_c(self, direction, order=1, acc=None):\n",
" direction_idx = dfu.axesdict[direction]\n",
"\n",
" # If there are no neighbouring cells in the specified direction, zero\n",
" # field is returned.\n",
"\n",
" # ASSUME that needs more cells than curvature!\n",
" if self.mesh.n[direction_idx] <= order:\n",
" return self.zero\n",
"\n",
" # Create FinDiff object.\n",
" derivative_array = fd_diff_edit_c(\n",
" self.array, self.mesh.cell[direction_idx], order, direction_idx, acc=acc\n",
" )\n",
"\n",
" return self.__class__(\n",
" self.mesh,\n",
" nvdim=self.nvdim,\n",
" value=derivative_array,\n",
" vdims=self.vdims,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[34mTesting diff_c...\u001b[0m\n",
"Test 1 cell(s) in z - order 1 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 2 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 4 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 3 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 5 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 3 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 4 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 3 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 4 \u001b[32mSuccess\u001b[0m\n",
"\u001b[34mSuccess rate: 24/24\u001b[0m\n"
]
}
],
"source": [
"test_function(diff_c)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mixed implementation"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def diff_d(self, direction, order=1, acc=None):\n",
" direction_idx = dfu.axesdict[direction]\n",
"\n",
" # If there are no neighbouring cells in the specified direction, zero\n",
" # field is returned.\n",
" if self.mesh.n[direction_idx] <= order:\n",
" return self.zero\n",
"\n",
" padded_array = self.array\n",
"\n",
" # Only use our implimentation if we have to\n",
" if order in (1, 2) and acc is None:\n",
" if order == 1:\n",
" if self.nvdim == 1:\n",
" derivative_array = np.gradient(\n",
" padded_array[..., 0],\n",
" self.mesh.cell[direction_idx],\n",
" axis=direction_idx,\n",
" )[..., np.newaxis]\n",
" else:\n",
" derivative_array = np.gradient(\n",
" padded_array, self.mesh.cell[direction_idx], axis=direction_idx\n",
" )\n",
"\n",
" elif order == 2:\n",
" derivative_array = np.zeros_like(padded_array)\n",
" for i in range(padded_array.shape[direction_idx]):\n",
" if i == 0:\n",
" i1, i2, i3 = i + 2, i + 1, i\n",
" elif i == padded_array.shape[direction_idx] - 1:\n",
" i1, i2, i3 = i, i - 1, i - 2\n",
" else:\n",
" i1, i2, i3 = i + 1, i, i - 1\n",
" index1 = dfu.assemble_index(slice(None), 4, {direction_idx: i1})\n",
" index2 = dfu.assemble_index(slice(None), 4, {direction_idx: i2})\n",
" index3 = dfu.assemble_index(slice(None), 4, {direction_idx: i3})\n",
" index = dfu.assemble_index(slice(None), 4, {direction_idx: i})\n",
" derivative_array[index] = (\n",
" padded_array[index1]\n",
" - 2 * padded_array[index2]\n",
" + padded_array[index3]\n",
" ) / self.mesh.cell[direction_idx] ** 2\n",
" else:\n",
" if acc is None:\n",
" acc = 2\n",
" coeffs = fd.coefficients(order, acc=acc)\n",
" stencil_len_back = len(coeffs[\"backward\"][\"offsets\"]) - 1\n",
" stencil_max_cent = max(coeffs[\"center\"][\"offsets\"]) - 1\n",
"\n",
" len_needed = stencil_len_back + stencil_max_cent\n",
" if self.mesh.n[direction_idx] <= len_needed:\n",
" raise ValueError(\n",
" f\"The minimum size of the mesh in the direction {direction} \"\n",
" f\"is {len_needed+1} for an order order {order} and acc {acc}.\"\n",
" )\n",
" diff_fd = fd.FinDiff(\n",
" direction_idx, self.mesh.cell[direction_idx], order, acc=acc\n",
" )\n",
" derivative_array = diff_fd(padded_array)\n",
"\n",
" return self.__class__(\n",
" self.mesh, dim=self.nvdim, value=derivative_array, vdims=self.vdims\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[34mTesting diff_d...\u001b[0m\n",
"Test 1 cell(s) in z - order 1 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 1 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 2 cell(s) in z - order 2 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 2 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 3 cell(s) in z - order 3 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 3 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 4 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 4 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m The minimum size of the mesh in the direction z is 6 for an order order 3 and acc 2.\n",
"Test 4 cell(s) in z - order 4 \u001b[33mSuccess with assumption\u001b[0m\n",
"Test 5 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 5 cell(s) in z - order 3 \u001b[31mFailed\u001b[0m The minimum size of the mesh in the direction z is 6 for an order order 3 and acc 2.\n",
"Test 5 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m The minimum size of the mesh in the direction z is 7 for an order order 4 and acc 2.\n",
"Test 6 cell(s) in z - order 1 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 2 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 3 \u001b[32mSuccess\u001b[0m\n",
"Test 6 cell(s) in z - order 4 \u001b[31mFailed\u001b[0m The minimum size of the mesh in the direction z is 7 for an order order 4 and acc 2.\n",
"\u001b[34mSuccess rate: 20/24\u001b[0m\n"
]
}
],
"source": [
"test_function(diff_d)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[34mTesting diff_d...\u001b[0m\n",
"Test 1 cell(s) in x - order 1 \u001b[31mFailed\u001b[0m Expected 0.22360679774997896, got 0.0\n",
"Test 1 cell(s) in x - order 2 \u001b[31mFailed\u001b[0m Expected -0.022360679774997897, got 0.0\n",
"Test 1 cell(s) in x - order 3 \u001b[31mFailed\u001b[0m Expected 0.006708203932499369, got 0.0\n",
"Test 1 cell(s) in x - order 4 \u001b[31mFailed\u001b[0m Expected -0.0033541019662496844, got 0.0\n",
"Test 2 cell(s) in x - order 1 \u001b[32mSuccess\u001b[0m Expected 0.22360679774997896, got 0.23149479148832816\n",
"Test 2 cell(s) in x - order 2 \u001b[31mFailed\u001b[0m Expected -0.022360679774997897, got 0.0\n",
"Test 2 cell(s) in x - order 3 \u001b[31mFailed\u001b[0m Expected 0.006708203932499369, got 0.0\n",
"Test 2 cell(s) in x - order 4 \u001b[31mFailed\u001b[0m Expected -0.0033541019662496844, got 0.0\n",
"Test 3 cell(s) in x - order 1 \u001b[31mFailed\u001b[0m Expected 0.22360679774997896, got 0.2393635345818485\n",
"Test 3 cell(s) in x - order 2 \u001b[31mFailed\u001b[0m Expected -0.022360679774997897, got -0.02649511442840804\n",
"Test 3 cell(s) in x - order 3 \u001b[31mFailed\u001b[0m Expected 0.006708203932499369, got 0.0\n",
"Test 3 cell(s) in x - order 4 \u001b[31mFailed\u001b[0m Expected -0.0033541019662496844, got 0.0\n",
"Test 4 cell(s) in x - order 1 \u001b[31mFailed\u001b[0m Expected 0.22360679774997896, got 0.20430964368921992\n",
"Test 4 cell(s) in x - order 2 \u001b[31mFailed\u001b[0m Expected -0.022360679774997897, got -0.016874949655437347\n",
"Test 4 cell(s) in x - order 3 \n",
"Test 4 cell(s) in x - order 4 \u001b[31mFailed\u001b[0m Expected -0.0033541019662496844, got 0.0\n",
"\u001b[34mSuccess rate: 1/16\u001b[0m\n"
]
}
],
"source": [
"test_value(diff_d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ndimage implementation"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"from scipy import ndimage"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"p1 = (0.0, 0.0, 0.0)\n",
"p2 = (10.0, 10.0, 10.0)\n",
"\n",
"\n",
"def f(pos):\n",
" x, y, z = pos\n",
" return (x**4, y**4, z**4)\n",
"\n",
"\n",
"mesh = df.Mesh(p1=p1, p2=p2, n=(10, 10, 20), bc=\"xyz\")\n",
"field = df.Field(mesh, dim=3, value=f)"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-4070., 0., 0.])"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"field.diff(\"x\").array[0, 5, 5]"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [],
"source": [
"k = np.zeros((3, 1, 1, 1))\n",
"k[0, 0, 0, 0] = 0.5\n",
"k[2, 0, 0, 0] = -0.5"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-4070., 0., 0.])"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(\n",
" ndimage.convolve(field.array, k, mode=\"wrap\", origin=(0, 0, 0, 0))\n",
" / field.mesh.cell[0]\n",
")[0, 5, 5]"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.8.13 ('ubermagdev')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.13"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "3ca9094a041cc4b421461943860d4006a32a83ae6de9d82790d0cc597a0ee9b5"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}