Coverage for pySDC/implementations/hooks/AllenCahn_monitor.py: 100%
60 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +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 super().__init__()
34 self.init_radius = None
35 self.ndim = None
36 self.phase_thresh = self.default_phase_thresh
38 @staticmethod
39 def get_real_space(L, u):
40 """Undo the transform if the problem carries its solution in spectral space."""
41 return L.prob.fft.backward(u) if getattr(L.prob, 'spectral', False) else u[:]
43 def get_volume(self, L, u):
44 """Volume of the high phase, summed over the space communicator if there is one."""
45 count = float(np.count_nonzero(self.get_real_space(L, u) > self.phase_thresh))
47 comm = getattr(L.prob, 'comm', None)
48 if comm is not None:
49 from mpi4py import MPI
51 count = comm.allreduce(sendobj=count, op=MPI.SUM)
53 return count * L.prob.dx**self.ndim
55 def radius_from_volume(self, vol):
56 """Radius of the ball of this volume."""
57 if self.ndim == 2:
58 return np.sqrt(vol / np.pi)
59 elif self.ndim == 3:
60 return (vol / (np.pi * 4.0 / 3.0)) ** (1.0 / 3.0)
61 raise NotImplementedError(f'Can only monitor 2D and 3D problems, got {self.ndim}D')
63 def exact_radius_squared(self, t):
64 r"""Mean curvature flow shrinks the blob as :math:`R(t)^2 = R_0^2 - 2 (d - 1) t`."""
65 return max(self.init_radius**2 - 2.0 * (self.ndim - 1) * t, 0)
67 def exact_radius(self, t):
68 return np.sqrt(self.exact_radius_squared(t))
70 def exact_volume(self, t):
71 r2 = self.exact_radius_squared(t)
72 return np.pi * r2 if self.ndim == 2 else np.pi * 4.0 / 3.0 * r2**1.5
74 def measures_interface_width(self, L):
75 """Only a 2D field held whole on this rank can be cut across to measure the interface."""
76 comm = getattr(L.prob, 'comm', None)
77 return self.ndim == 2 and (comm is None or comm.Get_size() == 1)
79 def get_interface_width(self, L, u):
80 """Width of the transition, in units of epsilon, along a cut through the middle."""
81 n = L.prob.init[0][0]
82 rows1 = np.where(u[n // 2, : n // 2] > 0.005)
83 rows2 = np.where(u[n // 2, : n // 2] < 0.995)
85 return (rows2[0][-1] - rows1[0][0]) * L.prob.dx / L.prob.eps
87 def get_diagnostics(self, L, u, t):
88 """Everything worth recording about ``u``, as a dict of stats entries."""
89 vol = self.get_volume(L, u)
91 diagnostics = {
92 'computed_radius': self.radius_from_volume(vol),
93 'exact_radius': self.exact_radius(t),
94 'computed_volume': vol,
95 'exact_volume': self.exact_volume(t),
96 }
98 if self.measures_interface_width(L):
99 diagnostics['interface_width'] = self.get_interface_width(L, self.get_real_space(L, u))
101 return diagnostics
103 def record(self, step, L, t, diagnostics):
104 for key, value in diagnostics.items():
105 self.add_to_stats(
106 process=step.status.slot,
107 time=t,
108 level=-1,
109 iter=step.status.iter,
110 sweep=L.status.sweep,
111 type=key,
112 value=value,
113 )
115 def pre_run(self, step, level_number):
116 super().pre_run(step, level_number)
117 L = step.levels[0]
119 self.init_radius = L.prob.radius
120 self.phase_thresh = getattr(L.prob, 'phase_thresh', self.default_phase_thresh)
121 self.ndim = len(self.get_real_space(L, L.u[0]).shape)
123 if L.time == 0.0:
124 self.record(step, L, L.time, self.get_diagnostics(L, L.u[0], 0.0))
126 def post_step(self, step, level_number):
127 super().post_step(step, level_number)
128 L = step.levels[0]
130 self.record(step, L, L.time + L.dt, self.get_diagnostics(L, L.uend, L.time + L.dt))