From 33630e3817ee077abb6e3636bb89e5f7f62b509a Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Thu, 13 Jun 2024 16:32:07 -0400 Subject: [PATCH 1/2] allow continuous and discontinuous compositions --- include/aspect/fe_variable_collection.h | 9 +- include/aspect/introspection.h | 2 +- source/fe_variable_collection.cc | 21 ++- source/particle/world.cc | 3 + source/simulator/core.cc | 116 +++++++----- source/simulator/introspection.cc | 99 +++++++---- source/simulator/parameters.cc | 2 +- tests/composition_cg_dg_fem.prm | 96 ++++++++++ tests/composition_cg_dg_fem/screen-output | 204 ++++++++++++++++++++++ tests/composition_cg_dg_fem/statistics | 40 +++++ unit_tests/introspection.cc | 59 ++++++- 11 files changed, 566 insertions(+), 85 deletions(-) create mode 100644 tests/composition_cg_dg_fem.prm create mode 100644 tests/composition_cg_dg_fem/screen-output create mode 100644 tests/composition_cg_dg_fem/statistics 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 &variable(const std::string &name) const; + /** + * Return a vector of pointers of all variables with name @p name. + */ + std::vector*> 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 &fevs); + ComponentMasks (const FEVariableCollection &fevs, const Introspection::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 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 + std::vector*> + FEVariableCollection::variables_with_name(const std::string &name) const + { + std::vector*> result; + for (unsigned int i=0; i bool FEVariableCollection::variable_exists(const std::string &name) const diff --git a/source/particle/world.cc b/source/particle/world.cc index 918b6ea4920..a970700c594 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 d9af44c7436..72bc66fe319 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -880,6 +880,30 @@ namespace aspect + template + bool compositional_field_needs_matrix_block(const Introspection &introspection, const unsigned int composition_index) + { + const typename Simulator::AdvectionField adv_field (Simulator::AdvectionField::composition(composition_index)); + switch (adv_field.advection_method(introspection)) + { + case Parameters::AdvectionFieldMethod::fem_field: + case Parameters::AdvectionFieldMethod::fem_melt_field: + case Parameters::AdvectionFieldMethod::fem_darcy_field: + case Parameters::AdvectionFieldMethod::prescribed_field_with_diffusion: + return true; + case Parameters::AdvectionFieldMethod::particles: + case Parameters::AdvectionFieldMethod::volume_of_fluid: + case Parameters::AdvectionFieldMethod::static_field: + case Parameters::AdvectionFieldMethod::prescribed_field: + break; + default: + Assert (false, ExcNotImplemented()); + } + return false; + } + + + template bool compositional_fields_need_matrix_block(const Introspection &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::AdvectionField adv_field (Simulator::AdvectionField::composition(c)); - switch (adv_field.advection_method(introspection)) - { - case Parameters::AdvectionFieldMethod::fem_field: - case Parameters::AdvectionFieldMethod::fem_melt_field: - case Parameters::AdvectionFieldMethod::fem_darcy_field: - case Parameters::AdvectionFieldMethod::prescribed_field_with_diffusion: - return true; - case Parameters::AdvectionFieldMethod::particles: - case Parameters::AdvectionFieldMethod::volume_of_fluid: - case Parameters::AdvectionFieldMethod::static_field: - case Parameters::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::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; imultiplicity; ++i) + { + ci.compositional_fields.push_back(fe->first_component_index+i); + } + } + } return ci; } @@ -56,20 +64,24 @@ namespace aspect */ template typename Introspection::BlockIndices - setup_blocks (FEVariableCollection &fevs, const unsigned int n_compositional_fields) + setup_blocks (FEVariableCollection &fevs) { typename Introspection::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; imultiplicity; ++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 typename Introspection::BaseElements - setup_base_elements (FEVariableCollection &fevs, const unsigned int n_compositional_fields) + setup_base_elements (FEVariableCollection &fevs) { typename Introspection::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({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; imultiplicity; ++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( - "compositions", - internal::new_FE_Q_or_DGQ(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; c0 && 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> fe = internal::new_FE_Q_or_DGQ(parameters.use_discontinuous_composition_discretization[c], + parameters.composition_degree); + variables.push_back(VariableDeclaration("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(*this)), - n_blocks(FEVariableCollection::n_blocks()), - block_indices (internal::setup_blocks(*this, parameters.n_compositional_fields)), + n_blocks (FEVariableCollection::n_blocks()), + block_indices (internal::setup_blocks(*this)), extractors (component_indices), - base_elements (internal::setup_base_elements(*this, parameters.n_compositional_fields)), + base_elements (internal::setup_base_elements(*this)), polynomial_degree (internal::setup_polynomial_degree(parameters)), quadratures (internal::setup_quadratures(parameters, ReferenceCells::get_hypercube())), face_quadratures (internal::setup_face_quadratures(parameters, ReferenceCells::get_hypercube())), - 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 std::vector - make_component_mask_sequence(const FEVariable &variable) + make_composition_component_mask_sequence (const FEVariableCollection &fevs, const typename Introspection::ComponentIndices &indices) { std::vector result; - for (unsigned int i=0; i - Introspection::ComponentMasks::ComponentMasks (FEVariableCollection &fevs) + Introspection::ComponentMasks::ComponentMasks (const FEVariableCollection &fevs, const Introspection::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 22f1c02e003..a328a31649c 100644 --- a/source/simulator/parameters.cc +++ b/source/simulator/parameters.cc @@ -1011,7 +1011,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 #include -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({3})); CHECK(introspection.get_compositional_field_indices_with_base_element(3) == std::vector({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> 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({3,4,4})); + CHECK(introspection.get_composition_base_element_indices() == std::vector({3,4})); + 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})); +} From ad6c5010981f8bfcf1dd16ae264e9363e05121eb Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Tue, 18 Jun 2024 12:42:13 -0400 Subject: [PATCH 2/2] changelog --- doc/modules/changes/20240618_tjhei | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 doc/modules/changes/20240618_tjhei 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. +
+(Timo Heister, 2024/06/18)