diff --git a/doc/modules/changes/20250113_gassmoeller b/doc/modules/changes/20250113_gassmoeller
new file mode 100644
index 00000000000..5d7159fca39
--- /dev/null
+++ b/doc/modules/changes/20250113_gassmoeller
@@ -0,0 +1,5 @@
+Changed: Prescribed compositional fields are now copied from the material
+model output in one combined operation (instead of being copied field by field).
+This change improves the efficiency of the copy process.
+<br>
+(Rene Gassmoeller, 2025/01/13)
diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h
index 7345e29efdc..27848c48831 100644
--- a/include/aspect/simulator.h
+++ b/include/aspect/simulator.h
@@ -1370,6 +1370,25 @@ namespace aspect
        */
       void apply_limiter_to_dg_solutions (const AdvectionField &advection_field);
 
+      /**
+       * Compute the unique support points for the advection fields @p advection_fields.
+       * The support points are collected by taking the union of the support points
+       * of all given advection fields and filtering out duplicate points. The resulting
+       * set of points is written into @p unique_support_points. @p support_point_index_by_field
+       * is a vector of vectors that contains the indices of the support points for each field.
+       * I.e. support_point_index_by_field[i][j] contains the j-th support point
+       * for the i-th field in @p advection_fields and
+       * unique_support_points[support_point_index_by_field[i][j]] is its location.
+       *
+       * For the common case that all fields have the same support points, each vector in
+       * @p support_point_index_by_field will contain the same indices for all fields.
+       *
+       * Note that existing content of @p unique_support_points and @p support_point_index_by_field
+       * will be overwritten in this function.
+       */
+      void compute_unique_advection_support_points (const std::vector<AdvectionField> &advection_fields,
+                                                    std::vector<Point<dim>> &unique_support_points,
+                                                    std::vector<std::vector<unsigned int>> &support_point_index_by_field) const;
 
       /**
        * Compute the reactions in case of operator splitting:
@@ -1439,7 +1458,7 @@ namespace aspect
        * This function is implemented in
        * <code>source/simulator/helper_functions.cc</code>.
        */
