diff --git a/doc/modules/changes/20240618_tjhei b/doc/modules/changes/20240618_tjhei
new file mode 100644
index 00000000000..15ff52c82de
--- /dev/null
+++ b/doc/modules/changes/20240618_tjhei
@@ -0,0 +1,5 @@
+New: ASPECT now supports compositional fields with
+different discretizations (continuous or discontinuous)
+at the same time.
+<br>
+(Timo Heister, 2024/06/18)
diff --git a/include/aspect/fe_variable_collection.h b/include/aspect/fe_variable_collection.h
index c0acba8189a..f98de32db12 100644
--- a/include/aspect/fe_variable_collection.h
+++ b/include/aspect/fe_variable_collection.h
@@ -194,10 +194,17 @@ namespace aspect
 
       /**
        * Return the variable with name @p name. Throws an exception if this
-       * variable does not exist.
+       * variable does not exist. If more than one variable with the same name
+       * exists, return the first one. Use variables_with_name() if you want
+       * to access all of them.
        */
       const FEVariable<dim> &variable(const std::string &name) const;
 
+      /**
+       * Return a vector of pointers of all variables with name @p name.
+       */
+      std::vector<const FEVariable<dim>*> variables_with_name(const std::string &name) const;
+
       /**
        * Returns true if the variable with @p name exists in the list of
        * variables.
diff --git a/include/aspect/introspection.h b/include/aspect/introspection.h
index f8f4e47789b..9ce0167be3a 100644
--- a/include/aspect/introspection.h
+++ b/include/aspect/introspection.h
@@ -350,7 +350,7 @@ namespace aspect
        */
       struct ComponentMasks
       {
-        ComponentMasks (FEVariableCollection<dim> &fevs);
+        ComponentMasks (const FEVariableCollection<dim> &fevs, const Introspection<dim>::ComponentIndices &indices);
 
         ComponentMask              velocities;
         ComponentMask              pressure;
diff --git a/source/fe_variable_collection.cc b/source/fe_variable_collection.cc
index 7e6750c2f7b..ef8b9030d77 100644
--- a/source/fe_variable_collection.cc
+++ b/source/fe_variable_collection.cc
@@ -117,9 +117,6 @@ namespace aspect
     variables.clear();
     variables.reserve(variable_definitions.size());
 
-    // store names temporarily to make sure they are unique
-    std::set<std::string> names;
-
     unsigned int component_index = 0;
     unsigned int block_index = 0;
 
@@ -129,10 +126,6 @@ namespace aspect
                                             component_index, block_index, i));
         component_index+= variables[i].n_components();
         block_index += variables[i].n_blocks;
-
-        Assert(names.find(variables[i].name)==names.end(),
-               ExcMessage("Can not add two variables with the same name."));
-        names.insert(variables[i].name);
       }
 
     Assert(variables.back().n_blocks != 0
@@ -182,6 +175,20 @@ namespace aspect
 
 
 
+  template <int dim>
+  std::vector<const FEVariable<dim>*>
+  FEVariableCollection<dim>::variables_with_name(const std::string &name) const
+  {
+    std::vector<const FEVariable<dim>*> result;
+    for (unsigned int i=0; i<variables.size(); ++i)
+      if (variables[i].name == name)
+        result.emplace_back(&variables[i]);
+
+    return result;
+  }
+
+
+
   template <int dim>
   bool
   FEVariableCollection<dim>::variable_exists(const std::string &name) const
diff --git a/source/particle/world.cc b/source/particle/world.cc
index 965901f3df5..06d37152878 100644
--- a/source/particle/world.cc
+++ b/source/particle/world.cc
@@ -79,6 +79,9 @@ namespace aspect
                                          property_manager->get_n_property_components());
 
       connect_to_signals(this->get_signals());
+
+      AssertThrow(this->introspection().get_composition_base_element_indices().size()<=1,
+                  ExcNotImplemented("Particles are not supported in computations with compositional fields with different finite element types."));
     }
 
 
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 138f456293a..aaa31fbd920 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -880,6 +880,30 @@ namespace aspect
 
 
 
