diff --git a/src/rightHandSide/CMakeLists.txt b/src/rightHandSide/CMakeLists.txt
index e6a4334..7056179 100644
--- a/src/rightHandSide/CMakeLists.txt
+++ b/src/rightHandSide/CMakeLists.txt
@@ -1,6 +1,7 @@
 target_sources(astra
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/advection.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/burgers.hpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/boussinesq.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/hydro.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/mhd.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/rightHandSide.hpp
diff --git a/src/rightHandSide/advection.hpp b/src/rightHandSide/advection.hpp
index 0ffa458..1788b9d 100644
--- a/src/rightHandSide/advection.hpp
+++ b/src/rightHandSide/advection.hpp
@@ -11,6 +11,7 @@
 #ifndef RIGHTHANDSIDE_ADVECTION_HPP_
 #define RIGHTHANDSIDE_ADVECTION_HPP_
 
+#include <memory>
 #include <string>
 #include <vector>
 #include "rightHandSide.hpp"
@@ -21,7 +22,7 @@ class Grid;
 template <typename Shear>
 class Advection : public RightHandSide<Array3D<complex>, Shear> {
  public:
-  Advection(Input &input, Grid *grid);
+  Advection(Input &input, Grid *grid, std::vector<std::unique_ptr<RightHandSideConcept<Array3D<complex>>>> &rhsVector);
 
   ~Advection();
 
@@ -45,8 +46,11 @@ class Advection : public RightHandSide<Array3D<complex>, Shear> {
 #include "global.hpp"
 #include "fft.hpp"
 
+using RhsPtr = std::unique_ptr<RightHandSideConcept<Array3D<complex>>>;
+
 template <typename Shear>
-Advection<Shear>::Advection(Input &input, Grid *grid) : RightHandSide<Array3D<complex>, Shear>(input, grid) {
+Advection<Shear>::Advection(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector) :
+                            RightHandSide<Array3D<complex>, Shear>(input, grid, rhsVector) {
   direction = input.Get<int>("Physics","direction",0);
   velocity = input.GetOrSet<real>("Physics","velocity", 0,1.0);
   if(direction < 0 || direction > 2) {
diff --git a/src/rightHandSide/boussinesq.hpp b/src/rightHandSide/boussinesq.hpp
new file mode 100644
index 0000000..824ac1e
--- /dev/null
+++ b/src/rightHandSide/boussinesq.hpp
@@ -0,0 +1,262 @@
+// ***********************************************************************************
+// ASTRA spectral code
+// Accelerated Spectral code for TuRbulent plasmA
+// Copyright(C) Geoffroy R. J. Lesur <geoffroy.lesur@univ-grenoble-alpes.fr>
+// and other code contributors
+// Licensed under CeCILL 2.1 License, see COPYING for more information
+// ***********************************************************************************
+
+#ifndef RIGHTHANDSIDE_BOUSSINESQ_HPP_
+#define RIGHTHANDSIDE_BOUSSINESQ_HPP_
+
+#include <memory>
+#include <vector>
+#include <string>
+#include "rightHandSide.hpp"
+#include "input.hpp"
+#include "arrays.hpp"
+#include "shear.hpp"
+#include "hydro.hpp"
+#include "mhd.hpp"
+
+using RhsPtr = std::unique_ptr<RightHandSideConcept<Array3D<complex>>>;
+
+class Grid;
+
+// A class for the boussinesq right hand side
+template <typename Shear>
+class Boussinesq : public RightHandSide<Array3D<complex>, Shear> {
+ public:
+  Boussinesq(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector);
+  ~Boussinesq();
+
+  void ExplicitStep(Field<Array3D<complex>>& fldin, Field<Array3D<complex>>& dfld, real t) override;
+  void ImplicitStep(Field<Array3D<complex>>& fldin, real t, real dt) override;
+  void PostStage(Field<Array3D<complex>>& fldin, real t) override;
+
+  real GetInvDt() override;
+  std::vector<std::string> GetVariables() override;
+
+ private:
+  real kappa;
+  real g1, g2, g3;
+  real N2;
+  int kappaOrder{1};
+
+  Array3D<real> vr1, vr2, vr3;
+  Array3D<real> thr;
+  Array3D<real> wr1, wr2, wr3;
+  Array3D<complex> wf1, wf2, wf3;
+  std::array<int,3> npr, npf;
+};
+
+
+// Implementation
+#include "grid.hpp"
+#include "loop.hpp"
+#include "reduce.hpp"
+#include "global.hpp"
+#include "fft.hpp"
+
+template <typename Shear>
+Boussinesq<Shear>::Boussinesq(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector) :
+                            RightHandSide<Array3D<complex>, Shear>(input, grid, rhsVector) {
+  // Allocate all of the temporary arrays
+  wr1 = astra::makeArray<Array3D<real>>("Boussinesq::wr1", grid->npr_t);
+  wf1 = astra::makeArray<Array3D<complex>>("Boussinesq::wf1", grid->npf);
+  wr2 = astra::makeArray<Array3D<real>>("Boussinesq::wr2", grid->npr_t);
+  wf2 = astra::makeArray<Array3D<complex>>("Boussinesq::wf2", grid->npf);
+  wr3 = astra::makeArray<Array3D<real>>("Boussinesq::wr3", grid->npr_t);
+  wf3 = astra::makeArray<Array3D<complex>>("Boussinesq::wf3", grid->npf);
+  thr = astra::makeArray<Array3D<real>>("Boussinesq::thr", grid->npr_t);
+
+  // Set the velocity field to be used by the boussinesq class
+  bool found = false;
+
+  // walk through the rhs vector and find the hydro or mhd class, and use their velocity field as the velocity field for the boussinesq class
+  bool foundHydro = false;
+  for(auto &rhs : rhsVector) {
+    if(dynamic_cast<Hydro<Shear>*>(rhs.get())) {
+        auto hydro = dynamic_cast<Hydro<Shear>*>(rhs.get());
+        this->vr1 = hydro->vr1;
+        this->vr2 = hydro->vr2;
+        this->vr3 = hydro->vr3;
+        foundHydro = true;
+        break;
+    } else if(dynamic_cast<Mhd<Shear>*>(rhs.get())) {
+        auto mhd = dynamic_cast<Mhd<Shear>*>(rhs.get());
+        this->vr1 = mhd->vr1;
+        this->vr2 = mhd->vr2;
+        this->vr3 = mhd->vr3;
+        foundHydro = true;
+        break;
+    }
+  }
+
+  if (!foundHydro) {
+    throw std::runtime_error("Boussinesq right hand side requires Hydro or Mhd right hand side to be present in the rhs");
+  }
+
+  npr=grid->npr_t;
+  npf=grid->npf;
+
+  this->kappa = input.GetOrSet<real>("Physics","kappa",0,1e-3);
+  this->kappaOrder = input.GetOrSet<int>("Physics","kappa",1,1);
+  this->g1 = input.GetOrSet<real>("Physics","gstrat",0,0.0);
+  this->g2 = input.GetOrSet<real>("Physics","gstrat",1,0.0);
+  this->g3 = input.GetOrSet<real>("Physics","gstrat",2,1.0);
+  // normalize the gravity vector
+  real gnorm = std::sqrt(g1*g1+g2*g2+g3*g3);
+  if(gnorm==0.0) {
+    throw std::runtime_error("Boussinesq right hand side requires a non-zero gravity vector");
+  }
+  this->g1 /= gnorm;
+  this->g2 /= gnorm;
+  this->g3 /= gnorm;
+  this->N2 = input.Get<real>("Physics","N2",0);
+
+  astra::cout << "Boussinesq: diffusivity kappa=" << kappa << std::endl;
+  astra::cout << "Boussinesq: gravity g=[" << g1 << "," << g2 << "," << g3 << "]" << std::endl;
+  astra::cout << "Boussinesq: N2=" << N2 << std::endl;
+}
+
+template <typename Shear>
+Boussinesq<Shear>::~Boussinesq() {
+}
+
+template <typename Shear>
+void Boussinesq<Shear>::ExplicitStep(Field<Array3D<complex>>& fldin, Field<Array3D<complex>>& dfld, real t) {
+  astra::pushRegion("Boussinesq::ExplicitStep");
+  // Fourier transform the velocity field to real space (use transposed arrays)
+  auto th = fldin["th"];
+  // Normally vr has already been computed by the Hydro or Mhd class
+  auto vr1 = this->vr1;
+  auto vr2 = this->vr2;
+  auto vr3 = this->vr3;
+
+  this->grid->fft->C2R(th, thr, false);
+  auto wr1 = this->wr1;
+  auto wr2 = this->wr2;
+  auto wr3 = this->wr3;
+
+  auto thr = this->thr;
+
+  astra_for("Boussinesq::HydroTransport",0,npr[IDIR],0,npr[JDIR],0,npr[KDIR],
+    KOKKOS_LAMBDA(int64_t i, int64_t j, int64_t k) {
+      wr1(i,j,k) = thr(i,j,k)*vr1(i,j,k);
+      wr2(i,j,k) = thr(i,j,k)*vr2(i,j,k);
+      wr3(i,j,k) = thr(i,j,k)*vr3(i,j,k);
+  });
+
+  // Fourier transform back to spectral space
+  this->grid->fft->R2C(wr1, wf1, false);
+  this->grid->fft->R2C(wr2, wf2, false);
+  this->grid->fft->R2C(wr3, wf3, false);
+
+  // Compute the cross-correlation
+  auto wf1 = this->wf1;
+  auto wf2 = this->wf2;
+  auto wf3 = this->wf3;
+
+  auto kx1 = this->grid->kx[IDIR];
+  auto kx2 = this->grid->kx[JDIR];
+  auto kx3 = this->grid->kx[KDIR];
+
+  auto dvx1 = dfld["vx1"];
+  auto dvx2 = dfld["vx2"];
+  auto dvx3 = dfld["vx3"];
+  auto dth = dfld["th"];
+  auto v1 = fldin["vx1"];
+  auto v2 = fldin["vx2"];
+  auto v3 = fldin["vx3"];
+
+  real kx1max = this->grid->kmax[IDIR];
+  real kx2max = this->grid->kmax[JDIR];
+  real kx3max = this->grid->kmax[KDIR];
+
+  real g1 = this->g1;
+  real g2 = this->g2;
+  real g3 = this->g3;
+  real N2 = this->N2;
+  Shear &shear = this->shear;
+  shear.Refresh(t);
+
+  // Compute the nonlinear term in spectral space
+  astra_for("hydro_nonlinear", 0,npf[IDIR],0,npf[JDIR],0,npf[KDIR],
+    KOKKOS_LAMBDA(int64_t i, int64_t j, int64_t k) {
+      // 2/3 de-aliasing rule
+      complex mask = (std::fabs(kx1(i))< 2./3*kx1max
+                   && std::fabs(kx2(j))< 2./3*kx2max
+                   && std::fabs(kx3(k))< 2./3*kx3max) ? 1.0 : 0.0;
+      const real kx1t = shear.kx1t(kx1(i),kx2(j),kx3(k));
+      const real kx2t = shear.kx2t(kx1(i),kx2(j),kx3(k));
+      const real kx3t = shear.kx3t(kx1(i),kx2(j),kx3(k));
+
+      // dv = -N^2 [g-k(k.g)/k^2] th
+      const real kdotg = kx1t*g1+kx2t*g2+kx3t*g3;
+      real k2 = kx1t*kx1t+kx2t*kx2t+kx3t*kx3t;
+      if(k2==0.0) k2=1.0; // avoid division by zero
+      dvx1(i,j,k) -= N2*(g1-kdotg*kx1t/k2)*th(i,j,k);
+      dvx2(i,j,k) -= N2*(g2-kdotg*kx2t/k2)*th(i,j,k);
+      dvx3(i,j,k) -= N2*(g3-kdotg*kx3t/k2)*th(i,j,k);
+
+      // dth = -nabla.(v th)+g.v
+      dth(i,j,k) = -mask * Kokkos::complex(0.0,1.0) * (kx1t*wf1(i,j,k) + kx2t*wf2(i,j,k) + kx3t*wf3(i,j,k));
+      dth(i,j,k) += g1*v1(i,j,k) + g2*v2(i,j,k) + g3*v3(i,j,k);
+  });
+
+  astra::popRegion();
+}
+
+template <typename Shear>
+void Boussinesq<Shear>::ImplicitStep(Field<Array3D<complex>>& fldin, real t, real dt) {
+  astra::pushRegion("Boussinesq::ImplicitStep");
+  auto kx1 = this->grid->kx[IDIR];
+  auto kx2 = this->grid->kx[JDIR];
+  auto kx3 = this->grid->kx[KDIR];
+
+  auto th = fldin["th"];
+
+  real kappa = this->kappa;
+  int n = this->kappaOrder;
+  Shear shear = this->shear;
+  shear.Refresh(t);
+  astra_for("boussinesq_diffusivity", fldin,
+    KOKKOS_LAMBDA(int64_t i, int64_t j, int64_t k) {
+      const real kx1t = shear.kx1t(kx1(i),kx2(j),kx3(k));
+      const real kx2t = shear.kx2t(kx1(i),kx2(j),kx3(k));
+      const real kx3t = shear.kx3t(kx1(i),kx2(j),kx3(k));
+      const real k2t = kx1t*kx1t + kx2t*kx2t + kx3t*kx3t;
+
+      //real factor = (1-0.5*dt*nu*k2)/(1+0.5*dt*nu*k2); // Crank-Nicholson
+      real factor = std::exp(-dt * pow(kappa*k2t, n)); // Exact integration
+      th(i,j,k) *= factor;
+    });
+  astra::popRegion();
+}
+
+template <typename Shear> real Boussinesq<Shear>::GetInvDt() {
+  astra::pushRegion("Boussinesq::GetInvDt");
+  real invdt = std::sqrt(std::fabs(this->N2));
+
+  astra::popRegion();
+  return invdt;
+}
+
+template <typename Shear>
+std::vector<std::string> Boussinesq<Shear>::GetVariables() {
+  return {"th"};
+}
+
+template <typename Shear>
+void Boussinesq<Shear>::PostStage(Field<Array3D<complex>>& fldin, real t) {
+    if constexpr(Shear::isEnabled) {
+    if(this->shear.NeedRemap(t)) {
+      this->shear.Remap(t, fldin["th"]);
+
+      this->shear.SetTinit(t);
+    }
+  }
+}
+
+#endif // RIGHTHANDSIDE_BOUSSINESQ_HPP_
diff --git a/src/rightHandSide/burgers.hpp b/src/rightHandSide/burgers.hpp
index 42270a6..6e9f163 100644
--- a/src/rightHandSide/burgers.hpp
+++ b/src/rightHandSide/burgers.hpp
@@ -9,6 +9,7 @@
 #ifndef RIGHTHANDSIDE_BURGERS_HPP_
 #define RIGHTHANDSIDE_BURGERS_HPP_
 
+#include <memory>
 #include <vector>
 #include <string>
 #include "rightHandSide.hpp"
@@ -16,13 +17,15 @@
 #include "arrays.hpp"
 #include "shear.hpp"
 
+using RhsPtr = std::unique_ptr<RightHandSideConcept<Array3D<complex>>>;
+
 class Grid;
 
 // A class for the burgers right hand side
 template <typename Shear>
 class Burgers : public RightHandSide<Array3D<complex>, Shear> {
  public:
-  Burgers(Input &input, Grid *grid);
+  Burgers(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector);
   ~Burgers();
 
   void ExplicitStep(Field<Array3D<complex>>& fldin, Field<Array3D<complex>>& dfld, real t) override;
@@ -55,7 +58,8 @@ class Burgers : public RightHandSide<Array3D<complex>, Shear> {
 #include "fft.hpp"
 
 template <typename Shear>
-Burgers<Shear>::Burgers(Input &input, Grid *grid) : RightHandSide<Array3D<complex>, Shear>(input, grid) {
+Burgers<Shear>::Burgers(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector) :
+                        RightHandSide<Array3D<complex>, Shear>(input, grid, rhsVector) {
   // Allocate all of the temporary arrays
   vr = astra::makeArray<Array3D<real>>("Burgers::vr1", grid->npr_t);
   wr11 = astra::makeArray<Array3D<real>>("Burgers::wr11", grid->npr_t);
diff --git a/src/rightHandSide/compressible_hydro.hpp b/src/rightHandSide/compressible_hydro.hpp
index e45b726..b723b02 100644
--- a/src/rightHandSide/compressible_hydro.hpp
+++ b/src/rightHandSide/compressible_hydro.hpp
@@ -9,6 +9,7 @@
 #ifndef RIGHTHANDSIDE_COMPRESSIBLE_HYDRO_HPP_
 #define RIGHTHANDSIDE_COMPRESSIBLE_HYDRO_HPP_
 
+#include <memory>
 #include <vector>
 #include <string>
 #include "rightHandSide.hpp"
@@ -16,13 +17,15 @@
 #include "arrays.hpp"
 #include "shear.hpp"
 
+using RhsPtr = std::unique_ptr<RightHandSideConcept<Array3D<complex>>>;
+
 class Grid;
 
 // A class for the hydrodynamics right hand side
 template <typename Shear>
 class CompressibleHydro : public RightHandSide<Array3D<complex>, Shear> {
  public:
-  CompressibleHydro(Input &input, Grid *grid);
+  CompressibleHydro(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector);
   ~CompressibleHydro();
 
   void ExplicitStep(Field<Array3D<complex>>& fldin, Field<Array3D<complex>>& dfld, real t) override;
@@ -57,7 +60,8 @@ class CompressibleHydro : public RightHandSide<Array3D<complex>, Shear> {
 #include "fft.hpp"
 
 template <typename Shear>
-CompressibleHydro<Shear>::CompressibleHydro(Input &input, Grid *grid) : RightHandSide<Array3D<complex>, Shear>(input, grid) {
+CompressibleHydro<Shear>::CompressibleHydro(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector) :
+                                            RightHandSide<Array3D<complex>, Shear>(input, grid, rhsVector) {
   // Allocate all of the temporary arrays
   //vr1 = astra::makeArray<Array3D<real>>("CompressibleHydro::vr1", grid->npr_t);
   //vr2 = astra::makeArray<Array3D<real>>("CompressibleHydro::vr2", grid->npr_t);
diff --git a/src/rightHandSide/hydro.hpp b/src/rightHandSide/hydro.hpp
index 8ba94b8..0431e24 100644
--- a/src/rightHandSide/hydro.hpp
+++ b/src/rightHandSide/hydro.hpp
@@ -9,6 +9,7 @@
 #ifndef RIGHTHANDSIDE_HYDRO_HPP_
 #define RIGHTHANDSIDE_HYDRO_HPP_
 
+#include <memory>
 #include <vector>
 #include <string>
 #include "rightHandSide.hpp"
@@ -16,13 +17,18 @@
 #include "arrays.hpp"
 #include "shear.hpp"
 
+
+
+using RhsPtr = std::unique_ptr<RightHandSideConcept<Array3D<complex>>>;
+
 class Grid;
+template <typename S> class Boussinesq;
 
 // A class for the hydrodynamics right hand side
 template <typename Shear>
 class Hydro : public RightHandSide<Array3D<complex>, Shear> {
  public:
-  Hydro(Input &input, Grid *grid);
+  Hydro(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector);
   ~Hydro();
 
   void ExplicitStep(Field<Array3D<complex>>& fldin, Field<Array3D<complex>>& dfld, real t) override;
@@ -42,6 +48,8 @@ class Hydro : public RightHandSide<Array3D<complex>, Shear> {
   Array3D<real> wr11,wr12,wr13,wr22,wr23,wr33;
   Array3D<complex> wf11,wf12,wf13,wf22,wf23,wf33;
   std::array<int,3> npr, npf;
+
+  friend class Boussinesq<Shear>;
 };
 
 
@@ -53,7 +61,8 @@ class Hydro : public RightHandSide<Array3D<complex>, Shear> {
 #include "fft.hpp"
 
 template <typename Shear>
-Hydro<Shear>::Hydro(Input &input, Grid *grid) : RightHandSide<Array3D<complex>, Shear>(input, grid) {
+Hydro<Shear>::Hydro(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector) :
+                    RightHandSide<Array3D<complex>, Shear>(input, grid, rhsVector) {
   // Allocate all of the temporary arrays
   vr1 = astra::makeArray<Array3D<real>>("Hydro::vr1", grid->npr_t);
   vr2 = astra::makeArray<Array3D<real>>("Hydro::vr2", grid->npr_t);
diff --git a/src/rightHandSide/mhd.hpp b/src/rightHandSide/mhd.hpp
index 84deab6..4b3b4fb 100644
--- a/src/rightHandSide/mhd.hpp
+++ b/src/rightHandSide/mhd.hpp
@@ -9,6 +9,7 @@
 #ifndef RIGHTHANDSIDE_MHD_HPP_
 #define RIGHTHANDSIDE_MHD_HPP_
 
+#include <memory>
 #include <string>
 #include <vector>
 #include "field.hpp"
@@ -17,13 +18,16 @@
 #include "arrays.hpp"
 #include "shear.hpp"
 
+using RhsPtr = std::unique_ptr<RightHandSideConcept<Array3D<complex>>>;
+
 class Grid;
+template <typename S> class Boussinesq;
 
 // A class for the hydrodynamics right hand side
 template <typename Shear>
 class Mhd : public RightHandSide<Array3D<complex>, Shear> {
  public:
-  Mhd(Input &input, Grid *grid);
+  Mhd(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector);
   ~Mhd();
 
   void ExplicitStep(Field<Array3D<complex>>& fldin, Field<Array3D<complex>>& dfld, real t) override;
@@ -46,6 +50,8 @@ class Mhd : public RightHandSide<Array3D<complex>, Shear> {
   Array3D<real> wr11,wr12,wr13,wr22,wr23,wr33;
   Array3D<complex> wf11,wf12,wf13,wf22,wf23,wf33;
   std::array<int,3> npr, npf;
+
+  friend class Boussinesq<Shear>;
 };
 
 
@@ -57,7 +63,8 @@ class Mhd : public RightHandSide<Array3D<complex>, Shear> {
 #include "fft.hpp"
 
 template <typename Shear>
-Mhd<Shear>::Mhd(Input &input, Grid *grid) : RightHandSide<Array3D<complex>, Shear>(input, grid) {
+Mhd<Shear>::Mhd(Input &input, Grid *grid, std::vector<RhsPtr> &rhsVector) :
+                RightHandSide<Array3D<complex>, Shear>(input, grid, rhsVector) {
   // Allocate all of the temporary arrays
   vr1 = astra::makeArray<Array3D<real>>("Mhd::vr1", grid->npr_t);
   vr2 = astra::makeArray<Array3D<real>>("Mhd::vr2", grid->npr_t);
diff --git a/src/rightHandSide/rightHandSide.hpp b/src/rightHandSide/rightHandSide.hpp
index 6db7303..d286ebe 100644
--- a/src/rightHandSide/rightHandSide.hpp
+++ b/src/rightHandSide/rightHandSide.hpp
@@ -33,7 +33,7 @@ class RightHandSideConcept {
 template <typename T, typename Shear>
 class RightHandSide : public RightHandSideConcept<T> {
  public:
-  RightHandSide(Input &input, Grid *grid) : grid(grid), shear(input, grid) {}
+  RightHandSide(Input &input, Grid *grid, std::vector<std::unique_ptr<RightHandSideConcept<T>>> &rhsVector) : grid(grid), shear(input, grid) {}
   virtual ~RightHandSide() {}
 
   virtual void ExplicitStep(Field<T>& fldin, Field<T>& dfld, real t) = 0;
diff --git a/src/rightHandSide/rightHandSideFactory.hpp b/src/rightHandSide/rightHandSideFactory.hpp
index 78bd581..a2f7967 100644
--- a/src/rightHandSide/rightHandSideFactory.hpp
+++ b/src/rightHandSide/rightHandSideFactory.hpp
@@ -18,6 +18,7 @@
 #include "grid.hpp"
 #include "advection.hpp"
 #include "burgers.hpp"
+#include "boussinesq.hpp"
 #include "hydro.hpp"
 #include "compressible_hydro.hpp"
 #include "mhd.hpp"
@@ -50,33 +51,39 @@ class RightHandSideFactory {
       std::string method = input.Get<std::string>("Physics","rhs",irhs);
       if(method == "advection") {
         if(shearType == ShearType::NoShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<Advection<NoShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<Advection<NoShear>>(input, grid, rhsVector)));
         } else if(shearType == ShearType::LinearShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<Advection<LinearShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<Advection<LinearShear>>(input, grid, rhsVector)));
         }
       } else if(method == "burgers") {
         if(shearType == ShearType::NoShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<Burgers<NoShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<Burgers<NoShear>>(input, grid, rhsVector)));
         } else if(shearType == ShearType::LinearShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<Burgers<LinearShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<Burgers<LinearShear>>(input, grid, rhsVector)));
         }
       } else if(method == "hydro") {
         if(shearType == ShearType::NoShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<Hydro<NoShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<Hydro<NoShear>>(input, grid, rhsVector)));
         } else if(shearType == ShearType::LinearShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<Hydro<LinearShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<Hydro<LinearShear>>(input, grid, rhsVector)));
         }
       } else if(method == "compressible_hydro") {
         if(shearType == ShearType::NoShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<CompressibleHydro<NoShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<CompressibleHydro<NoShear>>(input, grid, rhsVector)));
         } else if(shearType == ShearType::LinearShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<CompressibleHydro<LinearShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<CompressibleHydro<LinearShear>>(input, grid, rhsVector)));
         }
       } else if(method == "mhd") {
         if(shearType == ShearType::NoShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<Mhd<NoShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<Mhd<NoShear>>(input, grid, rhsVector)));
         } else if(shearType == ShearType::LinearShear) {
-          rhsVector.emplace_back(std::move(std::make_unique<Mhd<LinearShear>>(input, grid)));
+          rhsVector.emplace_back(std::move(std::make_unique<Mhd<LinearShear>>(input, grid, rhsVector)));
+        }
+      } else if(method == "boussinesq") {
+        if(shearType == ShearType::NoShear) {
+          rhsVector.emplace_back(std::move(std::make_unique<Boussinesq<NoShear>>(input, grid, rhsVector)));
+        } else if(shearType == ShearType::LinearShear) {
+          rhsVector.emplace_back(std::move(std::make_unique<Boussinesq<LinearShear>>(input, grid, rhsVector)));
         }
       } else {
         throw std::runtime_error("Unknown right hand side method: " + method);