From 38f3e2bc3ad710afc3249afb10420f836ee4c4a0 Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Sat, 27 Apr 2024 19:55:48 +0200 Subject: [PATCH] Fixed some issues with transfer classes (#422) * Fixes for transfer classes * Cleanup of `mesh_to_mesh` * Bug fixes * Fixed `fft_to_fft` --- pySDC/helpers/transfer_helper.py | 2 +- .../check_iteration_estimator.py | 2 +- .../implementations/datatype_classes/mesh.py | 2 +- .../problem_classes/Brusselator.py | 1 + .../problem_classes/GrayScott_MPIFFT.py | 1 + .../sweeper_classes/imex_1st_order.py | 8 +- .../transfer_classes/TransferMesh.py | 132 ++++++--------- .../transfer_classes/TransferMesh_FFT.py | 50 +++--- .../transfer_classes/TransferMesh_MPIFFT.py | 103 ++++++------ .../test_mesh_to_mesh.py} | 1 - .../test_transfer_order.py | 155 ++++++++++++++++++ 11 files changed, 282 insertions(+), 175 deletions(-) rename pySDC/tests/{test_spatial_transfer.py => test_transfer_classes/test_mesh_to_mesh.py} (99%) create mode 100644 pySDC/tests/test_transfer_classes/test_transfer_order.py diff --git a/pySDC/helpers/transfer_helper.py b/pySDC/helpers/transfer_helper.py index 4f5cd441f2..196d50f0a5 100644 --- a/pySDC/helpers/transfer_helper.py +++ b/pySDC/helpers/transfer_helper.py @@ -138,7 +138,7 @@ def restriction_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, pad=1): def interpolation_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, pad=1, equidist_nested=True): """ - Function to contruct the restriction matrix in 1d using barycentric interpolation + Function to construct the restriction matrix in 1d using barycentric interpolation Args: fine_grid (np.ndarray): a one dimensional 1d array containing the nodes of the fine grid diff --git a/pySDC/implementations/convergence_controller_classes/check_iteration_estimator.py b/pySDC/implementations/convergence_controller_classes/check_iteration_estimator.py index 06ee05f142..d2245ba99f 100644 --- a/pySDC/implementations/convergence_controller_classes/check_iteration_estimator.py +++ b/pySDC/implementations/convergence_controller_classes/check_iteration_estimator.py @@ -6,7 +6,7 @@ class CheckIterationEstimatorNonMPI(ConvergenceController): def __init__(self, controller, params, description, **kwargs): """ - Initalization routine + Initialization routine Args: controller (pySDC.Controller): The controller diff --git a/pySDC/implementations/datatype_classes/mesh.py b/pySDC/implementations/datatype_classes/mesh.py index ce99ebc78a..41d77b1540 100644 --- a/pySDC/implementations/datatype_classes/mesh.py +++ b/pySDC/implementations/datatype_classes/mesh.py @@ -199,7 +199,7 @@ def __new__(cls, init, *args, **kwargs): def __getattr__(self, name): if name in self.components: if self.shape[0] == len(self.components): - return self[self.components.index(name)] + return self[self.components.index(name)].view(mesh) else: raise AttributeError(f'Cannot access {name!r} in {type(self)!r} because the shape is unexpected.') else: diff --git a/pySDC/implementations/problem_classes/Brusselator.py b/pySDC/implementations/problem_classes/Brusselator.py index 946e676781..746747d4c8 100644 --- a/pySDC/implementations/problem_classes/Brusselator.py +++ b/pySDC/implementations/problem_classes/Brusselator.py @@ -30,6 +30,7 @@ def __init__(self, alpha=0.1, **kwargs): shape = (2,) + (self.init[0]) self.iU = 0 self.iV = 1 + self.ncomp = 2 # needed for transfer class self.init = (shape, self.comm, np.dtype('float')) def _eval_explicit_part(self, u, t, f_expl): diff --git a/pySDC/implementations/problem_classes/GrayScott_MPIFFT.py b/pySDC/implementations/problem_classes/GrayScott_MPIFFT.py index 293b1c846e..ef597addae 100644 --- a/pySDC/implementations/problem_classes/GrayScott_MPIFFT.py +++ b/pySDC/implementations/problem_classes/GrayScott_MPIFFT.py @@ -79,6 +79,7 @@ def __init__(self, Du=1.0, Dv=0.01, A=0.09, B=0.086, **kwargs): shape = (2,) + (self.init[0]) self.iU = 0 self.iV = 1 + self.ncomp = 2 # needed for transfer class self.init = (shape, self.comm, self.xp.dtype('float')) self._makeAttributeAndRegister('Du', 'Dv', 'A', 'B', localVars=locals(), readOnly=True) diff --git a/pySDC/implementations/sweeper_classes/imex_1st_order.py b/pySDC/implementations/sweeper_classes/imex_1st_order.py index bc9805119d..c86b4facab 100644 --- a/pySDC/implementations/sweeper_classes/imex_1st_order.py +++ b/pySDC/implementations/sweeper_classes/imex_1st_order.py @@ -42,16 +42,14 @@ def integrate(self): list of dtype_u: containing the integral as values """ - # get current level and problem description L = self.level + P = L.prob me = [] - # integrate RHS over all collocation nodes for m in range(1, self.coll.num_nodes + 1): - me.append(L.dt * self.coll.Qmat[m, 1] * (L.f[1].impl + L.f[1].expl)) - # new instance of dtype_u, initialize values with 0 - for j in range(2, self.coll.num_nodes + 1): + me.append(P.dtype_u(P.init, val=0.0)) + for j in range(1, self.coll.num_nodes + 1): me[m - 1] += L.dt * self.coll.Qmat[m, j] * (L.f[j].impl + L.f[j].expl) return me diff --git a/pySDC/implementations/transfer_classes/TransferMesh.py b/pySDC/implementations/transfer_classes/TransferMesh.py index b6a0c8dee8..6f89c70e7f 100644 --- a/pySDC/implementations/transfer_classes/TransferMesh.py +++ b/pySDC/implementations/transfer_classes/TransferMesh.py @@ -4,15 +4,14 @@ import pySDC.helpers.transfer_helper as th from pySDC.core.Errors import TransferError from pySDC.core.SpaceTransfer import space_transfer -from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh, comp2_mesh class mesh_to_mesh(space_transfer): """ - Custon base_transfer class, implements Transfer.py + Custom base_transfer class, implements Transfer.py This implementation can restrict and prolong between nd meshes with dirichlet-0 or periodic boundaries - via matrix-vector products + via matrix-vector products. Attributes: Rspace: spatial restriction matrix, dim. Nf x Nc @@ -30,7 +29,7 @@ def __init__(self, fine_prob, coarse_prob, params): """ # invoke super initialization - super(mesh_to_mesh, self).__init__(fine_prob, coarse_prob, params) + super().__init__(fine_prob, coarse_prob, params) if self.params.rorder % 2 != 0: raise TransferError('Need even order for restriction') @@ -153,51 +152,31 @@ def restrict(self, F): Args: F: the fine level data (easier to access than via the fine attribute) """ - if isinstance(F, mesh): - G = self.coarse_prob.dtype_u(self.coarse_prob.init) - if hasattr(self.fine_prob, 'ncomp'): - for i in range(self.fine_prob.ncomp): - tmpF = F[..., i].flatten() - tmpG = self.Rspace.dot(tmpF) - G[..., i] = tmpG.reshape(self.coarse_prob.nvars) - else: - tmpF = F.flatten() - tmpG = self.Rspace.dot(tmpF) - G[:] = tmpG.reshape(self.coarse_prob.nvars) - elif isinstance(F, imex_mesh): - G = self.coarse_prob.dtype_f(self.coarse_prob.init) - if hasattr(self.fine_prob, 'ncomp'): - for i in range(self.fine_prob.ncomp): - tmpF = F.impl[..., i].flatten() - tmpG = self.Rspace.dot(tmpF) - G.impl[..., i] = tmpG.reshape(self.coarse_prob.nvars) - tmpF = F.expl[..., i].flatten() - tmpG = self.Rspace.dot(tmpF) - G.expl[..., i] = tmpG.reshape(self.coarse_prob.nvars) - else: - tmpF = F.impl.flatten() - tmpG = self.Rspace.dot(tmpF) - G.impl[:] = tmpG.reshape(self.coarse_prob.nvars) - tmpF = F.expl.flatten() - tmpG = self.Rspace.dot(tmpF) - G.expl[:] = tmpG.reshape(self.coarse_prob.nvars) - elif isinstance(F, comp2_mesh): - G = self.coarse_prob.dtype_f(self.coarse_prob.init) + G = type(F)(self.coarse_prob.init) + + def _restrict(fine, coarse): if hasattr(self.fine_prob, 'ncomp'): for i in range(self.fine_prob.ncomp): - tmpF = F.comp1[..., i].flatten() - tmpG = self.Rspace.dot(tmpF) - G.comp1[..., i] = tmpG.reshape(self.coarse_prob.nvars) - tmpF = F.comp2[..., i].flatten() - tmpG = self.Rspace.dot(tmpF) - G.comp2[..., i] = tmpG.reshape(self.coarse_prob.nvars) + if fine.shape[-1] == self.fine_prob.ncomp: + tmpF = fine[..., i].flatten() + tmpG = self.Rspace.dot(tmpF) + coarse[..., i] = tmpG.reshape(self.coarse_prob.nvars) + elif fine.shape[0] == self.fine_prob.ncomp: + tmpF = fine[i, ...].flatten() + tmpG = self.Rspace.dot(tmpF) + coarse[i, ...] = tmpG.reshape(self.coarse_prob.nvars) + else: + raise TransferError('Don\'t know how to restrict for this problem with multiple components') else: - tmpF = F.comp1.flatten() - tmpG = self.Rspace.dot(tmpF) - G.comp1[:] = tmpG.reshape(self.coarse_prob.nvars) - tmpF = F.comp2.flatten() + tmpF = fine.flatten() tmpG = self.Rspace.dot(tmpF) - G.comp2[:] = tmpG.reshape(self.coarse_prob.nvars) + coarse[:] = tmpG.reshape(self.coarse_prob.nvars) + + if hasattr(type(F), 'components'): + for comp in F.components: + _restrict(F.__getattr__(comp), G.__getattr__(comp)) + elif type(F).__name__ == 'mesh': + _restrict(F, G) else: raise TransferError('Wrong data type for restriction, got %s' % type(F)) return G @@ -208,51 +187,32 @@ def prolong(self, G): Args: G: the coarse level data (easier to access than via the coarse attribute) """ - if isinstance(G, mesh): - F = self.fine_prob.dtype_u(self.fine_prob.init) - if hasattr(self.fine_prob, 'ncomp'): - for i in range(self.fine_prob.ncomp): - tmpG = G[..., i].flatten() - tmpF = self.Pspace.dot(tmpG) - F[..., i] = tmpF.reshape(self.fine_prob.nvars) - else: - tmpG = G.flatten() - tmpF = self.Pspace.dot(tmpG) - F[:] = tmpF.reshape(self.fine_prob.nvars) - elif isinstance(G, imex_mesh): - F = self.fine_prob.dtype_f(self.fine_prob.init) - if hasattr(self.fine_prob, 'ncomp'): - for i in range(self.fine_prob.ncomp): - tmpG = G.impl[..., i].flatten() - tmpF = self.Pspace.dot(tmpG) - F.impl[..., i] = tmpF.reshape(self.fine_prob.nvars) - tmpG = G.expl[..., i].flatten() - tmpF = self.Rspace.dot(tmpG) - F.expl[..., i] = tmpF.reshape(self.fine_prob.nvars) - else: - tmpG = G.impl.flatten() - tmpF = self.Pspace.dot(tmpG) - F.impl[:] = tmpF.reshape(self.fine_prob.nvars) - tmpG = G.expl.flatten() - tmpF = self.Pspace.dot(tmpG) - F.expl[:] = tmpF.reshape(self.fine_prob.nvars) - elif isinstance(G, comp2_mesh): - F = self.fine_prob.dtype_f(self.fine_prob.init) + F = type(G)(self.fine_prob.init) + + def _prolong(coarse, fine): if hasattr(self.fine_prob, 'ncomp'): for i in range(self.fine_prob.ncomp): - tmpG = G.comp1[..., i].flatten() - tmpF = self.Pspace.dot(tmpG) - F.comp1[..., i] = tmpF.reshape(self.fine_prob.nvars) - tmpG = G.comp2[..., i].flatten() - tmpF = self.Rspace.dot(tmpG) - F.comp2[..., i] = tmpF.reshape(self.fine_prob.nvars) + if coarse.shape[-1] == self.fine_prob.ncomp: + tmpG = coarse[..., i].flatten() + tmpF = self.Pspace.dot(tmpG) + fine[..., i] = tmpF.reshape(self.fine_prob.nvars) + elif coarse.shape[0] == self.fine_prob.ncomp: + tmpG = coarse[i, ...].flatten() + tmpF = self.Pspace.dot(tmpG) + fine[i, ...] = tmpF.reshape(self.fine_prob.nvars) + else: + raise TransferError('Don\'t know how to prolong for this problem with multiple components') else: - tmpG = G.comp1.flatten() - tmpF = self.Pspace.dot(tmpG) - F.comp1[:] = tmpF.reshape(self.fine_prob.nvars) - tmpG = G.comp2.flatten() + tmpG = coarse.flatten() tmpF = self.Pspace.dot(tmpG) - F.comp2[:] = tmpF.reshape(self.fine_prob.nvars) + fine[:] = tmpF.reshape(self.fine_prob.nvars) + return fine + + if hasattr(type(F), 'components'): + for comp in G.components: + _prolong(G.__getattr__(comp), F.__getattr__(comp)) + elif type(G).__name__ == 'mesh': + F[:] = _prolong(G, F) else: raise TransferError('Wrong data type for prolongation, got %s' % type(G)) return F diff --git a/pySDC/implementations/transfer_classes/TransferMesh_FFT.py b/pySDC/implementations/transfer_classes/TransferMesh_FFT.py index 18277a3c8c..2664a063a6 100644 --- a/pySDC/implementations/transfer_classes/TransferMesh_FFT.py +++ b/pySDC/implementations/transfer_classes/TransferMesh_FFT.py @@ -2,7 +2,6 @@ from pySDC.core.Errors import TransferError from pySDC.core.SpaceTransfer import space_transfer -from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh class mesh_to_mesh_fft(space_transfer): @@ -26,9 +25,9 @@ def __init__(self, fine_prob, coarse_prob, params): params: parameters for the transfer operators """ # invoke super initialization - super(mesh_to_mesh_fft, self).__init__(fine_prob, coarse_prob, params) + super().__init__(fine_prob, coarse_prob, params) - self.ratio = int(self.fine_prob.params.nvars / self.coarse_prob.params.nvars) + self.ratio = int(self.fine_prob.nvars / self.coarse_prob.nvars) def restrict(self, F): """ @@ -37,11 +36,11 @@ def restrict(self, F): Args: F: the fine level data (easier to access than via the fine attribute) """ - if isinstance(F, mesh): - G = mesh(self.coarse_prob.init, val=0.0) + G = type(F)(self.coarse_prob.init, val=0.0) + + if type(F).__name__ == 'mesh': G[:] = F[:: self.ratio] - elif isinstance(F, imex_mesh): - G = imex_mesh(self.coarse_prob.init, val=0.0) + elif type(F).__name__ == 'imex_mesh': G.impl[:] = F.impl[:: self.ratio] G.expl[:] = F.expl[:: self.ratio] else: @@ -55,28 +54,21 @@ def prolong(self, G): Args: G: the coarse level data (easier to access than via the coarse attribute) """ - if isinstance(G, mesh): - F = mesh(self.fine_prob.init, val=0.0) - tmpG = np.fft.rfft(G) - tmpF = np.zeros(self.fine_prob.init[0] // 2 + 1, dtype=np.complex128) - halfG = int(self.coarse_prob.init[0] / 2) - tmpF[0:halfG] = tmpG[0:halfG] - tmpF[-1] = tmpG[-1] - F[:] = np.fft.irfft(tmpF) * self.ratio - elif isinstance(G, imex_mesh): - F = imex_mesh(G) - tmpG_impl = np.fft.rfft(G.impl) - tmpF_impl = np.zeros(self.fine_prob.init[0] // 2 + 1, dtype=np.complex128) - halfG = int(self.coarse_prob.init[0] / 2) - tmpF_impl[0:halfG] = tmpG_impl[0:halfG] - tmpF_impl[-1] = tmpG_impl[-1] - F.impl[:] = np.fft.irfft(tmpF_impl) * self.ratio - tmpG_expl = np.fft.rfft(G.expl) - tmpF_expl = np.zeros(self.fine_prob.init[0] // 2 + 1, dtype=np.complex128) - halfG = int(self.coarse_prob.init[0] / 2) - tmpF_expl[0:halfG] = tmpG_expl[0:halfG] - tmpF_expl[-1] = tmpG_expl[-1] - F.expl[:] = np.fft.irfft(tmpF_expl) * self.ratio + F = type(G)(self.fine_prob.init, val=0.0) + + def _prolong(coarse): + coarse_hat = np.fft.rfft(coarse) + fine_hat = np.zeros(self.fine_prob.init[0] // 2 + 1, dtype=np.complex128) + half_idx = self.coarse_prob.init[0] // 2 + fine_hat[0:half_idx] = coarse_hat[0:half_idx] + fine_hat[-1] = coarse_hat[-1] + return np.fft.irfft(fine_hat) * self.ratio + + if type(G).__name__ == 'mesh': + F[:] = _prolong(G) + elif type(G).__name__ == 'imex_mesh': + F.impl[:] = _prolong(G.impl) + F.expl[:] = _prolong(G.expl) else: raise TransferError('Unknown data type, got %s' % type(G)) return F diff --git a/pySDC/implementations/transfer_classes/TransferMesh_MPIFFT.py b/pySDC/implementations/transfer_classes/TransferMesh_MPIFFT.py index 1e55ccbd49..35072e1a43 100644 --- a/pySDC/implementations/transfer_classes/TransferMesh_MPIFFT.py +++ b/pySDC/implementations/transfer_classes/TransferMesh_MPIFFT.py @@ -6,7 +6,7 @@ class fft_to_fft(space_transfer): """ - Custon base_transfer class, implements Transfer.py + Custom base_transfer class, implements Transfer.py This implementation can restrict and prolong between PMESH datatypes meshes with FFT for periodic boundaries @@ -22,7 +22,7 @@ def __init__(self, fine_prob, coarse_prob, params): params: parameters for the transfer operators """ # invoke super initialization - super(fft_to_fft, self).__init__(fine_prob, coarse_prob, params) + super().__init__(fine_prob, coarse_prob, params) assert self.fine_prob.spectral == self.coarse_prob.spectral @@ -49,24 +49,38 @@ def restrict(self, F): Args: F: the fine level data (easier to access than via the fine attribute) """ - if isinstance(F, mesh): + G = type(F)(self.coarse_prob.init) + + def _restrict(fine, coarse): if self.spectral: - G = self.coarse_prob.dtype_u(self.coarse_prob.init) if hasattr(self.fine_prob, 'ncomp'): for i in range(self.fine_prob.ncomp): - tmpF = newDistArray(self.fine_prob.fft, False) - tmpF = self.fine_prob.fft.backward(F[..., i], tmpF) - tmpG = tmpF[:: int(self.ratio[0]), :: int(self.ratio[1])] - G[..., i] = self.coarse_prob.fft.forward(tmpG, G[..., i]) + if fine.shape[-1] == self.fine_prob.ncomp: + tmpF = newDistArray(self.fine_prob.fft, False) + tmpF = self.fine_prob.fft.backward(fine[..., i], tmpF) + tmpG = tmpF[:: int(self.ratio[0]), :: int(self.ratio[1])] + coarse[..., i] = self.coarse_prob.fft.forward(tmpG, coarse[..., i]) + elif fine.shape[0] == self.fine_prob.ncomp: + tmpF = newDistArray(self.fine_prob.fft, False) + tmpF = self.fine_prob.fft.backward(fine[i, ...], tmpF) + tmpG = tmpF[:: int(self.ratio[0]), :: int(self.ratio[1])] + coarse[i, ...] = self.coarse_prob.fft.forward(tmpG, coarse[i, ...]) + else: + raise TransferError('Don\'t know how to restrict for this problem with multiple components') else: - tmpF = self.fine_prob.fft.backward(F) + tmpF = self.fine_prob.fft.backward(fine) tmpG = tmpF[:: int(self.ratio[0]), :: int(self.ratio[1])] - G[:] = self.coarse_prob.fft.forward(tmpG, G) + coarse[:] = self.coarse_prob.fft.forward(tmpG, coarse) else: - G = self.coarse_prob.dtype_u(self.coarse_prob.init) - G[:] = F[:: int(self.ratio[0]), :: int(self.ratio[1])] + coarse[:] = fine[:: int(self.ratio[0]), :: int(self.ratio[1])] + + if hasattr(type(F), 'components'): + for comp in F.components: + _restrict(F.__getattr__(comp), G.__getattr__(comp)) + elif type(F).__name__ == 'mesh': + _restrict(F, G) else: - raise TransferError('Unknown data type, got %s' % type(F)) + raise TransferError('Wrong data type for restriction, got %s' % type(F)) return G @@ -77,52 +91,39 @@ def prolong(self, G): Args: G: the coarse level data (easier to access than via the coarse attribute) """ - if isinstance(G, mesh): - if self.spectral: - F = self.fine_prob.dtype_u(self.fine_prob.init) - if hasattr(self.fine_prob, 'ncomp'): - for i in range(self.fine_prob.ncomp): - tmpF = self.fft_pad.backward(G[..., i]) - F[..., i] = self.fine_prob.fft.forward(tmpF, F[..., i]) - else: - tmpF = self.fft_pad.backward(G) - F[:] = self.fine_prob.fft.forward(tmpF, F) - else: - F = self.fine_prob.dtype_u(self.fine_prob.init) - if hasattr(self.fine_prob, 'ncomp'): - for i in range(self.fine_prob.ncomp): - G_hat = self.coarse_prob.fft.forward(G[..., i]) - F[..., i] = self.fft_pad.backward(G_hat, F[..., i]) - else: - G_hat = self.coarse_prob.fft.forward(G) - F[:] = self.fft_pad.backward(G_hat, F) - elif isinstance(G, imex_mesh): + F = type(G)(self.fine_prob.init) + + def _prolong(coarse, fine): if self.spectral: - F = self.fine_prob.dtype_f(self.fine_prob.init) if hasattr(self.fine_prob, 'ncomp'): for i in range(self.fine_prob.ncomp): - tmpF = self.fft_pad.backward(G.impl[..., i]) - F.impl[..., i] = self.fine_prob.fft.forward(tmpF, F.impl[..., i]) - tmpF = self.fft_pad.backward(G.expl[..., i]) - F.expl[..., i] = self.fine_prob.fft.forward(tmpF, F.expl[..., i]) + if coarse.shape[-1] == self.fine_prob.ncomp: + tmpF = self.fft_pad.backward(coarse[..., i]) + fine[..., i] = self.fine_prob.fft.forward(tmpF, fine[..., i]) + elif coarse.shape[0] == self.fine_prob.ncomp: + tmpF = self.fft_pad.backward(coarse[i, ...]) + fine[i, ...] = self.fine_prob.fft.forward(tmpF, fine[i, ...]) + else: + raise TransferError('Don\'t know how to prolong for this problem with multiple components') + else: - tmpF = self.fft_pad.backward(G.impl) - F.impl[:] = self.fine_prob.fft.forward(tmpF, F.impl) - tmpF = self.fft_pad.backward(G.expl) - F.expl[:] = self.fine_prob.fft.forward(tmpF, F.expl) + tmpF = self.fft_pad.backward(coarse) + fine[:] = self.fine_prob.fft.forward(tmpF, fine) else: - F = self.fine_prob.dtype_f(self.fine_prob.init) if hasattr(self.fine_prob, 'ncomp'): for i in range(self.fine_prob.ncomp): - G_hat = self.coarse_prob.fft.forward(G.impl[..., i]) - F.impl[..., i] = self.fft_pad.backward(G_hat, F.impl[..., i]) - G_hat = self.coarse_prob.fft.forward(G.expl[..., i]) - F.expl[..., i] = self.fft_pad.backward(G_hat, F.expl[..., i]) + G_hat = self.coarse_prob.fft.forward(coarse[..., i]) + fine[..., i] = self.fft_pad.backward(G_hat, fine[..., i]) else: - G_hat = self.coarse_prob.fft.forward(G.impl) - F.impl[:] = self.fft_pad.backward(G_hat, F.impl) - G_hat = self.coarse_prob.fft.forward(G.expl) - F.expl[:] = self.fft_pad.backward(G_hat, F.expl) + G_hat = self.coarse_prob.fft.forward(coarse) + fine[:] = self.fft_pad.backward(G_hat, fine) + + if hasattr(type(F), 'components'): + for comp in F.components: + _prolong(G.__getattr__(comp), F.__getattr__(comp)) + elif type(G).__name__ == 'mesh': + _prolong(G, F) + else: raise TransferError('Unknown data type, got %s' % type(G)) diff --git a/pySDC/tests/test_spatial_transfer.py b/pySDC/tests/test_transfer_classes/test_mesh_to_mesh.py similarity index 99% rename from pySDC/tests/test_spatial_transfer.py rename to pySDC/tests/test_transfer_classes/test_mesh_to_mesh.py index aecc3d11bf..ecb76ba69c 100644 --- a/pySDC/tests/test_spatial_transfer.py +++ b/pySDC/tests/test_transfer_classes/test_mesh_to_mesh.py @@ -250,4 +250,3 @@ def test_mesh_to_mesh_2d_periodic(): if __name__ == "__main__": test_mesh_to_mesh_1d_dirichlet() - pass diff --git a/pySDC/tests/test_transfer_classes/test_transfer_order.py b/pySDC/tests/test_transfer_classes/test_transfer_order.py new file mode 100644 index 0000000000..f6ac427015 --- /dev/null +++ b/pySDC/tests/test_transfer_classes/test_transfer_order.py @@ -0,0 +1,155 @@ +import pytest + + +def get_problem(nvars, xp, L, mpifft=False, spectral=False, x0=0): + """ + Get a dummy problem to test interpolation + + Args: + nvars (int): Number of DoF + xp: Numerical library, i.e. Numpy or CuPy + L (float): Length of space domain + + Returns: + Instance of pySDC problem class + """ + from pySDC.core.Problem import ptype + from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh + from pySDC.helpers.problem_helper import get_1d_grid + + class DummyProblem(ptype): + dtype_u = mesh + dtype_f = imex_mesh + + def __init__(self, nvars): + super().__init__(init=(nvars, None, xp.dtype('float64'))) + self._makeAttributeAndRegister('nvars', localVars=locals()) + self.dx = L / (nvars + 1) + self.x = xp.array([x0 + me * self.dx for me in range(nvars)]) + self.dx, self.x = get_1d_grid(nvars, 'periodic', right_boundary=L) + + def eval_f(self, *args, **kwargs): + me = self.f_init + me.impl[:] = xp.sin(2 * xp.pi / L * self.x) + me.expl[:] = xp.cos(2 * xp.pi / L * self.x) + return me + + def u_exact(self, *args, **kwargs): + me = self.u_init + me[:] = xp.sin(2 * xp.pi / L * self.x) + return me + + if mpifft: + from pySDC.implementations.problem_classes.generic_MPIFFT_Laplacian import IMEX_Laplacian_MPIFFT + + class DummyProblemMPIFFT(IMEX_Laplacian_MPIFFT): + def _fill(self): + me = xp.sin(2 * xp.pi / L * self.X[0]) + if self.spectral: + return self.fft.forward(me) + else: + return me + + def eval_f(self, *args, **kwargs): + me = self.f_init + me.impl[:] = self._fill() + me.expl[:] = self._fill() + return me + + def u_exact(self, *args, **kwargs): + me = self.u_init + me[:] = self._fill() + return me + + return DummyProblemMPIFFT(nvars=(nvars, nvars), spectral=spectral, x0=x0, L=L) + + return DummyProblem(nvars) + + +def get_transfer_class(name): + """ + Imports the problem class of the respective name and returns it. + """ + if name == 'mesh_to_mesh_fft': + from pySDC.implementations.transfer_classes.TransferMesh_FFT import mesh_to_mesh_fft as transfer_class + elif name == 'mesh_to_mesh': + from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh as transfer_class + elif name == 'fft_to_fft': + from pySDC.implementations.transfer_classes.TransferMesh_MPIFFT import fft_to_fft as transfer_class + else: + raise NotImplementedError(f'Don\'t know where to import transfer class {name!r}') + return transfer_class + + +@pytest.mark.base +@pytest.mark.parametrize('L', [1.0, 6.283185307179586]) +@pytest.mark.parametrize('x0', [0.0, -3.141592653589793]) +def test_mesh_to_mesh_fft(L, x0): + single_test('mesh_to_mesh_fft', -1.0, L, mpifft=False, x0=x0) + + +@pytest.mark.base +@pytest.mark.parametrize('order', [2, 4, 6, 8]) +@pytest.mark.parametrize('L', [1.0]) +@pytest.mark.parametrize('x0', [0.0, -3.141592653589793]) +def test_mesh_to_mesh(order, L, x0): + single_test('mesh_to_mesh', order, L, mpifft=False, x0=x0) + + +@pytest.mark.mpi4py +@pytest.mark.parametrize('L', [1.0, 6.283185307179586]) +@pytest.mark.parametrize('spectral', [True, False]) +@pytest.mark.parametrize('x0', [0.0, -3.141592653589793]) +def test_fft_to_fft(L, spectral, x0): + single_test('fft_to_fft', -1, L, mpifft=True, spectral=spectral, x0=x0) + + +def single_test(name, order, L, mpifft, spectral=False, x0=0): + import numpy as xp + + params = { + 'iorder': order, + 'rorder': order, + 'periodic': True, + } + resolutions = [2**8, 2**7, 2**6, 2**5] + + errors = { + 'restriction_u_exact': [], + 'restriction_eval_f': [], + 'prolongation_u_exact': [], + 'prolongation_eval_f': [], + } + + for function_name in ['u_exact', 'eval_f']: + + for res in resolutions[::-1]: + fine = get_problem(res, xp, L, mpifft=mpifft, spectral=spectral, x0=x0) + coarse = get_problem(res // 2, xp, L, mpifft=mpifft, spectral=spectral, x0=x0) + transfer = get_transfer_class(name)(fine, coarse, params) + + fine_exact = fine.__getattribute__(function_name)(t=0) + coarse_exact = coarse.__getattribute__(function_name)(t=0) + + errors[f'prolongation_{function_name}'] += [abs(fine_exact - transfer.prolong(coarse_exact))] + errors[f'restriction_{function_name}'] += [abs(coarse_exact - transfer.restrict(fine_exact))] + + resolutions = xp.array(resolutions) + for key, value in errors.items(): + value = xp.array(value) + + mask = value > 1e-14 + + if mask.any(): + num_order = abs( + xp.log(value[mask][1:] / value[mask][:-1]) / xp.log(resolutions[mask][1:] / resolutions[mask][:-1]) + ) + assert xp.isclose( + xp.median(num_order), order, atol=0.3 + ), f'Got unexpected order {xp.median(num_order):.2f} when expecting {order} in {key}. Errors: {value}' + + +if __name__ == '__main__': + import numpy as np + + test_fft_to_fft('fft_to_fft', 2, 2 * np.pi, True, 0.0)