-      void interpolate_material_output_into_advection_field (const AdvectionField &adv_field);
+      void interpolate_material_output_into_advection_field (const std::vector<AdvectionField> &adv_field);
 
 
       /**
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index bbd63f18b31..c4803ce814d 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -1701,6 +1701,53 @@ namespace aspect
 
 
 
+  template <int dim>
+  void Simulator<dim>::compute_unique_advection_support_points(const std::vector<AdvectionField> &advection_fields,
+                                                               std::vector<Point<dim>> &unique_support_points,
+                                                               std::vector<std::vector<unsigned int>> &support_point_index_by_field) const
+  {
+    const unsigned int n_fields = advection_fields.size();
+
+    unique_support_points.clear();
+    support_point_index_by_field.clear();
+    support_point_index_by_field.resize(n_fields);
+
+    // Reserve space to avoid most memory reallocations
+    if (n_fields > 0)
+      {
+        const unsigned int likely_number_of_support_points = dof_handler.get_fe().base_element(advection_fields[0].base_element(introspection)).get_unit_support_points().size();
+        unique_support_points.reserve(likely_number_of_support_points * n_fields);
+        for (unsigned int i=0; i<n_fields; ++i)
+          support_point_index_by_field[i].reserve(likely_number_of_support_points);
+      }
+
+    // Loop through all support points and fill the output data structures
+    for (unsigned int i=0; i<advection_fields.size(); ++i)
+      {
+        const std::vector<Point<dim>> &support_points = dof_handler.get_fe().base_element(advection_fields[i].base_element(introspection)).get_unit_support_points();
+
+        for (const auto &support_point: support_points)
+          {
+            // Naive n^2 algorithm to merge all points into the data structure, speed is likely irrelevant, because
+            // n is small and this is only done rarely.
+            const auto it = std::find(unique_support_points.begin(), unique_support_points.end(), support_point);
+            if (it != unique_support_points.end())
+              {
+                // We already have this support point, record its number:
+                support_point_index_by_field[i].push_back(std::distance(unique_support_points.begin(), it));
+              }
+            else
+              {
+                // This is a new point that needs to be added:
+                unique_support_points.push_back(support_point);
+                support_point_index_by_field[i].push_back(unique_support_points.size()-1);
+              }
+          }
+      }
+  }
+
+
+
   template <int dim>
   void Simulator<dim>::compute_reactions ()
   {
@@ -1742,44 +1789,19 @@ namespace aspect
     // write back all values that correspond to support points of that particular field. Field values we computed that are
     // not actually support points, we just throw away.
 
-    // First compute the union of all support points (temperature and compositions):
+    // First compute all unique support points (temperature and compositions):
     std::vector<Point<dim>> unique_support_points;
-    // For each field (T or composition) store where each support point is found in the list of indices of support points
-    // unique_support_points. This means index=support_point_index_by_field[0][3] is the index of the third support point
-    // of the temperature field and unique_support_points[index] is its location.
-    const unsigned int n_fields = 1 + introspection.n_compositional_fields;
-    std::vector<std::vector<unsigned int>> support_point_index_by_field(n_fields);
+    std::vector<std::vector<unsigned int>> support_point_index_by_field;
+    std::vector<AdvectionField> advection_fields;
 
-    // Fill in the support_point_index_by_field data structure. This is a bit complicated...
-    {
-      // first fill in temperature:
-      unique_support_points = dof_handler.get_fe().base_element(introspection.base_elements.temperature).get_unit_support_points();
-      for (unsigned int i=0; i<unique_support_points.size(); ++i)
-        support_point_index_by_field[0].push_back(i);
+    // First add the temperature field
+    advection_fields.push_back(Simulator<dim>::AdvectionField::temperature());
+    // Then add all compositional fields
+    for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
+      advection_fields.push_back(Simulator<dim>::AdvectionField::composition(c));
 
-      // Naive n^2 algorithm to merge all compositional fields into the data structure.
-      for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
-        {
-          const auto composition = Simulator<dim>::AdvectionField::composition(c);
-          const std::vector<Point<dim>> &points = dof_handler.get_fe().base_element(composition.base_element(introspection)).get_unit_support_points();
-
-          for (unsigned int i=0; i<points.size(); ++i)
-            {
-              const auto it = std::find(unique_support_points.begin(), unique_support_points.end(), points[i]);
-              if (it != unique_support_points.end())
-                {
-                  // We already have this support point, record its number:
-                  support_point_index_by_field[1+c].push_back(std::distance(unique_support_points.begin(), it));
-                }
-              else
-                {
-                  // This is a new point that needs to be added:
-                  unique_support_points.push_back(points[i]);
-                  support_point_index_by_field[1+c].push_back(unique_support_points.size()-1);
-                }
-            }
-        }
-    }
+    const unsigned int n_fields = advection_fields.size();
+    compute_unique_advection_support_points(advection_fields, unique_support_points, support_point_index_by_field);
 
     const Quadrature<dim> combined_support_points(unique_support_points);
     FEValues<dim> fe_values (*mapping,
@@ -2068,29 +2090,49 @@ namespace aspect
 
 
   template <int dim>
-  void Simulator<dim>::interpolate_material_output_into_advection_field (const AdvectionField &adv_field)
+  void Simulator<dim>::interpolate_material_output_into_advection_field (const std::vector<AdvectionField> &adv_fields)
   {
-    // we need some temporary vectors to store our updates to composition in
+    // we need a temporary vector to store our updates to the advection fields in
     // before we copy them over to the solution vector in the end
     LinearAlgebra::BlockVector distributed_vector (introspection.index_sets.system_partitioning,
                                                    mpi_communicator);
 
-    if (adv_field.is_temperature())
-      pcout << "   Copying properties into prescribed temperature field... "
-            << std::flush;
-    else
+    if (adv_fields.size() == 0)
+      return;
+    else if (adv_fields.size() == 1)
       {
-        const std::string name_of_field = introspection.name_for_compositional_index(adv_field.compositional_variable);
+        const AdvectionField &adv_field = adv_fields[0];
+        if (adv_field.is_temperature())
+          pcout << "   Copying properties into prescribed temperature field... "
+                << std::flush;
+        else
+          {
+            const std::string name_of_field = introspection.name_for_compositional_index(adv_field.compositional_variable);
 
-        pcout << "   Copying properties into prescribed compositional field " + name_of_field + "... "
+            pcout << "   Copying properties into prescribed compositional field " + name_of_field + "... "
+                  << std::flush;
+          }
+      }
+    else
+      {
+        pcout << "   Copying properties into prescribed compositional fields... "
               << std::flush;
       }
 
+    // Advection fields can have different Finite Element discretizations (degree, continuous/discontinuous)
+    // and will have different support points. We solve this by computing the union of all support points and evaluating
+    // the material output for all these support points for every cell.
+
+    /// First compute all unique support points (temperature and compositions):
+    const unsigned int n_fields = adv_fields.size();
+    std::vector<Point<dim>> unique_support_points;
+    std::vector<std::vector<unsigned int>> support_point_index_by_field;
+    compute_unique_advection_support_points(adv_fields, unique_support_points, support_point_index_by_field);
+
     // Create an FEValues object that allows us to interpolate onto the solution
     // vector. To make this happen, we need to have a quadrature formula that
-    // consists of the support points of the advection field finite element
-    const Quadrature<dim> quadrature(dof_handler.get_fe().base_element(adv_field.base_element(introspection))
-                                     .get_unit_support_points());
+    // consists of the support points of all advection field finite elements
+    const Quadrature<dim> quadrature(unique_support_points);
 
     FEValues<dim> fe_values (*mapping,
                              dof_handler.get_fe(),
@@ -2109,89 +2151,85 @@ namespace aspect
     MaterialModel::PrescribedTemperatureOutputs<dim> *prescribed_temperature_out
       = out.template get_additional_output<MaterialModel::PrescribedTemperatureOutputs<dim>>();
 
-    // check if the material model computes the correct prescribed field outputs
-    if (adv_field.is_temperature())
-      {
-        AssertThrow(prescribed_temperature_out != nullptr,
-                    ExcMessage("You are trying to use a prescribed temperature field, "
-                               "but the material model you use does not support interpolating properties "
-                               "(it does not create PrescribedTemperatureOutputs, which is required for this "
-                               "temperature field type)."));
-      }
-    else
-      {
-        AssertThrow(prescribed_field_out != nullptr,
-                    ExcMessage("You are trying to use a prescribed advection field, "
-                               "but the material model you use does not support interpolating properties "
-                               "(it does not create PrescribedFieldOutputs, which is required for this "
-                               "advection field type)."));
-      }
-
     // Make a loop first over all cells, and then over all degrees of freedom in each element
     // to interpolate material properties onto a solution vector.
-
-    // Note that the values for some degrees of freedom are set more than once in the loop
-    // below where we assign the new values to distributed_vector (if they are located on the
-    // interface between cells), as we loop over all cells, and then over all degrees of freedom
-    // on each cell. But even though we touch some DoF more than once, we always compute the same value,
-    // and then overwrite the same value in distributed_vector.
     for (const auto &cell : dof_handler.active_cell_iterators())
       if (cell->is_locally_owned())
         {
           fe_values.reinit (cell);
           cell->get_dof_indices (local_dof_indices);
           in.reinit(fe_values, cell, introspection, solution);
-
           material_model->evaluate(in, out);
 
-          // Interpolate material properties onto the advection fields
-          const unsigned int advection_dofs_per_cell =
-            dof_handler.get_fe().base_element(adv_field.base_element(introspection)).dofs_per_cell;
-
-          for (unsigned int j=0; j<advection_dofs_per_cell; ++j)
+          for (unsigned int i=0; i<n_fields; ++i)
             {
-              const unsigned int dof_idx
-                = dof_handler.get_fe().component_to_system_index(adv_field.component_index(introspection),
-                                                                 /*dof index within component=*/ j);
+              const AdvectionField &adv_field = adv_fields[i];
+
+              // Interpolate material properties onto the advection fields
+              const unsigned int advection_dofs_per_cell =
+                dof_handler.get_fe().base_element(adv_field.base_element(introspection)).dofs_per_cell;
 
