From 9c7da1a73b75cedc003835b442c7c3e5692c28a8 Mon Sep 17 00:00:00 2001 From: Jorgen Schartum Dokken Date: Wed, 11 Nov 2026 22:29:46 +0000 Subject: [PATCH] Improve time-stepping by having an initialize, restart and "continue" phase for a simulation, which means that we can use BDF2 and BDF3 with adaptive steps. --- docs/tutorials/04_time_stepping.py | 62 +++++-- src/sumfem/solvers/simulation.py | 254 ++++++++++++++++++--------- src/sumfem/solvers/time_stepping.py | 113 +++++++++--- tests/discrete/test_simulation.py | 35 ++++ tests/discrete/test_time_stepping.py | 138 ++++++++++++++- 5 files changed, 470 insertions(+), 132 deletions(-) diff --git a/docs/tutorials/04_time_stepping.py b/docs/tutorials/04_time_stepping.py index 1c6feb6..7838131 100644 --- a/docs/tutorials/04_time_stepping.py +++ b/docs/tutorials/04_time_stepping.py @@ -180,30 +180,65 @@ def run(scheme: solvers.TimeScheme, dt: float, t_end: float = 1.0) -> float: print(f"Crank-Nicolson: errors {errors}, rates {rates.round(2)}") assert np.all(rates > 2 - 0.15) +# ## Runs in parts +# +# A simulation may run in parts, e.g. to change a load or the step between +# them. The first run starts the history from the current values, and later +# runs continue it, so BDF2 keeps its second order. Its weights follow the +# sizes of the latest steps, so the step may change between the runs: here +# it is halved halfway, and the errors still fall at second order. +# {py:meth}`restart ` starts afresh +# instead, e.g. from new initial values: + + +# + +def run_in_parts(dt: float) -> tuple[float, int]: + disc = discrete.discretize([model, clock], discrete.DomainMaps({domain: marked})) + disc.set({clock.dt: dt}) + simulation = solvers.Simulation( + disc, model, bcs, clock=clock, scheme=solvers.BDF(2) + ) + disc.set({u: spatial}) + simulation.run(t_end=0.5) + disc.set({clock.dt: dt / 2}) + simulation.run(t_end=1.0) + e = error(disc[u], lambda x: (1 + x[0] ** 2 + 2 * x[1] ** 2) * (1 + np.sin(1.0))) + return e, simulation.steps_taken + + +results = [run_in_parts(dt) for dt in (0.1, 0.05, 0.025)] +errors = np.array([e for e, _ in results]) +rates = np.log2(errors[:-1] / errors[1:]) +print(f"BDF2 in two parts: steps {[n for _, n in results]}, rates {rates.round(2)}") +assert np.all(rates > 2 - 0.15) +# - + # ## Output and metadata # # With an {py:class}`Output `, a simulation writes its # fields with ADIOS2 and a `metadata.json` that records the versions, the # parameters with their units, the model, the time stepping and the meshes. -# A callback is called after every step, e.g. to post-process. The output -# goes to a temporary directory here, removed at the end: +# The output stays open across the runs of a simulation, until it is closed +# with {py:meth}`close ` or at the end of a +# `with` block. A callback is called after every step, e.g. to post-process. +# The output goes to a temporary directory here, removed at the end: with tempfile.TemporaryDirectory() as temporary: directory = pathlib.Path(MPI.COMM_WORLD.bcast(temporary, root=0)) disc = discrete.discretize([model], discrete.DomainMaps({domain: marked})) - simulation = solvers.Simulation( + with solvers.Simulation( disc, model, bcs, clock=clock, scheme=solvers.BDF(2), output=solvers.Output(directory, every=5), - ) - disc.set({u: spatial}) - simulation.run( - t_end=1.0, - callback=lambda t, step: step % 5 or print(f"t = {t:.1f}, step {step}"), - ) + ) as simulation: + disc.set({u: spatial}) + simulation.run( + t_end=1.0, + callback=lambda t, step: step % 5 or print(f"t = {t:.1f}, step {step}"), + ) MPI.COMM_WORLD.barrier() files = sorted(p.name for p in directory.iterdir()) metadata = json.loads((directory / "metadata.json").read_text()) @@ -219,16 +254,17 @@ def run(scheme: solvers.TimeScheme, dt: float, t_end: float = 1.0) -> float: # grows the step while Newton's method converges easily and cuts it when it # struggles, retrying failed steps, by the rules of an # {py:class}`AdaptiveTimeStep `. Steps end -# at milestones, e.g. output times, and at the end time. This problem is -# linear, so every step takes one iteration, and the steps grow up to the -# largest one allowed: +# at milestones, e.g. output times, and at the end time. BDF2 adapts its +# weights to the changing steps, as above. This problem is linear, so every +# step takes one iteration, and the steps grow up to the largest one +# allowed: disc = discrete.discretize([model], discrete.DomainMaps({domain: marked})) rules = solvers.AdaptiveTimeStep( initial=0.01, growth_factor=1.5, target_iterations=4, max_step=0.2, milestones=[0.5] ) adaptive = solvers.AdaptiveSimulation( - disc, model, rules, clock=clock, scheme=solvers.BDF(1), bcs=bcs + disc, model, rules, clock=clock, scheme=solvers.BDF(2), bcs=bcs ) disc.set({u: spatial}) times: list[float] = [] diff --git a/src/sumfem/solvers/simulation.py b/src/sumfem/solvers/simulation.py index c7709ef..d4ecf03 100644 --- a/src/sumfem/solvers/simulation.py +++ b/src/sumfem/solvers/simulation.py @@ -11,6 +11,7 @@ import pathlib import platform import typing +import weakref from collections.abc import Iterable, Mapping import petsc4py @@ -32,7 +33,7 @@ from ..symbolic.terminals import Parameter from ..symbolic.time import Clock, time_derivative_expressions, time_derivative_fields from .solver import DEFAULT_PETSC_OPTIONS, FACTORIZED_PETSC_OPTIONS, Solver -from .time_stepping import BDF, AdaptiveTimeStep, TimeScheme +from .time_stepping import AdaptiveTimeStep, TimeScheme Callback = typing.Callable[[float, int], None] """Called after a step with the time and the step number.""" @@ -64,6 +65,18 @@ class Simulation: the history. A steady problem needs no simulation: solve the model with a :class:`~sumfem.solvers.solver.Solver`. + A simulation may run in parts, e.g. to change a load or the step size + between them. The first run starts the history from the current values; + later runs continue from where the last one ended, with the history, the + numbering of the steps and the output of the earlier runs, so a + multistep scheme keeps its order across them. :meth:`restart` starts + afresh, e.g. from new initial values. + + Output stays open across the runs. Close it with :meth:`close`, or by + using the simulation in a ``with`` block. Every process must close it at + the same point, as writing is collective; an output left open is closed + when the simulation is garbage collected or the program ends. + Args: discretization: The meshes and the functions of all terminals. model: The model, solved monolithically. @@ -122,6 +135,18 @@ class Simulation: _time_scale: float """Reference time of the time derivatives, in seconds.""" + _started: bool + """Whether a run has started the history.""" + + _steps_taken: int + """Number of steps since the start.""" + + _step_sizes: list[float] + """Sizes of the latest steps, most recent first.""" + + _writer: _Writer | None + """The open output, if any.""" + def __init__( self, discretization: Discretization, @@ -147,6 +172,10 @@ def __init__( self._scheme = scheme self._petsc_options = dict(petsc_options or {}) self._output = output + self._started = False + self._steps_taken = 0 + self._step_sizes = [] + self._writer = None bcs = list(bcs) discretization.add(clock) discretization.add(*(bc for bc in bcs if isinstance(bc, symbolic.DirichletBC))) @@ -186,20 +215,50 @@ def solver(self) -> Solver: """The solver of one step.""" return self._solver + @property + def steps_taken(self) -> int: + """Number of steps since the start, over all runs.""" + return self._steps_taken + + def restart(self) -> None: + """Start afresh at the next run. + + The next run starts the history from the current values, e.g. new + initial values, numbers its steps from one, and opens the output + again, which overwrites the earlier files. + """ + self.close() + self._started = False + self._steps_taken = 0 + self._step_sizes = [] + + def close(self) -> None: + """Close the output. A later run opens it again, overwriting it.""" + if self._writer is not None: + self._writer.close() + self._writer = None + + def __enter__(self) -> Simulation: + return self + + def __exit__(self, *exc: typing.Any) -> None: + self.close() + def run( self, num_steps: int | None = None, t_end: float | None = None, callback: Callback | None = None, ) -> None: - """Take time steps of the clock's step size. + """Take time steps of the clock's step size, from the clock's time. Args: num_steps: Number of time steps. t_end: End time, a whole number of steps after the start. Give exactly one of ``num_steps`` and ``t_end``. callback: Called after each step with the time and the step - number, e.g. to post-process. + number, counted over all runs since the start, e.g. to + post-process. Raises: ValueError: For inconsistent time arguments, or without a time @@ -208,39 +267,59 @@ def run( disc, clock = self._discretization, self._clock if not disc.has_value(clock.dt): raise ValueError(f"Set the time step {clock.dt.name!r} before running.") - t0 = self._start() dt = float(disc.value(clock.dt)) + t0 = float(disc.value(clock.t)) num_steps = _number_of_steps(t0, dt, num_steps, t_end) - with _Writer(self) as writer: - writer.write(t0) - for step in range(1, num_steps + 1): - t = t0 + step * dt - self._set_time(step, t, dt) - self._solver.solve() - self._finish_step(t, step, writer, callback) + self._start() + for k in range(1, num_steps + 1): + t = t0 + k * dt + self._set_time(t, dt) + self._solver.solve() + self._finish_step(t, dt, callback) + self._end() def _start(self) -> float: - """Start a run: set up the history, and return the time.""" - disc = self._discretization - if self._scheme is None or not self._scheme.given_history: - disc.initialize_history() - return float(disc.value(self._clock.t)) + """Start a run, and return the time. - def _set_time(self, step: int, t: float, dt: float) -> None: - """Set the time and step size of a step, and the scheme's weights.""" + The first run after the start sets up the history, unless the scheme + has a given history, and opens the output with the initial values. + """ + disc = self._discretization + t = float(disc.value(self._clock.t)) + if not self._started: + if self._scheme is None or not self._scheme.given_history: + disc.initialize_history() + self._started = True + if self._output is not None: + self._writer = _Writer(self._output, disc) + self._writer.open(self.metadata()) + self._writer.write(t) + return t + + def _set_time(self, t: float, dt: float) -> None: + """Set the time and step size of the next step, and the scheme's weights.""" disc = self._discretization disc.set({self._clock.t: t, self._clock.dt: dt}) if self._scheme is not None: - self._scheme.set_step(step, dt, disc) + self._scheme.set_step( + self._steps_taken + 1, dt, disc, tuple(self._step_sizes) + ) - def _finish_step( - self, t: float, step: int, writer: _Writer, callback: Callback | None - ) -> None: + def _finish_step(self, t: float, dt: float, callback: Callback | None) -> None: """Advance the history, write, and call back after a step.""" self._discretization.advance_history() - writer.write(t, step) + self._steps_taken += 1 + # BDF3, the longest history, needs the two steps before the next one + self._step_sizes = [dt, *self._step_sizes][:2] + if self._writer is not None: + self._writer.write(t, self._steps_taken) if callback is not None: - callback(t, step) + callback(t, self._steps_taken) + + def _end(self) -> None: + """End a run: write the files that are rewritten after every run.""" + if self._writer is not None: + self._writer.flush() def metadata(self) -> dict[str, typing.Any]: """A description of the simulation, for archiving with its results.""" @@ -304,14 +383,17 @@ class AdaptiveSimulation(Simulation): model: See :class:`Simulation`. time_step: The rules of the steps. clock: The time and time step. - scheme: A one-step scheme: BDF of order one, or a theta method. + scheme: A time-stepping scheme. BDF of order two or three adapts + its weights to the changing steps; it stays stable while the + steps grow slowly, see :class:`~sumfem.solvers.time_stepping.BDF`. **kwargs: The other arguments of :class:`Simulation`. The ``time_scale`` defaults to the largest step, if it is a number, so that the absolute tolerance of Newton's method is no looser for any step than for the largest; otherwise to the first step. - Raises: - ValueError: With a scheme of more than one step. + Like a :class:`Simulation`, it may run in parts. Each run starts with + the initial step of the rules, as a run usually starts at a change, e.g. + of a load, that a long step could miss. """ _time_step: AdaptiveTimeStep @@ -326,11 +408,6 @@ def __init__( scheme: TimeScheme | None = None, **kwargs: typing.Any, ): - if isinstance(scheme, BDF) and scheme.order > 1: - raise ValueError( - f"BDF of order {scheme.order} assumes equal steps. Adapt the steps " - "with BDF(1) or a theta method." - ) if kwargs.get("time_scale") is None: largest = time_step.max_step kwargs["time_scale"] = ( @@ -344,12 +421,13 @@ def __init__( def run( # type: ignore[override] self, t_end: float, callback: Callback | None = None ) -> None: - """Step to the end time. + """Step from the clock's time to the end time. Args: t_end: End time. callback: Called after each step with the time and the step - number, e.g. to post-process. + number, counted over all runs since the start, e.g. to + post-process. Raises: RuntimeError: If a step fails after its retries. @@ -357,30 +435,28 @@ def run( # type: ignore[override] rules = self._time_step t = self._start() dt = rules.initial - step = 0 - with _Writer(self) as writer: - writer.write(t) - while t < t_end and not math.isclose( - t, t_end, rel_tol=rules.milestone_tolerance - ): - dt = rules.limit(dt, t, t_end) - before = [owned(f).copy() for f in self._solver.system.solution] - for attempt in range(rules.max_retries + 1): - self._set_time(step + 1, t + dt, dt) - if self._try_step(): - break - self._restore(before) - retried = rules.retry(dt) - if attempt == rules.max_retries or retried is None: - raise RuntimeError( - f"The step from t = {t} failed with dt = {dt}, after " - f"{attempt} retries." - ) - dt = retried - t += dt - step += 1 - self._finish_step(t, step, writer, callback) - dt = rules.next(dt, self._solver.iterations) + while t < t_end and not math.isclose( + t, t_end, rel_tol=rules.milestone_tolerance + ): + dt = rules.limit(dt, t, t_end) + before = [owned(f).copy() for f in self._solver.system.solution] + for attempt in range(rules.max_retries + 1): + self._set_time(t + dt, dt) + if self._try_step(): + break + self._restore(before) + retried = rules.retry(dt) + if attempt == rules.max_retries or retried is None: + self._end() + raise RuntimeError( + f"The step from t = {t} failed with dt = {dt}, after " + f"{attempt} retries." + ) + dt = retried + t += dt + self._finish_step(t, dt, callback) + dt = rules.next(dt, self._solver.iterations) + self._end() def _try_step(self) -> bool: """Solve a step, and whether the solve converged.""" @@ -453,25 +529,23 @@ def _number_of_steps( class _Writer: - """Writes fields and metadata for a simulation, if it has an output.""" + """Writes the fields of a discretization, and metadata, to an output.""" - def __init__(self, simulation: Simulation): - self._simulation = simulation - self._output = simulation._output + def __init__(self, output: Output, discretization: Discretization): + self._output = output + self._discretization = discretization + self._comm = next(iter(discretization.fields.values())).function_space.mesh.comm self._vtx: list[dolfinx.io.VTXWriter] = [] self._real: dict[str, dolfinx.fem.Function] = {} self._real_rows: list[list[float]] = [] + self._finalizer: weakref.finalize | None = None - def __enter__(self) -> _Writer: - output = self._output - if output is None: - return self - disc = self._simulation.discretization - comm = next(iter(disc.fields.values())).function_space.mesh.comm + def open(self, metadata: dict[str, typing.Any]) -> None: + """Create the directory, write the metadata and open the writers.""" + output, disc, comm = self._output, self._discretization, self._comm directory = pathlib.Path(output.directory) if comm.rank == 0: directory.mkdir(parents=True, exist_ok=True) - metadata = self._simulation.metadata() (directory / "metadata.json").write_text(json.dumps(metadata, indent=2)) comm.barrier() names = sorted(disc.fields) if output.fields is None else list(output.fields) @@ -484,10 +558,12 @@ def __enter__(self) -> _Writer: self._vtx.append( dolfinx.io.VTXWriter(comm, directory / f"{name}.bp", [function]) ) - return self + # Close the writers if the simulation is dropped without closing them + self._finalizer = weakref.finalize(self, _close_writers, self._vtx) def write(self, t: float, step: int = 0) -> None: - if self._output is None or step % self._output.every != 0: + """Write the fields at a time, at every chosen step.""" + if step % self._output.every != 0: return for writer in self._vtx: writer.write(t) @@ -496,19 +572,23 @@ def write(self, t: float, step: int = 0) -> None: [t, *(float(f.x.array[0]) for f in self._real.values())] ) - def __exit__(self, *exc: typing.Any) -> None: - for writer in self._vtx: - writer.close() - if self._output is not None and self._real: - disc = self._simulation.discretization - comm = next(iter(disc.fields.values())).function_space.mesh.comm - if comm.rank == 0: - path = pathlib.Path(self._output.directory) / "real_fields.csv" - header = ",".join(["t", *self._real]) - np.savetxt( - path, - np.array(self._real_rows), - delimiter=",", - header=header, - comments="", - ) + def flush(self) -> None: + """Write the values of the real fields so far.""" + if not self._real or self._comm.rank != 0: + return + path = pathlib.Path(self._output.directory) / "real_fields.csv" + header = ",".join(["t", *self._real]) + np.savetxt( + path, np.array(self._real_rows), delimiter=",", header=header, comments="" + ) + + def close(self) -> None: + """Write the real fields and close the writers.""" + self.flush() + if self._finalizer is not None: + self._finalizer() + + +def _close_writers(writers: list[dolfinx.io.VTXWriter]) -> None: + for writer in writers: + writer.close() diff --git a/src/sumfem/solvers/time_stepping.py b/src/sumfem/solvers/time_stepping.py index 09726d4..24f09ce 100644 --- a/src/sumfem/solvers/time_stepping.py +++ b/src/sumfem/solvers/time_stepping.py @@ -14,6 +14,7 @@ import functools import math import typing +from collections.abc import Sequence import sympy.physics.units as syu import ufl @@ -37,12 +38,30 @@ ] -# Weights of u_{n+1}, u_n, u_{n-1}, ... for each order, to be divided by dt -_BDF_WEIGHTS: dict[int, tuple[float, ...]] = { - 1: (1.0, -1.0), - 2: (1.5, -2.0, 0.5), - 3: (11 / 6, -3.0, 1.5, -1 / 3), -} +_BDF_ORDERS = (1, 2, 3) +"""The orders of the BDF schemes.""" + + +def _bdf_weights(steps: Sequence[float]) -> tuple[float, ...]: + """Weights of ``u_{n+1}, u_n, u_{n-1}, ...`` in a BDF formula. + + The derivative at the new time of the polynomial that interpolates the + values at the new time and at the previous ones, so the formula is exact + for polynomials in time of degree up to the number of steps. + + Args: + steps: The sizes of the current step and of the previous ones, most + recent first; their number is the order. + """ + # The times of the values, relative to the new time + times = [0.0, *(-sum(steps[: j + 1]) for j in range(len(steps)))] + weights = [sum(1.0 / (0.0 - tm) for tm in times[1:])] + for j, tj in enumerate(times[1:], start=1): + others = [tm for m, tm in enumerate(times) if m != j] + numerator = math.prod(0.0 - tm for tm in others[1:]) + weights.append(numerator / math.prod(tj - tm for tm in others)) + return tuple(weights) + Substitutions = dict[ufl.Coefficient, ufl.core.expr.Expr] """Replacements of terminals in the forms of a system.""" @@ -76,13 +95,21 @@ def prepare(self, model: Residual, clock: Clock) -> tuple[Residual, Substitution """ raise NotImplementedError - def set_step(self, step: int, dt: float, discretization: Discretization) -> None: + def set_step( + self, + step: int, + dt: float, + discretization: Discretization, + previous: Sequence[float] = (), + ) -> None: """Set the scheme's parameters for a time step. Args: step: Number of the step, starting at 1. dt: Time step size. discretization: The discretization holding the parameter values. + previous: Sizes of the steps before, most recent first. Missing + ones are taken equal to ``dt``. """ raise NotImplementedError @@ -101,17 +128,28 @@ class BDF(TimeScheme): ``w0 * field + w1 * field.previous(1) + ...``, and the time derivative ``Dt(q)`` of an expression by ``w0 * q + w1 * q(previous(1)) + ...``, the expression of the fields' previous values. The weights are - parameters, set at every step from the time step size. The first steps - of a higher order scheme use the highest order the available history - allows, so one compiled system serves every step. - - BDF3 is not A-stable, so it suits damped problems better than waves. + parameters, set at every step from the sizes of the current and the + previous steps: they differentiate, at the new time, the polynomial + through the values at the times of the steps. The steps may thus change + in size, e.g. between the runs of a + :class:`~sumfem.solvers.simulation.Simulation` or in an + :class:`~sumfem.solvers.simulation.AdaptiveSimulation`; with equal steps + the weights are the classical ones. The first steps of a higher order + scheme use the highest order the available history allows, so one + compiled system serves every step. + + A step that grows fast relative to the previous one can make the higher + orders unstable: BDF2 stays zero-stable while each step is less than + ``1 + sqrt(2)`` times the one before, and BDF3 needs smaller ratios. + Shrinking steps are safe. BDF3 is not A-stable, so it suits damped + problems better than waves. Without a given history, the first steps use lower orders. The first, backward Euler step then limits BDF3 to second order in general. To keep the full order, set the previous values ``field.previous(k)`` before running, e.g. to an exact solution or a finer start-up run, and pass - ``given_history=True``. + ``given_history=True``. The given values are taken one step of the first + step's size apart. Args: order: 1, 2 or 3. @@ -145,10 +183,8 @@ def __init__( blend: float = 1.0, given_history: bool = False, ): - if order not in _BDF_WEIGHTS: - raise ValueError( - f"BDF order must be one of {sorted(_BDF_WEIGHTS)}, got {order}." - ) + if order not in _BDF_ORDERS: + raise ValueError(f"BDF order must be one of {_BDF_ORDERS}, got {order}.") if not 0.0 <= blend <= 1.0: raise ValueError(f"The blend must lie in [0, 1], got {blend}.") if order != 2 and blend != 1.0: @@ -207,19 +243,32 @@ def substitutions(self, model: Residual) -> Substitutions: ) return result - def set_step(self, step: int, dt: float, discretization: Discretization) -> None: + def set_step( + self, + step: int, + dt: float, + discretization: Discretization, + previous: Sequence[float] = (), + ) -> None: """Set the weights for a time step. Args: step: Number of the step, starting at 1. dt: Time step size. discretization: The discretization holding the weight values. + previous: Sizes of the steps before, most recent first. Missing + ones, e.g. of a given history or of a loop that steps with + one size, are taken equal to ``dt``. """ _check_step(step, dt) + if any(not size > 0 for size in previous): + raise ValueError(f"Step sizes must be positive, got {list(previous)}.") order = self._order if self._given_history else min(self._order, step) - values = [*_BDF_WEIGHTS[order], *[0.0] * (self._order - order)] + steps = [dt, *previous[: order - 1]] + steps += [dt] * (order - len(steps)) + values = [*_bdf_weights(steps), *[0.0] * (self._order - order)] if order == 2: - first = [*_BDF_WEIGHTS[1], *[0.0] * (self._order - 1)] + first = [*_bdf_weights([dt]), *[0.0] * (self._order - 1)] values = [ self._blend * a + (1 - self._blend) * b for a, b in zip(values, first, strict=True) @@ -228,9 +277,7 @@ def set_step(self, step: int, dt: float, discretization: Discretization) -> None # Only the weights of domains this discretization's systems use if parameters[0] not in discretization: continue - discretization.set( - {p: v / dt for p, v in zip(parameters, values, strict=True)} - ) + discretization.set(dict(zip(parameters, values, strict=True))) class Theta(TimeScheme): @@ -319,20 +366,28 @@ def prepare(self, model: Residual, clock: Clock) -> tuple[Residual, Substitution return steps[0], {} return ResidualSum(steps, name=model.name), {} - def set_step(self, step: int, dt: float, discretization: Discretization) -> None: + def set_step( + self, + step: int, + dt: float, + discretization: Discretization, + previous: Sequence[float] = (), + ) -> None: + """Set ``1 / dt`` and the previous time; a one-step method needs no + previous step sizes.""" _check_step(step, dt) for rate in self._rates.values(): if rate in discretization: discretization.set({rate: 1 / dt}) - previous = self._previous_time + previous_time = self._previous_time # Only a step whose forms use the time holds the previous time if ( self._clock is not None - and previous is not None - and previous in discretization + and previous_time is not None + and previous_time in discretization ): t = float(discretization.value(self._clock.t)) - discretization.set({previous: t - dt}) + discretization.set({previous_time: t - dt}) class CrankNicolson(Theta): @@ -415,7 +470,7 @@ def residual(self) -> FormOrZero: class AdaptiveTimeStep: - """A time step adapted to the Newton iterations, for one-step schemes. + """A time step adapted to the Newton iterations. After a solve with fewer Newton iterations than ``target_iterations`` the step grows by ``growth_factor``, and after one with more it is cut diff --git a/tests/discrete/test_simulation.py b/tests/discrete/test_simulation.py index 0098bd4..74cef1a 100644 --- a/tests/discrete/test_simulation.py +++ b/tests/discrete/test_simulation.py @@ -78,6 +78,41 @@ def residual(self): assert metadata["meshes"]["cell"]["cell_type"] == "triangle" +def test_output_stays_open_across_runs(toy, square, sides, tmp_path): + """Runs in parts, with a change of step, write one output until closed.""" + domain = Domain("cell", "triangle") + c = Field("c", domain, ("Lagrange", 1)) + mean = Field("mean", domain, "Real") + + class MeanOf(Residual): + """The mean of c, as a Real field.""" + + def residual(self): + return (mean - c) * mean.v * domain.dx() + + model = toy.Accumulation("dc", c) + toy.Diffusion("d", c) + MeanOf("m") + directory = shared_directory(tmp_path) / "parts" + disc = discretize([], DomainMaps({"cell": square()})) + clock = Clock(domain, dt=0.5) + bcs = [DirichletBC(c, sides, 1.0)] + with Simulation( + disc, model, bcs, clock=clock, scheme=BDF(2), output=Output(directory) + ) as simulation: + simulation.run(num_steps=2) + # The real fields so far are written after every run + MPI.COMM_WORLD.barrier() + rows = np.loadtxt(directory / "real_fields.csv", delimiter=",", skiprows=1) + np.testing.assert_allclose(rows[:, 0], [0.0, 0.5, 1.0]) + disc.set({clock.dt: 0.25}) + simulation.run(num_steps=2) + MPI.COMM_WORLD.barrier() + rows = np.loadtxt(directory / "real_fields.csv", delimiter=",", skiprows=1) + np.testing.assert_allclose(rows[:, 0], [0.0, 0.5, 1.0, 1.25, 1.5]) + adios2 = pytest.importorskip("adios2") + with adios2.FileReader(str(directory / "c.bp")) as reader: + assert reader.num_steps() == 5 + + def test_output_writes_only_the_chosen_fields(toy, square, sides, tmp_path): domain = Domain("cell", "triangle") a = Field("a", domain, ("Lagrange", 1)) diff --git a/tests/discrete/test_time_stepping.py b/tests/discrete/test_time_stepping.py index 0f5d49c..0d7c5df 100644 --- a/tests/discrete/test_time_stepping.py +++ b/tests/discrete/test_time_stepping.py @@ -20,6 +20,7 @@ Solver, Theta, ) +from sumfem.solvers.time_stepping import _bdf_weights from sumfem.symbolic import ( Clock, DirichletBC, @@ -371,9 +372,140 @@ def failing_solve(self): ) -def test_adapted_steps_need_a_one_step_scheme(toy, square, sides): - with pytest.raises(ValueError, match="equal steps"): - adaptive_heat(toy, square, sides, AdaptiveTimeStep(0.1), scheme=BDF(2)) +def exact_heat(toy, square, sides, simulation): + """The problem of :func:`solve_heat`, run by ``simulation``. + + Returns the simulation, the discretization, the clock and the L2 error + at the clock's time. + """ + domain = Domain("cell", "triangle") + u = Field("u", domain, ("Lagrange", 2)) + clock = Clock(domain, dt=0.1) + x = ufl.SpatialCoordinate(domain) + exact = spatial(x) * (1 + ufl.sin(clock.t)) + source = spatial(x) * ufl.cos(clock.t) - 6 * (1 + ufl.sin(clock.t)) + model = ( + toy.Accumulation("dudt", u) + + toy.Diffusion("diffusion", u) + + toy.Source("source", u, f=source) + ) + disc = discretize([], DomainMaps({"cell": square()})) + run = simulation(disc, model, [DirichletBC(u, sides, exact)], clock) + disc.set({u: spatial}) + + def error(): + t = float(disc.value(clock.t)) + solution = disc.fields["u"] + mesh = solution.function_space.mesh + xm = ufl.SpatialCoordinate(mesh) + local = dolfinx.fem.assemble_scalar( + dolfinx.fem.form((solution - spatial(xm) * (1 + np.sin(t))) ** 2 * ufl.dx) + ) + return np.sqrt(mesh.comm.allreduce(local, op=MPI.SUM)) + + return run, disc, clock, error + + +@pytest.mark.parametrize("order", [1, 2, 3]) +def test_variable_step_bdf_weights_differentiate_polynomials_exactly(order): + """The weights of steps of any sizes are exact for degree up to the order.""" + steps = [0.3, 0.1, 0.7][:order] + times = -np.cumsum([0.0, *steps]) + weights = _bdf_weights(steps) + for degree in range(order + 1): + values = times**degree + derivative = 1.0 if degree == 1 else 0.0 + assert np.dot(weights, values) == pytest.approx(derivative, abs=1e-12) + + +def test_equal_steps_give_the_classical_weights(): + np.testing.assert_allclose(_bdf_weights([0.5, 0.5]), np.array([1.5, -2, 0.5]) / 0.5) + np.testing.assert_allclose(_bdf_weights([1.0] * 3), [11 / 6, -3, 1.5, -1 / 3]) + + +def test_runs_in_parts_continue_the_history(toy, square, sides): + """BDF2 run in two parts takes the steps of one run, numbered alike.""" + results, numbers = [], [] + for parts in ([8], [3, 5]): + simulation, disc, _, _ = exact_heat( + toy, + square, + sides, + lambda disc, model, bcs, clock: Simulation( + disc, model, bcs, clock=clock, scheme=BDF(2) + ), + ) + steps: list[int] = [] + for count in parts: + simulation.run(num_steps=count, callback=lambda t, n: steps.append(n)) + results.append(disc.fields["u"].x.array.copy()) + numbers.append(steps) + assert simulation.steps_taken == 8 + np.testing.assert_allclose(results[1], results[0], rtol=1e-12, atol=1e-14) + assert numbers[0] == numbers[1] == list(range(1, 9)) + + +def test_a_restart_starts_afresh(toy, square, sides): + """After a restart, a run from new values steps as a new simulation.""" + results, numbers = [], [] + for restart in (False, True): + simulation, disc, clock, _ = exact_heat( + toy, + square, + sides, + lambda disc, model, bcs, clock: Simulation( + disc, model, bcs, clock=clock, scheme=BDF(2) + ), + ) + if restart: + simulation.run(num_steps=3) + simulation.restart() + disc.set({clock.t: 0.0, disc.terminal("u"): spatial}) + steps: list[int] = [] + simulation.run(num_steps=4, callback=lambda t, n: steps.append(n)) + results.append(disc.fields["u"].x.array.copy()) + numbers.append(steps) + np.testing.assert_allclose(results[1], results[0], rtol=1e-12, atol=1e-14) + assert numbers[0] == numbers[1] == [1, 2, 3, 4] + + +def test_bdf2_keeps_its_order_when_the_step_changes_between_runs(toy, square, sides): + """Steps halved halfway, which constant-step weights would bring to order one.""" + errors = [] + for dt in (0.1, 0.05, 0.025): + simulation, disc, clock, error = exact_heat( + toy, + square, + sides, + lambda disc, model, bcs, clock: Simulation( + disc, model, bcs, clock=clock, scheme=BDF(2) + ), + ) + disc.set({clock.dt: dt}) + simulation.run(t_end=0.4) + disc.set({clock.dt: dt / 2}) + simulation.run(t_end=1.0) + errors.append(error()) + rates = np.log2(np.array(errors[:-1]) / np.array(errors[1:])) + assert np.all(rates > 1.85), rates + + +def test_adapted_steps_keep_the_accuracy_of_bdf2(toy, square, sides): + """BDF2 with growing steps is far more accurate than backward Euler.""" + rules = AdaptiveTimeStep(0.01, growth_factor=1.2, target_iterations=2, max_step=0.1) + errors = {} + for order in (1, 2): + simulation, _, _, error = exact_heat( + toy, + square, + sides, + lambda disc, model, bcs, clock, order=order: AdaptiveSimulation( + disc, model, rules, clock, scheme=BDF(order), bcs=bcs + ), + ) + simulation.run(t_end=1.0) + errors[order] = error() + assert errors[2] < 0.1 * errors[1], errors def test_steps_that_keep_failing_raise(toy, square, sides):