diff --git a/src/Makefile.am b/src/Makefile.am
index 7e2ceeb6067bb831620472b50ccd51fdfd50b494..4d6a67c7569464486a017d8306c3e54730ebb3b2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -56,7 +56,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
 		 kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
                  units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
-		 equation_of_state.h \
+		 dimension.h equation_of_state.h \
 		 gravity.h gravity_io.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \
diff --git a/src/const.h b/src/const.h
index 498b60b711e88e8054b5e0e73a18475d98277448..d9e4e09d2081361778df329660db3b6ad8f2e204 100644
--- a/src/const.h
+++ b/src/const.h
@@ -36,6 +36,11 @@
 /* Time integration constants. */
 #define const_max_u_change 0.1f
 
+/* Dimensionality of the problem */
+//#define HYDRO_DIMENSION_3D
+#define HYDRO_DIMENSION_2D
+//#define HYDRO_DIMENSION_1D
+
 /* Hydrodynamical adiabatic index. */
 #define HYDRO_GAMMA_5_3
 //#define HYDRO_GAMMA_4_3
diff --git a/src/dimension.h b/src/dimension.h
new file mode 100644
index 0000000000000000000000000000000000000000..92b56ac25a1027891cecaca696186c985e01a48b
--- /dev/null
+++ b/src/dimension.h
@@ -0,0 +1,158 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016   Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DIMENSION_H
+#define SWIFT_DIMENSION_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "inline.h"
+#include "vector.h"
+
+/* First define some constants */
+#if defined(HYDRO_DIMENSION_3D)
+
+#define hydro_dimension 3.f
+#define hydro_dimension_inv 0.3333333333f
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+#define hydro_dimension 2.f
+#define hydro_dimension_inv 0.5f
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+#define hydro_dimension 1.f
+#define hydro_dimension_inv 1.f
+
+#else
+
+#error "A problem dimensionality must be chosen in const.h !"
+
+#endif
+
+/**
+ * @brief Returns the argument to the power given by the dimension
+ *
+ * Computes \f$x^d\f$.
+ */
+__attribute__((always_inline)) INLINE static float pow_dimension(float x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+  return x * x * x;
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+  return x * x;
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+  return x;
+
+#else
+
+  error("The dimension is not defined !");
+  return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Returns the argument to the power given by the dimension plus one
+ *
+ * Computes \f$x^{d+1}\f$.
+ */
+__attribute__((always_inline)) INLINE static float pow_dimension_plus_one(
+    float x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+  const float x2 = x * x;
+  return x2 * x2;
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+  return x * x * x;
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+  return x * x;
+
+#else
+
+  error("The dimension is not defined !");
+  return 0.f;
+
+#endif
+}
+
+#ifdef WITH_VECTORIZATION
+
+__attribute__((always_inline)) INLINE static vector pow_dimension_vec(
+    vector x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+  return (vector)(x.v * x.v * x.v);
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+  return (vector)(x.v * x.v);
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+  return x;
+
+#else
+
+  error("The dimension is not defined !");
+  return vec_set(0.f);
+
+#endif
+}
+
+__attribute__((always_inline)) INLINE static vector pow_dimension_plus_one_vec(
+    vector x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+  const vector x2 = x.v * x.v;
+  return (vector)(x2.v * x2.v);
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+  return (vector)(x.v * x.v * x.v);
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+  return (vector)(x.v * x.v);
+
+#else
+
+  error("The dimension is not defined !");
+  return vec_set(0.f);
+
+#endif
+}
+#endif
+
+#endif /* SWIFT_DIMENSION_H */
diff --git a/src/drift.h b/src/drift.h
index a6e347cc337385f89bece53bc477bff7d163c202..880595dc59e3e5174ee5e888e595a9204ad383e2 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -25,6 +25,7 @@
 /* Local headers. */
 #include "const.h"
 #include "debug.h"
+#include "dimension.h"
 #include "hydro.h"
 #include "part.h"
 
@@ -85,7 +86,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
     p->h *= expf(w1);
 
   /* Predict density */
-  const float w2 = -3.0f * w1;
+  const float w2 = -hydro_dimension * w1;
   if (fabsf(w2) < 0.2f)
     p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
   else
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index a1448d8686a30de88e24629558f01661daca223a..dafec2939eebabcad809da8f75987618fae0c3f0 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -18,6 +18,7 @@
  ******************************************************************************/
 
 #include "adiabatic_index.h"
+#include "dimension.h"
 #include "equation_of_state.h"
 
 /**
@@ -169,33 +170,33 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Some smoothing length multiples. */
   const float h = p->h;
-  const float ih = 1.0f / h;
-  const float ih2 = ih * ih;
-  const float ih4 = ih2 * ih2;
+  const float h_inv = 1.0f / h;                       /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
+  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
 
   /* Final operation on the density (add self-contribution). */
   p->rho += p->mass * kernel_root;
-  p->rho_dh -= 3.0f * p->mass * kernel_root;
+  p->rho_dh -= hydro_dimension * p->mass * kernel_root;
   p->density.wcount += kernel_root;
 
   /* Finish the calculation by inserting the missing h-factors */
-  p->rho *= ih * ih2;
-  p->rho_dh *= ih4;
+  p->rho *= h_inv_dim;
+  p->rho_dh *= h_inv_dim_plus_one;
   p->density.wcount *= kernel_norm;
-  p->density.wcount_dh *= ih * kernel_gamma * kernel_norm;
+  p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
 
   const float irho = 1.f / p->rho;
 
   /* Compute the derivative term */
-  p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
+  p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
 
   /* Finish calculation of the velocity curl components */
-  p->density.rot_v[0] *= ih4 * irho;
-  p->density.rot_v[1] *= ih4 * irho;
-  p->density.rot_v[2] *= ih4 * irho;
+  p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
+  p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
+  p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
 
   /* Finish calculation of the velocity divergence */
-  p->density.div_v *= ih4 * irho;
+  p->density.div_v *= h_inv_dim_plus_one * irho;
 }
 
 /**
@@ -308,7 +309,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part *restrict p) {
 
-  p->force.h_dt *= p->h * 0.333333333f;
+  p->force.h_dt *= p->h * hydro_dimension_inv;
 
   p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
 }
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 323f066db6a047db7ae3401a3525afa941bfc047..d37ac491fe5d8ecdf127c217ca025080daf4bbfd 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -58,7 +58,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   /* Compute contribution to the density */
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * (3.f * wi + ui * wi_dx);
+  pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
 
   /* Compute contribution to the number of neighbours */
   pi->density.wcount += wi;
