diff --git a/src/diffusion/DiffusionTensorOp.H b/src/diffusion/DiffusionTensorOp.H
index cfba723d..1586764d 100644
--- a/src/diffusion/DiffusionTensorOp.H
+++ b/src/diffusion/DiffusionTensorOp.H
@@ -66,6 +66,7 @@ private:
     amrex::Real m_mg_rtol = amrex::Real(1.0e-11);
     amrex::Real m_mg_atol = amrex::Real(1.0e-14);
 #endif
+    bool m_mg_use_gauss_seidel = true;
     std::string m_bottom_solver = "bicgstab";
 };
 
diff --git a/src/diffusion/DiffusionTensorOp.cpp b/src/diffusion/DiffusionTensorOp.cpp
index 7aa030d7..883482a4 100644
--- a/src/diffusion/DiffusionTensorOp.cpp
+++ b/src/diffusion/DiffusionTensorOp.cpp
@@ -59,6 +59,7 @@ DiffusionTensorOp::DiffusionTensorOp (incflo* a_incflo)
                                                 m_incflo->boxArray(0,finest_level),
                                                 m_incflo->DistributionMap(0,finest_level),
                                                 info_solve);
+            m_reg_solve_op->setGaussSeidel(m_mg_use_gauss_seidel);
             m_reg_solve_op->setMaxOrder(m_mg_maxorder);
             m_reg_solve_op->setDomainBC(m_incflo->get_diffuse_tensor_bc(Orientation::low),
                                         m_incflo->get_diffuse_tensor_bc(Orientation::high));
@@ -96,6 +97,8 @@ DiffusionTensorOp::readParameters ()
 
     pp.query("num_pre_smooth", m_num_pre_smooth);
     pp.query("num_post_smooth", m_num_post_smooth);
+
+    pp.query("use_gauss_seidel", m_mg_use_gauss_seidel);
 }
 
 void