diff --git a/src/struphy/initial/perturbations.py b/src/struphy/initial/perturbations.py
index 4746b8bc3..60de28f11 100644
--- a/src/struphy/initial/perturbations.py
+++ b/src/struphy/initial/perturbations.py
@@ -1264,6 +1264,38 @@ def __call__(self, e1, e2, e3):
         return val
 
 
+class GaussianBlobEta1(Perturbation):
+    r"""Gaussian blob in eta1.
+
+    .. math::
+
+        u(\eta_1, \eta_2, \eta_3) = A \exp \left(- \frac{(\eta_1 - 0.5)^2}{2 \sigma^2} \right) \,.
+    """
+
+    def __init__(
+        self,
+        center: float = 0.5,
+        amp: float = 1e-1,
+        sigma: float = 0.1,
+        given_in_basis: LiteralOptions.GivenInBasis = None,
+        comp: int = 0,
+    ):
+        if given_in_basis is not None:
+            assert "physical" not in given_in_basis, f"Perturbation {self.__name__} can only be used in logical space."
+
+        self._center = center
+        self._amp = amp
+        self._sigma = sigma
+
+        # use the setters
+        self.given_in_basis = given_in_basis
+        self.comp = comp
+
+    def __call__(self, e1, e2, e3):
+        val = self._amp * xp.exp(-((e1 - self._center) ** 2) / (2.0 * self._sigma**2))
+        return val
+
+
 class Erf_z(Perturbation):
     r"""Shear layer in eta3 (-1 in lower regions, 1 in upper regions).
 
diff --git a/src/struphy/models/tests/verification/test_verif_ViscousEulerSPH.py b/src/struphy/models/tests/verification/test_verif_ViscousEulerSPH.py
index 3566e6b4c..cea5ff898 100644
--- a/src/struphy/models/tests/verification/test_verif_ViscousEulerSPH.py
+++ b/src/struphy/models/tests/verification/test_verif_ViscousEulerSPH.py
@@ -2,6 +2,7 @@
 import shutil
 
 import cunumpy as xp
+import h5py
 import pytest
 from feectools.ddm.mpi import mpi as MPI
 from matplotlib import pyplot as plt
@@ -161,5 +162,139 @@ def test_soundwave_1d(nx: int, plot_pts: int, do_plot: bool = False):
         shutil.rmtree(test_folder)
 
 
+@pytest.mark.parametrize("nx", [12, 24])
+@pytest.mark.parametrize("plot_pts", [11, 32])
+def test_viscosity_1d(nx: int, plot_pts: int, do_plot: bool = False):
+    """Verification test for SPH discretization of viscosity in Euler equations.
+    A Gaussian blob in vx diffuses in periodic boundary conditions.
+    """
+
+    # environment options
+    test_folder = os.path.join(os.getcwd(), "struphy_verification_tests")
+    out_folders = os.path.join(test_folder, "ViscousEulerSPH")
+    env = EnvironmentOptions(out_folders=out_folders, sim_folder="viscosity_1d")
+
+    # units
+    base_units = BaseUnits(kBT=1.0)
+
+    # time stepping
+    time_opts = Time(dt=0.01, Tend=1.0, split_algo="LieTrotter")
+
+    # geometry
+    r1 = 1.0
+    domain = domains.Cuboid(r1=r1)
+
+    # grid
+    grid = None
+
+    # derham options
+    derham_opts = None
+
+    # light-weight model instance
+    model = ViscousEulerSPH(with_B0=False, with_p=False, with_viscosity=True)
+
+    # species parameters
+    model.euler_fluid.set_species_properties()
+
+    loading_params = LoadingParameters(ppb=100, loading="tesselation")
+    weights_params = WeightsParameters()
+    boundary_params = BoundaryParameters(bc_sph=("periodic", "periodic", "periodic"))
+    model.euler_fluid.set_markers(
+        loading_params=loading_params,
+        weights_params=weights_params,
+        boundary_params=boundary_params,
+    )
+    model.euler_fluid.set_sorting_boxes(
+        boxes_per_dim=(nx, 1, 1),
+        dims_maks=(True, False, False),
+    )
+
+    bin_plot = BinningPlot(slice="e1", n_bins=(32,), ranges=(0.0, 1.0), output_quantity="current_1")
+    kd_plot = KernelDensityPlot(pts_e1=plot_pts, pts_e2=1)
+    model.euler_fluid.set_save_data(
+        binning_plots=(bin_plot,),
+        kernel_density_plots=(kd_plot,),
+    )
+
+    # propagator options
+    from struphy.ode.utils import ButcherTableau
+
+    butcher = ButcherTableau(algo="forward_euler")
+    # model.propagators.push_eta.options = model.propagators.push_eta.Options(butcher=butcher)
+    if model.with_viscosity:
+        model.propagators.push_viscous.options = model.propagators.push_viscous.Options(
+            kernel_type="gaussian_1d", mu=0.1
+        )
+
+    # background, perturbations and initial conditions
+    background = equils.ConstantVelocity()
+    model.euler_fluid.var.add_background(background)
+    # perturbation = perturbations.GaussianBlobEta1(center=0.5, amp=1.0, sigma=0.1)
+    perturbation = perturbations.ModesSin(ls=(1,), amps=(1.0,))
+    model.euler_fluid.var.add_perturbation(del_u1=perturbation)
+
+    # instance of simulation
+    sim = Simulation(
+        model=model,
+        env=env,
+        base_units=base_units,
+        time_opts=time_opts,
+        domain=domain,
+        grid=grid,
+        derham_opts=derham_opts,
+        verbose=True,
+    )
+
+    # run
+    sim.run(verbose=True)
+
+    # get scalar data
+    if MPI.COMM_WORLD.Get_rank() == 0:
+        pa_data = os.path.join(env.path_out, "data")
+        with h5py.File(os.path.join(pa_data, "data_proc0.hdf5"), "r") as f:
+            time = f["time"]["value"][()]
+            en_kin = f["scalar"]["en_kin"][()]
+
+        # plot
+        if do_plot:
+            plt.figure(figsize=(18, 12))
+            plt.plot(time, en_kin, label="numerical")
+            plt.legend()
+            plt.show()
+
+    # post processing
+    if MPI.COMM_WORLD.Get_rank() == 0:
+        sim.pproc(verbose=True)
+
+        # diagnostics
+        sim.load_plotting_data(env.path_out)
+
+        shp = sim.t_grid.size
+        grid_e1 = sim.f.euler_fluid.e1_current_1.grid_e1
+        f_binned = sim.f.euler_fluid.e1_current_1.f_binned
+        print(f_binned.shape)
+
+        if do_plot:
+            plt.figure(figsize=(20, 8))
+            plt.subplot(1, 4, 1)
+            plt.plot(grid_e1, f_binned[0, :], label=f"time {sim.t_grid[0]}")
+            plt.title(f"time {sim.t_grid[0]}")
+            plt.ylim([-1, 1])
+            plt.subplot(1, 4, 2)
+            plt.plot(grid_e1, f_binned[shp // 3, :], label=f"time {sim.t_grid[shp // 3]}")
+            plt.title(f"time {sim.t_grid[shp // 3]}")
+            plt.ylim([-1, 1])
+            plt.subplot(1, 4, 3)
+            plt.plot(grid_e1, f_binned[2 * shp // 3, :], label=f"time {sim.t_grid[2 * shp // 3]}")
+            plt.title(f"time {sim.t_grid[2 * shp // 3]}")
+            plt.ylim([-1, 1])
+            plt.subplot(1, 4, 4)
+            plt.plot(grid_e1, f_binned[-1, :], label=f"time {sim.t_grid[-1]}")
+            plt.title(f"time {sim.t_grid[-1]}")
+            plt.ylim([-1, 1])
+            plt.show()
+
+
 if __name__ == "__main__":
-    test_soundwave_1d(nx=12, plot_pts=11, do_plot=True)
+    # test_soundwave_1d(nx=12, plot_pts=11, do_plot=True)
+    test_viscosity_1d(nx=12, plot_pts=11, do_plot=True)
diff --git a/src/struphy/models/viscous_euler_sph.py b/src/struphy/models/viscous_euler_sph.py
index 07ef82b99..2cdc59e6c 100644
--- a/src/struphy/models/viscous_euler_sph.py
+++ b/src/struphy/models/viscous_euler_sph.py
@@ -75,7 +75,7 @@ def __init__(self):
 
     class Propagators:
         def __init__(self, with_B0: bool = True, with_p: bool = True, with_viscosity: bool = True):
-            self.push_eta = propagators_markers.PushEta()
+            # self.push_eta = propagators_markers.PushEta()
             if with_B0:
                 self.push_vxb = propagators_markers.PushVxB()
             if with_p:
@@ -98,7 +98,7 @@ def __init__(self, with_B0: bool = True, with_p: bool = True, with_viscosity: bo
         self.propagators = self.Propagators(with_B0=with_B0, with_p=with_p, with_viscosity=with_viscosity)
 
         # 3. assign variables to propagators
-        self.propagators.push_eta.variables.var = self.euler_fluid.var
+        # self.propagators.push_eta.variables.var = self.euler_fluid.var
         if with_B0:
             self.propagators.push_vxb.variables.ions = self.euler_fluid.var
         if with_p:
diff --git a/src/struphy/pic/base.py b/src/struphy/pic/base.py
index 36bb6b18e..77e79004e 100644
--- a/src/struphy/pic/base.py
+++ b/src/struphy/pic/base.py
@@ -1291,8 +1291,10 @@ def _set_initial_condition(self):
         else:
             assert isinstance(self.f0, FluidEquilibrium)
 
-            # get vector-field representation of the fluid velocity
-            self._u_init = self.f0.uv
+            _del_n = None
+            _del_u1 = None
+            _del_u2 = None
+            _del_u3 = None
 
             if self.perturbations is not None:
                 for (
@@ -1308,7 +1310,7 @@ def _set_initial_condition(self):
                         if pert.given_in_basis is None:
                             pert.given_in_basis = "0"
 
-                        _fun = TransformedPformComponent(
+                        _del_n = TransformedPformComponent(
                             pert,
                             pert.given_in_basis,
                             "0",
@@ -1318,24 +1320,24 @@ def _set_initial_condition(self):
                     elif moment == "u1":
                         if pert.given_in_basis is None:
                             pert.given_in_basis = "v"
-                        _fun = TransformedPformComponent(
+                        _del_u1 = TransformedPformComponent(
                             pert,
                             pert.given_in_basis,
                             "v",
                             comp=pert.comp,
                             domain=self.domain,
                         )
-                        self._u_init = lambda e1, e2, e3: self.f0.uv(e1, e2, e3) + _fun(e1, e2, e3)
-                        # TODO: add other velocity components
+                        # self._u_init = lambda e1, e2, e3: self.f0.uv(e1, e2, e3) + _del_u1(e1, e2, e3)
+                        # # TODO: add other velocity components
             else:
-                _fun = None
+                _del_n = None
 
             def _f_init(*etas, flat_eval=False):
                 if len(etas) == 1:
-                    if _fun is None:
+                    if _del_n is None:
                         out = self.f0.n0(etas[0])
                     else:
-                        out = self.f0.n0(etas[0]) + _fun(*etas[0].T)
+                        out = self.f0.n0(etas[0]) + _del_n(*etas[0].T)
                 else:
                     assert len(etas) == 3
                     E1, E2, E3, is_sparse_meshgrid = Domain.prepare_eval_pts(
@@ -1347,10 +1349,42 @@ def _f_init(*etas, flat_eval=False):
 
                     out0 = self.f0.n0(E1, E2, E3)
 
-                    if _fun is None:
+                    if _del_n is None:
                         out = out0
                     else:
-                        out1 = _fun(E1, E2, E3)
+                        out1 = _del_n(E1, E2, E3)
+                        assert out0.shape == out1.shape
+                        out = out0 + out1
+
+                    if flat_eval:
+                        out = xp.squeeze(out)
+
+                return out
+
+            def _u_init(*etas, flat_eval=False):
+                if len(etas) == 1:
+                    out = self.f0.uv(etas[0])
+                    if _del_u1 is not None:
+                        out[0] += _del_u1(*etas[0].T)
+                    if _del_u2 is not None:
+                        out[1] += _del_u2(*etas[0].T)
+                    if _del_u3 is not None:
+                        out[2] += _del_u3(*etas[0].T)
+                else:
+                    assert len(etas) == 3
+                    E1, E2, E3, is_sparse_meshgrid = Domain.prepare_eval_pts(
+                        etas[0],
+                        etas[1],
+                        etas[2],
+                        flat_eval=flat_eval,
+                    )
+
+                    out0 = self.f0.uv(E1, E2, E3)
+
+                    if _del_u1 is None:
+                        out = out0
+                    else:
+                        out1 = _del_u1(E1, E2, E3)
                         assert out0.shape == out1.shape
                         out = out0 + out1
 
@@ -1360,6 +1394,7 @@ def _f_init(*etas, flat_eval=False):
                 return out
 
             self._f_init = _f_init
+            self._u_init = _u_init
 
     def _load_external(
         self,
@@ -1568,6 +1603,7 @@ def draw_markers(
             self._load_tesselation()
             if self.type == "sph":
                 self._set_initial_condition()
+                print()
                 self.velocities = xp.array(self.u_init(self.positions)).T
             # set markers ID in last column
             self.marker_ids = _first_marker_id + xp.arange(n_mks_load_loc, dtype=float)
diff --git a/src/struphy/pic/pushing/pusher_kernels.py b/src/struphy/pic/pushing/pusher_kernels.py
index b770b1151..a32c6e0f3 100644
--- a/src/struphy/pic/pushing/pusher_kernels.py
+++ b/src/struphy/pic/pushing/pusher_kernels.py
@@ -3127,6 +3127,7 @@ def push_v_viscosity(
         # n_at_eta = markers[ip, first_free_idx]
         loc_box = int(markers[ip, n_cols - 2])
 
+        f_visc[:] = 0.0
         for j in range(3):  # row of viscosity tensor
             for k in range(3):  # column = derivative direction
                 coeff_idx = first_free_idx + 3 * (j + 1) + k
diff --git a/src/struphy/pic/tests/test_sph.py b/src/struphy/pic/tests/test_sph.py
index 52c350244..cfcb3f769 100644
--- a/src/struphy/pic/tests/test_sph.py
+++ b/src/struphy/pic/tests/test_sph.py
@@ -1959,11 +1959,10 @@ def u_xyz(x, y, z):
 
 
 if __name__ == "__main__":
-    # test_sph_no_slip_boundary_1d(
-    #     tesselation=False,
-    #     direction="x",
-    #     show_plot=True,
-    # )
-    test_sph_velocity_evaluation_2d(
-        (12, 12, 1), "gaussian_2d", 1, "periodic", "periodic", 11, tesselation=False, show_plot=True
+    test_sph_no_slip_boundary_1d(
+        (1, 1, 12),
+        "gaussian_3d",
+        tesselation=False,
+        direction="z",
+        show_plot=True,
     )