diff --git a/doc/modules/changes/20240618_tjhei b/doc/modules/changes/20240618_tjhei
index 15ff52c82de..4c9269de89b 100644
--- 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.
(Timo Heister, 2024/06/18)
diff --git a/include/aspect/introspection.h b/include/aspect/introspection.h
index 9ce0167be3a..7dc66827e73 100644
--- 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 compositional_fields;
+ unsigned int max_compositional_field;
};
/**
@@ -310,7 +312,8 @@ namespace aspect
Quadrature velocities;
Quadrature pressure;
Quadrature temperature;
- Quadrature compositional_fields;
+ Quadrature compositional_field_max;
+ std::vector> compositional_fields;
Quadrature system;
};
diff --git a/include/aspect/parameters.h b/include/aspect/parameters.h
index 201fb425efa..3edcaecf291 100644
--- a/include/aspect/parameters.h
+++ b/include/aspect/parameters.h
@@ -689,7 +689,8 @@ namespace aspect
std::vector use_discontinuous_composition_discretization;
bool have_discontinuous_composition_discretization;
unsigned int temperature_degree;
- unsigned int composition_degree;
+ std::vector 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
index dd4022c4650..57ec63a1f2b 100644
--- 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(quadrature_positions),
this->get_mapping(),
in.requested_properties,
out_copy);
diff --git a/source/postprocess/composition_statistics.cc b/source/postprocess/composition_statistics.cc
index d0073fcf11c..b9ba83e4690 100644
--- 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 &quadrature_formula = this->introspection().quadratures.compositional_fields;
+ const Quadrature &quadrature_formula = this->introspection().quadratures.compositional_field_max;
const unsigned int n_q_points = quadrature_formula.size();
FEValues fe_values (this->get_mapping(),
diff --git a/source/postprocess/max_depth_field.cc b/source/postprocess/max_depth_field.cc
index d0cba1d2a56..9c1b73aecbe 100644
--- 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 &quadrature_formula = this->introspection().quadratures.compositional_fields;
+ const Quadrature &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
index 935afb8cb7e..a8eff1ef296 100644
--- 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 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
index 9ddae9f4ff6..ed47b059337 100644
--- a/source/simulator/entropy_viscosity.cc
+++ b/source/simulator/entropy_viscosity.cc
@@ -44,10 +44,10 @@ namespace aspect
return numbers::signaling_nan();
// record maximal entropy on Gauss quadrature points
- const Quadrature &quadrature_formula
- = (advection_field.is_temperature() ?
- introspection.quadratures.temperature :
- introspection.quadratures.compositional_fields);
+ const Quadrature &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
index 845c11c0c8c..f0a760024cf 100644
--- 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 &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
index c52c4e91181..22c564a0f82 100644
--- 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(parameters.temperature_degree+1);
- quadratures.compositional_fields = reference_cell.get_gauss_type_quadrature(parameters.composition_degree+1);
+ quadratures.compositional_field_max = reference_cell.get_gauss_type_quadrature(parameters.max_composition_degree+1);
+ for (const auto °ree: parameters.composition_degrees)
+ quadratures.compositional_fields.emplace_back(reference_cell.get_gauss_type_quadrature(degree+1));
quadratures.system = reference_cell.get_gauss_type_quadrature(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(parameters.temperature_degree+1);
- quadratures.compositional_fields = reference_cell.face_reference_cell(0).get_gauss_type_quadrature(parameters.composition_degree+1);
- quadratures.system = reference_cell.face_reference_cell(0).get_gauss_type_quadrature(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(parameters.max_composition_degree+1);
+ quadratures.system = reference_cell.face_reference_cell(0).get_gauss_type_quadrature(
+ 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; c0 && 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> fe = internal::new_FE_Q_or_DGQ(parameters.use_discontinuous_composition_discretization[c],
- parameters.composition_degree);
+ parameters.composition_degrees[c]);
variables.push_back(VariableDeclaration("compositions", fe, 1, 1));
}
}
diff --git a/source/simulator/lateral_averaging.cc b/source/simulator/lateral_averaging.cc
index 3b1218c80ec..cdb57d88d14 100644
--- 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(std::ceil((max_fe_degree+1.0)/2.0));
+ const unsigned int lateral_quadrature_degree = static_cast(std::ceil((this->introspection().polynomial_degree.max_degree+1.0)/2.0));
std::unique_ptr> quadrature_formula;
if (geometry_unique_depth_direction != numbers::invalid_unsigned_int)
diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
index 630301c8fce..b983c88de84 100644
--- 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 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
index 4959bbadef3..0484b116e9f 100644
--- 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
new file mode 100644
index 00000000000..2132bdbec85
--- /dev/null
+++ 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
diff --git a/tests/composition_cg_dg_fem2/screen-output b/tests/composition_cg_dg_fem2/screen-output
new file mode 100644
index 00000000000..a226728bd69
--- /dev/null
+++ b/tests/composition_cg_dg_fem2/screen-output
@@ -0,0 +1,113 @@
+
+Number of active cells: 16 (on 3 levels)
+Number of degrees of freedom: 535 (162+25+81+16+81+81+25+64)
+
+*** 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 C_3 system ... 0 iterations.
+ Solving C_4 system ... 0 iterations.
+ Solving C_5 system ... 0 iterations.
+ Solving Stokes system... 10+0 iterations.
+
+Number of active cells: 40 (on 4 levels)
+Number of degrees of freedom: 1,275 (386+55+193+40+193+193+55+160)
+
+*** 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 C_3 system ... 0 iterations.
+ Solving C_4 system ... 0 iterations.
+ Solving C_5 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.5 // 0/1/0.4583 // 0/1/0.4583 // 0/1/0.375 // 0/1/0.375
+
+Total system matrix memory consumption: 0.21 MB.
+Total system matrix nnz: 8,489
+system matrix nnz by block:
+ 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0
+ 0 0 2,157 0 0 0 0 0
+ 0 0 0 180 0 0 0 0
+ 0 0 0 0 2,845 0 0 0
+ 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 427 0
+ 0 0 0 0 0 0 0 2,880
+
+Total system preconditioner matrix memory consumption: 0.00 MB.
+Total system preconditioner matrix nnz: 0
+system preconditioner matrix nnz by block:
+
+
+*** 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 C_3 system ... 10 iterations.
+ Solving C_4 system ... 6 iterations.
+ Solving C_5 system ... 3 iterations.
+ Solving Stokes system... 8+0 iterations.
+
+ Postprocessing:
+ Temperature min/avg/max: 0 K, 0.5 K, 1 K
+ Compositions min/max/mass: 1.106e-10/0.9992/0.5 // -0.0158/1.028/0.4583 // -0.0158/1.028/0.4583 // -0.07232/1.07/0.375 // -0.1145/1.119/0.375
+
+Total system matrix memory consumption: 0.21 MB.
+Total system matrix nnz: 8,489
+system matrix nnz by block:
+ 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0
+ 0 0 2,157 0 0 0 0 0
+ 0 0 0 180 0 0 0 0
+ 0 0 0 0 2,845 0 0 0
+ 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 427 0
+ 0 0 0 0 0 0 0 2,880
+
+Total system preconditioner matrix memory consumption: 0.00 MB.
+Total system preconditioner matrix nnz: 0
+system preconditioner matrix nnz by block:
+
+
+*** Timestep 2: t=0.1 seconds, dt=0.0375 seconds
+ Solving temperature system... 10 iterations.
+ Solving C_1 system ... 3 iterations.
+ Solving C_2 system ... 10 iterations.
+ Solving C_3 system ... 10 iterations.
+ Solving C_4 system ... 6 iterations.
+ Solving C_5 system ... 3 iterations.
+ Solving Stokes system... 13+0 iterations.
+
+ Postprocessing:
+ Temperature min/avg/max: 0 K, 0.5 K, 1 K
+ Compositions min/max/mass: 2.282e-10/0.9985/0.5 // -0.02062/1.033/0.4583 // -0.02062/1.033/0.4583 // -0.1054/1.098/0.375 // -0.1495/1.157/0.3749
+
+Total system matrix memory consumption: 0.21 MB.
+Total system matrix nnz: 8,489
+system matrix nnz by block:
+ 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 0 0
+ 0 0 2,157 0 0 0 0 0
+ 0 0 0 180 0 0 0 0
+ 0 0 0 0 2,845 0 0 0
+ 0 0 0 0 0 0 0 0
+ 0 0 0 0 0 0 427 0
+ 0 0 0 0 0 0 0 2,880
+
+Total system preconditioner matrix memory consumption: 0.00 MB.
+Total system preconditioner matrix nnz: 0
+system preconditioner matrix nnz by block:
+
+
+Number of active cells: 64 (on 4 levels)
+Number of degrees of freedom: 1,927 (578+81+289+64+289+289+81+256)
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/composition_cg_dg_fem2/statistics b/tests/composition_cg_dg_fem2/statistics
new file mode 100644
index 00000000000..31a3fcc7646
--- /dev/null
+++ b/tests/composition_cg_dg_fem2/statistics
@@ -0,0 +1,38 @@
+# 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 composition solver 4
+# 13: Iterations for composition solver 5
+# 14: Iterations for Stokes solver
+# 15: Velocity iterations in Stokes preconditioner
+# 16: Schur complement iterations in Stokes preconditioner
+# 17: Minimal temperature (K)
+# 18: Average temperature (K)
+# 19: Maximal temperature (K)
+# 20: Average nondimensional temperature (K)
+# 21: Minimal value for composition C_1
+# 22: Maximal value for composition C_1
+# 23: Global mass for composition C_1
+# 24: Minimal value for composition C_2
+# 25: Maximal value for composition C_2
+# 26: Global mass for composition C_2
+# 27: Minimal value for composition C_3
+# 28: Maximal value for composition C_3
+# 29: Global mass for composition C_3
+# 30: Minimal value for composition C_4
+# 31: Maximal value for composition C_4
+# 32: Global mass for composition C_4
+# 33: Minimal value for composition C_5
+# 34: Maximal value for composition C_5
+# 35: Global mass for composition C_5
+0 0.000000000000e+00 0.000000000000e+00 40 441 193 200 0 0 0 0 0 0 12 14 14 0.00000000e+00 5.00000000e-01 1.00000000e+00 5.00000000e-01 0.00000000e+00 1.00000000e+00 5.00000000e-01 0.00000000e+00 1.00000000e+00 4.58333333e-01 0.00000000e+00 1.00000000e+00 4.58333333e-01 0.00000000e+00 1.00000000e+00 3.75000000e-01 0.00000000e+00 1.00000000e+00 3.75000000e-01
+1 6.250000000000e-02 6.250000000000e-02 40 441 193 200 10 3 10 10 6 3 7 9 9 0.00000000e+00 5.00021358e-01 1.00000000e+00 5.00021358e-01 1.10561158e-10 9.99186176e-01 5.00020085e-01 -1.57952550e-02 1.02777854e+00 4.58314382e-01 -1.57952550e-02 1.02777854e+00 4.58314382e-01 -7.23249290e-02 1.06953675e+00 3.75000000e-01 -1.14531625e-01 1.11905441e+00 3.74972012e-01
+2 1.000000000000e-01 3.750000000000e-02 40 441 193 200 10 3 10 10 6 3 12 14 14 0.00000000e+00 5.00024722e-01 1.00000000e+00 5.00024722e-01 2.28231019e-10 9.98489448e-01 5.00036945e-01 -2.06184196e-02 1.03261208e+00 4.58295584e-01 -2.06184196e-02 1.03261208e+00 4.58295584e-01 -1.05388605e-01 1.09844303e+00 3.75000000e-01 -1.49536172e-01 1.15735921e+00 3.74948780e-01
diff --git a/unit_tests/introspection.cc b/unit_tests/introspection.cc
index a74ea3834fc..0f86b2dc714 100644
--- a/unit_tests/introspection.cc
+++ b/unit_tests/introspection.cc
@@ -59,16 +59,9 @@ TEST_CASE("Introspection::basic")
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);
+ CHECK(introspection.block_indices.compositional_fields == std::vector({3,4,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);
+ CHECK(introspection.block_indices.compositional_field_sparsity_pattern == std::vector({3,3,3}));
CHECK(introspection.base_elements.velocities == 0);
CHECK(introspection.base_elements.pressure == 1);
@@ -118,14 +111,9 @@ TEST_CASE("Introspection::different-composition-types")
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);
+ CHECK(introspection.block_indices.compositional_fields == std::vector({3,4,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);
+ CHECK(introspection.block_indices.compositional_field_sparsity_pattern == std::vector({3,4,4}));
// base elements:
CHECK(introspection.base_elements.compositional_fields == std::vector({3,4,4}));
@@ -133,3 +121,57 @@ TEST_CASE("Introspection::different-composition-types")
CHECK(introspection.get_compositional_field_indices_with_base_element(3) == std::vector({0}));
CHECK(introspection.get_compositional_field_indices_with_base_element(4) == std::vector({1,2}));
}
+
+TEST_CASE("Introspection::different-composition-degrees")
+{
+ 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", "5");
+ prm.set("Names of fields", "c1,c2,c3,c4,c5");
+ prm.leave_subsection();
+ prm.enter_subsection("Discretization");
+ prm.set("Use discontinuous temperature discretization", "false");
+ prm.set("Use discontinuous composition discretization", "false,true,false,false,false");
+ prm.set("Composition polynomial degree", "1,0,2,2,3");
+ prm.leave_subsection();
+
+ Parameters<2> parameters(prm, MPI_COMM_WORLD);
+
+ const std::vector> 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_Q<2>(1)-FE_DGQ<2>(0)-FE_Q<2>(2)^2-FE_Q<2>(3)]");
+
+ CHECK(introspection.n_components == 9);
+ CHECK(introspection.n_compositional_fields == 5);
+ CHECK(introspection.n_blocks == 8);
+
+ 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 == std::vector({4,5,6,7,8}));
+
+ CHECK(introspection.block_indices.velocities == 0);
+ CHECK(introspection.block_indices.pressure == 1);
+ CHECK(introspection.block_indices.temperature == 2);
+ CHECK(introspection.block_indices.compositional_fields == std::vector({3,4,5,6,7}));
+
+ // sparsity pattern block index:
+ CHECK(introspection.block_indices.compositional_field_sparsity_pattern == std::vector({3,4,5,5,7}));
+
+ // base elements:
+ CHECK(introspection.base_elements.compositional_fields == std::vector({3,4,5,5,6}));
+ CHECK(introspection.get_composition_base_element_indices() == std::vector({3,4,5,6}));
+ CHECK(introspection.get_compositional_field_indices_with_base_element(3) == std::vector({0}));
+ CHECK(introspection.get_compositional_field_indices_with_base_element(4) == std::vector({1}));
+ CHECK(introspection.get_compositional_field_indices_with_base_element(5) == std::vector({2,3}));
+ CHECK(introspection.get_compositional_field_indices_with_base_element(6) == std::vector({4}));
+}