From 2a663861d842772f41843eddcc8a51f43b4acfad Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Wed, 22 May 2024 11:14:55 +0200 Subject: [PATCH 1/3] Removed nested explicit loop in Allen Cahn 2D FD exact solution --- pySDC/implementations/problem_classes/AllenCahn_2D_FD.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py b/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py index 04260f30ec..ca4ec22037 100644 --- a/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py +++ b/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py @@ -286,10 +286,9 @@ def eval_rhs(t, u): me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init) else: - for i in range(self.nvars[0]): - for j in range(self.nvars[1]): - r2 = self.xvalues[i] ** 2 + self.xvalues[j] ** 2 - me[i, j] = np.tanh((self.radius - np.sqrt(r2)) / (np.sqrt(2) * self.eps)) + X, Y = np.meshgrid(self.xvalues, self.xvalues) + r2 = X**2 + Y**2 + me[:] = np.tanh((self.radius - np.sqrt(r2)) / (np.sqrt(2) * self.eps)) return me From 9c967baf9abfcb231754275f68a7babcf9333f8b Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Wed, 22 May 2024 12:31:50 +0200 Subject: [PATCH 2/3] Simplified hook for computing the radius for Allen-Cahn in TOMS project --- pySDC/projects/TOMS/AllenCahn_monitor.py | 51 +++++++++++++----------- 1 file changed, 28 insertions(+), 23 deletions(-) diff --git a/pySDC/projects/TOMS/AllenCahn_monitor.py b/pySDC/projects/TOMS/AllenCahn_monitor.py index 9755340f81..dca4e457a2 100644 --- a/pySDC/projects/TOMS/AllenCahn_monitor.py +++ b/pySDC/projects/TOMS/AllenCahn_monitor.py @@ -4,37 +4,45 @@ class monitor(hooks): + phase_thresh = 0.0 # count everything above this threshold to the high phase. + def __init__(self): """ Initialization of Allen-Cahn monitoring """ - super(monitor, self).__init__() + super().__init__() self.init_radius = None + def get_exact_radius(self, t): + return np.sqrt(max(self.init_radius**2 - 2.0 * t, 0)) + + @classmethod + def get_radius(cls, u, dx): + c = np.count_nonzero(u > cls.phase_thresh) + return np.sqrt(c / np.pi) * dx + + @staticmethod + def get_interface_width(u, L): + # TODO: How does this generalize to different phase transitions? + rows1 = np.where(u[L.prob.init[0][0] // 2, : L.prob.init[0][0] // 2] > -0.99) + rows2 = np.where(u[L.prob.init[0][0] // 2, : L.prob.init[0][0] // 2] < 0.99) + + return (rows2[0][-1] - rows1[0][0]) * L.prob.dx / L.prob.eps + def pre_run(self, step, level_number): """ - Overwrite standard pre run hook + Record radius of the blob, exact radius and interface width. Args: step (pySDC.Step.step): the current step level_number (int): the current level number """ - super(monitor, self).pre_run(step, level_number) + super().pre_run(step, level_number) L = step.levels[0] - c = np.count_nonzero(L.u[0] > 0.0) - radius = np.sqrt(c / np.pi) * L.prob.dx - - radius1 = 0 - rows, cols = np.where(L.u[0] > 0.0) - for r in rows: - radius1 = max(radius1, abs(L.prob.xvalues[r])) - - rows1 = np.where(L.u[0][int((L.prob.init[0][0]) / 2), : int((L.prob.init[0][0]) / 2)] > -0.99) - rows2 = np.where(L.u[0][int((L.prob.init[0][0]) / 2), : int((L.prob.init[0][0]) / 2)] < 0.99) - interface_width = (rows2[0][-1] - rows1[0][0]) * L.prob.dx / L.prob.eps - + radius = self.get_radius(L.u[0], L.prob.dx) + interface_width = self.get_interface_width(L.u[0], L) self.init_radius = L.prob.radius if L.time == 0.0: @@ -68,24 +76,21 @@ def pre_run(self, step, level_number): def post_step(self, step, level_number): """ - Overwrite standard post step hook + Record radius of the blob, exact radius and interface width. Args: step (pySDC.Step.step): the current step level_number (int): the current level number """ - super(monitor, self).post_step(step, level_number) + super().post_step(step, level_number) # some abbreviations L = step.levels[0] - c = np.count_nonzero(L.uend >= 0.0) - radius = np.sqrt(c / np.pi) * L.prob.dx + radius = self.get_radius(L.uend, L.prob.dx) + interface_width = self.get_interface_width(L.uend, L) - exact_radius = np.sqrt(max(self.init_radius**2 - 2.0 * (L.time + L.dt), 0)) - rows1 = np.where(L.uend[int((L.prob.init[0][0]) / 2), : int((L.prob.init[0][0]) / 2)] > -0.99) - rows2 = np.where(L.uend[int((L.prob.init[0][0]) / 2), : int((L.prob.init[0][0]) / 2)] < 0.99) - interface_width = (rows2[0][-1] - rows1[0][0]) * L.prob.dx / L.prob.eps + exact_radius = self.get_exact_radius(L.time + L.dt) self.add_to_stats( process=step.status.slot, From e306b7b7e03a35fc3d5b91d296d84e6242dded6c Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Wed, 22 May 2024 13:41:37 +0200 Subject: [PATCH 3/3] Removed unused legacy function --- .../problem_classes/AllenCahn_2D_FD.py | 35 ------------------- 1 file changed, 35 deletions(-) diff --git a/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py b/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py index ca4ec22037..32275368fc 100644 --- a/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py +++ b/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py @@ -134,41 +134,6 @@ def __init__( self.work_counters['rhs'] = WorkCounter() self.work_counters['linear'] = WorkCounter() - @staticmethod - def __get_A(N, dx): - """ - Helper function to assemble FD matrix A in sparse format. - - Parameters - ---------- - N : list - Number of degrees of freedom. - dx : float - Distance between two spatial nodes. - - Returns - ------- - A : scipy.sparse.csc_matrix - Matrix in CSC format. - """ - - stencil = [1, -2, 1] - zero_pos = 2 - - dstencil = np.concatenate((stencil, np.delete(stencil, zero_pos - 1))) - offsets = np.concatenate( - ( - [N[0] - i - 1 for i in reversed(range(zero_pos - 1))], - [i - zero_pos + 1 for i in range(zero_pos - 1, len(stencil))], - ) - ) - doffsets = np.concatenate((offsets, np.delete(offsets, zero_pos - 1) - N[0])) - - A = sp.diags(dstencil, doffsets, shape=(N[0], N[0]), format='csc') - A = sp.kron(A, sp.eye(N[0])) + sp.kron(sp.eye(N[1]), A) - A *= 1.0 / (dx**2) - return A - # noinspection PyTypeChecker def solve_system(self, rhs, factor, u0, t): """