-              // Skip degrees of freedom that are not locally owned. These
-              // will eventually be handled by one of the other processors.
-              if (dof_handler.locally_owned_dofs().is_element(local_dof_indices[dof_idx]))
+              // Make sure data structures have the expected size
+              Assert(advection_dofs_per_cell == support_point_index_by_field[i].size(), ExcInternalError());
+
+              for (unsigned int j=0; j<advection_dofs_per_cell; ++j)
                 {
-                  if (adv_field.is_temperature())
-                    {
-                      Assert(numbers::is_finite(prescribed_temperature_out->prescribed_temperature_outputs[j]),
-                             ExcMessage("You are trying to use a prescribed advection field, "
-                                        "but the material model you use does not fill the PrescribedFieldOutputs "
-                                        "for your prescribed field, which is required for this method."));
+                  const unsigned int dof_idx
+                    = dof_handler.get_fe().component_to_system_index(adv_field.component_index(introspection),
+                                                                     /*dof index within component=*/ j);
 
-                      distributed_vector(local_dof_indices[dof_idx])
-                        = prescribed_temperature_out->prescribed_temperature_outputs[j];
-                    }
-                  else
+                  // Skip degrees of freedom that are not locally owned. These
+                  // will eventually be handled by one of the other processors.
+                  if (dof_handler.locally_owned_dofs().is_element(local_dof_indices[dof_idx]))
                     {
-                      Assert(numbers::is_finite(prescribed_field_out->prescribed_field_outputs[j][adv_field.compositional_variable]),
-                             ExcMessage("You are trying to use a prescribed advection field, "
-                                        "but the material model you use does not fill the PrescribedFieldOutputs "
-                                        "for your prescribed field, which is required for this method."));
-
-                      distributed_vector(local_dof_indices[dof_idx])
-                        = prescribed_field_out->prescribed_field_outputs[j][adv_field.compositional_variable];
+                      if (adv_field.is_temperature())
+                        {
+                          Assert(prescribed_temperature_out != nullptr,
+                                 ExcMessage("You are trying to use a prescribed temperature field, "
+                                            "but the material model you use does not support interpolating properties "
+                                            "(it does not create PrescribedTemperatureOutputs, which is required for this "
+                                            "temperature field type)."));
+                          Assert(numbers::is_finite(prescribed_temperature_out->prescribed_temperature_outputs[j]),
+                                 ExcMessage("You are trying to use a prescribed advection field, "
+                                            "but the material model you use does not fill the PrescribedFieldOutputs "
+                                            "for your prescribed field, which is required for this method."));
+
+                          distributed_vector(local_dof_indices[dof_idx])
+                            = prescribed_temperature_out->prescribed_temperature_outputs[support_point_index_by_field[i][j]];
+                        }
+                      else
+                        {
+                          Assert(prescribed_field_out != nullptr,
+                                 ExcMessage("You are trying to use a prescribed advection field, "
+                                            "but the material model you use does not support interpolating properties "
+                                            "(it does not create PrescribedFieldOutputs, which is required for this "
+                                            "advection field type)."));
+                          Assert(numbers::is_finite(prescribed_field_out->prescribed_field_outputs[j][adv_field.compositional_variable]),
+                                 ExcMessage("You are trying to use a prescribed advection field, "
+                                            "but the material model you use does not fill the PrescribedFieldOutputs "
+                                            "for your prescribed field, which is required for this method."));
+
+                          distributed_vector(local_dof_indices[dof_idx])
+                            = prescribed_field_out->prescribed_field_outputs[support_point_index_by_field[i][j]][adv_field.compositional_variable];
+                        }
                     }
                 }
             }
         }
 
-    // Put the final values into the solution vector, also
-    // updating the ghost elements of the 'solution' vector.
-    const unsigned int advection_block = adv_field.block_index(introspection);
-    distributed_vector.block(advection_block).compress(VectorOperation::insert);
+    for (const auto &adv_field: adv_fields)
+      {
+        // Put the final values into the solution vector, also
+        // updating the ghost elements of the solution vector.
+        const unsigned int advection_block = adv_field.block_index(introspection);
+        distributed_vector.block(advection_block).compress(VectorOperation::insert);
 
-    if (adv_field.is_temperature() ||
-        adv_field.compositional_variable != introspection.find_composition_type(CompositionalFieldDescription::density))
-      current_constraints.distribute (distributed_vector);
+        if (adv_field.is_temperature() ||
+            adv_field.compositional_variable != introspection.find_composition_type(CompositionalFieldDescription::density))
+          current_constraints.distribute (distributed_vector);
 
-    solution.block(advection_block) = distributed_vector.block(advection_block);
+        solution.block(advection_block) = distributed_vector.block(advection_block);
+      }
 
     pcout << "done." << std::endl;
   }
@@ -2781,9 +2819,12 @@ namespace aspect
   template bool Simulator<dim>::stokes_A_block_is_symmetric() const; \
   template void Simulator<dim>::interpolate_onto_velocity_system(const TensorFunction<1,dim> &func, LinearAlgebra::Vector &vec) const;\
   template void Simulator<dim>::apply_limiter_to_dg_solutions(const AdvectionField &advection_field); \
+  template void Simulator<dim>::compute_unique_advection_support_points(const std::vector<AdvectionField> &advection_fields, \
+                                                                        std::vector<Point<dim>> &support_points, \
+                                                                        std::vector<std::vector<unsigned int>> &support_point_index_by_field) const; \
   template void Simulator<dim>::compute_reactions(); \
   template void Simulator<dim>::initialize_current_linearization_point (); \
-  template void Simulator<dim>::interpolate_material_output_into_advection_field(const AdvectionField &adv_field); \
+  template void Simulator<dim>::interpolate_material_output_into_advection_field(const std::vector<AdvectionField> &adv_field); \
   template void Simulator<dim>::check_consistency_of_formulation(); \
   template void Simulator<dim>::replace_outflow_boundary_ids(const unsigned int boundary_id_offset); \
   template void Simulator<dim>::restore_outflow_boundary_ids(const unsigned int boundary_id_offset); \
diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc
index 1b377c0e7c2..0cede57e5cb 100644
--- a/source/simulator/solver_schemes.cc
+++ b/source/simulator/solver_schemes.cc
@@ -157,7 +157,7 @@ namespace aspect
             {
               TimerOutput::Scope timer (computing_timer, "Interpolate prescribed temperature");
 
-              interpolate_material_output_into_advection_field(adv_field);
+              interpolate_material_output_into_advection_field({adv_field});
 
               // Also set the old_solution block to the prescribed field. The old
               // solution is the one that is used to assemble the diffusion system in
@@ -180,7 +180,7 @@ namespace aspect
 
           TimerOutput::Scope timer (computing_timer, "Interpolate prescribed temperature");
 
-          interpolate_material_output_into_advection_field(adv_field);
+          interpolate_material_output_into_advection_field({adv_field});
 
           // Call the signal in case the user wants to do something with the variable:
           SolverControl dummy;
@@ -188,6 +188,7 @@ namespace aspect
                                         adv_field.is_temperature(),
                                         adv_field.compositional_variable,
                                         dummy);
+
           break;
         }
 
@@ -245,6 +246,7 @@ namespace aspect
       Assert(residual->size() == introspection.n_compositional_fields, ExcInternalError());
 
     std::vector<AdvectionField> fields_advected_by_particles;
+    std::vector<AdvectionField> fields_interpolated_from_material_output;
 
     for (unsigned int c=0; c < introspection.n_compositional_fields; ++c)
       {
@@ -263,7 +265,7 @@ namespace aspect
                 {
                   TimerOutput::Scope timer (computing_timer, "Interpolate prescribed composition");
 
-                  interpolate_material_output_into_advection_field(adv_field);
+                  interpolate_material_output_into_advection_field({adv_field});
 
                   // Also set the old_solution block to the prescribed field. The old
                   // solution is the one that is used to assemble the diffusion system in
@@ -298,16 +300,8 @@ namespace aspect
 
             case Parameters<dim>::AdvectionFieldMethod::prescribed_field:
             {
-              TimerOutput::Scope timer (computing_timer, "Interpolate prescribed composition");
-
-              interpolate_material_output_into_advection_field(adv_field);
-
-              // Call the signal in case the user wants to do something with the variable:
-              SolverControl dummy;
-              signals.post_advection_solver(*this,
-                                            adv_field.is_temperature(),
-                                            adv_field.compositional_variable,
-                                            dummy);
+              // handle all prescribed fields together to increase efficiency
+              fields_interpolated_from_material_output.push_back(adv_field);
               break;
             }
 
@@ -334,6 +328,23 @@ namespace aspect
           }
       }
 
+    if (fields_interpolated_from_material_output.size() > 0)
+      {
+        TimerOutput::Scope timer (computing_timer, "Interpolate prescribed composition");
+
+        interpolate_material_output_into_advection_field(fields_interpolated_from_material_output);
+
+        for (const auto &adv_field: fields_interpolated_from_material_output)
+          {
+            // Call the signal in case the user wants to do something with the variable:
+            SolverControl dummy;
+            signals.post_advection_solver(*this,
+                                          adv_field.is_temperature(),
+                                          adv_field.compositional_variable,
+                                          dummy);
+          }
+      }
+
     if (fields_advected_by_particles.size() > 0)
       interpolate_particle_properties(fields_advected_by_particles);
 
diff --git a/tests/entropy_half_space/statistics b/tests/entropy_half_space/statistics
index 0bebe1880af..6f25adf508c 100644
--- a/tests/entropy_half_space/statistics
+++ b/tests/entropy_half_space/statistics
@@ -8,8 +8,8 @@
 # 8: Number of nonlinear iterations
 # 9: Iterations for temperature solver
 # 10: Iterations for composition solver 1
-# 11: Iterations for composition solver 2
-# 12: Iterations for composition solver 3
+# 11: Iterations for composition solver 3
+# 12: Iterations for composition solver 2
 # 13: Iterations for Stokes solver
 # 14: Velocity iterations in Stokes preconditioner
 # 15: Schur complement iterations in Stokes preconditioner
