diff --git a/include/aspect/newton.h b/include/aspect/newton.h
index 7615a346575..1a11e4717c5 100644
--- a/include/aspect/newton.h
+++ b/include/aspect/newton.h
@@ -307,74 +307,6 @@ namespace aspect
         execute(internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
                 internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
     };
-
-    /**
-     * This class assembles the right-hand-side term of the Newton Stokes system
-     * that is caused by the compressibility in the mass conservation equation.
-     * This function approximates this term as
-     * $- \nabla \mathbf{u} = \frac{1}{\rho} \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u}$
-     */
-    template <int dim>
-    class NewtonStokesReferenceDensityCompressibilityTerm : public Assemblers::Interface<dim>,
-      public SimulatorAccess<dim>
-    {
-      public:
-        void
-        execute(internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
-                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
-    };
-
-    /**
-     * This class assembles the compressibility term of the Newton Stokes system
-     * that is caused by the compressibility in the mass conservation equation.
-     * It includes this term implicitly in the matrix,
-     * which is therefore not longer symmetric.
-     * This function approximates this term as
-     * $ - \nabla \mathbf{u} - \frac{1}{\rho} \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$
-     */
-    template <int dim>
-    class NewtonStokesImplicitReferenceDensityCompressibilityTerm : public Assemblers::Interface<dim>,
-      public SimulatorAccess<dim>
-    {
-      public:
-        void
-        execute(internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
-                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
-    };
-
-    /**
-     * This class assembles the right-hand-side term of the Newton Stokes system
-     * that is caused by the compressibility in the mass conservation equation.
-     * This function approximates this term as
-     * $ - \nabla \mathbf{u} = \frac{1}{\rho} \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$
-     */
-    template <int dim>
-    class NewtonStokesIsentropicCompressionTerm : public Assemblers::Interface<dim>,
-      public SimulatorAccess<dim>
-    {
-      public:
-        void
-        execute(internal::Assembly::Scratch::ScratchBase<dim>  &scratch_base,
-                internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
-    };
-
-    /**
-     * This class assembles the right-hand-side term of the Stokes equation
-     * that is caused by the variable density in the mass conservation equation.
-     * This class approximates this term as
-     * $ - \nabla \cdot \mathbf{u} = \frac{1}{\rho} \frac{\partial \rho}{\partial t} + \frac{1}{\rho} \nabla \rho \cdot \mathbf{u}$
-     * where the right-hand side velocity is explicitly taken from the last timestep,
-     * and the density is taken from a compositional field of the type 'density'.
-     */
-    template <int dim>
-    class NewtonStokesProjectedDensityFieldTerm : public Assemblers::Interface<dim>,
-      public SimulatorAccess<dim>
-    {
-      public:
-        void
-        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch,
-                internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
-    };
   }
 }
 
diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc
index 47eee3d1df7..19f186f7330 100644
--- a/source/simulator/assemblers/newton_stokes.cc
+++ b/source/simulator/assemblers/newton_stokes.cc
@@ -992,259 +992,6 @@ namespace aspect
             }
         }
     }