+    template <int dim>
+    bool compositional_field_needs_matrix_block(const Introspection<dim> &introspection, const unsigned int composition_index)
+    {
+      const typename Simulator<dim>::AdvectionField adv_field (Simulator<dim>::AdvectionField::composition(composition_index));
+      switch (adv_field.advection_method(introspection))
+        {
+          case Parameters<dim>::AdvectionFieldMethod::fem_field:
+          case Parameters<dim>::AdvectionFieldMethod::fem_melt_field:
+          case Parameters<dim>::AdvectionFieldMethod::fem_darcy_field:
+          case Parameters<dim>::AdvectionFieldMethod::prescribed_field_with_diffusion:
+            return true;
+          case Parameters<dim>::AdvectionFieldMethod::particles:
+          case Parameters<dim>::AdvectionFieldMethod::volume_of_fluid:
+          case Parameters<dim>::AdvectionFieldMethod::static_field:
+          case Parameters<dim>::AdvectionFieldMethod::prescribed_field:
+            break;
+          default:
+            Assert (false, ExcNotImplemented());
+        }
+      return false;
+    }
+
+
+
     template <int dim>
     bool compositional_fields_need_matrix_block(const Introspection<dim> &introspection)
     {
@@ -887,22 +911,8 @@ namespace aspect
       // (as opposed to all are advected by other means or prescribed fields)
       for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
         {
-          const typename Simulator<dim>::AdvectionField adv_field (Simulator<dim>::AdvectionField::composition(c));
-          switch (adv_field.advection_method(introspection))
-            {
-              case Parameters<dim>::AdvectionFieldMethod::fem_field:
-              case Parameters<dim>::AdvectionFieldMethod::fem_melt_field:
-              case Parameters<dim>::AdvectionFieldMethod::fem_darcy_field:
-              case Parameters<dim>::AdvectionFieldMethod::prescribed_field_with_diffusion:
-                return true;
-              case Parameters<dim>::AdvectionFieldMethod::particles:
-              case Parameters<dim>::AdvectionFieldMethod::volume_of_fluid:
-              case Parameters<dim>::AdvectionFieldMethod::static_field:
-              case Parameters<dim>::AdvectionFieldMethod::prescribed_field:
-                break;
-              default:
-                Assert (false, ExcNotImplemented());
-            }
+          if (compositional_field_needs_matrix_block(introspection, c))
+            return true;
         }
       return false;
     }
@@ -1010,11 +1020,21 @@ namespace aspect
         &&
         compositional_fields_need_matrix_block(introspection))
       {
-        // If we need at least one compositional field block, we
-        // create a matrix block in the first compositional block. Its sparsity
-        // pattern will later be used to allocate composition matrices as
-        // needed. All other matrix blocks are left empty to save memory.
-        coupling[x.compositional_fields[0]][x.compositional_fields[0]] = DoFTools::always;
+        // We reuse matrix blocks for compositional fields, but we need different blocks for each base_element.
+        // All other matrix blocks are left empty to save memory.
+        for (const unsigned int base_element_index : introspection.get_composition_base_element_indices())
+          {
+            bool block_needed = false;
+            for (const unsigned int c : introspection.get_compositional_field_indices_with_base_element(base_element_index))
+              if (compositional_field_needs_matrix_block(introspection, c))
+                block_needed = true;
+
+            if (block_needed)
+              {
+                const unsigned int first_c = introspection.get_compositional_field_indices_with_base_element(base_element_index).front();
+                coupling[x.compositional_fields[first_c]][x.compositional_fields[first_c]] = DoFTools::always;
+              }
+          }
       }
 
     // If we are using volume of fluid interface tracking, create a matrix block in the
@@ -1062,10 +1082,20 @@ namespace aspect
             parameters.temperature_method != Parameters<dim>::AdvectionFieldMethod::static_field)
           face_coupling[x.temperature][x.temperature] = DoFTools::always;
 
-        if (parameters.have_discontinuous_composition_discretization &&
-            solver_scheme_solves_advection_equations(parameters) &&
-            compositional_fields_need_matrix_block(introspection))
-          face_coupling[x.compositional_fields[0]][x.compositional_fields[0]] = DoFTools::always;
+        for (const unsigned int base_element_index : introspection.get_composition_base_element_indices())
+          {
+            bool block_needed = false;
+            for (const unsigned int c : introspection.get_compositional_field_indices_with_base_element(base_element_index))
+              if (compositional_field_needs_matrix_block(introspection, c))
+                block_needed = true;
+
+            const unsigned int first_c = introspection.get_compositional_field_indices_with_base_element(base_element_index).front();
+
+            if (parameters.use_discontinuous_composition_discretization[first_c]
+                && solver_scheme_solves_advection_equations(parameters)
+                && block_needed)
+              face_coupling[x.compositional_fields[first_c]][x.compositional_fields[first_c]] = DoFTools::always;
+          }
 
         if (parameters.volume_of_fluid_tracking_enabled)
           {
@@ -1082,9 +1112,7 @@ namespace aspect
                                               Utilities::MPI::
                                               this_mpi_process(mpi_communicator));
 
-        if (solver_scheme_solves_advection_equations(parameters)
-            &&
-            compositional_fields_need_matrix_block(introspection))
+        if (solver_scheme_solves_advection_equations(parameters))
           {
             // If we solve for more than one compositional field make sure we keep constrained entries
             // to allow different boundary conditions for different fields. In order to keep constrained
@@ -1094,20 +1122,24 @@ namespace aspect
                                                              introspection.n_components);
             composition_coupling.fill (DoFTools::none);
 
