Step-2: Data structures and my first sweeper¶
In this step, we will dig a bit deeper into the data structures of
pySDC and work with our first SDC sweeper.
Part A: Step data structure¶
We start with creating the basic data structure pySDC is built from:
the step. It represents a single time step with whatever hierarchy
we prescribe (more on that later). The relevant data (e.g. the solution
and right-hand side vectors as well as the problem instances and so on)
are all part of a level and a set of levels make a step
(together with transfer operators, see the following tutorials). In this
first example, we simply create a step and test the problem instance of
its only level. This is the same test we ran
here.
Important things to note:
This part is for demonstration purpose only. Normally, users do not have to deal with the internal data structures, see Part C below.
The
descriptiondictionary is the central steering tool forpySDC.Happily passing around classes instead of instances make life much easier, although it is not the best way in terms of programming…
Full code: pySDC/tutorial/step_2/A_step_data_structure.py
from pathlib import Path
from pySDC.core.Step import step
from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
from pySDC.tutorial.step_1.A_spatial_problem_setup import run_accuracy_check
def main():
"""
A simple test program to setup a full step instance
"""
# initialize level parameters
level_params = dict()
level_params['restol'] = 1e-10
level_params['dt'] = 0.1
# initialize sweeper parameters
sweeper_params = dict()
sweeper_params['quad_type'] = 'RADAU-RIGHT'
sweeper_params['num_nodes'] = 3
sweeper_params['QI'] = 'LU'
# initialize problem parameters
problem_params = dict()
problem_params['nu'] = 0.1 # diffusion coefficient
problem_params['freq'] = 4 # frequency for the test value
problem_params['nvars'] = 1023 # number of degrees of freedom
problem_params['bc'] = 'dirichlet-zero' # boundary conditions
# initialize step parameters
step_params = dict()
step_params['maxiter'] = 20
# fill description dictionary for easy step instantiation
description = dict()
description['problem_class'] = heatNd_unforced
description['problem_params'] = problem_params
description['sweeper_class'] = generic_implicit
description['sweeper_params'] = sweeper_params
description['level_params'] = level_params
description['step_params'] = step_params
# now the description contains more or less everything we need to create a step
S = step(description=description)
# we only have a single level, make a shortcut
L = S.levels[0]
# one of the integral parts of each level is the problem class, make a shortcut
P = L.prob
# now we can do e.g. what we did before with the problem
err = run_accuracy_check(P)
Path("data").mkdir(parents=True, exist_ok=True)
f = open('data/step_2_A_out.txt', 'w')
out = 'Error of the spatial accuracy test: %8.6e' % err
f.write(out)
print(out)
f.close()
assert err <= 2e-04, "ERROR: the spatial accuracy is higher than expected, got %s" % err
if __name__ == "__main__":
main()
Results:
Error of the spatial accuracy test: 1.981783e-04
Part B: My first sweeper¶
Since we know how to create a step, we now create our first SDC
iteration (by hand, this time). We make use of the IMEX SDC sweeper
imex_1st_order and of the problem class
HeatEquation_1D_FD_forced. Also, the data structure for the
right-hand side is now rhs_imex_mesh, since we need implicit and
explicit parts of the right-hand side. The rest is rather
straightforward: we set initial values and times, start by spreading the
data and then do the iteration until the maximum number of iterations
is reached or until the residual is small enough. Yet, this example
sheds light on the key functionalities of the sweeper:
compute_residual, update_nodes and compute_end_point. Also,
the step and level structures are explored a bit more deeply,
since we make use of parameters and status objects here.
Important things to note:
Again, this part is for demonstration purpose only. Normally, users do not have to deal with the internal data structures, see Part C below.
Note the difference between status and parameter objects: parameters are user-defined flags created using the dictionaries (e.g.
maxiteras part of thestep_paramsin this example), while status objects are internal control objects which reflect the current status of a level or a step (e.g.iterortime).The logic required to implement an SDC iteration is simple but also rather tedious. This will get worse if you want to deal with more than one step or, behold, multiple parallel steps each with a space-time hierarchy. Therefore, please proceed quickly to part C!
Full code: pySDC/tutorial/step_2/B_my_first_sweeper.py
from pathlib import Path
from pySDC.core.Step import step
from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced
from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
def main():
"""
A simple test program to run IMEX SDC for a single time step
"""
# initialize level parameters
level_params = dict()
level_params['restol'] = 1e-10
level_params['dt'] = 0.1
# initialize sweeper parameters
sweeper_params = dict()
sweeper_params['quad_type'] = 'RADAU-RIGHT'
sweeper_params['num_nodes'] = 3
# initialize problem parameters
problem_params = dict()
problem_params['nu'] = 0.1 # diffusion coefficient
problem_params['freq'] = 4 # frequency for the test value
problem_params['nvars'] = 1023 # number of degrees of freedom
problem_params['bc'] = 'dirichlet-zero' # boundary conditions
# initialize step parameters
step_params = dict()
step_params['maxiter'] = 20
# Fill description dictionary for easy hierarchy creation
description = dict()
description['problem_class'] = heatNd_forced
description['problem_params'] = problem_params
description['sweeper_class'] = imex_1st_order
description['sweeper_params'] = sweeper_params
description['level_params'] = level_params
description['step_params'] = step_params
# instantiate the step we are going to work on
S = step(description=description)
# run IMEX SDC test and check error, residual and number of iterations
err, res, niter = run_imex_sdc(S)
print('Error and residual: %12.8e -- %12.8e' % (err, res))
assert err <= 1e-5, "ERROR: IMEX SDC iteration did not reduce the error enough, got %s" % err
assert res <= level_params['restol'], "ERROR: IMEX SDC iteration did not reduce the residual enough, got %s" % res
assert niter <= 12, "ERROR: IMEX SDC took too many iterations, got %s" % niter
def run_imex_sdc(S):
"""
Routine to run IMEX SDC on a single time step
Args:
S: an instance of a step representing the time step
Returns:
the error of SDC vs. exact solution
the residual after the SDC sweeps
the number of iterations
"""
# make shortcuts for the level and the problem
L = S.levels[0]
P = L.prob
# set initial time in the status of the level
L.status.time = 0.1
# compute initial value (using the exact function here)
L.u[0] = P.u_exact(L.time)
# access the sweeper's predict routine to get things started
# if we don't do this, the values at the nodes are not initialized
L.sweep.predict()
# compute the residual (we may be done already!)
L.sweep.compute_residual()
# reset iteration counter
S.status.iter = 0
# run the SDC iteration until either the maximum number of iterations is reached or the residual is small enough
Path("data").mkdir(parents=True, exist_ok=True)
f = open('data/step_2_B_out.txt', 'w')
while S.status.iter < S.params.maxiter and L.status.residual > L.params.restol:
# this is where the nodes are actually updated according to the SDC formulas
L.sweep.update_nodes()
# compute/update the residual
L.sweep.compute_residual()
# increment the iteration counter
S.status.iter += 1
out = 'Time %4.2f of %s -- Iteration: %2i -- Residual: %12.8e' % (
L.time,
L.level_index,
S.status.iter,
L.status.residual,
)
f.write(out + '\n')
print(out)
f.close()
# compute the interval's endpoint: this (and only this) will set uend, depending on the collocation nodes
L.sweep.compute_end_point()
# update the simulation time
L.status.time += L.dt
# compute exact solution and compare
uex = P.u_exact(L.status.time)
err = abs(uex - L.uend)
return err, L.status.residual, S.status.iter
if __name__ == "__main__":
main()
Results:
Time 0.10 of 0 -- Iteration: 1 -- Residual: 4.11190756e-03
Time 0.10 of 0 -- Iteration: 2 -- Residual: 6.68442667e-04
Time 0.10 of 0 -- Iteration: 3 -- Residual: 8.80377591e-05
Time 0.10 of 0 -- Iteration: 4 -- Residual: 1.21707909e-05
Time 0.10 of 0 -- Iteration: 5 -- Residual: 1.38272147e-06
Time 0.10 of 0 -- Iteration: 6 -- Residual: 6.36445413e-07
Time 0.10 of 0 -- Iteration: 7 -- Residual: 1.68953216e-07
Time 0.10 of 0 -- Iteration: 8 -- Residual: 3.52601840e-08
Time 0.10 of 0 -- Iteration: 9 -- Residual: 6.07249025e-09
Time 0.10 of 0 -- Iteration: 10 -- Residual: 8.27343378e-10
Time 0.10 of 0 -- Iteration: 11 -- Residual: 1.18931339e-10
Time 0.10 of 0 -- Iteration: 12 -- Residual: 1.48499772e-11
Part C: Using pySDC’s frontend¶
Finally, we arrive at the user-friendliest interface pySDC has to offer.
We use one of the controller implementations to do the whole
iteration logic for us, namely controller_nonMPI. It can
do SDC, multi-level SDC, multi-step SDC and PFASST, depending on the
input dictionary (more on this in later tutorials). This is the default
controller and does not require mpi4py to work.
Important things to note:
By using one of the controllers, the whole code relevant for the user is reduced to setting up the
descriptiondictionary, some pre- and some post-processing.During initialization, the parameters used for the run are printed out. Also, user-defined/-changed parameters are indicated. This can be suppressed by setting the controller parameter
dump_setupto False.We make use of
controller_parametersin order to provide logging to file capabilities.In contrast to Part B, we do not have direct access to residuals or iteration counts for now. We will deal with these later.
This example is the prototype for a user to work with pySDC. Most of the logic and most of the data structures are hidden, but all relevant parameters are accessible using the
description.
Full code: pySDC/tutorial/step_2/C_using_pySDCs_frontend.py
from pathlib import Path
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced
from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
def main():
"""
A simple test program to run IMEX SDC for a single time step
"""
# initialize level parameters
level_params = dict()
level_params['restol'] = 1e-10
level_params['dt'] = 0.1
# initialize sweeper parameters
sweeper_params = dict()
sweeper_params['quad_type'] = 'RADAU-RIGHT'
sweeper_params['num_nodes'] = 3
# initialize problem parameters
problem_params = dict()
problem_params['nu'] = 0.1 # diffusion coefficient
problem_params['freq'] = 4 # frequency for the test value
problem_params['nvars'] = 1023 # number of degrees of freedom
problem_params['bc'] = 'dirichlet-zero' # boundary conditions
# initialize step parameters
step_params = dict()
step_params['maxiter'] = 20
# initialize controller parameters
controller_params = dict()
controller_params['log_to_file'] = True
controller_params['fname'] = 'data/step_2_C_out.txt'
# Fill description dictionary for easy hierarchy creation
description = dict()
description['problem_class'] = heatNd_forced
description['problem_params'] = problem_params
description['sweeper_class'] = imex_1st_order
description['sweeper_params'] = sweeper_params
description['level_params'] = level_params
description['step_params'] = step_params
Path("data").mkdir(parents=True, exist_ok=True)
# instantiate the controller
controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
# set time parameters
t0 = 0.1
Tend = 0.3 # note that we are requesting 2 time steps here (dt is 0.1)
# get initial values on finest level
P = controller.MS[0].levels[0].prob
uinit = P.u_exact(t0)
# call main function to get things done...
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
# compute exact solution and compare
uex = P.u_exact(Tend)
err = abs(uex - uend)
f = open('data/step_2_C_out.txt', 'a')
out = 'Error after SDC iterations: %8.6e' % err
f.write(out)
print(out)
f.close()
assert err <= 2e-5, "ERROR: controller doing IMEX SDC iteration did not reduce the error enough, got %s" % err
if __name__ == "__main__":
main()
Results:
2024-04-03 16:19:19,227 - controller - Controller - welcome_message - 146 - INFO: Welcome to the one and only, really very astonishing and 87.3% bug free
_____ _____ _____
/ ____| __ \ / ____|
_ __ _ _| (___ | | | | |
| '_ \| | | |\___ \| | | | |
| |_) | |_| |____) | |__| | |____
| .__/ \__, |_____/|_____/ \_____|
| | __/ |
|_| |___/
2024-04-03 16:19:19,227 - controller - Controller - dump_setup - 160 - INFO: Setup overview (--> user-defined, -> dependency) -- BEGIN
2024-04-03 16:19:19,227 - controller - Controller - dump_setup - 227 - INFO: ----------------------------------------------------------------------------------------------------
Controller: <class 'pySDC.implementations.controller_classes.controller_nonMPI.controller_nonMPI'>
all_to_done = False
dump_setup = True
--> fname = data/step_2_C_out.txt
--> hook_class = [<class 'pySDC.implementations.hooks.default_hook.DefaultHooks'>]
--> log_to_file = True
logger_level = 20
mssdc_jac = True
predict_type = None
use_iteration_estimator = False
Step: <class 'pySDC.core.Step.step'>
--> maxiter = 20
Number of steps: None
Level: <class 'pySDC.core.Level.level'>
Level 0
--> dt = 0.1
dt_initial = 0.1
nsweeps = 1
residual_type = full_abs
--> restol = 1e-10
--> Problem: <class 'pySDC.implementations.problem_classes.HeatEquation_ND_FD.heatNd_forced'>
--> bc = dirichlet-zero
--> freq = (4,)
liniter = 10000
lintol = 1e-12
--> nu = 0.1
--> nvars = (1023,)
order = 2
sigma = 0.06
solver_type = direct
stencil_type = center
--> Data type u: <class 'pySDC.implementations.datatype_classes.mesh.mesh'>
--> Data type f: <class 'pySDC.implementations.datatype_classes.mesh.imex_mesh'>
--> Sweeper: <class 'pySDC.implementations.sweeper_classes.imex_1st_order.imex_1st_order'>
QE = EE
QI = IE
do_coll_update = False
initial_guess = spread
--> num_nodes = 3
--> quad_type = RADAU-RIGHT
skip_residual_computation = ()
--> Collocation: <class 'pySDC.core.Collocation.CollBase'>
Active convergence controllers:
| # | order | convergence controller
----+----+-------+---------------------------------------------------------------------------------------
| 0 | 95 | BasicRestartingNonMPI
-> | 1 | 100 | SpreadStepSizesBlockwiseNonMPI
| 2 | 200 | CheckConvergence
2024-04-03 16:19:19,227 - controller - Controller - dump_setup - 230 - INFO: ----------------------------------------------------------------------------------------------------
2024-04-03 16:19:19,227 - controller - Controller - dump_setup - 232 - INFO: Setup overview (--> user-defined, -> dependency) -- END
2024-04-03 16:19:19,231 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 1 -- Sweep: 1 -- residual: 4.11190756e-03
2024-04-03 16:19:19,235 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 2 -- Sweep: 1 -- residual: 6.68442667e-04
2024-04-03 16:19:19,238 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 3 -- Sweep: 1 -- residual: 8.80377591e-05
2024-04-03 16:19:19,241 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 4 -- Sweep: 1 -- residual: 1.21707909e-05
2024-04-03 16:19:19,245 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 5 -- Sweep: 1 -- residual: 1.38272147e-06
2024-04-03 16:19:19,248 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 6 -- Sweep: 1 -- residual: 6.36445413e-07
2024-04-03 16:19:19,251 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 7 -- Sweep: 1 -- residual: 1.68953216e-07
2024-04-03 16:19:19,255 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 8 -- Sweep: 1 -- residual: 3.52601840e-08
2024-04-03 16:19:19,258 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 9 -- Sweep: 1 -- residual: 6.07249025e-09
2024-04-03 16:19:19,262 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 10 -- Sweep: 1 -- residual: 8.27343378e-10
2024-04-03 16:19:19,265 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 11 -- Sweep: 1 -- residual: 1.18931339e-10
2024-04-03 16:19:19,268 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.100000 at stage IT_FINE: Level: 0 -- Iteration: 12 -- Sweep: 1 -- residual: 1.48499772e-11
2024-04-03 16:19:19,273 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 1 -- Sweep: 1 -- residual: 6.69984764e-03
2024-04-03 16:19:19,276 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 2 -- Sweep: 1 -- residual: 1.05518433e-03
2024-04-03 16:19:19,279 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 3 -- Sweep: 1 -- residual: 1.40642621e-04
2024-04-03 16:19:19,283 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 4 -- Sweep: 1 -- residual: 1.85982063e-05
2024-04-03 16:19:19,286 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 5 -- Sweep: 1 -- residual: 2.79216702e-06
2024-04-03 16:19:19,289 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 6 -- Sweep: 1 -- residual: 1.12278839e-06
2024-04-03 16:19:19,293 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 7 -- Sweep: 1 -- residual: 2.85495353e-07
2024-04-03 16:19:19,296 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 8 -- Sweep: 1 -- residual: 5.78947003e-08
2024-04-03 16:19:19,299 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 9 -- Sweep: 1 -- residual: 9.68230621e-09
2024-04-03 16:19:19,303 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 10 -- Sweep: 1 -- residual: 1.26313315e-09
2024-04-03 16:19:19,306 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 11 -- Sweep: 1 -- residual: 1.82951499e-10
2024-04-03 16:19:19,309 - hooks - default_hook - post_sweep - 170 - INFO: Process 0 on time 0.200000 at stage IT_FINE: Level: 0 -- Iteration: 12 -- Sweep: 1 -- residual: 1.99691114e-11
2024-04-03 16:19:19,310 - hooks - default_hook - post_run - 340 - INFO: Finished run after 0.08s
Error after SDC iterations: 1.166689e-05