-
-
-
-    template <int dim>
-    void
-    NewtonStokesReferenceDensityCompressibilityTerm<dim>::
-    execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
-             internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
-    {
-      internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
-      internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);
-
-      // assemble RHS of:
-      //  - div u = 1/rho * drho/dz g/||g||* u
-      Assert(this->get_parameters().formulation_mass_conservation ==
-             Parameters<dim>::Formulation::MassConservation::reference_density_profile,
-             ExcInternalError());
-
-      const Introspection<dim> &introspection = this->introspection();
-      const FiniteElement<dim> &fe = this->get_fe();
-      const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
-      const unsigned int n_q_points    = scratch.finite_element_values.n_quadrature_points;
-      const double pressure_scaling = this->get_pressure_scaling();
-
-      for (unsigned int q=0; q<n_q_points; ++q)
-        {
-          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
-            {
-              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
-                {
-                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
-                  ++i_stokes;
-                }
-              ++i;
-            }
-
-          const Tensor<1,dim>
-          gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
-          const double drho_dz_u = scratch.reference_densities_depth_derivative[q]
-                                   * (gravity * scratch.velocity_values[q]) / gravity.norm();
-          const double one_over_rho = 1.0/scratch.reference_densities[q];
-          const double JxW = scratch.finite_element_values.JxW(q);
-
-          for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
-            data.local_rhs(i) += (pressure_scaling *
-                                  one_over_rho * drho_dz_u * scratch.phi_p[i])
-                                 * JxW;
-        }
-    }
-
-
-
-    template <int dim>
-    void
-    NewtonStokesImplicitReferenceDensityCompressibilityTerm<dim>::
-    execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
-             internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
-    {
-      internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
-      internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);
-
-      // assemble compressibility term of:
-      //  - div u - 1/rho * drho/dz g/||g||* u = 0
-      Assert(this->get_parameters().formulation_mass_conservation ==
-             Parameters<dim>::Formulation::MassConservation::implicit_reference_density_profile,
-             ExcInternalError());
-
-      if (!scratch.rebuild_stokes_matrix)
-        return;
-
-      const Introspection<dim> &introspection = this->introspection();
-      const FiniteElement<dim> &fe = this->get_fe();
-      const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
-      const unsigned int n_q_points    = scratch.finite_element_values.n_quadrature_points;
-      const double pressure_scaling = this->get_pressure_scaling();
-
-      for (unsigned int q=0; q<n_q_points; ++q)
-        {
-          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
-            {
-              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
-                {
-                  scratch.phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].value (i,q);
-                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i,q);
-                  ++i_stokes;
-                }
-              ++i;
-            }
-
-          const Tensor<1,dim>
-          gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
-          const Tensor<1,dim> drho_dz = scratch.reference_densities_depth_derivative[q]
-                                        * gravity / gravity.norm();
-          const double one_over_rho = 1.0/scratch.reference_densities[q];
-          const double JxW = scratch.finite_element_values.JxW(q);
-
-          for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
-            for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
-              data.local_matrix(i,j) += (pressure_scaling *
-                                         one_over_rho * drho_dz * scratch.phi_u[j] * scratch.phi_p[i])
-                                        * JxW;
-        }
-    }
-
-
-
-    template <int dim>
-    void
-    NewtonStokesIsentropicCompressionTerm<dim>::
-    execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
-             internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
-    {
-      internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
-      internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);
-
-      // assemble RHS of:
-      //  - div u = 1/rho * drho/dp rho * g * u
-      Assert(this->get_parameters().formulation_mass_conservation ==
-             Parameters<dim>::Formulation::MassConservation::isentropic_compression,
-             ExcInternalError());
-
-      const Introspection<dim> &introspection = this->introspection();
-      const FiniteElement<dim> &fe = this->get_fe();
-      const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
-      const unsigned int n_q_points    = scratch.finite_element_values.n_quadrature_points;
-      const double pressure_scaling = this->get_pressure_scaling();
-
-      for (unsigned int q=0; q<n_q_points; ++q)
-        {
-          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
-            {
-              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
-                {
-                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
-                  ++i_stokes;
-                }
-              ++i;
-            }
-
-          const Tensor<1,dim>
-          gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
-
-          const double compressibility
-            = scratch.material_model_outputs.compressibilities[q];
-
-          const double density = scratch.material_model_outputs.densities[q];
-          const double JxW = scratch.finite_element_values.JxW(q);
-
-          for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
-            data.local_rhs(i) += (
-                                   // add the term that results from the compressibility. compared
-                                   // to the manual, this term seems to have the wrong sign, but this
-                                   // is because we negate the entire equation to make sure we get
-                                   // -div(u) as the adjoint operator of grad(p)
-                                   (pressure_scaling *
-                                    compressibility * density *
-                                    (scratch.velocity_values[q] * gravity) *
-                                    scratch.phi_p[i])
-                                 )
-                                 * JxW;
-        }
-    }
-
-
-
-    template <int dim>
-    void
-    NewtonStokesProjectedDensityFieldTerm<dim>::
-    execute (internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
-             internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
-    {
-      internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
-      internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);
-
-      // assemble RHS of:
-      // $ - \nabla \cdot \mathbf{u} = \frac{1}{\rho} \frac{\partial \rho}{\partial t} + \frac{1}{\rho} \nabla \rho \cdot \mathbf{u}$
-
-      // Compared to the manual, this term seems to have the wrong sign, but
-      // this is because we negate the entire equation to make sure we get
-      // -div(u) as the adjoint operator of grad(p)
-
-      Assert(this->get_parameters().formulation_mass_conservation ==
-             Parameters<dim>::Formulation::MassConservation::projected_density_field,
-             ExcInternalError());
-
-      const Introspection<dim> &introspection = this->introspection();
-      const FiniteElement<dim> &fe = this->get_fe();
-      const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
-      const unsigned int n_q_points    = scratch.finite_element_values.n_quadrature_points;
-      const double pressure_scaling = this->get_pressure_scaling();
-      const unsigned int density_idx = this->introspection().find_composition_type(CompositionalFieldDescription::density);
-
-      const double time_step = this->get_timestep();
-      const double old_time_step = this->get_old_timestep();
-
-      std::vector<Tensor<1,dim>> density_gradients(n_q_points);
-      std::vector<double> density(n_q_points);
-      std::vector<double> density_old(n_q_points);
-      std::vector<double> density_old_old(n_q_points);
-
-      scratch.finite_element_values[introspection.extractors.compositional_fields[density_idx]].get_function_gradients (this->get_current_linearization_point(),
-          density_gradients);
-      scratch.finite_element_values[introspection.extractors.compositional_fields[density_idx]].get_function_values (this->get_current_linearization_point(),
-          density);
-      scratch.finite_element_values[introspection.extractors.compositional_fields[density_idx]].get_function_values (this->get_old_solution(),
-          density_old);
-      scratch.finite_element_values[introspection.extractors.compositional_fields[density_idx]].get_function_values (this->get_old_old_solution(),
-          density_old_old);
-
-      for (unsigned int q=0; q<n_q_points; ++q)
-        {
-          for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
-            {
-              if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
-                {
-                  scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
-                  ++i_stokes;
-                }
-              ++i;
-            }
-
-          double drho_dt;
-
-          if (this->get_timestep_number() > 1)
-            drho_dt = (1.0/time_step) *
-                      (density[q] *
-                       (2*time_step + old_time_step) / (time_step + old_time_step)
-                       -
-                       density_old[q] *
-                       (1 + time_step/old_time_step)
-                       +
-                       density_old_old[q] *
-                       (time_step * time_step) / (old_time_step * (time_step + old_time_step)));
-          else if (this->get_timestep_number() == 1)
-            drho_dt =
-              (density[q] - density_old[q]) / time_step;
-          else
-            drho_dt = 0.0;
-
-          const double JxW = scratch.finite_element_values.JxW(q);
-
-          for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
-            data.local_rhs(i) += (
-                                   (pressure_scaling *
-                                    (1.0 / density[q]) *
-                                    (density_gradients[q] *
-                                     scratch.velocity_values[q]) *
-                                    scratch.phi_p[i])
-                                   + pressure_scaling * (1.0 / density[q]) * drho_dt * scratch.phi_p[i]
-                                 )
-                                 * JxW;
-        }
-    }
   }
 } // namespace aspect
 