-            const unsigned int component = introspection.component_indices.compositional_fields[0];
-            composition_coupling[component][component] = coupling[component][component];
-
-            const unsigned int block = introspection.get_components_to_blocks()[component];
-            sp.block(block,block).reinit(sp.block(block,block).locally_owned_range_indices(),
-                                         sp.block(block,block).locally_owned_domain_indices());
-
-            DoFTools::make_flux_sparsity_pattern (dof_handler,
-                                                  sp,
-                                                  current_constraints, true,
-                                                  composition_coupling,
-                                                  face_coupling,
-                                                  Utilities::MPI::
-                                                  this_mpi_process(mpi_communicator));
+            for (const unsigned int base_element_index : introspection.get_composition_base_element_indices())
+              {
+                const unsigned int first_c = introspection.get_compositional_field_indices_with_base_element(base_element_index).front();
+                const unsigned int component = x.compositional_fields[first_c];
+                composition_coupling[component][component] = coupling[component][component];
+
+                const unsigned int block = introspection.block_indices.compositional_field_sparsity_pattern[first_c];
+                sp.block(block,block).reinit(sp.block(block,block).locally_owned_range_indices(),
+                                             sp.block(block,block).locally_owned_domain_indices());
+
+                DoFTools::make_flux_sparsity_pattern (dof_handler,
+                                                      sp,
+                                                      current_constraints, true,
+                                                      composition_coupling,
+                                                      face_coupling,
+                                                      Utilities::MPI::
+                                                      this_mpi_process(mpi_communicator));
+              }
           }
       }
     else
diff --git a/source/simulator/introspection.cc b/source/simulator/introspection.cc
index 7a1963fe230..c52c4e91181 100644
--- a/source/simulator/introspection.cc
+++ b/source/simulator/introspection.cc
@@ -44,9 +44,17 @@ namespace aspect
 
       ci.pressure = fevs.variable("pressure").first_component_index;
       ci.temperature = fevs.variable("temperature").first_component_index;
-      unsigned int n_compositional_fields = fevs.variable("compositions").n_components();
-      for (unsigned int i=0; i<n_compositional_fields; ++i)
-        ci.compositional_fields.push_back(fevs.variable("compositions").first_component_index+i);
+
+      {
+        const auto composition_variables = fevs.variables_with_name("compositions");
+        for (const auto &fe : composition_variables)
+          {
+            for (unsigned int i=0; i<fe->multiplicity; ++i)
+              {
+                ci.compositional_fields.push_back(fe->first_component_index+i);
+              }
+          }
+      }
 
       return ci;
     }
@@ -56,20 +64,24 @@ namespace aspect
      */
     template <int dim>
     typename Introspection<dim>::BlockIndices
-    setup_blocks (FEVariableCollection<dim> &fevs, const unsigned int n_compositional_fields)
+    setup_blocks (FEVariableCollection<dim> &fevs)
     {
       typename Introspection<dim>::BlockIndices b;
       b.velocities = fevs.variable("velocity").block_index;
       b.pressure = fevs.variable("pressure").block_index;
       b.temperature = fevs.variable("temperature").block_index;
 
-      const unsigned int first_composition_block_index = fevs.variable("compositions").block_index;
-      for (unsigned int i=0; i<n_compositional_fields; ++i)
-        {
-          b.compositional_fields.push_back(fevs.variable("compositions").block_index+i);
-          b.compositional_field_sparsity_pattern.push_back(first_composition_block_index);
-        }
-
+      {
+        const auto variables = fevs.variables_with_name("compositions");
+        for (const auto &fe : variables)
+          {
+            for (unsigned int i=0; i<fe->multiplicity; ++i)
+              {
+                b.compositional_fields.push_back(fe->block_index + i);
+                b.compositional_field_sparsity_pattern.push_back(fe->block_index);
+              }
+          }
+      }
       return b;
     }
 
@@ -77,15 +89,27 @@ namespace aspect
 
     template <int dim>
     typename Introspection<dim>::BaseElements
-    setup_base_elements (FEVariableCollection<dim> &fevs, const unsigned int n_compositional_fields)
+    setup_base_elements (FEVariableCollection<dim> &fevs)
     {
       typename Introspection<dim>::BaseElements base_elements;
 
       base_elements.velocities = fevs.variable("velocity").base_index;
       base_elements.pressure = fevs.variable("pressure").base_index;
       base_elements.temperature = fevs.variable("temperature").base_index;
-      base_elements.compositional_fields = Utilities::possibly_extend_from_1_to_N(std::vector<unsigned int>({fevs.variable("compositions").base_index}),
-                                                                                  n_compositional_fields, "");
+
+      {
+        // TODO: We would ideally reuse base elements also when they are
+        // not consecutively encountered. For now, just ignore that fact.
+        const auto variables = fevs.variables_with_name("compositions");
+        for (const auto &fe : variables)
+          {
+            for (unsigned int i=0; i<fe->multiplicity; ++i)
+              base_elements.compositional_fields.push_back(fe->base_index);
+          }
+      }
+
+
+
 
       return base_elements;
     }
@@ -223,14 +247,25 @@ namespace aspect
         1,
         1));
 
