Coverage for pySDC/implementations/hooks/AllenCahn_monitor.py: 100%
60 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
3from pySDC.core.hooks import Hooks
6class AllenCahnMonitor(Hooks):
7 r"""
8 Track the shrinking circle (or sphere) of an Allen-Cahn run.
10 Under mean curvature flow a blob of initial radius :math:`R_0` obeys
11 :math:`R(t)^2 = R_0^2 - 2 (d - 1) t`, so comparing the measured radius against that is the
12 standard diagnostic for these problems. This hook records ``computed_radius``,
13 ``exact_radius``, ``computed_volume`` and ``exact_volume`` at :math:`t = 0` and after every
14 step, plus ``interface_width`` where that is well defined.
16 The volume is the number of cells in the high phase times the cell volume. Counting is used
17 rather than integrating the field, even though pySDC's phase fields run from :math:`0` to
18 :math:`1` and so integrate straight to a volume: the integral also picks up the diffuse
19 interface, which biases it by :math:`O(\varepsilon)` no matter how fine the mesh, whereas
20 counting is consistent and its :math:`O(\Delta x)` bias refines away.
22 Attributes
23 ----------
24 phase_thresh : float
25 Cells above this count towards the high phase. Taken from the problem when it says, and
26 :math:`0.5` -- the midpoint of the two wells -- otherwise.
27 """
29 default_phase_thresh = 0.5
31 def __init__(self):
32 """Start with no initial radius or dimension; both are filled in by :meth:`pre_run`."""
33 super().__init__()
35 self.init_radius = None
36 self.ndim = None
37 self.phase_thresh = self.default_phase_thresh
39 @staticmethod
40 def get_real_space(L, u):
41 """Undo the transform if the problem carries its solution in spectral space."""
42 return L.prob.fft.backward(u) if getattr(L.prob, 'spectral', False) else u[:]
44 def get_volume(self, L, u):
45 """Volume of the high phase, summed over the space communicator if there is one."""
46 count = float(np.count_nonzero(self.get_real_space(L, u) > self.phase_thresh))
48 comm = getattr(L.prob, 'comm', None)
49 if comm is not None:
50 from mpi4py import MPI
52 count = comm.allreduce(sendobj=count, op=MPI.SUM)
54 return count * L.prob.dx**self.ndim
56 def radius_from_volume(self, vol):
57 """Radius of the ball of this volume."""
58 if self.ndim == 2:
59 return np.sqrt(vol / np.pi)
60 elif self.ndim == 3:
61 return (vol / (np.pi * 4.0 / 3.0)) ** (1.0 / 3.0)
62 raise NotImplementedError(f'Can only monitor 2D and 3D problems, got {self.ndim}D')
64 def exact_radius_squared(self, t):
65 r"""Mean curvature flow shrinks the blob as :math:`R(t)^2 = R_0^2 - 2 (d - 1) t`."""
66 return max(self.init_radius**2 - 2.0 * (self.ndim - 1) * t, 0)
68 def exact_radius(self, t):
69 """Exact radius at time ``t`` under mean curvature flow, clipped at zero once the blob has vanished."""
70 return np.sqrt(self.exact_radius_squared(t))
72 def exact_volume(self, t):
73 """Exact area (2D) or volume (3D) of the blob at time ``t``."""
74 r2 = self.exact_radius_squared(t)
75 return np.pi * r2 if self.ndim == 2 else np.pi * 4.0 / 3.0 * r2**1.5
77 def measures_interface_width(self, L):
78 """Only a 2D field held whole on this rank can be cut across to measure the interface."""
79 comm = getattr(L.prob, 'comm', None)
80 return self.ndim == 2 and (comm is None or comm.Get_size() == 1)
82 def get_interface_width(self, L, u):
83 """Width of the transition, in units of epsilon, along a cut through the middle."""
84 n = L.prob.init[0][0]
85 rows1 = np.where(u[n // 2, : n // 2] > 0.005)
86 rows2 = np.where(u[n // 2, : n // 2] < 0.995)
88 return (rows2[0][-1] - rows1[0][0]) * L.prob.dx / L.prob.eps
90 def get_diagnostics(self, L, u, t):
91 """Everything worth recording about ``u``, as a dict of stats entries."""
92 vol = self.get_volume(L, u)
94 diagnostics = {
95 'computed_radius': self.radius_from_volume(vol),
96 'exact_radius': self.exact_radius(t),
97 'computed_volume': vol,
98 'exact_volume': self.exact_volume(t),
99 }
101 if self.measures_interface_width(L):
102 diagnostics['interface_width'] = self.get_interface_width(L, self.get_real_space(L, u))
104 return diagnostics
106 def record(self, step, L, t, diagnostics):
107 """
108 Add each diagnostic to the stats, with its key as ``type``, on level ``-1``.
110 Args:
111 step (pySDC.Step.step): the current step
112 L (pySDC.Level.level): the finest level
113 t (float): the time the diagnostics refer to
114 diagnostics (dict): stats entries as returned by :meth:`get_diagnostics`
115 """
116 for key, value in diagnostics.items():
117 self.add_to_stats(
118 process=step.status.slot,
119 time=t,
120 level=-1,
121 iter=step.status.iter,
122 sweep=L.status.sweep,
123 type=key,
124 value=value,
125 )
127 def pre_run(self, step, level_number):
128 """
129 Read the initial radius, phase threshold and dimension off the problem, and record the diagnostics of the
130 initial condition if the run starts at :math:`t = 0`.
132 Args:
133 step (pySDC.Step.step): the current step
134 level_number (int): the current level number
135 """
136 super().pre_run(step, level_number)
137 L = step.levels[0]
139 self.init_radius = L.prob.radius
140 self.phase_thresh = getattr(L.prob, 'phase_thresh', self.default_phase_thresh)
141 self.ndim = len(self.get_real_space(L, L.u[0]).shape)
143 if L.time == 0.0:
144 self.record(step, L, L.time, self.get_diagnostics(L, L.u[0], 0.0))
146 def post_step(self, step, level_number):
147 """
148 Record the diagnostics of the solution at the end of the step, at time ``L.time + L.dt``.
150 Args:
151 step (pySDC.Step.step): the current step
152 level_number (int): the current level number
153 """
154 super().post_step(step, level_number)
155 L = step.levels[0]
157 self.record(step, L, L.time + L.dt, self.get_diagnostics(L, L.uend, L.time + L.dt))