Coverage for pySDC/tutorial/step_3/D_writing_solutions_to_file.py: 100%
50 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-29 12:50 +0000
1# ---
2# jupyter:
3# jupytext:
4# formats: py:percent
5# kernelspec:
6# display_name: Python 3
7# name: python3
8# ---
10# %% [markdown]
11# # Part D: Writing solutions to file
12#
13# The statistics of Parts A to C live in memory, which is fine for a few numbers per step. The solution of a PDE on
14# a fine grid, recorded for a long run, is better written to disk. pySDC has a hook for that: `LogToFile` writes the
15# solution at regular intervals to a single binary file, using the file handlers of
16# {py:mod}`pySDC.helpers.fieldsIO`. The same file can be read back for analysis, or to restart a run where it
17# stopped.
18#
19# The problem class decides what goes into the file, with two methods: `getOutputFile` creates the file and writes
20# its header, and `processSolutionForOutput` turns a solution into the array that is stored. The problems that
21# implement them are `testequation0d` and those built on `GenericSpectralLinear`, such as the heat, Burgers and
22# Rayleigh-Bénard equations.
23#
24# ## A heat equation
25#
26# We take `Heat2DChebychev`, the heat equation on a grid that is periodic in $x$ and has Dirichlet boundaries in
27# $y$. The boundary values are $0$ at the bottom and $1$ at the top, so a sine mode decays towards the linear
28# profile in between. The problem is written in first-order form, with the derivatives $u_x$ and $u_y$ as extra
29# components, which makes three variables per grid point.
31# %%
32import matplotlib.pyplot as plt
33import numpy as np
34from pathlib import Path
36from pySDC.helpers.fieldsIO import FieldsIO
37from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
38from pySDC.implementations.hooks.log_solution import LogToFile
39from pySDC.implementations.problem_classes.HeatEquation_Chebychev import Heat2DChebychev
40from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
42# initialize level parameters
43level_params = {'restol': 1e-10, 'dt': 0.05}
45# initialize sweeper parameters
46sweeper_params = {'quad_type': 'RADAU-RIGHT', 'num_nodes': 3, 'QI': 'LU'}
48# initialize problem parameters: 32 x 33 points, u = 0 at the bottom and u = 1 at the top
49problem_params = {'nx': 32, 'ny': 33, 'a': 0, 'b': 0, 'c': 1, 'nu': 0.1}
51# initialize step parameters
52step_params = {'maxiter': 20}
54# Fill description dictionary for easy hierarchy creation
55description = {
56 'problem_class': Heat2DChebychev,
57 'problem_params': problem_params,
58 'sweeper_class': generic_implicit,
59 'sweeper_params': sweeper_params,
60 'level_params': level_params,
61 'step_params': step_params,
62}
64Path("data").mkdir(parents=True, exist_ok=True)
66# %% [markdown]
67# ## Configuring the hook
68#
69# `LogToFile` is configured with class attributes. We set them on a subclass of our own, so that other runs in the
70# same Python process keep the defaults:
71#
72# - `filename`: the file to write,
73# - `time_increment`: the time between two solutions in the file, here two steps,
74# - `allow_overwriting`: whether an existing file may be replaced. We allow it, so that this page can run again.
75#
76# Besides the solutions every `time_increment`, the hook writes the initial conditions and the solution at the end.
79# %%
80class MyLogToFile(LogToFile):
81 filename = 'data/step_3_D_heat.pysdc'
82 time_increment = 0.1
83 allow_overwriting = True
86# initialize controller parameters
87controller_params = {
88 'hook_class': MyLogToFile,
89 'logger_level': 30, # reduce verbosity, the hook tells what it writes on level 20
90}
92# %% [markdown]
93# We run up to $t=0.5$ for now:
95# %%
96controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
97P = controller.MS[0].levels[0].prob
98uend, stats = controller.run(u0=P.u_exact(0), t0=0.0, Tend=0.5)
100# %% [markdown]
101# ## Reading the file
102#
103# `FieldsIO.fromFile` reads any file written this way. From the header it finds out what kind of file it is (here a
104# `Rectilinear` grid) and returns the matching handler:
106# %%
107file = FieldsIO.fromFile(MyLogToFile.filename)
108print(file)
109print('variables:', file.nVar, ' grid:', file.gridSizes)
110print('times:', np.round(file.times, 12))
112# %% [markdown]
113# `readField` returns the time and the solution at an index in the file. The solution has one array per variable,
114# on the grid whose coordinates are stored in the header:
116# %%
117t, u = file.readField(-1)
118print(f't = {t:.2f}, u has shape {u.shape}')
119x, y = file.header['coords']
120print(f'x from {x.min():.2f} to {x.max():.2f}, y from {y.min():.2f} to {y.max():.2f}')
122# %% [markdown]
123# ## Restarting from the file
124#
125# Now we carry on to $t=1$. A new run that starts at $t_0 > 0$ appends to an existing file instead of replacing it,
126# whatever `allow_overwriting` says. `load` gets the solution at an index of the file back. What was written is the
127# output of `processSolutionForOutput`, the real solution on the grid, which is what this problem computes with anyway,
128# so it can go straight into a new initial condition:
130# %%
131restart = MyLogToFile.load(-1)
132u0 = P.u_init
133u0[:] = restart['u']
135controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
136uend, stats = controller.run(u0=u0, t0=restart['t'], Tend=1.0)
138file = FieldsIO.fromFile(MyLogToFile.filename)
139print('times:', np.round(file.times, 12))
141# %% [markdown]
142# The file now holds the whole run. Here is $u$, the first variable, at three of its times:
144# %% tags=["hide-input"]
145fig, axs = plt.subplots(1, 3, figsize=(9, 3), sharey=True)
146for ax, idx in zip(axs, [0, 5, 10], strict=True):
147 t, u = file.readField(idx)
148 im = ax.pcolormesh(x, y, u[0].T, vmin=-0.5, vmax=1.5, cmap='plasma', shading='gouraud')
149 ax.set_title(f'$t={t:.1f}$')
150 ax.set_xlabel('$x$')
151axs[0].set_ylabel('$y$')
152fig.colorbar(im, ax=axs, label='$u$')
154# %% [markdown]
155# The sine mode decays and leaves the linear profile between the boundaries. For this problem the exact solution is
156# known, so we can check the last solution in the file:
158# %%
159t, u = file.readField(-1)
160err = np.max(np.abs(u[0] - P.u_exact(t)[0]))
161print(f'error at t={t:.1f}: {err:.3e}')
163# %% [markdown]
164# ## The file format
165#
166# A file starts with a header: the kind of file and the data type of the values, then, for a `Rectilinear` grid, the
167# number of variables, the dimension, the number of points and the coordinates in each direction. After that come
168# the solutions, each one a time followed by the values of all variables on the whole grid. As each solution takes
169# the same number of bytes, `readField` jumps straight to the one it needs, and appending is cheap.
170#
171# The handlers are not tied to pySDC's controller. {py:mod}`pySDC.helpers.fieldsIO` shows how to write and read a
172# file by hand, and `Rectilinear.toVTR` converts the solutions of a 3D grid into VTR files for ParaView.
173#
174# ## With MPI
175#
176# When the problem is distributed over several processes in space, all of them write into the same file, each its
177# own part of the grid, with collective MPI-IO. The problem sets this up in `setUpFieldsIO`, which tells
178# `Rectilinear.setupMPI` the part of the grid of each process. The header always holds the global grid, so a file
179# written by any number of processes can be read by any other number, in serial as well. The spectral problems do all
180# this by themselves, and the hook is used exactly as above.
181#
182# :::{admonition} Important things to note
183# - Configure `LogToFile` with class attributes on a subclass, so that other runs keep the defaults.
184# - Without `allow_overwriting`, the hook will neither replace an existing file nor write a solution at a time
185# the file already has.
186# - To write another problem to file, implement `getOutputFile` and `processSolutionForOutput` in its class, and,
187# for MPI, `setUpFieldsIO`.
188# :::
189#
190# The checks the tests run:
192# %%
193assert np.allclose(file.times, np.arange(11) * 0.1), f'ERROR: unexpected times in the file: {file.times}'
194assert err < 1e-8, f'ERROR: solution is not as exact as expected, got {err}'