diff --git a/evoxels/boundary_conditions.py b/evoxels/boundary_conditions.py index ed4bf8a..3f05151 100644 --- a/evoxels/boundary_conditions.py +++ b/evoxels/boundary_conditions.py @@ -30,14 +30,32 @@ def pad_zero_flux_periodic(self, field): padded = self.vg.set(padded, (__,-1,__,__), padded[:,-2,:,:]) return padded - def pad_zero_flux(self, field): - padded = self.vg.pad_zeros(field) - padded = self.vg.set(padded, (__, 0,__,__), padded[:, 1,:,:]) - padded = self.vg.set(padded, (__,-1,__,__), padded[:,-2,:,:]) - padded = self.vg.set(padded, (__,__, 0,__), padded[:,:, 1,:]) - padded = self.vg.set(padded, (__,__,-1,__), padded[:,:,-2,:]) - padded = self.vg.set(padded, (__,__,__, 0), padded[:,:,:, 1]) - padded = self.vg.set(padded, (__,__,__,-1), padded[:,:,:,-2]) + def pad_bc(self, field, bc): + padded = self.vg.pad_periodic(field) + axis_slices = [ + ((__, 0, __, __), (__, 1, __, __), (__, -1, __, __), (__, -2, __, __)), + ((__, __, 0, __), (__, __, 1, __), (__, __, -1, __), (__, __, -2, __)), + ((__, __, __, 0), (__, __, __, 1), (__, __, __, -1), (__, __, __, -2)), + ] + + for axis_bc, (left_ghost, left_inner, right_ghost, right_inner) in zip(bc, axis_slices): + kind, values = axis_bc + + if kind == 'periodic': + continue + + if kind == 'neumann': + padded = self.vg.set(padded, left_ghost, padded[left_inner]) + padded = self.vg.set(padded, right_ghost, padded[right_inner]) + continue + + if kind == 'dirichlet': + padded = self.vg.set(padded, left_ghost, 2.0 * values[0] - padded[left_inner]) + padded = self.vg.set(padded, right_ghost, 2.0 * values[1] - padded[right_inner]) + continue + + raise ValueError(f"Unsupported BC type: {kind}") + return padded def pad_fft_periodic(self, field): @@ -101,16 +119,16 @@ def pad_zero_flux_periodic(self, field): padded = self.vg.set(padded, (__,-1,__,__), fac1*padded[:,-2,:,:] - fac2*padded[:,-3,:,:]) return padded - def pad_zero_flux(self, field): + def pad_bc(self, field, bc): raise NotImplementedError - + def pad_fft_periodic(self, field): """ If field is fully periodic it should be in cell center convention! """ raise NotImplementedError - + def pad_fft_dirichlet_periodic(self, field): """Pad with inverse of flipped field in x direction.""" bc = self.vg.lib.zeros_like(field[:,0:1]) @@ -135,4 +153,4 @@ def trim_ghost_nodes(self, field): else: raise ValueError( f"The provided field has the wrong shape {self.vg.shape}." - ) \ No newline at end of file + ) diff --git a/evoxels/problem_definition.py b/evoxels/problem_definition.py index 373adb3..5fb6a27 100644 --- a/evoxels/problem_definition.py +++ b/evoxels/problem_definition.py @@ -111,12 +111,8 @@ def initialize_boundary_conditions(self): ) elif self.bc_type == ('neumann','periodic','periodic'): self._pad_bc = self.vg.bc.pad_zero_flux_periodic - elif self.bc_type == ('neumann','neumann','neumann'): - self._pad_bc = self.vg.bc.pad_zero_flux else: - raise NotImplementedError( - f"Finite-difference BC combination {self.bc} is not implemented." - ) + self._pad_bc = lambda field: self.vg.bc.pad_bc(field, self.bc) @property def bc_type(self): @@ -140,11 +136,22 @@ def fourier_symbol(self): pass def verify_fft_bc_config(self): - x_bc, y_bc, z_bc = self.bc_type - if y_bc != 'periodic' or z_bc != 'periodic': + x_bc, _, _ = self.bc_type + nonperiodic_axes = tuple( + axis for axis, kind in zip(('x', 'y', 'z'), self.bc_type) + if kind != 'periodic' + ) + + if len(nonperiodic_axes) > 1: raise ValueError( - "FFT-based timesteppers require periodic boundary conditions " - f"in y and z, got {self.bc_type}." + "FFT-based timesteppers currently support at most one non-periodic axis, " + f"got {self.bc_type}." + ) + + if len(nonperiodic_axes) == 1 and nonperiodic_axes[0] != 'x': + raise NotImplementedError( + "FFT-based timesteppers currently only implement the single non-periodic axis " + f"case for x; got {self.bc_type}. Axis permutation is not implemented yet." ) if x_bc == 'periodic': diff --git a/tests/test_solvers.py b/tests/test_solvers.py index f7259e1..878020a 100644 --- a/tests/test_solvers.py +++ b/tests/test_solvers.py @@ -85,11 +85,15 @@ def test_1D_analytical_tanh_profile(): def test_reaction_diffusion_normalizes_bc(): vf = evo.VoxelFields((4, 4, 4)) vg = VoxelGridTorch(vf.grid_info(), device="cpu") - problem = ReactionDiffusion( - vg, - D=1.0, - bc=(('dirichlet', (1, -1)), 'periodic', 'periodic'), - ) + with pytest.warns( + UserWarning, + match="Applying Dirichlet BCs on a cell_center grid reduces the spatial order of convergence to 0.5!", + ): + problem = ReactionDiffusion( + vg, + D=1.0, + bc=(('dirichlet', (1, -1)), 'periodic', 'periodic'), + ) assert problem.bc == ( ("dirichlet", (1, -1)), @@ -98,10 +102,70 @@ def test_reaction_diffusion_normalizes_bc(): ) +def test_reaction_diffusion_mixed_bc_uses_generic_padding_fallback(): + vf = evo.VoxelFields((2, 2, 2)) + vg = VoxelGridTorch(vf.grid_info(), device="cpu") + with pytest.warns( + UserWarning, + match="Applying Dirichlet BCs on a cell_center grid reduces the spatial order of convergence to 0.5!", + ): + problem = ReactionDiffusion( + vg, + D=1.0, + bc=(('dirichlet', (10.0, 20.0)), 'neumann', 'periodic'), + ) + field = vg.init_scalar_field(np.arange(1, 9, dtype=np.float32).reshape(2, 2, 2)) + + padded = vg.to_numpy(problem.pad_bc(field))[0] + expected = np.pad(np.arange(1, 9, dtype=np.float32).reshape(2, 2, 2), 1, mode='wrap') + expected[0, :, :] = 2.0 * 10.0 - expected[1, :, :] + expected[-1, :, :] = 2.0 * 20.0 - expected[-2, :, :] + expected[:, 0, :] = expected[:, 1, :] + expected[:, -1, :] = expected[:, -2, :] + + assert np.allclose(padded, expected) + + +def test_reaction_diffusion_dirichlet_periodic_keeps_specialized_padding(): + vf = evo.VoxelFields((2, 2, 2)) + vg = VoxelGridTorch(vf.grid_info(), device="cpu") + with pytest.warns( + UserWarning, + match="Applying Dirichlet BCs on a cell_center grid reduces the spatial order of convergence to 0.5!", + ): + problem = ReactionDiffusion( + vg, + D=1.0, + bc=(('dirichlet', (1.0, -1.0)), 'periodic', 'periodic'), + ) + field = vg.init_scalar_field(np.arange(1, 9, dtype=np.float32).reshape(2, 2, 2)) + + padded = problem.pad_bc(field) + expected = vg.bc.pad_dirichlet_periodic(field, 1.0, -1.0) + + assert np.allclose(vg.to_numpy(padded), vg.to_numpy(expected)) + + +def test_reaction_diffusion_neumann_periodic_keeps_specialized_padding(): + vf = evo.VoxelFields((2, 2, 2)) + vg = VoxelGridTorch(vf.grid_info(), device="cpu") + problem = ReactionDiffusion( + vg, + D=1.0, + bc=('neumann', 'periodic', 'periodic'), + ) + field = vg.init_scalar_field(np.arange(1, 9, dtype=np.float32).reshape(2, 2, 2)) + + padded = problem.pad_bc(field) + expected = vg.bc.pad_zero_flux_periodic(field) + + assert np.allclose(vg.to_numpy(padded), vg.to_numpy(expected)) + + def test_exponential_euler_rejects_full_neumann_semilinear_problem(): vf = evo.VoxelFields((4, 4, 4)) vg = VoxelGridTorch(vf.grid_info(), device="cpu") problem = TwoPhaseAllenCahn(vg) - with pytest.raises(ValueError, match="periodic boundary conditions in y and z"): + with pytest.raises(ValueError, match="support at most one non-periodic axis"): ExponentialEuler(problem, 0.1)