|
49 | 49 | "N = 128\n", |
50 | 50 | "\n", |
51 | 51 | "# setup grid\n", |
52 | | - "dx, grid = get_1d_grid(size=N, bc='periodic', left_boundary=0, right_boundary=2*np.pi)\n", |
| 52 | + "dx, grid = get_1d_grid(size=N, bc='periodic', left_boundary=0, right_boundary=2 * np.pi)\n", |
53 | 53 | "\n", |
54 | 54 | "# setup finite difference matrices\n", |
55 | 55 | "fd_params = {'order': 4, 'stencil_type': 'center', 'dx': dx, 'size': N, 'dim': 1, 'bc': 'periodic'}\n", |
56 | 56 | "advection, _ = get_finite_difference_matrix(derivative=1, **fd_params)\n", |
57 | 57 | "diffusion, _ = get_finite_difference_matrix(derivative=2, **fd_params)\n", |
58 | 58 | "\n", |
| 59 | + "\n", |
59 | 60 | "def eval_f(u):\n", |
60 | | - " return (c*advection + nu*diffusion) @ u" |
| 61 | + " return (c * advection + nu * diffusion) @ u" |
61 | 62 | ] |
62 | 63 | }, |
63 | 64 | { |
|
76 | 77 | "outputs": [], |
77 | 78 | "source": [ |
78 | 79 | "u = np.sin(grid)\n", |
79 | | - "assert np.allclose(diffusion@u, -np.sin(grid)) # check diffusion part\n", |
80 | | - "assert np.allclose(advection@u, np.cos(grid)) # check advection part\n", |
81 | | - "assert np.allclose(eval_f(u), -nu*np.sin(grid)+c*np.cos(grid)) # check f evaluation" |
| 80 | + "assert np.allclose(diffusion @ u, -np.sin(grid)) # check diffusion part\n", |
| 81 | + "assert np.allclose(advection @ u, np.cos(grid)) # check advection part\n", |
| 82 | + "assert np.allclose(eval_f(u), -nu * np.sin(grid) + c * np.cos(grid)) # check f evaluation" |
82 | 83 | ] |
83 | 84 | }, |
84 | 85 | { |
|
101 | 102 | "from scipy.sparse.linalg import spsolve\n", |
102 | 103 | "from scipy.sparse import eye\n", |
103 | 104 | "\n", |
| 105 | + "\n", |
104 | 106 | "def implicit_euler(rhs, dt):\n", |
105 | | - " A = eye(N) - dt*(c*advection + nu*diffusion)\n", |
| 107 | + " A = eye(N) - dt * (c * advection + nu * diffusion)\n", |
106 | 108 | " return spsolve(A, rhs)" |
107 | 109 | ] |
108 | 110 | }, |
|
127 | 129 | "# check advection\n", |
128 | 130 | "nu = 0\n", |
129 | 131 | "c = 1e-1\n", |
130 | | - "assert np.allclose(implicit_euler(u, dt), np.sin(grid+c*dt))\n", |
| 132 | + "assert np.allclose(implicit_euler(u, dt), np.sin(grid + c * dt))\n", |
131 | 133 | "\n", |
132 | 134 | "# check diffusion\n", |
133 | 135 | "nu = 1e-1\n", |
134 | 136 | "c = 0\n", |
135 | | - "assert np.allclose(implicit_euler(u, dt), np.sin(grid)*np.exp(-nu*dt))" |
| 137 | + "assert np.allclose(implicit_euler(u, dt), np.sin(grid) * np.exp(-nu * dt))" |
136 | 138 | ] |
137 | 139 | }, |
138 | 140 | { |
|
161 | 163 | "import scipy.sparse as sp\n", |
162 | 164 | "import numpy as np\n", |
163 | 165 | "\n", |
| 166 | + "\n", |
164 | 167 | "class AdvectionDiffusion(Problem):\n", |
165 | 168 | " dtype_u = mesh # wraps numpy ndarray\n", |
166 | 169 | " dtype_f = mesh\n", |
167 | | - " \n", |
| 170 | + "\n", |
168 | 171 | " def __init__(self, N, c, nu):\n", |
169 | 172 | " init = (N, None, np.dtype('float64')) # describes how to initialize data\n", |
170 | 173 | " super().__init__(init=init)\n", |
171 | | - " \n", |
| 174 | + "\n", |
172 | 175 | " # setup grid\n", |
173 | | - " dx, self.grid = get_1d_grid(size=N, bc='periodic', left_boundary=0, right_boundary=2*np.pi)\n", |
| 176 | + " dx, self.grid = get_1d_grid(size=N, bc='periodic', left_boundary=0, right_boundary=2 * np.pi)\n", |
174 | 177 | "\n", |
175 | 178 | " # setup finite difference matrices\n", |
176 | 179 | " fd_params = {'order': 4, 'stencil_type': 'center', 'dx': dx, 'size': N, 'dim': 1, 'bc': 'periodic'}\n", |
|
201 | 204 | " # prepare a new work counter for this factor / step size. Seems weird now, but bear with me.\n", |
202 | 205 | " key = f'gmres-{dt:.2e}'\n", |
203 | 206 | " self.work_counters[f'gmres-{dt:.2e}'] = self.work_counters.get(key, WorkCounter())\n", |
204 | | - " \n", |
205 | | - " me[...], _ = sp.linalg.gmres(self.Id - dt*self.A, rhs, x0=u0, callback=self.work_counters[key], rtol=1e-10, callback_type='pr_norm')\n", |
| 207 | + "\n", |
| 208 | + " me[...], _ = sp.linalg.gmres(\n", |
| 209 | + " self.Id - dt * self.A, rhs, x0=u0, callback=self.work_counters[key], rtol=1e-10, callback_type='pr_norm'\n", |
| 210 | + " )\n", |
206 | 211 | " return me\n", |
207 | 212 | "\n", |
208 | 213 | " def u0(self):\n", |
|
229 | 234 | "source": [ |
230 | 235 | "prob = AdvectionDiffusion(128, 1e-1, 1e-3)\n", |
231 | 236 | "u = prob.u0()\n", |
232 | | - "assert np.allclose(prob.eval_f(u, 0), -prob.nu*np.sin(prob.grid)+prob.c*np.cos(prob.grid))\n", |
| 237 | + "assert np.allclose(prob.eval_f(u, 0), -prob.nu * np.sin(prob.grid) + prob.c * np.cos(prob.grid))\n", |
233 | 238 | "\n", |
234 | 239 | "dt = 1e-2\n", |
235 | 240 | "advection_prob = AdvectionDiffusion(N=128, c=1e-1, nu=0)\n", |
236 | 241 | "u = advection_prob.u0()\n", |
237 | | - "assert abs(advection_prob.solve_system(u, dt, u, 0)- np.sin(advection_prob.grid+advection_prob.c*dt)) < 1e-6\n", |
| 242 | + "assert abs(advection_prob.solve_system(u, dt, u, 0) - np.sin(advection_prob.grid + advection_prob.c * dt)) < 1e-6\n", |
238 | 243 | "\n", |
239 | 244 | "diffusion_prob = AdvectionDiffusion(N=128, c=0, nu=1e-1)\n", |
240 | 245 | "u = diffusion_prob.u0()\n", |
241 | | - "assert abs(diffusion_prob.solve_system(u, dt, u, 0)- np.sin(diffusion_prob.grid)*np.exp(-diffusion_prob.nu*dt)) < 1e-6" |
| 246 | + "assert (\n", |
| 247 | + " abs(diffusion_prob.solve_system(u, dt, u, 0) - np.sin(diffusion_prob.grid) * np.exp(-diffusion_prob.nu * dt)) < 1e-6\n", |
| 248 | + ")" |
242 | 249 | ] |
243 | 250 | }, |
244 | 251 | { |
|
283 | 290 | "Controller: <class 'pySDC.implementations.controller_classes.controller_nonMPI.controller_nonMPI'>\n", |
284 | 291 | " all_to_done = False\n", |
285 | 292 | " dump_setup = True\n", |
286 | | - " fname = run_pid77606.log\n", |
| 293 | + " fname = run_pid94051.log\n", |
287 | 294 | "--> hook_class = [<class 'pySDC.implementations.hooks.default_hook.DefaultHooks'>, <class 'pySDC.implementations.hooks.log_timings.CPUTimings'>, <class 'pySDC.implementations.hooks.log_solution.LogSolution'>, <class 'pySDC.implementations.hooks.log_work.LogWork'>]\n", |
288 | 295 | " log_to_file = False\n", |
289 | 296 | "--> logger_level = 15\n", |
|
349 | 356 | "sweeper_params['num_nodes'] = 3\n", |
350 | 357 | "sweeper_params['QI'] = 'MIN-SR-S'\n", |
351 | 358 | "\n", |
352 | | - "problem_params = {'N': 128, 'c': 8e-1, 'nu':1e-1,}\n", |
| 359 | + "problem_params = {\n", |
| 360 | + " 'N': 128,\n", |
| 361 | + " 'c': 8e-1,\n", |
| 362 | + " 'nu': 1e-1,\n", |
| 363 | + "}\n", |
353 | 364 | "\n", |
354 | 365 | "# gather all parameters in one dictionary and add problem and sweeper class\n", |
355 | 366 | "description = {}\n", |
|
430 | 441 | "hooks - INFO: Process 0 on time 0.900000 at stage IT_FINE: Level: 0 -- Iteration: 2 -- Sweep: 1 -- residual: 5.78007793e-06\n", |
431 | 442 | "hooks - INFO: Process 0 on time 0.900000 at stage IT_FINE: Level: 0 -- Iteration: 3 -- Sweep: 1 -- residual: 1.36492229e-07\n", |
432 | 443 | "hooks - INFO: Process 0 on time 0.900000 at stage IT_FINE: Level: 0 -- Iteration: 4 -- Sweep: 1 -- residual: 3.15861094e-09\n", |
433 | | - "hooks - INFO: Finished run after 8.13e-02s\n" |
| 444 | + "hooks - INFO: Finished run after 8.51e-02s\n" |
434 | 445 | ] |
435 | 446 | } |
436 | 447 | ], |
|
563 | 574 | "source": [ |
564 | 575 | "# get values on the diagonal of the preconditioner\n", |
565 | 576 | "sweeper = controller.MS[0].levels[0].sweep\n", |
566 | | - "diag_vals = [sweeper.QI[i+1, i+1] for i in range(sweeper.coll.num_nodes)]\n", |
| 577 | + "diag_vals = [sweeper.QI[i + 1, i + 1] for i in range(sweeper.coll.num_nodes)]\n", |
567 | 578 | "\n", |
568 | 579 | "# get the keys for gmres iterations at various step sizes from this\n", |
569 | 580 | "work_keys = [f'work_gmres-{sweeper.level.dt * me:.2e}' for me in diag_vals]\n", |
|
573 | 584 | "total_gmres_iterations = sum(gmres_iter_per_task)\n", |
574 | 585 | "\n", |
575 | 586 | "# compute speedup in implicit solves\n", |
576 | | - "max_rel_gmres_per_task = max(gmres_iter_per_task)/total_gmres_iterations\n", |
577 | | - "speedup = 1/max_rel_gmres_per_task\n", |
| 587 | + "max_rel_gmres_per_task = max(gmres_iter_per_task) / total_gmres_iterations\n", |
| 588 | + "speedup = 1 / max_rel_gmres_per_task\n", |
578 | 589 | "parallel_efficiency = speedup / sweeper.coll.num_nodes\n", |
579 | 590 | "\n", |
580 | 591 | "# visualize and print the results\n", |
581 | 592 | "import matplotlib.pyplot as plt\n", |
| 593 | + "\n", |
582 | 594 | "plt.bar(np.arange(sweeper.coll.num_nodes), gmres_iter_per_task)\n", |
583 | 595 | "plt.ylabel('GMRES iterations')\n", |
584 | 596 | "plt.xlabel('Collocation node / task')\n", |
585 | 597 | "\n", |
586 | | - "print(f'Performed {max_rel_gmres_per_task*100:.0f}% of GMRES iterations on task {np.argmax(gmres_iter_per_task)}. This translates to speedup of {speedup:.2f} and {parallel_efficiency*100:.0f}% parallel efficiency in the implicit solves.')\n" |
| 598 | + "print(\n", |
| 599 | + " f'Performed {max_rel_gmres_per_task*100:.0f}% of GMRES iterations on task {np.argmax(gmres_iter_per_task)}. This translates to speedup of {speedup:.2f} and {parallel_efficiency*100:.0f}% parallel efficiency in the implicit solves.'\n", |
| 600 | + ")" |
587 | 601 | ] |
588 | 602 | }, |
589 | 603 | { |
|
0 commit comments