@@ -71,7 +71,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   /* Compute contribution to the density */
   pj->rho += mi * wj;
-  pj->rho_dh -= mi * (3.f * wj + uj * wj_dx);
+  pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
 
   /* Compute contribution to the number of neighbours */
   pj->density.wcount += wj;
@@ -191,7 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
 
   /* Compute density of pi. */
   rhoi.v = mj.v * wi.v;
-  rhoi_dh.v = mj.v * (vec_set1(3.0f) * wi.v + xi.v * wi_dx.v);
+  rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
   wcounti.v = wi.v;
   wcounti_dh.v = xi.v * wi_dx.v;
   div_vi.v = mj.v * dvdr.v * wi_dx.v;
@@ -199,7 +199,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
 
   /* Compute density of pj. */
   rhoj.v = mi.v * wj.v;
-  rhoj_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v);
+  rhoj_dh.v = mi.v * (vec_set1(hydro_dimension) * wj.v + xj.v * wj_dx.v);
   wcountj.v = wj.v;
   wcountj_dh.v = xj.v * wj_dx.v;
   div_vj.v = mi.v * dvdr.v * wj_dx.v;
@@ -253,7 +253,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
 
   /* Compute contribution to the density */
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * (3.f * wi + u * wi_dx);
+  pi->rho_dh -= mj * (hydro_dimension * wi + u * wi_dx);
 
   /* Compute contribution to the number of neighbours */
   pi->density.wcount += wi;
