diff --git a/okl/elliptic/ellipticLinAlg.okl b/okl/elliptic/ellipticLinAlg.okl
index e999cdf80..410f3f9d8 100644
--- a/okl/elliptic/ellipticLinAlg.okl
+++ b/okl/elliptic/ellipticLinAlg.okl
@@ -37,6 +37,18 @@ SOFTWARE.
   }
 }
 
+@kernel void fusedCopyDfloatToPfloat(const dlong N,
+                      @restrict const dfloat * x_dfloat,
+                      @restrict const dfloat * y_dfloat,
+                      @restrict pfloat * x_pfloat,
+                      @restrict pfloat * y_pfloat){
+  for(dlong n=0;n<N;++n;@tile(p_blockSize,@outer,@inner)){
+    if(n<N){
+      x_pfloat[n] = x_dfloat[n];
+      y_pfloat[n] = y_dfloat[n];
+    }
+  }
+}
 @kernel void copyDfloatToPfloat(const dlong N,
                       @restrict pfloat *  y,
                       @restrict const dfloat * x){
diff --git a/src/elliptic/elliptic.h b/src/elliptic/elliptic.h
index df57e6b81..f11efd14d 100644
--- a/src/elliptic/elliptic.h
+++ b/src/elliptic/elliptic.h
@@ -127,6 +127,7 @@ struct elliptic_t
   occa::kernel scaledAddPfloatKernel;
   occa::kernel dotMultiplyPfloatKernel;
   occa::kernel copyDfloatToPfloatKernel;
+  occa::kernel fusedCopyDfloatToPfloatKernel;
   occa::kernel copyPfloatToDPfloatKernel;
   occa::kernel updateSmoothedSolutionVecKernel;
   occa::kernel updateChebyshevSolutionVecKernel;
diff --git a/src/elliptic/ellipticMultiGrid.h b/src/elliptic/ellipticMultiGrid.h
index 97380a4ae..0aa10ea5e 100644
--- a/src/elliptic/ellipticMultiGrid.h
+++ b/src/elliptic/ellipticMultiGrid.h
@@ -27,11 +27,16 @@
 #ifndef ELLIPTIC_MGLEVEL_HPP
 #define ELLIPTIC_MGLEVEL_HPP
 
-typedef enum {RICHARDSON = 1,
-              CHEBYSHEV = 2,
-              SCHWARZ = 3} SmoothType;
-typedef enum {JACOBI = 1,
-              SCHWARZ_SMOOTH = 3} SmootherType;
+enum class SmootherType
+{
+  CHEBYSHEV,
+  SCHWARZ,
+};
+enum class SecondarySmootherType
+{  
+  JACOBI,
+  SCHWARZ,
+};
 
 class MGLevel : public parAlmond::multigridLevel
 {
@@ -50,9 +55,9 @@ class MGLevel : public parAlmond::multigridLevel
   occa::memory o_invDegree;
 
   //smoothing params
-  SmoothType stype;
-  SmootherType smtypeUp;
-  SmootherType smtypeDown;
+  SmootherType stype;
+  SecondarySmootherType smtypeUp;
+  SecondarySmootherType smtypeDown;
 
   dfloat lambda1, lambda0;
   int ChebyshevIterations;
@@ -128,7 +133,6 @@ class MGLevel : public parAlmond::multigridLevel
 
   void smoother(occa::memory o_x, occa::memory o_Sx, bool xIsZero);
 
-  void smoothRichardson(occa::memory &o_r, occa::memory &o_x, bool xIsZero);
   void smoothChebyshev (occa::memory &o_r, occa::memory &o_x, bool xIsZero);
   void smoothSchwarz (occa::memory &o_r, occa::memory &o_x, bool xIsZero);
 
@@ -139,7 +143,6 @@ class MGLevel : public parAlmond::multigridLevel
   void setupSmoother(elliptic_t* base);
   dfloat maxEigSmoothAx();
 
-  void buildCoarsenerTriTet(mesh_t** meshLevels, int Nf, int Nc);
   void buildCoarsenerQuadHex(mesh_t** meshLevels, int Nf, int Nc);
 private:
   void smoothChebyshevOneIteration (occa::memory &o_r, occa::memory &o_x, bool xIsZero);
diff --git a/src/elliptic/ellipticMultiGridLevel.cpp b/src/elliptic/ellipticMultiGridLevel.cpp
index 202022e50..1df450b7f 100644
--- a/src/elliptic/ellipticMultiGridLevel.cpp
+++ b/src/elliptic/ellipticMultiGridLevel.cpp
@@ -34,7 +34,7 @@ void MGLevel::Ax(occa::memory o_x, occa::memory o_Ax)
 
 void MGLevel::residual(occa::memory o_rhs, occa::memory o_x, occa::memory o_res)
 {
-  if(stype != SCHWARZ) {
+  if(stype != SmootherType::SCHWARZ) {
     ellipticOperator(elliptic,o_x,o_res, dfloatString);
     // subtract r = b - A*x
     const dlong Nlocal = mesh->Np * mesh->Nelements;
@@ -73,24 +73,18 @@ void MGLevel::prolongate(occa::memory o_x, occa::memory o_Px)
 
 void MGLevel::smooth(occa::memory o_rhs, occa::memory o_x, bool x_is_zero)
 {
-  if(!x_is_zero && stype == SCHWARZ) return;
+  if(!x_is_zero && stype == SmootherType::SCHWARZ) return;
   if(!strstr(pfloatString,dfloatString)) {
-    elliptic->copyDfloatToPfloatKernel(Nrows, o_xPfloat, o_x);
-    elliptic->copyDfloatToPfloatKernel(Nrows, o_rhsPfloat, o_rhs);
-    if (stype == RICHARDSON)
-      this->smoothRichardson(o_rhsPfloat, o_xPfloat, x_is_zero);
-    else if (stype == CHEBYSHEV)
+    elliptic->fusedCopyDfloatToPfloatKernel(Nrows, o_x, o_rhs, o_xPfloat, o_rhsPfloat);
+    if (stype == SmootherType::CHEBYSHEV)
       this->smoothChebyshev(o_rhsPfloat, o_xPfloat, x_is_zero);
-    else if (stype == SCHWARZ)
+    else
       this->smoothSchwarz(o_rhsPfloat, o_xPfloat, x_is_zero);
     elliptic->copyPfloatToDPfloatKernel(Nrows, o_xPfloat, o_x);
-    elliptic->copyPfloatToDPfloatKernel(Nrows, o_rhsPfloat, o_rhs);
   } else {
-    if (stype == RICHARDSON)
-      this->smoothRichardson(o_rhs, o_x, x_is_zero);
-    else if (stype == CHEBYSHEV)
+    if (stype == SmootherType::CHEBYSHEV)
       this->smoothChebyshev(o_rhs, o_x, x_is_zero);
-    else if (stype == SCHWARZ)
+    else
       this->smoothSchwarz(o_rhs, o_x, x_is_zero);
   }
 }
@@ -99,41 +93,18 @@ void MGLevel::smoother(occa::memory o_x, occa::memory o_Sx, bool x_is_zero)
 {
   // x_is_zero = true <-> downward leg
   if(x_is_zero) {
-    if (smtypeDown == JACOBI)
+    if (smtypeDown == SecondarySmootherType::JACOBI)
       this->smootherJacobi(o_x, o_Sx);
-    else if (smtypeDown == SCHWARZ_SMOOTH)
-      //this->smootherSchwarz(o_x, o_Sx);
+    else
       this->smoothSchwarz(o_x, o_Sx, true); // no-op if false
   } else {
-    if (smtypeUp == JACOBI)
+    if (smtypeUp == SecondarySmootherType::JACOBI)
       this->smootherJacobi(o_x, o_Sx);
-    else if (smtypeUp == SCHWARZ_SMOOTH)
-      //this->smootherSchwarz(o_x, o_Sx);
+    else
       this->smoothSchwarz(o_x, o_Sx, true); // no-op if false
   }
 }
 
-void MGLevel::smoothRichardson(occa::memory &o_r, occa::memory &o_x, bool xIsZero)
-{
-  occa::memory o_res = o_smootherResidual;
-
-  if (xIsZero) {
-    this->smoother(o_r, o_x, xIsZero);
-    return;
-  }
-
-  pfloat one = 1.;
-  pfloat mone = -1.;
-
-  //res = r-Ax
-  this->Ax(o_x,o_res);
-  elliptic->scaledAddPfloatKernel(Nrows, one, o_r, mone, o_res);
-
-  //smooth the fine problem x = x + S(r-Ax)
-  this->smoother(o_res, o_res, xIsZero);
-  elliptic->scaledAddPfloatKernel(Nrows, one, o_res, one, o_x);
-}
-
 void MGLevel::smoothChebyshevOneIteration (occa::memory &o_r, occa::memory &o_x, bool xIsZero)
 {
   const pfloat theta = 0.5 * (lambda1 + lambda0);
diff --git a/src/elliptic/ellipticMultiGridLevelSetup.cpp b/src/elliptic/ellipticMultiGridLevelSetup.cpp
index 2559f091f..a39a2545d 100644
--- a/src/elliptic/ellipticMultiGridLevelSetup.cpp
+++ b/src/elliptic/ellipticMultiGridLevelSetup.cpp
@@ -113,15 +113,15 @@ void MGLevel::setupSmoother(elliptic_t* ellipticBase)
 
   if (options.compareArgs("MULTIGRID SMOOTHER","ASM") ||
       options.compareArgs("MULTIGRID SMOOTHER","RAS")) {
-    stype = SCHWARZ;
-    smtypeUp = JACOBI;
-    smtypeDown = JACOBI;
+    stype = SmootherType::SCHWARZ;
+    smtypeUp = SecondarySmootherType::JACOBI;
+    smtypeDown = SecondarySmootherType::JACOBI;
     build(ellipticBase);
 
     if(options.compareArgs("MULTIGRID SMOOTHER","CHEBYSHEV")) {
-      smtypeUp = SCHWARZ_SMOOTH;
-      smtypeDown = SCHWARZ_SMOOTH;
-      stype = CHEBYSHEV;
+      smtypeUp = SecondarySmootherType::SCHWARZ;
+      smtypeDown = SecondarySmootherType::SCHWARZ;
+      stype = SmootherType::CHEBYSHEV;
       if (!options.getArgs("MULTIGRID CHEBYSHEV DEGREE", ChebyshevIterations))
         ChebyshevIterations = 2;   //default to degree 2
       //estimate the max eigenvalue of S*A
@@ -138,13 +138,13 @@ void MGLevel::setupSmoother(elliptic_t* ellipticBase)
         casted_invDiagA[i] = static_cast<pfloat>(invDiagA[i]);
       o_invDiagA = platform->device.malloc(mesh->Np * mesh->Nelements * sizeof(pfloat), casted_invDiagA.data());
       if(options.compareArgs("MULTIGRID UPWARD SMOOTHER","JACOBI"))
-        smtypeUp = JACOBI;
+        smtypeUp = SecondarySmootherType::JACOBI;
       if(options.compareArgs("MULTIGRID DOWNWARD SMOOTHER","JACOBI"))
-        smtypeDown = JACOBI;
+        smtypeDown = SecondarySmootherType::JACOBI;
     }
   } else if (options.compareArgs("MULTIGRID SMOOTHER","DAMPEDJACOBI")) { //default to damped jacobi
-    smtypeUp = JACOBI;
-    smtypeDown = JACOBI;
+    smtypeUp = SecondarySmootherType::JACOBI;
+    smtypeDown = SecondarySmootherType::JACOBI;
     dfloat* invDiagA;
     ellipticBuildJacobi(elliptic,&invDiagA);
     std::vector<pfloat> casted_invDiagA(mesh->Np * mesh->Nelements, 0.0);
@@ -154,7 +154,7 @@ void MGLevel::setupSmoother(elliptic_t* ellipticBase)
     o_invDiagA = platform->device.malloc(mesh->Np * mesh->Nelements * sizeof(pfloat), casted_invDiagA.data());
 
     if (options.compareArgs("MULTIGRID SMOOTHER","CHEBYSHEV")) {
-      stype = CHEBYSHEV;
+      stype = SmootherType::CHEBYSHEV;
 
       if (!options.getArgs("MULTIGRID CHEBYSHEV DEGREE", ChebyshevIterations))
         ChebyshevIterations = 2; //default to degree 2
@@ -165,19 +165,11 @@ void MGLevel::setupSmoother(elliptic_t* ellipticBase)
       lambda1 = 1.1 * rho;
       lambda0 = rho / 10.;
     }else {
-      stype = RICHARDSON;
 
-      //estimate the max eigenvalue of S*A
-      dfloat rho = this->maxEigSmoothAx();
-
-      //set the stabilty weight (jacobi-type interation)
-      lambda0 = (4. / 3.) / rho;
-
-      for (dlong n = 0; n < mesh->Np * mesh->Nelements; n++)
-        casted_invDiagA[n] *= static_cast<pfloat>(lambda0);
-
-      //update diagonal with weight
-      o_invDiagA.copyFrom(casted_invDiagA.data());
+      std::string invalidSmootherName;
+      options.getArgs("MULTIGRID SMOOTHER", invalidSmootherName);
+      if(platform->comm.mpiRank == 0) printf("Smoother %s is not supported!\n", invalidSmootherName.c_str());
+      ABORT(EXIT_FAILURE);
     }
     free(invDiagA);
   }
@@ -201,13 +193,11 @@ void MGLevel::Report()
 
   char smootherString[BUFSIZ];
   if (degree != 1) {
-    if (stype == RICHARDSON && smtypeDown == JACOBI)
-      strcpy(smootherString, "Damped Jacobi    ");
-    else if (stype == CHEBYSHEV && smtypeDown == JACOBI)
+    if (stype == SmootherType::CHEBYSHEV && smtypeDown == SecondarySmootherType::JACOBI)
       strcpy(smootherString, "Chebyshev+Jacobi ");
-    else if (stype == SCHWARZ)
+    else if (stype == SmootherType::SCHWARZ)
       strcpy(smootherString, "Schwarz          ");
-    else if (stype == CHEBYSHEV && smtypeDown == SCHWARZ_SMOOTH)
+    else if (stype == SmootherType::CHEBYSHEV && smtypeDown == SecondarySmootherType::SCHWARZ)
       strcpy(smootherString, "Chebyshev+Schwarz");
   }
 
diff --git a/src/elliptic/ellipticSolveSetup.cpp b/src/elliptic/ellipticSolveSetup.cpp
index 21f897b41..eef68d86e 100644
--- a/src/elliptic/ellipticSolveSetup.cpp
+++ b/src/elliptic/ellipticSolveSetup.cpp
@@ -271,6 +271,10 @@ void ellipticSolveSetup(elliptic_t* elliptic, occa::properties kernelInfo)
                                  "mask",
                                  pfloatKernelInfo);
         filename = install_dir + "/okl/elliptic/ellipticLinAlg.okl";
+        elliptic->fusedCopyDfloatToPfloatKernel =
+          platform->device.buildKernel(filename,
+                                   "fusedCopyDfloatToPfloat",
+                                   kernelInfo);
         elliptic->copyDfloatToPfloatKernel =
           platform->device.buildKernel(filename,
                                    "copyDfloatToPfloat",