Coverage for pySDC/tutorial/step_9/D_adaptive_alpha.py: 100%
77 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +0000
1"""
2This script shows how to let ParaDiag choose its alpha by itself.
4ParaDiag replaces the time-stepping matrix by an alpha-circulant approximation. Alpha trades two error
5sources against each other: a small value approximates the original problem better and converges in
6fewer iterations, but conditions the diagonalization worse, so round-off and inexact inner solves get
7amplified. A single fixed alpha therefore has to be a compromise for the whole run, even though the
8right balance shifts as the residual falls.
10The `AdaptiveAlpha` convergence controller updates alpha after every iteration instead, following
11`Caklovic et al. <https://doi.org/10.2140/camcos.2023.18.55>`_:
13 gamma = L * (3 * eps + tau)
14 alpha_k = sqrt(gamma * r_k / e_k)
15 e_{k+1} = 2 * sqrt(gamma * e_k * r_k)
17with L the block size, eps machine precision, tau the inner solver tolerance, r_k the residual and
18e_k a running bound on the error. Gamma is an accuracy floor: there is no point pushing alpha below
19the level at which round-off and the inner solver dominate anyway.
21We compare a few fixed alphas against the adaptive one on the advection problem from Part C.
22The interesting result is not that adaptive wins on iteration count -- it ties with the best fixed
23value -- but that it gets there without being told, and while keeping alpha orders of magnitude
24larger, which is exactly the margin that protects you once the inner solves are inexact.
26Everything here runs with the "virtually parallel" controller, which keeps all steps in one process.
27Alpha is a property of the method, not of the parallelization, so this is the right place to pin it
28down; Part E then runs the same thing across MPI ranks and checks it comes out the same.
29"""
31from pathlib import Path
33# we always do this many time-steps in total, no matter how many of them run in parallel
34num_steps_total = 4
37def get_description():
38 """
39 Set up the same advection problem as in Part C.
41 Returns:
42 dict: the description for the ParaDiag controller
43 """
44 from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd
45 from pySDC.implementations.sweeper_classes.ParaDiagSweepers import QDiagonalization
47 level_params = {}
48 level_params['dt'] = 0.1
49 level_params['restol'] = 1e-6
51 sweeper_params = {}
52 sweeper_params['quad_type'] = 'RADAU-RIGHT'
53 sweeper_params['num_nodes'] = 3
54 sweeper_params['initial_guess'] = 'copy'
56 # Part C uses GMRES here to count linear solver work. We only care about the parallelism, and the
57 # complex shifted systems ParaDiag produces are hard for GMRES, so we solve them directly instead.
58 problem_params = {'nvars': 64, 'order': 8, 'c': 1, 'solver_type': 'direct'}
60 step_params = {}
61 step_params['maxiter'] = 99
63 description = {}
64 description['problem_class'] = advectionNd
65 description['problem_params'] = problem_params
66 description['sweeper_class'] = QDiagonalization
67 description['sweeper_params'] = sweeper_params
68 description['level_params'] = level_params
69 description['step_params'] = step_params
71 return description
74# the fixed values we compare against, plus the adaptive strategy
75alpha_settings = [1e-2, 1e-4, 1e-8, 'adaptive']
78def get_controller_params(alpha):
79 """
80 Controller parameters for one alpha setting.
82 Args:
83 alpha: a number, or the string 'adaptive'
85 Returns:
86 tuple: the controller parameters and the extra description entries
87 """
88 from pySDC.implementations.convergence_controller_classes.adaptive_alpha import AdaptiveAlpha
90 controller_params = {}
91 controller_params['logger_level'] = 30
92 controller_params['average_jacobian'] = False
94 extra_description = {}
95 if alpha == 'adaptive':
96 # the adaptive controller overwrites this from the first iteration onwards, but ParaDiag needs
97 # some alpha to build its first transform with
98 controller_params['alpha'] = 1e-4
99 extra_description['convergence_controllers'] = {AdaptiveAlpha: {}}
100 else:
101 controller_params['alpha'] = alpha
103 return controller_params, extra_description
106def format_result(mode, alpha, niter, error, final_alpha):
107 """
108 One line of output, in the same shape for both controllers so they can be compared.
110 Args:
111 mode (str): which controller produced it, e.g. 'virtual' or 'MPI on 4'
112 alpha: the alpha setting used
113 niter (int): number of iterations needed
114 error (float): error against the exact solution
115 final_alpha (float): the alpha in use when the run finished
117 Returns:
118 str: the formatted line
119 """
120 return (
121 f'{mode:>11s}: alpha {str(alpha):>9s} -> {niter:2d} iterations, '
122 f'error {error:.4e}, final alpha {final_alpha:.3e}'
123 )
126def run(alpha, block_size, comm=None):
127 """
128 Run the advection problem with one alpha setting.
130 Args:
131 alpha: a number, or the string 'adaptive'
132 block_size (int): number of time-steps in one block
133 comm: MPI communicator, or None for the virtually parallel controller
135 Returns:
136 tuple: the end value, the iteration count, the error and the final alpha
137 """
138 import numpy as np
139 from pySDC.helpers.stats_helper import get_sorted
141 controller_params, extra_description = get_controller_params(alpha)
142 description = {**get_description(), **extra_description}
144 if comm is None:
145 from pySDC.implementations.controller_classes.controller_ParaDiag_nonMPI import controller_ParaDiag_nonMPI
147 controller_params['mssdc_jac'] = False
148 controller = controller_ParaDiag_nonMPI(
149 controller_params=controller_params, description=description, num_procs=block_size
150 )
151 steps = controller.MS
152 else:
153 from pySDC.implementations.controller_classes.controller_ParaDiag_MPI import controller_ParaDiag_MPI
155 controller = controller_ParaDiag_MPI(controller_params=controller_params, description=description, comm=comm)
156 steps = [controller.S]
158 # ParaDiag diagonalizes in time, so the solution becomes complex
159 for S in steps:
160 S.levels[0].prob.init = tuple([*S.levels[0].prob.init[:2]] + [np.dtype('complex128')])
162 P = steps[0].levels[0].prob
163 dt = steps[0].levels[0].params.dt
164 Tend = num_steps_total * dt
166 uend, stats = controller.run(u0=P.u_exact(0.0), t0=0.0, Tend=Tend)
167 niter = max(int(me[1]) for me in get_sorted(stats, type='niter', sortby='time', comm=comm))
169 return uend, niter, abs(uend - P.u_exact(Tend)), controller.params.alpha
172def main(fname='step_9_D_out.txt'):
173 """
174 Compare fixed and adaptive alpha with the virtually parallel controller.
176 Args:
177 fname (str): file under ``data/`` to write the results to
178 """
180 import numpy as np
182 # one block holding every step; Part E runs the same settings across MPI ranks
183 block_size = num_steps_total
185 results = {}
186 lines = []
187 for alpha in alpha_settings:
188 uend, niter, error, final_alpha = run(alpha, block_size)
189 results[alpha] = (uend, niter)
190 lines.append(format_result('virtual', alpha, niter, error, final_alpha))
192 Path("data").mkdir(parents=True, exist_ok=True)
193 with open('data/' + fname, 'w') as f:
194 for line in lines:
195 f.write(line + '\n')
196 print(line)
198 # the adaptive strategy should need no more iterations than the best fixed alpha we tried
199 best_fixed = min(results[a][1] for a in alpha_settings if a != 'adaptive')
200 assert (
201 results['adaptive'][1] <= best_fixed
202 ), 'ERROR: adaptive alpha needed %s iterations, the best fixed alpha only %s' % (results['adaptive'][1], best_fixed)
204 # alpha changes the iteration, not the problem, so all settings solve the same thing
205 reference = results[alpha_settings[0]][0]
206 for alpha in alpha_settings[1:]:
207 assert np.allclose(results[alpha][0], reference, atol=1e-5), (
208 'ERROR: alpha %s gives a different solution' % alpha
209 )
212if __name__ == "__main__":
213 main()