@@ -356,7 +356,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
 
   /* Compute density of pi. */
   rhoi.v = mj.v * wi.v;
-  rhoi_dh.v = mj.v * (vec_set1(3.0f) * wi.v + xi.v * wi_dx.v);
+  rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
   wcounti.v = wi.v;
   wcounti_dh.v = xi.v * wi_dx.v;
   div_vi.v = mj.v * dvdr.v * wi_dx.v;
@@ -402,17 +402,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Get the kernel for hi. */
   const float hi_inv = 1.0f / hi;
-  const float hi2_inv = hi_inv * hi_inv;
+  const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
   const float ui = r * hi_inv;
   kernel_deval(ui, &wi, &wi_dx);
-  const float wi_dr = hi2_inv * hi2_inv * wi_dx;
+  const float wi_dr = hid_inv * wi_dx;
 
   /* Get the kernel for hj. */
   const float hj_inv = 1.0f / hj;
-  const float hj2_inv = hj_inv * hj_inv;
+  const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
   const float xj = r * hj_inv;
   kernel_deval(xj, &wj, &wj_dx);
-  const float wj_dr = hj2_inv * hj2_inv * wj_dx;
+  const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
   const float P_over_rho2_i = pi->force.P_over_rho2;
@@ -485,7 +485,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
   vector r, r2, ri;
   vector xi, xj;
   vector hi, hj, hi_inv, hj_inv;
-  vector hi2_inv, hj2_inv;
+  vector hid_inv, hjd_inv;
   vector wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
   vector piPOrho, pjPOrho, pirho, pjrho;
   vector mi, mj;
@@ -580,19 +580,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
   hi.v = vec_load(Hi);
   hi_inv.v = vec_rcp(hi.v);
   hi_inv.v = hi_inv.v - hi_inv.v * (hi.v * hi_inv.v - vec_set1(1.0f));
-  hi2_inv.v = hi_inv.v * hi_inv.v;
+  hid_inv = pow_dimension_plus_one_vec(hi_inv); /* 1/h^(d+1) */
   xi.v = r.v * hi_inv.v;
   kernel_deval_vec(&xi, &wi, &wi_dx);
-  wi_dr.v = hi2_inv.v * hi2_inv.v * wi_dx.v;
+  wi_dr.v = hid_inv.v * wi_dx.v;
 
   /* Get the kernel for hj. */
   hj.v = vec_load(Hj);
   hj_inv.v = vec_rcp(hj.v);
   hj_inv.v = hj_inv.v - hj_inv.v * (hj.v * hj_inv.v - vec_set1(1.0f));
-  hj2_inv.v = hj_inv.v * hj_inv.v;
+  hjd_inv = pow_dimension_plus_one_vec(hj_inv); /* 1/h^(d+1) */
   xj.v = r.v * hj_inv.v;
   kernel_deval_vec(&xj, &wj, &wj_dx);
-  wj_dr.v = hj2_inv.v * hj2_inv.v * wj_dx.v;
+  wj_dr.v = hjd_inv.v * wj_dx.v;
 
   /* Compute dv dot r. */
   dvdr.v = ((vi[0].v - vj[0].v) * dx[0].v) + ((vi[1].v - vj[1].v) * dx[1].v) +
@@ -677,17 +677,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Get the kernel for hi. */
   const float hi_inv = 1.0f / hi;
-  const float hi2_inv = hi_inv * hi_inv;
+  const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
   const float ui = r * hi_inv;
   kernel_deval(ui, &wi, &wi_dx);