-
-    variables.push_back(
-      VariableDeclaration<dim>(
-        "compositions",
-        internal::new_FE_Q_or_DGQ<dim>(parameters.have_discontinuous_composition_discretization, // TODO: this is of course incorrect for now.
-                                       parameters.composition_degree),
-        parameters.n_compositional_fields,
-        parameters.n_compositional_fields));
+    // We can not create a single composition with multiplicity n (this assumes all have the same FiniteElement) and
+    // we also don't want to create n different compositional fields for performance reasons. Instead, we combine
+    // consecutive compositions of the same type into the same base element. For example, if the user requests
+    // "Q2,Q2,DGQ1,Q2", we actually create "Q2^2, DGQ1, Q2":
+    for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+      {
+        if (c>0 && parameters.use_discontinuous_composition_discretization[c] == parameters.use_discontinuous_composition_discretization[c-1])
+          {
+            // reuse last one because it is the same
+            variables.back().multiplicity += 1;
+            variables.back().n_blocks += 1;
+          }
+        else
+          {
+            std::shared_ptr<FiniteElement<dim>> fe = internal::new_FE_Q_or_DGQ<dim>(parameters.use_discontinuous_composition_discretization[c],
+                                                                                     parameters.composition_degree);
+            variables.push_back(VariableDeclaration<dim>("compositions", fe, 1, 1));
+          }
+      }
 
     return variables;
   }
@@ -247,14 +282,14 @@ namespace aspect
     use_discontinuous_temperature_discretization (parameters.use_discontinuous_temperature_discretization),
     use_discontinuous_composition_discretization (parameters.use_discontinuous_composition_discretization),
     component_indices (internal::setup_component_indices<dim>(*this)),
-    n_blocks(FEVariableCollection<dim>::n_blocks()),
-    block_indices (internal::setup_blocks<dim>(*this, parameters.n_compositional_fields)),
+    n_blocks (FEVariableCollection<dim>::n_blocks()),
+    block_indices (internal::setup_blocks<dim>(*this)),
     extractors (component_indices),
-    base_elements (internal::setup_base_elements<dim>(*this, parameters.n_compositional_fields)),
+    base_elements (internal::setup_base_elements<dim>(*this)),
     polynomial_degree (internal::setup_polynomial_degree<dim>(parameters)),
     quadratures (internal::setup_quadratures<dim>(parameters, ReferenceCells::get_hypercube<dim>())),
     face_quadratures (internal::setup_face_quadratures<dim>(parameters, ReferenceCells::get_hypercube<dim>())),
-    component_masks (*this),
+    component_masks (*this, component_indices),
     system_dofs_per_block (n_blocks),
     temperature_method(parameters.temperature_method),
     compositional_field_methods(parameters.compositional_field_methods),
