BC handling

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}."
-            )
+            )

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':

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)