-  const float wi_dr = hi2_inv * hi2_inv * wi_dx;
+  const float wi_dr = hid_inv * wi_dx;
 
   /* Get the kernel for hj. */
   const float hj_inv = 1.0f / hj;
-  const float hj2_inv = hj_inv * hj_inv;
+  const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
   const float xj = r * hj_inv;
   kernel_deval(xj, &wj, &wj_dx);
-  const float wj_dr = hj2_inv * hj2_inv * wj_dx;
+  const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
   const float P_over_rho2_i = pi->force.P_over_rho2;
@@ -753,7 +753,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
   vector r, r2, ri;
   vector xi, xj;
   vector hi, hj, hi_inv, hj_inv;
-  vector hi2_inv, hj2_inv;
+  vector hid_inv, hjd_inv;
   vector wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
   vector piPOrho, pjPOrho, pirho, pjrho;
   vector mj;
@@ -845,19 +845,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
   hi.v = vec_load(Hi);
   hi_inv.v = vec_rcp(hi.v);
   hi_inv.v = hi_inv.v - hi_inv.v * (hi.v * hi_inv.v - vec_set1(1.0f));
-  hi2_inv.v = hi_inv.v * hi_inv.v;
+  hid_inv = pow_dimension_plus_one_vec(hi_inv);
   xi.v = r.v * hi_inv.v;
   kernel_deval_vec(&xi, &wi, &wi_dx);
-  wi_dr.v = hi2_inv.v * hi2_inv.v * wi_dx.v;
+  wi_dr.v = hid_inv.v * wi_dx.v;
 
   /* Get the kernel for hj. */
   hj.v = vec_load(Hj);
   hj_inv.v = vec_rcp(hj.v);
   hj_inv.v = hj_inv.v - hj_inv.v * (hj.v * hj_inv.v - vec_set1(1.0f));
-  hj2_inv.v = hj_inv.v * hj_inv.v;
+  hjd_inv = pow_dimension_plus_one_vec(hj_inv);
   xj.v = r.v * hj_inv.v;
   kernel_deval_vec(&xj, &wj, &wj_dx);
-  wj_dr.v = hj2_inv.v * hj2_inv.v * wj_dx.v;
+  wj_dr.v = hjd_inv.v * wj_dx.v;
 
   /* Compute dv dot r. */
   dvdr.v = ((vi[0].v - vj[0].v) * dx[0].v) + ((vi[1].v - vj[1].v) * dx[1].v) +
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 8dde4a23be3705b4a679c6d17eb8f28c1a8ff6b7..c425ecce6de38b32645b367886e039bb950df68d 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -27,6 +27,7 @@
 /* Local headers. */
 #include "adiabatic_index.h"
 #include "common_io.h"
+#include "dimension.h"
 #include "error.h"
 #include "hydro.h"
 #include "kernel_hydro.h"
@@ -39,8 +40,7 @@ void hydro_props_init(struct hydro_props *p,
 
   /* Kernel properties */
   p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
-  const float eta3 = p->eta_neighbours * p->eta_neighbours * p->eta_neighbours;
-  p->target_neighbours = eta3 * kernel_norm;
+  p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
   p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
 
   /* Ghost stuff */
@@ -51,21 +51,26 @@ void hydro_props_init(struct hydro_props *p,
   p->CFL_condition = parser_get_param_float(params, "SPH:CFL_condition");
   const float max_volume_change = parser_get_opt_param_float(
       params, "SPH:max_volume_change", hydro_props_default_volume_change);
-  p->log_max_h_change = logf(powf(max_volume_change, 0.33333333333f));
+  p->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
 }
 
 void hydro_props_print(const struct hydro_props *p) {
 
   message("Adiabatic index gamma: %f.", hydro_gamma);
-  message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION);
+
+  message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
+          (int)hydro_dimension);
+
   message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
           kernel_name, p->target_neighbours, p->delta_neighbours,
           p->eta_neighbours);