@@ -345,13 +380,13 @@ namespace aspect
   {
     template <int dim>
     std::vector<ComponentMask>
-    make_component_mask_sequence(const FEVariable<dim> &variable)
+    make_composition_component_mask_sequence (const FEVariableCollection<dim> &fevs, const typename Introspection<dim>::ComponentIndices &indices)
     {
       std::vector<ComponentMask> result;
-      for (unsigned int i=0; i<variable.multiplicity; ++i)
+      for (const unsigned int idx : indices.compositional_fields)
         {
-          result.push_back(ComponentMask(variable.component_mask.size(), false));
-          result.back().set(variable.first_component_index+i, true);
+          result.push_back(ComponentMask(fevs.n_components(), false));
+          result.back().set(idx, true);
         }
       return result;
     }
@@ -360,12 +395,12 @@ namespace aspect
 
 
   template <int dim>
-  Introspection<dim>::ComponentMasks::ComponentMasks (FEVariableCollection<dim> &fevs)
+  Introspection<dim>::ComponentMasks::ComponentMasks (const FEVariableCollection<dim> &fevs, const Introspection<dim>::ComponentIndices &indices)
     :
     velocities (fevs.variable("velocity").component_mask),
     pressure (fevs.variable("pressure").component_mask),
     temperature (fevs.variable("temperature").component_mask),
-    compositional_fields (make_component_mask_sequence (fevs.variable("compositions")))
+    compositional_fields (make_composition_component_mask_sequence (fevs, indices))
   {}
 
 
diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
index 73c4f56a4ef..630301c8fce 100644
--- a/source/simulator/parameters.cc
+++ b/source/simulator/parameters.cc
@@ -1038,7 +1038,7 @@ namespace aspect
                          "between cells, and weak imposition of boundary terms for the temperature "
                          "field via the interior-penalty discontinuous Galerkin method.");
       prm.declare_entry ("Use discontinuous composition discretization", "false",
-                         Patterns::Bool (),
+                         Patterns::List(Patterns::Bool ()),
                          "Whether to use a composition discretization that is discontinuous "
                          "as opposed to continuous. This then requires the assembly of face terms "
                          "between cells, and weak imposition of boundary terms for the composition "
diff --git a/tests/composition_cg_dg_fem.prm b/tests/composition_cg_dg_fem.prm
new file mode 100644
index 00000000000..4959bbadef3
--- /dev/null
+++ b/tests/composition_cg_dg_fem.prm
@@ -0,0 +1,96 @@
+#
+#
+# copied from discontinuous_composition_1.prm
+
+
+set Dimension                              = 2
+set Start time                             = 0
+set End time                               = 1
+set Use years in output instead of seconds = false
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 2
+    set Y extent = 1
+  end
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators   = bottom, top
+  set List of model names = box
+
+  subsection Box
+    set Bottom temperature = 1
+    set Top temperature    = 0
+  end
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right, bottom
+  set Prescribed velocity boundary indicators = top: function
+
+  subsection Function
+    set Variable names      = x,z,t
+    set Function constants  = pi=3.1415926
+    set Function expression = if(x>1+sin(0.5*pi*t), 1, -1); 0
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,z
+    set Function expression = (1-z)
+  end
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Thermal conductivity          = 1e-6
+    set Thermal expansion coefficient = 1e-4
+    set Viscosity                     = 1
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 1
+  set Initial global refinement          = 2
+  set Time steps between mesh refinement = 2
+end
+
+subsection Discretization
+  set Use discontinuous temperature discretization = false
+  set Use discontinuous composition discretization = true,false
+end
+
+subsection Postprocess
+  set List of postprocessors = temperature statistics, composition statistics
+end
+
+# This is the new part: We declare that there will
+# be two compositional fields that will be
+# advected along. Their initial conditions are given by
+# a function that is one for the lowermost 0.2 height
+# units of the domain and zero otherwise in the first case,
+# and one in the top most 0.2 height units in the latter.
+subsection Compositional fields
+  set Number of fields = 2
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,y
+    set Function expression = if(y<0.2, 1, 0) ; if(y>0.8, 1, 0)
+  end
+end
diff --git a/tests/composition_cg_dg_fem/screen-output b/tests/composition_cg_dg_fem/screen-output
new file mode 100644
index 00000000000..7e6c5bdf812
--- /dev/null
+++ b/tests/composition_cg_dg_fem/screen-output
@@ -0,0 +1,204 @@
+
+Number of active cells: 16 (on 3 levels)
+Number of degrees of freedom: 493 (162+25+81+144+81)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Solving temperature system... 0 iterations.
+   Solving C_1 system ... 0 iterations.
+   Solving C_2 system ... 0 iterations.
+   Solving Stokes system... 10+0 iterations.
+
+Number of active cells: 40 (on 4 levels)
+Number of degrees of freedom: 1,187 (386+55+193+360+193)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Solving temperature system... 0 iterations.
+   Solving C_1 system ... 0 iterations.
+   Solving C_2 system ... 0 iterations.
+   Solving Stokes system... 13+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   0 K, 0.5 K, 1 K
+     Compositions min/max/mass: 0/1/0.4167 // 0/1/0.4583
+
+*** Timestep 1:  t=0.0625 seconds, dt=0.0625 seconds
+   Solving temperature system... 10 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 10 iterations.
+   Solving Stokes system... 8+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   0 K, 0.5 K, 1 K
+     Compositions min/max/mass: -0.1081/1.021/0.4167 // -0.0158/1.028/0.4583
+
+*** Timestep 2:  t=0.125 seconds, dt=0.0625 seconds
+   Solving temperature system... 11 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 10 iterations.
+   Solving Stokes system... 13+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   0 K, 0.5 K, 1 K
+     Compositions min/max/mass: -0.1914/1.043/0.4167 // -0.02143/1.03/0.4583
+
+Number of active cells: 64 (on 4 levels)
+Number of degrees of freedom: 1,813 (578+81+289+576+289)
+
+*** Timestep 3:  t=0.1875 seconds, dt=0.0625 seconds
+   Solving temperature system... 12 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 12 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   0 K, 0.5 K, 1 K
+     Compositions min/max/mass: -0.1976/1.126/0.4167 // -0.01355/1.028/0.4582
+
+*** Timestep 4:  t=0.25 seconds, dt=0.0625 seconds
+   Solving temperature system... 12 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.004362 K, 0.5 K, 1 K
+     Compositions min/max/mass: -0.1982/1.126/0.4167 // -0.004788/1.019/0.4583
+
+Skipping mesh refinement, because the mesh did not change.
+
+*** Timestep 5:  t=0.3125 seconds, dt=0.0625 seconds
+   Solving temperature system... 12 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 12 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.007754 K, 0.5 K, 1 K
+     Compositions min/max/mass: -0.1702/1.125/0.4167 // -0.004272/1.021/0.4582
+
+*** Timestep 6:  t=0.375 seconds, dt=0.0625 seconds
+   Solving temperature system... 11 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 11+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.01149 K, 0.5 K, 1 K
+     Compositions min/max/mass: -0.1629/1.125/0.4167 // -0.006268/1.014/0.4581
+
+Skipping mesh refinement, because the mesh did not change.
+
+*** Timestep 7:  t=0.4375 seconds, dt=0.0625 seconds
+   Solving temperature system... 12 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.01532 K, 0.4999 K, 1 K
+     Compositions min/max/mass: -0.1349/1.125/0.4167 // -0.009569/1.009/0.4581
+
+*** Timestep 8:  t=0.5 seconds, dt=0.0625 seconds
+   Solving temperature system... 12 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 12 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.0157 K, 0.4999 K, 1 K
+     Compositions min/max/mass: -0.119/1.125/0.4167 // -0.007126/1.013/0.4581
+
+Skipping mesh refinement, because the mesh did not change.
+
+*** Timestep 9:  t=0.5625 seconds, dt=0.0625 seconds
+   Solving temperature system... 11 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 11+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.01422 K, 0.4998 K, 1 K
+     Compositions min/max/mass: -0.1125/1.125/0.4167 // -0.003882/1.012/0.458
+
+*** Timestep 10:  t=0.625 seconds, dt=0.0625 seconds
+   Solving temperature system... 12 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 11+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.01272 K, 0.4997 K, 1 K
+     Compositions min/max/mass: -0.1223/1.125/0.4167 // -0.004574/1.01/0.458
+
+Skipping mesh refinement, because the mesh did not change.
+
+*** Timestep 11:  t=0.6875 seconds, dt=0.0625 seconds
+   Solving temperature system... 11 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.01329 K, 0.4996 K, 1 K
+     Compositions min/max/mass: -0.1302/1.125/0.4167 // -0.00472/1.01/0.4581
+
+*** Timestep 12:  t=0.75 seconds, dt=0.0625 seconds
+   Solving temperature system... 13 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 12 iterations.
+   Solving Stokes system... 12+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.0123 K, 0.4996 K, 1 K
+     Compositions min/max/mass: -0.1197/1.124/0.4167 // -0.00434/1.012/0.4581
+
+Skipping mesh refinement, because the mesh did not change.
+
+*** Timestep 13:  t=0.8125 seconds, dt=0.0625 seconds
+   Solving temperature system... 11 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 6+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.009901 K, 0.4996 K, 1 K
+     Compositions min/max/mass: -0.09665/1.125/0.4167 // -0.0038/1.01/0.4581
+
+*** Timestep 14:  t=0.875 seconds, dt=0.0625 seconds
+   Solving temperature system... 11 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 5+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.007248 K, 0.4996 K, 1 K
+     Compositions min/max/mass: -0.09368/1.125/0.4167 // -0.003175/1.007/0.4581
+
+Skipping mesh refinement, because the mesh did not change.
+
+*** Timestep 15:  t=0.9375 seconds, dt=0.0625 seconds
+   Solving temperature system... 11 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 5+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.006538 K, 0.4996 K, 1 K
+     Compositions min/max/mass: -0.1627/1.125/0.4167 // -0.002995/1.006/0.4581
+
+*** Timestep 16:  t=1 seconds, dt=0.0625 seconds
+   Solving temperature system... 12 iterations.
+   Solving C_1 system ... 3 iterations.
+   Solving C_2 system ... 11 iterations.
+   Solving Stokes system... 5+0 iterations.
+
+   Postprocessing:
+     Temperature min/avg/max:   -0.006603 K, 0.4996 K, 1 K
+     Compositions min/max/mass: -0.2114/1.149/0.4167 // -0.003116/1.006/0.4581
+
+Skipping mesh refinement, because the mesh did not change.
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/composition_cg_dg_fem/statistics b/tests/composition_cg_dg_fem/statistics
new file mode 100644
index 00000000000..806577a20a7
--- /dev/null
+++ b/tests/composition_cg_dg_fem/statistics
@@ -0,0 +1,40 @@
+# 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 Stokes solver
+# 12: Velocity iterations in Stokes preconditioner
+# 13: Schur complement iterations in Stokes preconditioner
+# 14: Minimal temperature (K)
+# 15: Average temperature (K)
+# 16: Maximal temperature (K)
+# 17: Average nondimensional temperature (K)
+# 18: Minimal value for composition C_1
+# 19: Maximal value for composition C_1
+# 20: Global mass for composition C_1
+# 21: Minimal value for composition C_2
+# 22: Maximal value for composition C_2
+# 23: Global mass for composition C_2
+ 0 0.000000000000e+00 0.000000000000e+00 40 441 193  720  0 0  0 12 14 14  0.00000000e+00 5.00000000e-01 1.00000000e+00 5.00000000e-01  0.00000000e+00 1.00000000e+00 4.16666667e-01  0.00000000e+00 1.00000000e+00 4.58333333e-01 
+ 1 6.250000000000e-02 6.250000000000e-02 40 441 193  720 10 3 10  7  9  9  0.00000000e+00 5.00021358e-01 1.00000000e+00 5.00021358e-01 -1.08142271e-01 1.02099034e+00 4.16676104e-01 -1.57952550e-02 1.02777854e+00 4.58314382e-01 
+ 2 1.250000000000e-01 6.250000000000e-02 40 441 193  720 11 3 10 12 14 14  0.00000000e+00 5.00023331e-01 1.00000000e+00 5.00023331e-01 -1.91381311e-01 1.04313258e+00 4.16691480e-01 -2.14326826e-02 1.03045963e+00 4.58276741e-01 
+ 3 1.875000000000e-01 6.250000000000e-02 64 659 289 1152 12 3 12 11 13 13  0.00000000e+00 5.00015279e-01 1.00000000e+00 5.00015279e-01 -1.97577121e-01 1.12631282e+00 4.16716449e-01 -1.35466918e-02 1.02789244e+00 4.58202929e-01 
+ 4 2.500000000000e-01 6.250000000000e-02 64 659 289 1152 12 3 11 11 13 13 -4.36194966e-03 4.99974296e-01 1.00000000e+00 5.04336246e-01 -1.98196716e-01 1.12603887e+00 4.16700097e-01 -4.78757254e-03 1.01907300e+00 4.58260855e-01 
+ 5 3.125000000000e-01 6.250000000000e-02 64 659 289 1152 12 3 12 11 13 13 -7.75404720e-03 4.99973164e-01 1.00000000e+00 5.07727212e-01 -1.70176772e-01 1.12528924e+00 4.16694932e-01 -4.27243422e-03 1.02127786e+00 4.58181536e-01 
+ 6 3.750000000000e-01 6.250000000000e-02 64 659 289 1152 11 3 11 10 12 12 -1.14918131e-02 4.99952303e-01 1.00000000e+00 5.11444116e-01 -1.62924178e-01 1.12481952e+00 4.16693908e-01 -6.26773689e-03 1.01441071e+00 4.58076795e-01 
+ 7 4.375000000000e-01 6.250000000000e-02 64 659 289 1152 12 3 11 11 13 13 -1.53184334e-02 4.99886877e-01 1.00000000e+00 5.15205310e-01 -1.34869926e-01 1.12500770e+00 4.16693633e-01 -9.56947843e-03 1.00940664e+00 4.58143056e-01 
+ 8 5.000000000000e-01 6.250000000000e-02 64 659 289 1152 12 3 12 11 13 13 -1.57000868e-02 4.99858770e-01 1.00000000e+00 5.15558856e-01 -1.19022547e-01 1.12485992e+00 4.16692755e-01 -7.12550187e-03 1.01343953e+00 4.58067091e-01 
+ 9 5.625000000000e-01 6.250000000000e-02 64 659 289 1152 11 3 11 10 12 12 -1.42223941e-02 4.99805947e-01 1.00000000e+00 5.14028341e-01 -1.12529822e-01 1.12461501e+00 4.16692325e-01 -3.88235609e-03 1.01243622e+00 4.57956650e-01 
+10 6.250000000000e-01 6.250000000000e-02 64 659 289 1152 12 3 11 10 12 12 -1.27182304e-02 4.99705542e-01 1.00000000e+00 5.12423772e-01 -1.22341396e-01 1.12477645e+00 4.16691691e-01 -4.57448318e-03 1.00956903e+00 4.58041114e-01 
+11 6.875000000000e-01 6.250000000000e-02 64 659 289 1152 11 3 11 11 13 13 -1.32874660e-02 4.99605244e-01 1.00000000e+00 5.12892710e-01 -1.30213764e-01 1.12466407e+00 4.16691406e-01 -4.71974132e-03 1.00967343e+00 4.58096559e-01 
+12 7.500000000000e-01 6.250000000000e-02 64 659 289 1152 13 3 12 11 13 13 -1.23038710e-02 4.99596773e-01 1.00000000e+00 5.11900644e-01 -1.19692682e-01 1.12447346e+00 4.16690576e-01 -4.34005789e-03 1.01208622e+00 4.58110921e-01 
+13 8.125000000000e-01 6.250000000000e-02 64 659 289 1152 11 3 11  5  7  7 -9.90118078e-03 4.99585050e-01 1.00000000e+00 5.09486231e-01 -9.66547894e-02 1.12487312e+00 4.16690023e-01 -3.79980720e-03 1.00983989e+00 4.58114652e-01 
+14 8.750000000000e-01 6.250000000000e-02 64 659 289 1152 11 3 11  4  6  6 -7.24766468e-03 4.99575434e-01 1.00000000e+00 5.06823098e-01 -9.36832178e-02 1.12509844e+00 4.16689591e-01 -3.17497789e-03 1.00703753e+00 4.58114974e-01 
+15 9.375000000000e-01 6.250000000000e-02 64 659 289 1152 11 3 11  4  6  6 -6.53794341e-03 4.99565727e-01 1.00000000e+00 5.06103670e-01 -1.62670326e-01 1.12504051e+00 4.16689228e-01 -2.99536198e-03 1.00612713e+00 4.58114092e-01 
+16 1.000000000000e+00 6.250000000000e-02 64 659 289 1152 12 3 11  4  6  6 -6.60290612e-03 4.99553058e-01 1.00000000e+00 5.06155964e-01 -2.11393866e-01 1.14894448e+00 4.16688918e-01 -3.11575787e-03 1.00627874e+00 4.58112412e-01 
diff --git a/unit_tests/introspection.cc b/unit_tests/introspection.cc
index f663efe382f..a74ea3834fc 100644
--- a/unit_tests/introspection.cc
+++ b/unit_tests/introspection.cc
@@ -24,7 +24,7 @@
 #include <aspect/simulator.h>
 #include <aspect/introspection.h>
 
-TEST_CASE("Introspection::1")
+TEST_CASE("Introspection::basic")
 {
   using namespace aspect;
   dealii::ParameterHandler prm;
@@ -59,11 +59,13 @@ TEST_CASE("Introspection::1")
   CHECK(introspection.block_indices.velocities == 0);
   CHECK(introspection.block_indices.pressure == 1);
   CHECK(introspection.block_indices.temperature == 2);
+  CHECK(introspection.block_indices.compositional_fields.size() == 3);
   CHECK(introspection.block_indices.compositional_fields[0] == 3);
   CHECK(introspection.block_indices.compositional_fields[1] == 4);
   CHECK(introspection.block_indices.compositional_fields[2] == 5);
 
   // sparsity pattern block index:
+  CHECK(introspection.block_indices.compositional_field_sparsity_pattern.size() == 3);
   CHECK(introspection.block_indices.compositional_field_sparsity_pattern[0] == 3);
   CHECK(introspection.block_indices.compositional_field_sparsity_pattern[1] == 3);
   CHECK(introspection.block_indices.compositional_field_sparsity_pattern[2] == 3);
@@ -76,3 +78,58 @@ TEST_CASE("Introspection::1")
   CHECK(introspection.get_composition_base_element_indices() == std::vector<unsigned int>({3}));
   CHECK(introspection.get_compositional_field_indices_with_base_element(3) == std::vector<unsigned int>({0,1,2}));
 }
+
+TEST_CASE("Introspection::different-composition-types")
+{
+  using namespace aspect;
+  dealii::ParameterHandler prm;
+  Simulator<2>::declare_parameters(prm);
+
+
+  prm.set("Output directory", "");
+
+  prm.enter_subsection("Compositional fields");
+  prm.set("Number of fields", "3");
+  prm.set("Names of fields", "a,b,c");
+  prm.leave_subsection();
+  prm.enter_subsection("Discretization");
+  prm.set("Use discontinuous temperature discretization", "false");
+  prm.set("Use discontinuous composition discretization", "true,false,false");
+  prm.leave_subsection();
+
+  Parameters<2> parameters(prm, MPI_COMM_WORLD);
+
+  const std::vector<VariableDeclaration<2>> variables = construct_default_variables (parameters);
+  Introspection<2> introspection (variables, parameters);
+
+  dealii::FESystem<2> fe(introspection.get_fes(), introspection.get_multiplicities());
+  CHECK(fe.get_name() == "FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(1)-FE_Q<2>(2)-FE_DGQ<2>(2)-FE_Q<2>(2)^2]");
+
+  CHECK(introspection.n_components == 7);
+  CHECK(introspection.n_compositional_fields == 3);
+  CHECK(introspection.n_blocks == 6);
+
+  CHECK(introspection.component_indices.velocities[0] == 0);
+  CHECK(introspection.component_indices.velocities[1] == 1);
+  CHECK(introspection.component_indices.pressure == 2);
+  CHECK(introspection.component_indices.temperature == 3);
+  CHECK(introspection.component_indices.compositional_fields[0] == 4);
+
+  CHECK(introspection.block_indices.velocities == 0);
+  CHECK(introspection.block_indices.pressure == 1);
+  CHECK(introspection.block_indices.temperature == 2);
+  CHECK(introspection.block_indices.compositional_fields[0] == 3);
+  CHECK(introspection.block_indices.compositional_fields[1] == 4);
+  CHECK(introspection.block_indices.compositional_fields[2] == 5);
+
+  // sparsity pattern block index:
+  CHECK(introspection.block_indices.compositional_field_sparsity_pattern[0] == 3);
+  CHECK(introspection.block_indices.compositional_field_sparsity_pattern[1] == 4);
+  CHECK(introspection.block_indices.compositional_field_sparsity_pattern[2] == 4);
+
+  // base elements:
+  CHECK(introspection.base_elements.compositional_fields == std::vector<unsigned int>({3,4,4}));
+  CHECK(introspection.get_composition_base_element_indices() == std::vector<unsigned int>({3,4}));
+  CHECK(introspection.get_compositional_field_indices_with_base_element(3) == std::vector<unsigned int>({0}));
+  CHECK(introspection.get_compositional_field_indices_with_base_element(4) == std::vector<unsigned int>({1,2}));
+}