@@ -1258,11 +1005,7 @@ namespace aspect
   template class NewtonStokesPreconditioner<dim>; \
   template class NewtonStokesCompressiblePreconditioner<dim>; \
   template class NewtonStokesIncompressibleTerms<dim>; \
-  template class NewtonStokesCompressibleStrainRateViscosityTerm<dim>; \
-  template class NewtonStokesReferenceDensityCompressibilityTerm<dim>; \
-  template class NewtonStokesImplicitReferenceDensityCompressibilityTerm<dim>; \
-  template class NewtonStokesIsentropicCompressionTerm<dim>; \
-  template class NewtonStokesProjectedDensityFieldTerm<dim>;
+  template class NewtonStokesCompressibleStrainRateViscosityTerm<dim>;
 
     ASPECT_INSTANTIATE(INSTANTIATE)
 
diff --git a/source/simulator/newton.cc b/source/simulator/newton.cc
index 408611fcacf..06a77d29ebb 100644
--- a/source/simulator/newton.cc
+++ b/source/simulator/newton.cc
@@ -67,13 +67,13 @@ namespace aspect
         Parameters<dim>::Formulation::MassConservation::implicit_reference_density_profile)
       {
         assemblers.stokes_system.push_back(
-          std::make_unique<aspect::Assemblers::NewtonStokesImplicitReferenceDensityCompressibilityTerm<dim>>());
+          std::make_unique<aspect::Assemblers::StokesImplicitReferenceDensityCompressibilityTerm<dim>>());
       }
     else if (this->get_parameters().formulation_mass_conservation ==
              Parameters<dim>::Formulation::MassConservation::reference_density_profile)
       {
         assemblers.stokes_system.push_back(
-          std::make_unique<aspect::Assemblers::NewtonStokesReferenceDensityCompressibilityTerm<dim>>());
+          std::make_unique<aspect::Assemblers::StokesReferenceDensityCompressibilityTerm<dim>>());
       }
     else if (this->get_parameters().formulation_mass_conservation ==
              Parameters<dim>::Formulation::MassConservation::incompressible)
@@ -84,14 +84,14 @@ namespace aspect
              Parameters<dim>::Formulation::MassConservation::isentropic_compression)
       {
         assemblers.stokes_system.push_back(
-          std::make_unique<aspect::Assemblers::NewtonStokesIsentropicCompressionTerm<dim>>());
+          std::make_unique<aspect::Assemblers::StokesIsentropicCompressionTerm<dim>>());
       }
     else if (this->get_parameters().formulation_mass_conservation ==
              Parameters<dim>::Formulation::MassConservation::projected_density_field)
       {
         CitationInfo::add("pda");
         assemblers.stokes_system.push_back(
-          std::make_unique<aspect::Assemblers::NewtonStokesProjectedDensityFieldTerm<dim>>());
+          std::make_unique<aspect::Assemblers::StokesProjectedDensityFieldTerm<dim>>());
       }
     else
       AssertThrow(false,