+
   message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
+
   message(
       "Hydrodynamic integration: Max change of volume: %.2f "
       "(max|dlog(h)/dt|=%f).",
-      powf(expf(p->log_max_h_change), 3.f), p->log_max_h_change);
+      powf(expf(p->log_max_h_change), hydro_dimension), p->log_max_h_change);
 
   if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
     message("Maximal iterations in ghost task set to %d (default is %d)",
@@ -76,6 +81,7 @@ void hydro_props_print(const struct hydro_props *p) {
 void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
 
   writeAttribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+  writeAttribute_i(h_grpsph, "Dimension", (int)hydro_dimension);
   writeAttribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION);
   writeAttribute_s(h_grpsph, "Kernel function", kernel_name);
   writeAttribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
diff --git a/src/kernel_hydro.h b/src/kernel_hydro.h
index 7a73bf4f3d51550c5bf9096c0121d0faab819c0b..be2815cae3fec8e99d2c59aca167d6f50ce3fd97 100644
--- a/src/kernel_hydro.h
+++ b/src/kernel_hydro.h
@@ -24,6 +24,7 @@
 
 /* Includes. */
 #include "const.h"
+#include "dimension.h"
 #include "error.h"
 #include "inline.h"
 #include "vector.h"
@@ -35,8 +36,16 @@
 #define kernel_name "Cubic spline (M4)"
 #define kernel_degree 3 /* Degree of the polynomial */
 #define kernel_ivals 2  /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
 #define kernel_gamma ((float)(1.825742))
 #define kernel_constant ((float)(16. * M_1_PI))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(1.778002))
+#define kernel_constant ((float)(80. * M_1_PI / 7.))
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_gamma ((float)(1.732051))
+#define kernel_constant ((float)(8. / 3.))
+#endif
 static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
     __attribute__((aligned(16))) = {3.f,  -3.f, 0.f,  0.5f, /* 0 < u < 0.5 */
                                     -1.f, 3.f,  -3.f, 1.f,  /* 0.5 < u < 1 */
@@ -47,10 +56,18 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
 
 /* Coefficients for the kernel. */
 #define kernel_name "Quartic spline (M5)"
-#define kernel_degree 4
-#define kernel_ivals 5
+#define kernel_degree 4 /* Degree of the polynomial */
+#define kernel_ivals 5  /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
 #define kernel_gamma ((float)(2.018932))
 #define kernel_constant ((float)(15625. * M_1_PI / 512.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(1.977173))
+#define kernel_constant ((float)(46875. * M_1_PI / 2398.))
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_gamma ((float)(1.936492))
+#define kernel_constant ((float)(3125. / 768.))
+#endif
 static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
     __attribute__((aligned(16))) = {
         6.f,  0.f,  -2.4f, 0.f,   0.368f, /* 0 < u < 0.2 */
@@ -65,10 +82,18 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
 
 /* Coefficients for the kernel. */
 #define kernel_name "Quintic spline (M6)"
-#define kernel_degree 5
-#define kernel_ivals 3
+#define kernel_degree 5 /* Degree of the polynomial */
+#define kernel_ivals 3  /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
 #define kernel_gamma ((float)(2.195775))
 #define kernel_constant ((float)(2187. * M_1_PI / 40.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(2.158131))
+#define kernel_constant ((float)(15309. * M_1_PI / 478.))
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_gamma ((float)(2.121321))
+#define kernel_constant ((float)(243. / 40.))
+#endif
 static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
     __attribute__((aligned(16))) = {
         -10.f,        10.f,      0.f,
@@ -85,10 +110,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
 
 /* Coefficients for the kernel. */
 #define kernel_name "Wendland C2"
-#define kernel_degree 5
-#define kernel_ivals 1
+#define kernel_degree 5 /* Degree of the polynomial */
+#define kernel_ivals 1  /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
 #define kernel_gamma ((float)(1.936492))
 #define kernel_constant ((float)(21. * M_1_PI / 2.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(1.897367))
+#define kernel_constant ((float)(7. * M_1_PI))
+#elif defined(HYDRO_DIMENSION_1D)
+#error "Wendland C2 kernel not defined in 1D."
+#endif
 static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
     __attribute__((aligned(16))) = {
         4.f, -15.f, 20.f, -10.f, 0.f, 1.f,  /* 0 < u < 1 */
@@ -99,10 +131,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
 
 /* Coefficients for the kernel. */
 #define kernel_name "Wendland C4"
-#define kernel_degree 8
+#define kernel_degree 8 /* Degree of the polynomial */
 #define kernel_ivals 1
+#if defined(HYDRO_DIMENSION_3D)
 #define kernel_gamma ((float)(2.207940))
 #define kernel_constant ((float)(495. * M_1_PI / 32.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(2.171239))
+#define kernel_constant ((float)(9. * M_1_PI))
+#elif defined(HYDRO_DIMENSION_1D)
+#error "Wendland C4 kernel not defined in 1D."
+#endif
 static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
     __attribute__((aligned(16))) = {
         11.666667f, -64.f,       140.f, -149.333333f, 70.f,
@@ -115,10 +154,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
 
 /* Coefficients for the kernel. */
 #define kernel_name "Wendland C6"
-#define kernel_degree 11
-#define kernel_ivals 1
+#define kernel_degree 11 /* Degree of the polynomial */
+#define kernel_ivals 1   /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
 #define kernel_gamma ((float)(2.449490))
 #define kernel_constant ((float)(1365. * M_1_PI / 64.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(2.415230))
+#define kernel_constant ((float)(78. * M_1_PI / 7.))
+#elif defined(HYDRO_DIMENSION_1D)
+#error "Wendland C6 kernel not defined in 1D."
+#endif
 static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
     __attribute__((aligned(16))) = {
         32.f, -231.f, 704.f, -1155.f, 1056.f, -462.f,
@@ -151,16 +197,32 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
 #define kernel_ivals_f ((float)(kernel_ivals))
 
 /* Kernel self contribution (i.e. W(0,h)) */
+#if defined(HYDRO_DIMENSION_3D)
 #define kernel_root \
   ((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_igamma3)
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_root \
+  ((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_igamma2)
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_root \
+  ((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_igamma)
+#endif
 
 /* Kernel normalisation constant (volume term) */
+#if defined(HYDRO_DIMENSION_3D)
 #define kernel_norm ((float)(4.0 * M_PI * kernel_gamma3 / 3.0))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_norm ((float)(M_PI * kernel_gamma2))
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_norm ((float)(2.0 * kernel_gamma))
+#endif
+
+/* ------------------------------------------------------------------------- */
 
 /**
  * @brief Computes the kernel function and its derivative.
  *
- * Return 0 if $u > \\gamma = H/h$
+ * Returns garbage if $u > \\gamma = H/h$
  *
  * @param u The ratio of the distance to the smoothing length $u = x/h$.
  * @param W (return) The value of the kernel function $W(x,h)$.
@@ -200,6 +262,8 @@ __attribute__((always_inline)) INLINE static void kernel_deval(
 /**
  * @brief Computes the kernel function.
  *
+ * Returns garbage if $u > \\gamma = H/h$
+ *
  * @param u The ratio of the distance to the smoothing length $u = x/h$.
  * @param W (return) The value of the kernel function $W(x,h)$.
  */
@@ -228,6 +292,8 @@ __attribute__((always_inline)) INLINE static void kernel_eval(
   *W = w * kernel_constant * kernel_igamma3;
 }
 
+/* ------------------------------------------------------------------------- */
+
 #ifdef WITH_VECTORIZATION
 
 static const vector kernel_igamma_vec = FILL_VEC((float)kernel_igamma);