diff --git a/tests/multiple_prescribed_fields.cc b/tests/multiple_prescribed_fields.cc
new file mode 100644
index 00000000000..4aa67f62cff
--- /dev/null
+++ b/tests/multiple_prescribed_fields.cc
@@ -0,0 +1,95 @@
+/*
+  Copyright (C) 2022 - 2024 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT 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 General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include <aspect/material_model/simple.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+  namespace MaterialModel
+  {
+    using namespace dealii;
+
+    template <int dim>
+    class PrescribedFieldsMaterial : public MaterialModel::Simple<dim>
+    {
+      public:
+
+        void evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
+                      MaterialModel::MaterialModelOutputs<dim> &out) const override;
+
+        void
+        create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const override;
+    };
+
+  }
+}
+
+namespace aspect
+{
+  namespace MaterialModel
+  {
+
+    template <int dim>
+    void
+    PrescribedFieldsMaterial<dim>::
+    evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
+             MaterialModel::MaterialModelOutputs<dim> &out) const
+    {
+      Simple<dim>::evaluate(in, out);
+
+      // set up variable to interpolate prescribed field outputs onto compositional fields
+      PrescribedFieldOutputs<dim> *prescribed_field_out = out.template get_additional_output<PrescribedFieldOutputs<dim>>();
+
+      if (prescribed_field_out != nullptr)
+        for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
+          {
+            const double y = in.position[i](1);
+            for (unsigned int j=0; j<this->n_compositional_fields(); ++j)
+              prescribed_field_out->prescribed_field_outputs[i][j] = y*j;
+          }
+    }
+
+
+    template <int dim>
+    void
+    PrescribedFieldsMaterial<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
+    {
+      if (out.template get_additional_output<PrescribedFieldOutputs<dim>>() == nullptr)
+        {
+          const unsigned int n_points = out.n_evaluation_points();
+          out.additional_outputs.push_back(
+            std::make_unique<MaterialModel::PrescribedFieldOutputs<dim>> (n_points,this->n_compositional_fields()));
+        }
+    }
+  }
+}
+
+// explicit instantiations
+namespace aspect
+{
+  namespace MaterialModel
+  {
+    ASPECT_REGISTER_MATERIAL_MODEL(PrescribedFieldsMaterial,
+                                   "prescribed fields material",
+                                   "A simple material model that is like the "
+                                   "'Simple' model, but creates multiple prescribed field outputs.")
+  }
+}
diff --git a/tests/multiple_prescribed_fields.prm b/tests/multiple_prescribed_fields.prm
new file mode 100644
index 00000000000..c33b23e8c21
--- /dev/null
+++ b/tests/multiple_prescribed_fields.prm
@@ -0,0 +1,83 @@
+# A testcase that demonstrates the named additional outputs postprocessor for multiple
+# named outputs. The test is based on prescribed_field_with_diffusion only
+# with more than one additional output.
+
+set Dimension = 2
+set End time                               = 0
+set Start time                             = 0
+set Adiabatic surface temperature          = 0
+set Surface pressure                       = 0
+set Use years in output instead of seconds = false
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 1.0
+  end
+end
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 1
+    set Y extent = 1
+  end
+end
+
+subsection Compositional fields
+  set Number of fields = 3
+  set Names of fields  = prescribed_1, prescribed_2, prescribed_3
+  set Compositional field methods = prescribed field, prescribed field, prescribed field
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.1
+  end
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.0; 0.2; 0.3
+  end
+end
+
+subsection Material model
+  set Model name = prescribed fields material
+
+  subsection Simple model
+    set Reference density             = 1
+    set Reference specific heat       = 1
+    set Reference temperature         = 0
+    set Thermal conductivity          = 1
+    set Thermal expansion coefficient = 1
+    set Viscosity                     = 1
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 1
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = 0, 1, 2, 3
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization, composition statistics
+
+  subsection Visualization
+    set Interpolate output = false
+    set List of output variables      = named additional outputs
+    set Number of grouped files       = 0
+    set Output format                 = gnuplot
+    set Time between graphical output = 0
+  end
+end
diff --git a/tests/multiple_prescribed_fields/screen-output b/tests/multiple_prescribed_fields/screen-output
new file mode 100644
index 00000000000..4f7622364c2
--- /dev/null
+++ b/tests/multiple_prescribed_fields/screen-output
@@ -0,0 +1,19 @@
+
+Loading shared library <./libmultiple_prescribed_fields.debug.so>
+
+Number of active cells: 4 (on 2 levels)
+Number of degrees of freedom: 159 (50+9+25+25+25+25)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Solving temperature system... 0 iterations.
+   Copying properties into prescribed compositional fields... done.
+   Solving Stokes system (GMG)... 6+0 iterations.
+
+   Postprocessing:
+     Writing graphical output:  output-multiple_prescribed_fields/solution/solution-00000
+     Compositions min/max/mass: 0/0/0 // 0/1/0.5 // 0/2/1
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/multiple_prescribed_fields/solution/solution-00000.0000.gnuplot b/tests/multiple_prescribed_fields/solution/solution-00000.0000.gnuplot
new file mode 100644
index 00000000000..49ed29ce85b
--- /dev/null
+++ b/tests/multiple_prescribed_fields/solution/solution-00000.0000.gnuplot
@@ -0,0 +1,35 @@
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <velocity> <velocity> <p> <T> <prescribed_1> <prescribed_2> <prescribed_3> <prescribed_field_output_C0> <prescribed_field_output_C1> <prescribed_field_output_C2> 
+0 0 0 0 0.9 0.1 0 0 0 0 0 0 
+0.5 0 -4.42931e-27 0 0.9 0.1 0 0 0 0 0 0 
+
+0 0.5 0 2.208e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+0.5 0.5 5.45949e-28 -3.03919e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+
+
+0.5 0 -4.42931e-27 0 0.9 0.1 0 0 0 0 0 0 
+1 0 0 0 0.9 0.1 0 0 0 0 0 0 
+
+0.5 0.5 5.45949e-28 -3.03919e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+1 0.5 0 2.208e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+
+
+0 0.5 0 2.208e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+0.5 0.5 5.45949e-28 -3.03919e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+
+0 1 0 0 -5.3257e-11 0.1 0 1 2 0 1 2 
+0.5 1 -5.39777e-27 0 5.3257e-11 0.1 0 1 2 0 1 2 
+
+
+0.5 0.5 5.45949e-28 -3.03919e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+1 0.5 0 2.208e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+
+0.5 1 -5.39777e-27 0 5.3257e-11 0.1 0 1 2 0 1 2 
+1 1 0 0 -5.3257e-11 0.1 0 1 2 0 1 2 
+
+
diff --git a/tests/multiple_prescribed_fields/statistics b/tests/multiple_prescribed_fields/statistics
new file mode 100644
index 00000000000..9e58a680ea8
--- /dev/null
+++ b/tests/multiple_prescribed_fields/statistics
@@ -0,0 +1,25 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Time step size (seconds)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for composition solver 1
+# 10: Iterations for composition solver 2
+# 11: Iterations for composition solver 3
+# 12: Iterations for Stokes solver
+# 13: Velocity iterations in Stokes preconditioner
+# 14: Schur complement iterations in Stokes preconditioner
+# 15: Visualization file name
+# 16: Minimal value for composition prescribed_1
+# 17: Maximal value for composition prescribed_1
+# 18: Global mass for composition prescribed_1
+# 19: Minimal value for composition prescribed_2
+# 20: Maximal value for composition prescribed_2
+# 21: Global mass for composition prescribed_2
+# 22: Minimal value for composition prescribed_3
+# 23: Maximal value for composition prescribed_3
+# 24: Global mass for composition prescribed_3
+0 0.000000000000e+00 0.000000000000e+00 4 59 25 75 0 4294967294 4294967294 4294967294 5 7 7 output-multiple_prescribed_fields/solution/solution-00000 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 5.00000000e-01 0.00000000e+00 2.00000000e+00 1.00000000e+00 
diff --git a/tests/multiple_prescribed_fields_dc.cc b/tests/multiple_prescribed_fields_dc.cc
new file mode 100644
index 00000000000..4aa67f62cff
--- /dev/null
+++ b/tests/multiple_prescribed_fields_dc.cc
@@ -0,0 +1,95 @@
+/*
+  Copyright (C) 2022 - 2024 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT 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 General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include <aspect/material_model/simple.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+  namespace MaterialModel
+  {
+    using namespace dealii;
+
+    template <int dim>
+    class PrescribedFieldsMaterial : public MaterialModel::Simple<dim>
+    {
+      public:
+
+        void evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
+                      MaterialModel::MaterialModelOutputs<dim> &out) const override;
+
+        void
+        create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const override;
+    };
+
+  }
+}
+
+namespace aspect
+{
+  namespace MaterialModel
+  {
+
+    template <int dim>
+    void
+    PrescribedFieldsMaterial<dim>::
+    evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
+             MaterialModel::MaterialModelOutputs<dim> &out) const
+    {
+      Simple<dim>::evaluate(in, out);
+
+      // set up variable to interpolate prescribed field outputs onto compositional fields
+      PrescribedFieldOutputs<dim> *prescribed_field_out = out.template get_additional_output<PrescribedFieldOutputs<dim>>();
+
+      if (prescribed_field_out != nullptr)
+        for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
+          {
+            const double y = in.position[i](1);
+            for (unsigned int j=0; j<this->n_compositional_fields(); ++j)
+              prescribed_field_out->prescribed_field_outputs[i][j] = y*j;
+          }
+    }
+
+
+    template <int dim>
+    void
+    PrescribedFieldsMaterial<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
+    {
+      if (out.template get_additional_output<PrescribedFieldOutputs<dim>>() == nullptr)
+        {
+          const unsigned int n_points = out.n_evaluation_points();
+          out.additional_outputs.push_back(
+            std::make_unique<MaterialModel::PrescribedFieldOutputs<dim>> (n_points,this->n_compositional_fields()));
+        }
+    }
+  }
+}
+
+// explicit instantiations
+namespace aspect
+{
+  namespace MaterialModel
+  {
+    ASPECT_REGISTER_MATERIAL_MODEL(PrescribedFieldsMaterial,
+                                   "prescribed fields material",
+                                   "A simple material model that is like the "
+                                   "'Simple' model, but creates multiple prescribed field outputs.")
+  }
+}
diff --git a/tests/multiple_prescribed_fields_dc.prm b/tests/multiple_prescribed_fields_dc.prm
new file mode 100644
index 00000000000..2150ba83011
--- /dev/null
+++ b/tests/multiple_prescribed_fields_dc.prm
@@ -0,0 +1,88 @@
+# A testcase that demonstrates the named additional outputs postprocessor for multiple
+# named outputs. The test is based on prescribed_field_with_diffusion only
+# with more than one additional output.
+
+set Dimension = 2
+set End time                               = 0
+set Start time                             = 0
+set Adiabatic surface temperature          = 0
+set Surface pressure                       = 0
+set Use years in output instead of seconds = false
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 1.0
+  end
+end
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 1
+    set Y extent = 1
+  end
+end
+
+subsection Discretization
+  set Composition polynomial degree = 2,2,3
+  set Use discontinuous composition discretization = false, true, true
+end
+
+subsection Compositional fields
+  set Number of fields = 3
+  set Names of fields  = prescribed_1, prescribed_2, prescribed_3
+  set Compositional field methods = prescribed field, prescribed field, prescribed field
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.1
+  end
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.0; 0.2; 0.3
+  end
+end
+
+subsection Material model
+  set Model name = prescribed fields material
+
+  subsection Simple model
+    set Reference density             = 1
+    set Reference specific heat       = 1
+    set Reference temperature         = 0
+    set Thermal conductivity          = 1
+    set Thermal expansion coefficient = 1
+    set Viscosity                     = 1
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 1
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = 0, 1, 2, 3
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization, composition statistics
+
+  subsection Visualization
+    set Interpolate output = false
+    set List of output variables      = named additional outputs
+    set Number of grouped files       = 0
+    set Output format                 = gnuplot
+    set Time between graphical output = 0
+  end
+end
diff --git a/tests/multiple_prescribed_fields_dc/screen-output b/tests/multiple_prescribed_fields_dc/screen-output
new file mode 100644
index 00000000000..b0e8a33b14d
--- /dev/null
+++ b/tests/multiple_prescribed_fields_dc/screen-output
@@ -0,0 +1,19 @@
+
+Loading shared library <./libmultiple_prescribed_fields_dc.debug.so>
+
+Number of active cells: 4 (on 2 levels)
+Number of degrees of freedom: 209 (50+9+25+25+36+64)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Solving temperature system... 0 iterations.
+   Copying properties into prescribed compositional fields... done.
+   Solving Stokes system (GMG)... 6+0 iterations.
+
+   Postprocessing:
+     Writing graphical output:  output-multiple_prescribed_fields_dc/solution/solution-00000
+     Compositions min/max/mass: 0/0/0 // 0/1/0.5 // 0/2/1
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/multiple_prescribed_fields_dc/solution/solution-00000.0000.gnuplot b/tests/multiple_prescribed_fields_dc/solution/solution-00000.0000.gnuplot
new file mode 100644
index 00000000000..49ed29ce85b
--- /dev/null
+++ b/tests/multiple_prescribed_fields_dc/solution/solution-00000.0000.gnuplot
@@ -0,0 +1,35 @@
+# This file was generated by the deal.II library.
+
+
+#
+# For a description of the GNUPLOT format see the GNUPLOT manual.
+#
+# <x> <y> <velocity> <velocity> <p> <T> <prescribed_1> <prescribed_2> <prescribed_3> <prescribed_field_output_C0> <prescribed_field_output_C1> <prescribed_field_output_C2> 
+0 0 0 0 0.9 0.1 0 0 0 0 0 0 
+0.5 0 -4.42931e-27 0 0.9 0.1 0 0 0 0 0 0 
+
+0 0.5 0 2.208e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+0.5 0.5 5.45949e-28 -3.03919e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+
+
+0.5 0 -4.42931e-27 0 0.9 0.1 0 0 0 0 0 0 
+1 0 0 0 0.9 0.1 0 0 0 0 0 0 
+
+0.5 0.5 5.45949e-28 -3.03919e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+1 0.5 0 2.208e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+
+
+0 0.5 0 2.208e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+0.5 0.5 5.45949e-28 -3.03919e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+
+0 1 0 0 -5.3257e-11 0.1 0 1 2 0 1 2 
+0.5 1 -5.39777e-27 0 5.3257e-11 0.1 0 1 2 0 1 2 
+
+
+0.5 0.5 5.45949e-28 -3.03919e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+1 0.5 0 2.208e-11 0.45 0.1 0 0.5 1 0 0.5 1 
+
+0.5 1 -5.39777e-27 0 5.3257e-11 0.1 0 1 2 0 1 2 
+1 1 0 0 -5.3257e-11 0.1 0 1 2 0 1 2 
+
+
diff --git a/tests/multiple_prescribed_fields_dc/statistics b/tests/multiple_prescribed_fields_dc/statistics
new file mode 100644
index 00000000000..bb50c08ae8f
--- /dev/null
+++ b/tests/multiple_prescribed_fields_dc/statistics
@@ -0,0 +1,25 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Time step size (seconds)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for composition solver 1
+# 10: Iterations for composition solver 2
+# 11: Iterations for composition solver 3
+# 12: Iterations for Stokes solver
+# 13: Velocity iterations in Stokes preconditioner
+# 14: Schur complement iterations in Stokes preconditioner
+# 15: Visualization file name
+# 16: Minimal value for composition prescribed_1
+# 17: Maximal value for composition prescribed_1
+# 18: Global mass for composition prescribed_1
+# 19: Minimal value for composition prescribed_2
+# 20: Maximal value for composition prescribed_2
+# 21: Global mass for composition prescribed_2
+# 22: Minimal value for composition prescribed_3
+# 23: Maximal value for composition prescribed_3
+# 24: Global mass for composition prescribed_3
+0 0.000000000000e+00 0.000000000000e+00 4 59 25 75 0 4294967295 4294967295 4294967295 5 7 7 output-multiple_prescribed_fields_dc/solution/solution-00000 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 5.00000000e-01 0.00000000e+00 2.00000000e+00 1.00000000e+00 
diff --git a/tests/steinberger_projected_density/screen-output b/tests/steinberger_projected_density/screen-output
index 6bb5c2b8853..7ecc08644e9 100644
--- a/tests/steinberger_projected_density/screen-output
+++ b/tests/steinberger_projected_density/screen-output
@@ -5,8 +5,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 *** Timestep 0:  t=0 years, dt=0 years
    Solving temperature system... 0 iterations.
    Solving composition system ... 0 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 0 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 15+0 iterations.
 
    Postprocessing:
@@ -19,8 +19,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 *** Timestep 1:  t=100000 years, dt=100000 years
    Solving temperature system... 13 iterations.
    Solving composition system ... 11 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 11 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 13+0 iterations.
 
    Postprocessing:
diff --git a/tests/steinberger_projected_density/statistics b/tests/steinberger_projected_density/statistics
index 35dc81177cc..677d83d8110 100644
--- a/tests/steinberger_projected_density/statistics
+++ b/tests/steinberger_projected_density/statistics
@@ -7,8 +7,8 @@
 # 7: Number of degrees of freedom for all compositions
 # 8: Iterations for temperature solver
 # 9: Iterations for composition solver 1
-# 10: Iterations for composition solver 2
-# 11: Iterations for composition solver 3
+# 10: Iterations for composition solver 3
+# 11: Iterations for composition solver 2
 # 12: Iterations for Stokes solver
 # 13: Velocity iterations in Stokes preconditioner
 # 14: Schur complement iterations in Stokes preconditioner
@@ -23,5 +23,5 @@
 # 23: Outward heat flux through boundary with indicator 1 ("top") (W)
 # 24: Outward heat flux through boundary with indicator 2 ("left") (W)
 # 25: Outward heat flux through boundary with indicator 3 ("right") (W)
-0 0.000000000000e+00 0.000000000000e+00 192 1891 833 2499  0  0 4294967295  0 14 16 16 output-steinberger_projected_density/solution/solution-00000 5.81342027e-01 1.11483000e+00 2.73000000e+02 2.43678322e+03 4.25000000e+03 5.44074231e-01  4.85706828e+05 -4.25073533e+06 0.00000000e+00 0.00000000e+00 
-1 1.000000000000e+05 1.000000000000e+05 192 1891 833 2499 13 11 4294967295 11 12 14 14 output-steinberger_projected_density/solution/solution-00001 5.78335271e-01 1.11468431e+00 2.73000000e+02 2.43540495e+03 4.25000000e+03 5.43727671e-01 -2.81356102e+06 -1.21533631e+06 0.00000000e+00 0.00000000e+00 
+0 0.000000000000e+00 0.000000000000e+00 192 1891 833 2499  0  0  0 4294967295 14 16 16 output-steinberger_projected_density/solution/solution-00000 5.81342027e-01 1.11483000e+00 2.73000000e+02 2.43678322e+03 4.25000000e+03 5.44074231e-01  4.85706828e+05 -4.25073533e+06 0.00000000e+00 0.00000000e+00 
+1 1.000000000000e+05 1.000000000000e+05 192 1891 833 2499 13 11 11 4294967295 12 14 14 output-steinberger_projected_density/solution/solution-00001 5.78335271e-01 1.11468431e+00 2.73000000e+02 2.43540495e+03 4.25000000e+03 5.43727671e-01 -2.81356102e+06 -1.21533631e+06 0.00000000e+00 0.00000000e+00 
diff --git a/tests/steinberger_projected_density_dcPicard/screen-output b/tests/steinberger_projected_density_dcPicard/screen-output
index f485c874bf1..b0d33000cef 100644
--- a/tests/steinberger_projected_density_dcPicard/screen-output
+++ b/tests/steinberger_projected_density_dcPicard/screen-output
@@ -5,8 +5,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 *** Timestep 0:  t=0 years, dt=0 years
    Solving temperature system... 0 iterations.
    Solving composition system ... 0 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 0 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Initial Newton Stokes residual = 1.22791e+14, v = 1.22791e+14, p = 0
 
    Solving Stokes system (GMG)... 15+0 iterations.
@@ -16,8 +16,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 
    Solving temperature system... 0 iterations.
    Solving composition system ... 0 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 0 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 11+0 iterations.
       Newton system information: Norm of the rhs: 2.23856e+12, Derivative scaling factor: 0
       Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.71633e-16, 1.61486e-16, 0, 1.56133e-16, 0.0182306
@@ -25,8 +25,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 
    Solving temperature system... 0 iterations.
    Solving composition system ... 0 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 0 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 11+0 iterations.
       Newton system information: Norm of the rhs: 6.56739e+10, Derivative scaling factor: 0
       Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.71633e-16, 1.61486e-16, 0, 1.56133e-16, 0.000534843
@@ -34,8 +34,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 
    Solving temperature system... 0 iterations.
    Solving composition system ... 0 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 0 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 11+0 iterations.
       Newton system information: Norm of the rhs: 1.62309e+09, Derivative scaling factor: 0
       Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.71633e-16, 1.61486e-16, 0, 1.56133e-16, 1.32183e-05
@@ -43,8 +43,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 
    Solving temperature system... 0 iterations.
    Solving composition system ... 0 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 0 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 11+0 iterations.
       Newton system information: Norm of the rhs: 4.47762e+07, Derivative scaling factor: 0
       Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.71633e-16, 1.61486e-16, 0, 1.56133e-16, 3.64653e-07
@@ -61,8 +61,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 *** Timestep 1:  t=100000 years, dt=100000 years
    Solving temperature system... 13 iterations.
    Solving composition system ... 11 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 11 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Initial Newton Stokes residual = 8.39912e+13, v = 8.31745e+13, p = 1.16846e+13
 
    Solving Stokes system (GMG)... 13+0 iterations.
@@ -72,8 +72,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 
    Solving temperature system... 11 iterations.
    Solving composition system ... 9 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 9 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 13+0 iterations.
       Newton system information: Norm of the rhs: 2.55278e+11, Derivative scaling factor: 0
       Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0.00012418, 0.000414973, 0, 0.000193229, 0.00303934
@@ -81,8 +81,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 
    Solving temperature system... 10 iterations.
    Solving composition system ... 8 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 8 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 13+0 iterations.
       Newton system information: Norm of the rhs: 2.2233e+10, Derivative scaling factor: 0
       Relative nonlinear residuals (temperature, compositional fields, Stokes system): 6.8685e-06, 3.11445e-05, 0, 1.43921e-05, 0.000264707
@@ -90,8 +90,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 
    Solving temperature system... 8 iterations.
    Solving composition system ... 7 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 7 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 12+0 iterations.
       Newton system information: Norm of the rhs: 1.00197e+09, Derivative scaling factor: 0
       Relative nonlinear residuals (temperature, compositional fields, Stokes system): 6.11361e-07, 2.60566e-06, 0, 1.38988e-06, 1.19295e-05
@@ -99,8 +99,8 @@ Number of degrees of freedom: 5,223 (1,666+225+833+833+833+833)
 
    Solving temperature system... 7 iterations.
    Solving composition system ... 6 iterations.
-   Copying properties into prescribed compositional field density_field... done.
    Solving another_field system ... 5 iterations.
+   Copying properties into prescribed compositional field density_field... done.
    Solving Stokes system (GMG)... 13+0 iterations.
       Newton system information: Norm of the rhs: 9.33525e+07, Derivative scaling factor: 0
       Relative nonlinear residuals (temperature, compositional fields, Stokes system): 4.31558e-08, 1.25138e-07, 0, 8.13074e-08, 1.11146e-06
diff --git a/tests/steinberger_projected_density_dcPicard/statistics b/tests/steinberger_projected_density_dcPicard/statistics
index 47b5ef6c887..93ef215908c 100644
--- a/tests/steinberger_projected_density_dcPicard/statistics
+++ b/tests/steinberger_projected_density_dcPicard/statistics
@@ -8,8 +8,8 @@
 # 8: Number of nonlinear iterations
 # 9: Iterations for temperature solver
 # 10: Iterations for composition solver 1
-# 11: Iterations for composition solver 2
-# 12: Iterations for composition solver 3
+# 11: Iterations for composition solver 3
+# 12: Iterations for composition solver 2
 # 13: Iterations for Stokes solver
 # 14: Velocity iterations in Stokes preconditioner
 # 15: Schur complement iterations in Stokes preconditioner
@@ -24,5 +24,5 @@
 # 24: Outward heat flux through boundary with indicator 1 ("top") (W)
 # 25: Outward heat flux through boundary with indicator 2 ("left") (W)
 # 26: Outward heat flux through boundary with indicator 3 ("right") (W)
-0 0.000000000000e+00 0.000000000000e+00 192 1891 833 2499 5  0  0 4294967291  0 54 64 64 output-steinberger_projected_density_dcPicard/solution/solution-00000 5.78314979e-01 1.09945267e+00 2.73000000e+02 2.43678322e+03 4.25000000e+03 5.44074231e-01  4.99327014e+05 -4.27418614e+06 0.00000000e+00 0.00000000e+00 
-1 1.000000000000e+05 1.000000000000e+05 192 1891 833 2499 5 49 41 4294967291 40 59 69 69 output-steinberger_projected_density_dcPicard/solution/solution-00001 5.77687184e-01 1.11532571e+00 2.73000000e+02 2.43583623e+03 4.25000000e+03 5.43836114e-01 -2.71821909e+06 -1.13647007e+06 0.00000000e+00 0.00000000e+00 
+0 0.000000000000e+00 0.000000000000e+00 192 1891 833 2499 5  0  0  0 4294967291 54 64 64 output-steinberger_projected_density_dcPicard/solution/solution-00000 5.78314979e-01 1.09945267e+00 2.73000000e+02 2.43678322e+03 4.25000000e+03 5.44074231e-01  4.99327014e+05 -4.27418614e+06 0.00000000e+00 0.00000000e+00 
+1 1.000000000000e+05 1.000000000000e+05 192 1891 833 2499 5 49 41 40 4294967291 59 69 69 output-steinberger_projected_density_dcPicard/solution/solution-00001 5.77687184e-01 1.11532571e+00 2.73000000e+02 2.43583623e+03 4.25000000e+03 5.43836114e-01 -2.71821909e+06 -1.13647007e+06 0.00000000e+00 0.00000000e+00