diff --git a/examples/PyMPDATA_examples/Jaruga_et_al_2015/fig19.ipynb b/examples/PyMPDATA_examples/Jaruga_et_al_2015/fig19.ipynb index 220d0798..3c3ce2c4 100644 --- a/examples/PyMPDATA_examples/Jaruga_et_al_2015/fig19.ipynb +++ b/examples/PyMPDATA_examples/Jaruga_et_al_2015/fig19.ipynb @@ -10,7 +10,14 @@ "### equations solved:\n", "\n", "$$\n", - "\\partial_t \\theta + ...$$" + "\\partial_t \\theta + \n", + "$$\n", + "$$\n", + "\\partial_t \\theta + \n", + "$$\n", + "$$\n", + "\\partial_t \\theta + \n", + "$$" ] }, { @@ -39,13 +46,13 @@ "from PyMPDATA import Options,ScalarField,VectorField,Solver,Stepper\n", "from matplotlib import colors, pyplot\n", "from PyMPDATA.boundary_conditions import Periodic\n", - "from PyMPDATA_examples.Jaruga_et_al_2015.temp import *\n", "import numpy as np\n", "from numba import jit\n", "import warnings\n", "import ipywidgets\n", "from numba.core.errors import NumbaExperimentalFeatureWarning\n", - "from PyMPDATA.impl.enumerations import IMPL_META_AND_DATA, META_AND_DATA_DATA, IMPL_BC" + "from PyMPDATA.impl.enumerations import IMPL_META_AND_DATA, META_AND_DATA_DATA, IMPL_BC\n", + "from IPython.display import display" ] }, { @@ -115,7 +122,7 @@ "def lap(Phi, dxy, err_init, lap_tmp, tmp_uvw = None):\n", " xchng_pres(Phi.impl)\n", " calc_grad([x.impl for x in lap_tmp], Phi, dxy)\n", - " if (err_init):\n", + " if err_init:\n", " for k in (0, 1):\n", " lap_tmp[k].data[:] -= data(tmp_uvw[k])\n", "\n", @@ -128,54 +135,54 @@ "def pressure_solver_loop_init(err,p_err,lap_p_err,dxy,lap_tmp,tmp_uvw):\n", " p_err[0].get()[:] = err.get()[:]\n", " lap_p_err[0][:] = lap(p_err[0], dxy, False, lap_tmp,tmp_uvw)\n", - " \n", - "##@njit\n", + "\n", + "#@njit\n", "def pressure_solver_loop_body(Phi,beta,converged,err,p_err,lap_p_err,dxy,\n", " k_iters,err_tol,lap_err,lap_tmp):\n", " tmp_den = [1.]*k_iters\n", " alpha = [1.]*k_iters\n", " for v in range(0,k_iters):\n", - " tmp_den[v] = np.sum(lap_p_err[v]**2)\n", - " if (tmp_den[v] != 0):\n", - " beta = - np.dot(\n", - " err.get().ravel(),\n", - " lap_p_err[v].ravel()\n", - " ) / tmp_den[v]\n", - " Phi.get()[:] += beta * p_err[v].get()\n", - " err.get()[:] += beta * lap_p_err[v]\n", - "\n", - " error = max(\n", - " abs(np.amax(err.get())),\n", - " abs(np.amin(err.get()))\n", - " )\n", - "\n", - " if (error <= err_tol): \n", - " converged = True\n", - "\n", - " lap_err[:] = lap(err, dxy, False, lap_tmp)\n", - " \n", - " for l in range(v):\n", - " if (tmp_den[l] != 0):\n", - " alpha[l] = - np.dot(lap_err.ravel(), lap_p_err[l].ravel()) / tmp_den[l]\n", - " if (v < (k_iters - 1)):\n", - " p_err[v + 1].get()[:] = err.get()\n", - " lap_p_err[v + 1][:] = lap_err[:]\n", - " for l in range(v):\n", - " p_err[v + 1].get()[:] += alpha[l] * p_err[l].get()\n", - " lap_p_err[v + 1][:] += alpha[l] * lap_p_err[l]\n", - " else:\n", - " p_err[0].get()[:] = err.get()[:] + alpha[0] * p_err[0].get()\n", - " lap_p_err[0][:] = lap_err[:] + alpha[0] * lap_p_err[0]\n", - " for l in range(1,v+1):\n", - " p_err[0].get()[:] += alpha[l] * p_err[l].get()\n", - " lap_p_err[0][:] += alpha[l] * lap_p_err[l]\n", + " tmp_den[v] = np.sum(lap_p_err[v]**2)\n", + " if tmp_den[v] != 0:\n", + " beta = - np.dot(\n", + " err.get().ravel(),\n", + " lap_p_err[v].ravel()\n", + " ) / tmp_den[v]\n", + " Phi.get()[:] += beta * p_err[v].get()\n", + " err.get()[:] += beta * lap_p_err[v]\n", + "\n", + " error = max(\n", + " abs(np.amax(err.get())),\n", + " abs(np.amin(err.get()))\n", + " )\n", + "\n", + " if error <= err_tol: \n", + " converged = True\n", + "\n", + " lap_err[:] = lap(err, dxy, False, lap_tmp)\n", + "\n", + " for l in range(v):\n", + " if tmp_den[l] != 0:\n", + " alpha[l] = - np.dot(lap_err.ravel(), lap_p_err[l].ravel()) / tmp_den[l]\n", + " if v < k_iters - 1:\n", + " p_err[v + 1].get()[:] = err.get()\n", + " lap_p_err[v + 1][:] = lap_err[:]\n", + " for l in range(v):\n", + " p_err[v + 1].get()[:] += alpha[l] * p_err[l].get()\n", + " lap_p_err[v + 1][:] += alpha[l] * lap_p_err[l]\n", + " else:\n", + " p_err[0].get()[:] = err.get()[:] + alpha[0] * p_err[0].get()\n", + " lap_p_err[0][:] = lap_err[:] + alpha[0] * lap_p_err[0]\n", + " for l in range(1,v+1):\n", + " p_err[0].get()[:] += alpha[l] * p_err[l].get()\n", + " lap_p_err[0][:] += alpha[l] * lap_p_err[l]\n", " return converged\n", "\n", "#@njit\n", - "def pressure_solver_update(solvers,Phi,beta,lap_tmp,tmp_uvw,err,p_err,lap_p_err,dxy,\n", - " k_iters,err_tol,lap_err,simple = False):\n", + "def pressure_solver_update(solvers,Phi,beta,lap_tmp,tmp_uvw,err,p_err,\n", + " lap_p_err,dxy,k_iters,err_tol,lap_err):\n", " for k in (0, 1):\n", - " data(tmp_uvw[k])[:] = solvers[k].advectee.data\n", + " data(tmp_uvw[k])[:] = solvers[k].advectee.data\n", " \n", " #initial error\n", " err.get()[:] = lap(Phi, dxy, True,lap_tmp,tmp_uvw)\n", @@ -186,16 +193,15 @@ " pressure_solver_loop_init(err,p_err,lap_p_err,dxy,lap_tmp,tmp_uvw)\n", " #pseudo-time loop\n", " while not converged: \n", - " converged = pressure_solver_loop_body(Phi,beta,converged,err,p_err,lap_p_err,\n", + " converged = pressure_solver_loop_body(Phi,beta,converged,err,p_err,lap_p_err,\n", " dxy,k_iters,err_tol,lap_err,lap_tmp)\n", - " iters += 1\n", - "\n", - " if (iters > 10000): # going beyond 10000 iters means something is really wrong,\n", - " # usually boundary conditions but not always !\n", - " raise Exception(\"stuck in pressure solver\")\n", + " iters += 1\n", "\n", + " if (iters > 10000): # going beyond 10000 iters means something is really wrong,\n", + " # usually boundary conditions but not always !\n", + " raise ArithmeticError(\"stuck in pressure solver\")\n", + " \n", " xchng_pres(Phi.impl)\n", - "\n", " calc_grad(tmp_uvw, Phi, dxy)\n", "\n", "#@njit\n", @@ -292,7 +298,7 @@ " return field\n", "\n", "def new_sfs(N, M, impl=False):\n", - " return tuple([new_sf(N,M).impl if impl else new_sf(N,M) for _ in ('u', 'w')])\n", + " return tuple(new_sf(N,M).impl if impl else new_sf(N,M) for _ in ('u', 'w'))\n", "\n", "solver_ctor_kwargs = {'stepper': stepper, 'advector': advector}\n", "solvers = {\n", @@ -300,9 +306,9 @@ " _W: Solver(advectee=new_sf(N,M), **solver_ctor_kwargs),\n", " _T: Solver(advectee=ScalarField(data=mesh, **field_ctor_kwargs), **solver_ctor_kwargs),\n", "}\n", - "advectee_data = tuple([\n", + "advectee_data = tuple(\n", " solvers[idx].advectee.data for idx in range(3)\n", - "])\n", + ")\n", "\n", "rhs_w = np.zeros((N,M))\n", "stash = new_sfs(N, M)\n", @@ -339,9 +345,9 @@ "\n", "@jit(**options.jit_flags)\n", "def vip_rhs_apply(dt, vip_rhs, advectee_data):\n", - " for k in (0, 1):\n", - " advectee_data[k][:] += 0.5 * dt * data(vip_rhs[k])\n", - " data(vip_rhs[k])[:] = 0\n", + " for k in (0, 1):\n", + " advectee_data[k][:] += 0.5 * dt * data(vip_rhs[k])\n", + " data(vip_rhs[k])[:] = 0\n", "\n", "@jit(**options.jit_flags)\n", "def apply_rhs(w, rhs_w : np.ndarray, dt : float):\n", @@ -374,7 +380,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "8b84b7919b314b9383b748fd32ec0d70", + "model_id": "a3a66e7d490249288d4b38f9cc3efa5c", "version_major": 2, "version_minor": 0 }, @@ -389,8 +395,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "CPU times: total: 1min 36s\n", - "Wall time: 1min 46s\n" + "CPU times: total: 1min 31s\n", + "Wall time: 1min 41s\n" ] } ], @@ -406,7 +412,7 @@ "\n", " Phi.get()[:] = 0\n", " pressure_solver_update(solvers,Phi,beta,lap_tmp,tmp_uvw,err,p_err,lap_p_err,dxy,\n", - " k_iters,err_tol,lap_err,simple = True)\n", + " k_iters,err_tol,lap_err)\n", " xchng_pres(Phi.impl)\n", " calc_grad(tmp_uvw, Phi, dxy)\n", " pressure_solver_apply(advectee_data, tmp_uvw)\n", @@ -470,7 +476,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "0a3762134cfa452cacece883690eabc7", + "model_id": "baa8f55f01f5441bb96110e9b0d523b6", "version_major": 2, "version_minor": 0 }, @@ -496,12 +502,12 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "56da990921454749b499d3d9f2e94776", + "model_id": "e1864a55f28a46f5850e0e2b10732fb9", "version_major": 2, "version_minor": 0 }, "text/plain": [ - "HTML(value=\".\\\\tmpgsxuio24.gif
\")" + "HTML(value=\".\\\\tmp531j16v4.gif
\")" ] }, "metadata": {}, @@ -533,6 +539,22 @@ "display(progress)\n", "show_anim(plot, frame_range=range(len(output)))" ] + }, + { + "cell_type": "markdown", + "id": "58f716ed-b17f-49de-9bdf-c655942e63bd", + "metadata": {}, + "source": [ + "add showplot\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e19fba6c-c5b1-43ba-8f99-a23222e0afd1", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": {