Merge pull request #5930 from tjhei/composition-degrees

support different composition degrees
geodynamics · Jun 27, 2024 · 9b91208 · 9b91208
2 parents 3451d1e + e91dfae
commit 9b91208
Show file tree

Hide file tree

Showing 17 changed files with 363 additions and 57 deletions.
diff --git a/doc/modules/changes/20240618_tjhei b/doc/modules/changes/20240618_tjhei
@@ -1,5 +1,5 @@
 New: ASPECT now supports compositional fields with
 different discretizations (continuous or discontinuous)
-at the same time.
+and different polynomial degrees at the same time.
 <br>
 (Timo Heister, 2024/06/18)
diff --git a/include/aspect/introspection.h b/include/aspect/introspection.h
@@ -275,9 +275,11 @@ namespace aspect
        */
       struct PolynomialDegree
       {
-        unsigned int       velocities;
-        unsigned int       temperature;
-        unsigned int       compositional_fields;
+        unsigned int              max_degree;
+        unsigned int              velocities;
+        unsigned int              temperature;
+        std::vector<unsigned int> compositional_fields;
+        unsigned int              max_compositional_field;
       };
 
       /**
@@ -310,7 +312,8 @@ namespace aspect
         Quadrature<dim>       velocities;
         Quadrature<dim>       pressure;
         Quadrature<dim>       temperature;
-        Quadrature<dim>       compositional_fields;
+        Quadrature<dim>       compositional_field_max;
+        std::vector<Quadrature<dim>> compositional_fields;
         Quadrature<dim>       system;
       };
 

diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h
@@ -689,7 +689,8 @@ namespace aspect
     std::vector<bool>              use_discontinuous_composition_discretization;
     bool                           have_discontinuous_composition_discretization;
     unsigned int                   temperature_degree;
-    unsigned int                   composition_degree;
+    std::vector<unsigned int>      composition_degrees;
+    unsigned int                   max_composition_degree;
     std::string                    pressure_normalization;
     MaterialModel::MaterialAveraging::AveragingOperation material_averaging;
 

diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc
@@ -376,7 +376,7 @@ namespace aspect
                   get_averaging_operation_for_viscosity(this->get_parameters().material_averaging);
                 MaterialAveraging::average(averaging_operation_for_viscosity,
                                            in.current_cell,
-                                           this->introspection().quadratures.compositional_fields,
+                                           Quadrature<dim>(quadrature_positions),
                                            this->get_mapping(),
                                            in.requested_properties,
                                            out_copy);

diff --git a/source/postprocess/composition_statistics.cc b/source/postprocess/composition_statistics.cc
@@ -37,7 +37,7 @@ namespace aspect
         return {"", ""};
 
       // create a quadrature formula based on the compositional element alone.
-      const Quadrature<dim> &quadrature_formula = this->introspection().quadratures.compositional_fields;
+      const Quadrature<dim> &quadrature_formula = this->introspection().quadratures.compositional_field_max;
       const unsigned int n_q_points = quadrature_formula.size();
 
       FEValues<dim> fe_values (this->get_mapping(),

diff --git a/source/postprocess/max_depth_field.cc b/source/postprocess/max_depth_field.cc
@@ -38,7 +38,7 @@ namespace aspect
       // be defensive about determining that a compositional field actually exists
       AssertThrow(this->n_compositional_fields() > 0,
                   ExcMessage("This postprocessor cannot be used without compositional fields."));
-      const Quadrature<dim> &quadrature_formula = this->introspection().quadratures.compositional_fields;
+      const Quadrature<dim> &quadrature_formula = this->introspection().quadratures.compositional_field_max;
 
       const unsigned int n_q_points = quadrature_formula.size();
 

diff --git a/source/simulator/checkpoint_restart.cc b/source/simulator/checkpoint_restart.cc
@@ -91,7 +91,7 @@ namespace aspect
       oa << parameters.use_discontinuous_temperature_discretization;
       oa << parameters.use_discontinuous_composition_discretization;
       oa << parameters.temperature_degree;
-      oa << parameters.composition_degree;
+      oa << parameters.composition_degrees;
       oa << parameters.pressure_normalization;
       oa << parameters.n_compositional_fields;
       oa << parameters.names_of_compositional_fields;
@@ -196,9 +196,9 @@ namespace aspect
                                "These need to be the same during restarting "
                                "from a checkpoint."));
 
-      unsigned int composition_degree;
-      ia >> composition_degree;
-      AssertThrow (composition_degree == parameters.composition_degree,
+      std::vector<unsigned int> composition_degrees;
+      ia >> composition_degrees;
+      AssertThrow (composition_degrees == parameters.composition_degrees,
                    ExcMessage ("The composition polynomial degree that was stored "
                                "in the checkpoint file is not the same as the one "
                                "you currently set in your input file. "

diff --git a/source/simulator/entropy_viscosity.cc b/source/simulator/entropy_viscosity.cc
@@ -44,10 +44,10 @@ namespace aspect
       return numbers::signaling_nan<double>();
 
     // record maximal entropy on Gauss quadrature points
-    const Quadrature<dim> &quadrature_formula
-      = (advection_field.is_temperature() ?
-         introspection.quadratures.temperature :
-         introspection.quadratures.compositional_fields);
+    const Quadrature<dim> &quadrature_formula =
+      (advection_field.is_temperature() ?
+       introspection.quadratures.temperature :
+       introspection.quadratures.compositional_fields[advection_field.compositional_variable]);
     const unsigned int n_q_points = quadrature_formula.size();
 
     const FEValuesExtractors::Scalar field = advection_field.scalar_extractor(introspection);

diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
@@ -201,7 +201,7 @@ namespace aspect
     if (this->is_temperature())
       return introspection.polynomial_degree.temperature;
     else
-      return introspection.polynomial_degree.compositional_fields;
+      return introspection.polynomial_degree.compositional_fields[compositional_variable];
   }
 
 
@@ -1507,7 +1507,7 @@ namespace aspect
     const Quadrature<dim> &quadrature_formula_0
       = (advection_field.is_temperature() ?
          introspection.quadratures.temperature :
-         introspection.quadratures.compositional_fields);
+         introspection.quadratures.compositional_fields[advection_field.compositional_variable]);
     const unsigned int n_q_points_0 = quadrature_formula_0.size();
 
     // fe values for points evaluation

diff --git a/source/simulator/introspection.cc b/source/simulator/introspection.cc
@@ -124,7 +124,13 @@ namespace aspect
 
       polynomial_degree.velocities = parameters.stokes_velocity_degree;
       polynomial_degree.temperature = parameters.temperature_degree;
-      polynomial_degree.compositional_fields = parameters.composition_degree;
+      polynomial_degree.compositional_fields = parameters.composition_degrees;
+      polynomial_degree.max_compositional_field = parameters.max_composition_degree;
+      polynomial_degree.max_degree = std::max({parameters.stokes_velocity_degree,
+                                               parameters.temperature_degree,
+                                               (parameters.n_compositional_fields>0 ? parameters.max_composition_degree : 0u)
+                                              }
+                                             );
 
       return polynomial_degree;
     }
@@ -145,10 +151,12 @@ namespace aspect
                                                                            :
                                                                            parameters.stokes_velocity_degree);
       quadratures.temperature = reference_cell.get_gauss_type_quadrature<dim>(parameters.temperature_degree+1);
-      quadratures.compositional_fields = reference_cell.get_gauss_type_quadrature<dim>(parameters.composition_degree+1);
+      quadratures.compositional_field_max = reference_cell.get_gauss_type_quadrature<dim>(parameters.max_composition_degree+1);
+      for (const auto &degree: parameters.composition_degrees)
+        quadratures.compositional_fields.emplace_back(reference_cell.get_gauss_type_quadrature<dim>(degree+1));
       quadratures.system = reference_cell.get_gauss_type_quadrature<dim>(std::max({parameters.stokes_velocity_degree,
                                                                                    parameters.temperature_degree,
-                                                                                   parameters.composition_degree
+                                                                                   parameters.max_composition_degree
                                                                                   }) + 1);
 
       return quadratures;
@@ -170,11 +178,9 @@ namespace aspect
                              :
                              parameters.stokes_velocity_degree);
       quadratures.temperature = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.temperature_degree+1);
-      quadratures.compositional_fields = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.composition_degree+1);
-      quadratures.system = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(std::max({parameters.stokes_velocity_degree,
-                           parameters.temperature_degree,
-                           parameters.composition_degree
-                                                                                                           }) + 1);
+      quadratures.compositional_fields = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.max_composition_degree+1);
+      quadratures.system = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(
+                             std::max({parameters.stokes_velocity_degree, parameters.temperature_degree, parameters.max_composition_degree}) + 1);
 
       return quadratures;
     }
@@ -253,7 +259,9 @@ namespace aspect
     // "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])
+        if (c>0
+            && parameters.use_discontinuous_composition_discretization[c] == parameters.use_discontinuous_composition_discretization[c-1]
+            && parameters.composition_degrees[c] == parameters.composition_degrees[c-1])
           {
             // reuse last one because it is the same
             variables.back().multiplicity += 1;
@@ -262,7 +270,7 @@ namespace aspect
         else
           {
             std::shared_ptr<FiniteElement<dim>> fe = internal::new_FE_Q_or_DGQ<dim>(parameters.use_discontinuous_composition_discretization[c],
-                                                                                     parameters.composition_degree);
+                                                                                     parameters.composition_degrees[c]);
             variables.push_back(VariableDeclaration<dim>("compositions", fe, 1, 1));
           }
       }

diff --git a/source/simulator/lateral_averaging.cc b/source/simulator/lateral_averaging.cc
@@ -589,13 +589,9 @@ namespace aspect
     else
       geometry_unique_depth_direction = numbers::invalid_unsigned_int;
 
-    const unsigned int max_fe_degree = std::max(this->introspection().polynomial_degree.velocities,
-                                                std::max(this->introspection().polynomial_degree.temperature,
-                                                         this->introspection().polynomial_degree.compositional_fields));
-
     // We want to integrate over a polynomial of degree p = max_fe_degree, for which we
     // need a quadrature of at least q, with p <= 2q-1 --> q >= (p+1)/2
-    const unsigned int lateral_quadrature_degree = static_cast<unsigned int>(std::ceil((max_fe_degree+1.0)/2.0));
+    const unsigned int lateral_quadrature_degree = static_cast<unsigned int>(std::ceil((this->introspection().polynomial_degree.max_degree+1.0)/2.0));
 
     std::unique_ptr<Quadrature<dim>> quadrature_formula;
     if (geometry_unique_depth_direction != numbers::invalid_unsigned_int)

diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
@@ -943,7 +943,7 @@ namespace aspect
                          "discontinuous field. "
                          "Units: None.");
       prm.declare_entry ("Composition polynomial degree", "2",
-                         Patterns::Integer (0),
+                         Patterns::List(Patterns::Integer (0)),
                          "The polynomial degree to use for the composition variable(s). "
                          "As an example, a value of 2 for this parameter will yield "
                          "either the element $Q_2$ or $DGQ_2$ for the compositional "
@@ -1752,7 +1752,14 @@ namespace aspect
     {
       stokes_velocity_degree = prm.get_integer ("Stokes velocity polynomial degree");
       temperature_degree     = prm.get_integer ("Temperature polynomial degree");
-      composition_degree     = prm.get_integer ("Composition polynomial degree");
+      composition_degrees    = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_unsigned_int(Utilities::split_string_list(prm.get("Composition polynomial degree"))),
+                                                                       n_compositional_fields,
+                                                                       "Composition polynomial degree");
+      if (n_compositional_fields > 0)
+        max_composition_degree = *std::max_element(composition_degrees.begin(), composition_degrees.end());
+      else
+        max_composition_degree = numbers::invalid_unsigned_int;
+
       use_locally_conservative_discretization
         = prm.get_bool ("Use locally conservative discretization");
       use_equal_order_interpolation_for_stokes
@@ -1767,10 +1774,10 @@ namespace aspect
         (std::find(use_discontinuous_composition_discretization.begin(), use_discontinuous_composition_discretization.end(), true)
          != use_discontinuous_composition_discretization.end());
 
-      // TODO: this needs to be modified once we lift the restriction that all compositions have the same degree.
-      AssertThrow(have_discontinuous_composition_discretization == true || composition_degree > 0,
-                  ExcMessage("Using a composition polynomial degree of 0 (cell-wise constant composition) "
-                             "is only supported if a discontinuous composition discretization is selected."));
+      for (unsigned int c=0; c<n_compositional_fields; ++c)
+        AssertThrow(use_discontinuous_composition_discretization[c] == true || composition_degrees[c] > 0,
+                    ExcMessage("Using a composition polynomial degree of 0 (cell-wise constant composition) "
+                               "is only supported if a discontinuous composition discretization is selected."));
 
       prm.enter_subsection ("Stabilization parameters");
       {

diff --git a/tests/composition_cg_dg_fem.prm b/tests/composition_cg_dg_fem.prm
@@ -1,4 +1,4 @@
-#
+# Test CG and DG compositional fields at the same time
 #
 # copied from discontinuous_composition_1.prm
 

diff --git a/tests/composition_cg_dg_fem2.prm b/tests/composition_cg_dg_fem2.prm
@@ -0,0 +1,98 @@
+# Test different composition degrees and cg/dg
+#
+# copied from discontinuous_composition_1.prm
+
+
+set Dimension                              = 2
+set Start time                             = 0
+set End time                               = 0.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
+  set Material averaging = harmonic average only viscosity
+
+  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,false,false,true
+  set Composition polynomial degree = 0,2,2,1,1
+end
+
+subsection Postprocess
+  set List of postprocessors = temperature statistics, composition statistics, matrix 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 = 5
+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); if(y>0.8, 1, 0); if(y>0.8, 1, 0); if(y>0.8, 1, 0)
+  end
+end