Skip to content

Commit

Permalink
bad indentation, superflous parens and various different errors fixed
Browse files Browse the repository at this point in the history
  • Loading branch information
Castlefox committed Jan 6, 2025
1 parent 57a029b commit 7fa6b92
Showing 1 changed file with 85 additions and 63 deletions.
148 changes: 85 additions & 63 deletions examples/PyMPDATA_examples/Jaruga_et_al_2015/fig19.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
"$$"
]
},
{
Expand Down Expand Up @@ -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"
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand Down Expand Up @@ -292,17 +298,17 @@
" 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",
" _U: Solver(advectee=new_sf(N,M), **solver_ctor_kwargs),\n",
" _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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -374,7 +380,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "8b84b7919b314b9383b748fd32ec0d70",
"model_id": "a3a66e7d490249288d4b38f9cc3efa5c",
"version_major": 2,
"version_minor": 0
},
Expand All @@ -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"
]
}
],
Expand All @@ -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",
Expand Down Expand Up @@ -470,7 +476,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "0a3762134cfa452cacece883690eabc7",
"model_id": "baa8f55f01f5441bb96110e9b0d523b6",
"version_major": 2,
"version_minor": 0
},
Expand All @@ -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=\"<a href='.\\\\tmpgsxuio24.gif' target='_blank'>.\\\\tmpgsxuio24.gif</a><br>\")"
"HTML(value=\"<a href='.\\\\tmp531j16v4.gif' target='_blank'>.\\\\tmp531j16v4.gif</a><br>\")"
]
},
"metadata": {},
Expand Down Expand Up @@ -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": {
Expand Down

0 comments on commit 7fa6b92

Please sign in to comment.