diff --git a/doc/modules/changes/20260311_lhy11009 b/doc/modules/changes/20260311_lhy11009 new file mode 100644 index 00000000000..62f8eda00b1 --- /dev/null +++ b/doc/modules/changes/20260311_lhy11009 @@ -0,0 +1,5 @@ +Added: a new instance of 'initial composition' in the prescribed solution plugin system. +This plugin prescribes compositional fields from the initial composition model within a +region defined by an indicator function. +
+(Haoyuan Li, 2026/03/11) diff --git a/doc/sphinx/parameters/Prescribed_20solution.md b/doc/sphinx/parameters/Prescribed_20solution.md index 4742bd37279..a40f8196b6e 100644 --- a/doc/sphinx/parameters/Prescribed_20solution.md +++ b/doc/sphinx/parameters/Prescribed_20solution.md @@ -9,12 +9,14 @@ :name: parameters:Prescribed_20solution/List_20of_20model_20names **Default value:** -**Pattern:** [MultipleSelection initial temperature|temperature function|velocity function ] +**Pattern:** [MultipleSelection initial composition|initial temperature|temperature function|velocity function ] **Documentation:** A comma-separated list of prescribed solution models that will be used to compute the solution in certain regions. These plugins are loaded in the order given, and are combined via the operators listed in ’List of model operators’. The following prescribed solution models are available: +‘initial composition’: Prescribe compositional fields from the initial composition model within a region defined by an indicator function. + ‘initial temperature’: Prescribe the temperature in a selected region using the active initial temperature model. The selected region is defined through an indicator function. At locations where the indicator value is greater than 0.5, the temperature is constrained to the initial temperature evaluated at that position. ‘temperature function’: Prescribe the temperature in terms of an explicit formula. The format of these functions follows the syntax understood by the muparser library. @@ -22,6 +24,50 @@ The following prescribed solution models are available: ‘velocity function’: Prescribe the velocity in terms of an explicit formula. The format of these functions follows the syntax understood by the muparser library, see {ref}`sec:run-aspect:parameters-overview:muparser-format`. :::: +(parameters:Prescribed_20solution/Initial_20composition)= +## **Subsection:** Prescribed solution / Initial composition +::::{dropdown} __Parameter:__ {ref}`Coordinate system` +:name: parameters:Prescribed_20solution/Initial_20composition/Coordinate_20system +**Default value:** cartesian + +**Pattern:** [Selection cartesian|spherical|depth ] + +**Documentation:** Coordinate system used for evaluating the indicator function. +:::: + +(parameters:Prescribed_20solution/Initial_20composition/Indicator_20function)= +## **Subsection:** Prescribed solution / Initial composition / Indicator function +::::{dropdown} __Parameter:__ {ref}`Function constants` +:name: parameters:Prescribed_20solution/Initial_20composition/Indicator_20function/Function_20constants +**Default value:** + +**Pattern:** [Anything] + +**Documentation:** Sometimes it is convenient to use symbolic constants in the expression that describes the function, rather than having to use its numeric value everywhere the constant appears. These values can be defined using this parameter, in the form ‘var1=value1, var2=value2, ...’. + +A typical example would be to set this runtime parameter to ‘pi=3.1415926536’ and then use ‘pi’ in the expression of the actual formula. (That said, for convenience this class actually defines both ‘pi’ and ‘Pi’ by default, but you get the idea.) +:::: + +::::{dropdown} __Parameter:__ {ref}`Function expression` +:name: parameters:Prescribed_20solution/Initial_20composition/Indicator_20function/Function_20expression +**Default value:** 0 + +**Pattern:** [Anything] + +**Documentation:** The formula that denotes the function you want to evaluate for particular values of the independent variables. This expression may contain any of the usual operations such as addition or multiplication, as well as all of the common functions such as ‘sin’ or ‘cos’. In addition, it may contain expressions like ‘if(x>0, 1, -1)’ where the expression evaluates to the second argument if the first argument is true, and to the third argument otherwise. For a full overview of possible expressions accepted see the documentation of the muparser library at http://muparser.beltoforion.de/. + +If the function you are describing represents a vector-valued function with multiple components, then separate the expressions for individual components by a semicolon. +:::: + +::::{dropdown} __Parameter:__ {ref}`Variable names` +:name: parameters:Prescribed_20solution/Initial_20composition/Indicator_20function/Variable_20names +**Default value:** x,y,t + +**Pattern:** [Anything] + +**Documentation:** The names of the variables as they will be used in the function, separated by commas. By default, the names of variables at which the function will be evaluated are ‘x’ (in 1d), ‘x,y’ (in 2d) or ‘x,y,z’ (in 3d) for spatial coordinates and ‘t’ for time. You can then use these variable names in your function expression and they will be replaced by the values of these variables at which the function is currently evaluated. However, you can also choose a different set of names for the independent variables at which to evaluate your function expression. For example, if you work in spherical coordinates, you may wish to set this input parameter to ‘r,phi,theta,t’ and then use these variable names in your function expression. +:::: + (parameters:Prescribed_20solution/Initial_20temperature)= ## **Subsection:** Prescribed solution / Initial temperature ::::{dropdown} __Parameter:__ {ref}`Coordinate system` diff --git a/include/aspect/prescribed_solution/initial_composition.h b/include/aspect/prescribed_solution/initial_composition.h new file mode 100644 index 00000000000..bbf57884f0d --- /dev/null +++ b/include/aspect/prescribed_solution/initial_composition.h @@ -0,0 +1,98 @@ +/* + Copyright (C) 2015 - 2022 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#ifndef _aspect_prescribed_solution_initial_composition_h +#define _aspect_prescribed_solution_initial_composition_h + +#include +#include +#include + +#include + +namespace aspect +{ + namespace PrescribedSolution + { + /** + * Prescribe compositional fields using the initial composition model. + * The region where the constraint applies is determined by an + * indicator function. + */ + template + class InitialComposition + : public Interface, + public SimulatorAccess + { + public: + + InitialComposition (); + + /** + * Store shared pointer to initial composition manager. + */ + void initialize () override; + + /** + * Update function time. + */ + void update () override; + + /** + * Declare parameters. + */ + static void declare_parameters (ParameterHandler &prm); + + /** + * Parse parameters. + */ + void parse_parameters (ParameterHandler &prm) override; + + /** + * Constrain compositional solution. + */ + void constrain_solution (const typename DoFHandler::active_cell_iterator &cell, + const std::vector> &positions, + const std::vector &component_indices, + std::vector &should_be_constrained, + std::vector &solution) override; + + private: + + /** + * Indicator function defining region where composition + * is prescribed. + */ + Functions::ParsedFunction indicator_function; + + /** + * Coordinate system used for evaluating indicator. + */ + Utilities::Coordinates::CoordinateSystem coordinate_system; + + /** + * Pointer to initial composition manager. + */ + std::shared_ptr> initial_composition_manager; + }; + } +} + +#endif diff --git a/source/prescribed_solution/initial_composition.cc b/source/prescribed_solution/initial_composition.cc new file mode 100644 index 00000000000..57b8b0391f9 --- /dev/null +++ b/source/prescribed_solution/initial_composition.cc @@ -0,0 +1,176 @@ +/* + Copyright (C) 2011 - 2022 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#include +#include +#include + +namespace aspect +{ + namespace PrescribedSolution + { + + template + InitialComposition::InitialComposition () + : + indicator_function(1) + {} + + + + template + void + InitialComposition::initialize () + { + initial_composition_manager = + this->get_initial_composition_manager_pointer(); + } + + + + template + void + InitialComposition::update () + { + if (this->convert_output_to_years()) + indicator_function.set_time(this->get_time() / year_in_seconds); + else + indicator_function.set_time(this->get_time()); + } + + + + template + void + InitialComposition::constrain_solution (const typename DoFHandler::active_cell_iterator &, + const std::vector> &positions, + const std::vector &component_indices, + std::vector &should_be_constrained, + std::vector &solution) + { + // Determine the component range corresponding to compositional fields. + // These component indices are originally mapped from local DoFs and + // determine which DoFs in the system belong to compositional fields. + const unsigned int first_comp = + this->introspection().component_indices.compositional_fields[0]; + + const unsigned int last_comp = + this->introspection().component_indices.compositional_fields[this->introspection().n_compositional_fields-1]; + + const auto &geometry_model = this->get_geometry_model(); + + for (unsigned int q=0; q last_comp) + continue; + + const auto point = geometry_model.cartesian_to_other_coordinates( + positions[q], coordinate_system); + + const double indicator = indicator_function.value( + Utilities::convert_array_to_point(point.get_coordinates())); + + if (indicator > 0.5) + { + const unsigned int field_index = component - first_comp; + + solution[q] = + initial_composition_manager->initial_composition( + positions[q], field_index); + + should_be_constrained[q] = true; + } + } + } + + + + template + void + InitialComposition::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Prescribed solution"); + { + prm.enter_subsection("Initial composition"); + { + + prm.declare_entry("Coordinate system", + "cartesian", + Patterns::Selection("cartesian|spherical|depth"), + "Coordinate system used for evaluating the indicator function."); + + prm.enter_subsection("Indicator function"); + { + Functions::ParsedFunction::declare_parameters(prm,1); + } + prm.leave_subsection(); + + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + + + template + void + InitialComposition::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Prescribed solution"); + { + prm.enter_subsection("Initial composition"); + { + + coordinate_system = + Utilities::Coordinates::string_to_coordinate_system( + prm.get("Coordinate system")); + + prm.enter_subsection("Indicator function"); + { + indicator_function.parse_parameters(prm); + } + prm.leave_subsection(); + + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + } +} + + + +namespace aspect +{ + namespace PrescribedSolution + { + + ASPECT_REGISTER_PRESCRIBED_SOLUTION(InitialComposition, + "initial composition", + "Prescribe compositional fields from the initial composition model " + "within a region defined by an indicator function.") + + } +} diff --git a/tests/prescribed_solution_composition_initial_condition.prm b/tests/prescribed_solution_composition_initial_condition.prm new file mode 100644 index 00000000000..7332729b385 --- /dev/null +++ b/tests/prescribed_solution_composition_initial_condition.prm @@ -0,0 +1,41 @@ +######################################################### +# This is modified from the composition_passive.prm +# parameter file by prescribing the compositional solution +# from initial compositional fields. + +include $ASPECT_SOURCE_DIR/cookbooks/composition_passive/composition_passive.prm + +# Modify the original cookbook to end earlier +set Start time = 0 +set End time = 1 +set Use years instead of seconds = false + +# Modify the original cookbook to have coarser meshes +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 3 + set Time steps between mesh refinement = 0 +end + +# The next section prescribes the composition inside the domain +# using a spatially and temporally varying function. An indicator +# function selects the region where the composition is constrained. +subsection Prescribed solution + set List of model names = initial composition + subsection Initial composition + subsection Indicator function + set Variable names = x, y + set Function expression = (x<1) ? 1:0 + end + end +end + +subsection Postprocess + set List of postprocessors = visualization + + subsection Visualization + set Interpolate output = false + set Time between graphical output = 0.1 + set Output format = gnuplot + end +end diff --git a/tests/prescribed_solution_composition_initial_condition/screen-output b/tests/prescribed_solution_composition_initial_condition/screen-output new file mode 100644 index 00000000000..7fa3cc4e629 --- /dev/null +++ b/tests/prescribed_solution_composition_initial_condition/screen-output @@ -0,0 +1,155 @@ + +Number of active cells: 64 (on 4 levels) +Number of degrees of freedom: 1,526 (578+81+289+289+289) + +*** 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 (GMG)... 11+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00000 + +*** Timestep 1: t=0.0625 seconds, dt=0.0625 seconds + Solving temperature system... 10 iterations. + Solving C_1 system ... 9 iterations. + Solving C_2 system ... 10 iterations. + Solving Stokes system (GMG)... 7+0 iterations. + + Postprocessing: + +*** Timestep 2: t=0.125 seconds, dt=0.0625 seconds + Solving temperature system... 11 iterations. + Solving C_1 system ... 12 iterations. + Solving C_2 system ... 11 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00001 + +*** Timestep 3: t=0.1875 seconds, dt=0.0625 seconds + Solving temperature system... 12 iterations. + Solving C_1 system ... 12 iterations. + Solving C_2 system ... 11 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + + Postprocessing: + +*** Timestep 4: t=0.25 seconds, dt=0.0625 seconds + Solving temperature system... 12 iterations. + Solving C_1 system ... 11 iterations. + Solving C_2 system ... 11 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00002 + +*** Timestep 5: t=0.3125 seconds, dt=0.0625 seconds + Solving temperature system... 12 iterations. + Solving C_1 system ... 10 iterations. + Solving C_2 system ... 11 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00003 + +*** Timestep 6: t=0.375 seconds, dt=0.0625 seconds + Solving temperature system... 11 iterations. + Solving C_1 system ... 10 iterations. + Solving C_2 system ... 10 iterations. + Solving Stokes system (GMG)... 11+0 iterations. + + Postprocessing: + +*** Timestep 7: t=0.4375 seconds, dt=0.0625 seconds + Solving temperature system... 12 iterations. + Solving C_1 system ... 10 iterations. + Solving C_2 system ... 10 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00004 + +*** Timestep 8: t=0.5 seconds, dt=0.0625 seconds + Solving temperature system... 12 iterations. + Solving C_1 system ... 10 iterations. + Solving C_2 system ... 11 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00005 + +*** Timestep 9: t=0.5625 seconds, dt=0.0625 seconds + Solving temperature system... 11 iterations. + Solving C_1 system ... 9 iterations. + Solving C_2 system ... 10 iterations. + Solving Stokes system (GMG)... 11+0 iterations. + + Postprocessing: + +*** Timestep 10: t=0.625 seconds, dt=0.0625 seconds + Solving temperature system... 12 iterations. + Solving C_1 system ... 9 iterations. + Solving C_2 system ... 10 iterations. + Solving Stokes system (GMG)... 11+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00006 + +*** Timestep 11: t=0.6875 seconds, dt=0.0625 seconds + Solving temperature system... 11 iterations. + Solving C_1 system ... 9 iterations. + Solving C_2 system ... 10 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + + Postprocessing: + +*** Timestep 12: t=0.75 seconds, dt=0.0625 seconds + Solving temperature system... 13 iterations. + Solving C_1 system ... 11 iterations. + Solving C_2 system ... 12 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00007 + +*** Timestep 13: t=0.8125 seconds, dt=0.0625 seconds + Solving temperature system... 11 iterations. + Solving C_1 system ... 9 iterations. + Solving C_2 system ... 10 iterations. + Solving Stokes system (GMG)... 6+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00008 + +*** Timestep 14: t=0.875 seconds, dt=0.0625 seconds + Solving temperature system... 11 iterations. + Solving C_1 system ... 9 iterations. + Solving C_2 system ... 11 iterations. + Solving Stokes system (GMG)... 5+0 iterations. + + Postprocessing: + +*** Timestep 15: t=0.9375 seconds, dt=0.0625 seconds + Solving temperature system... 11 iterations. + Solving C_1 system ... 9 iterations. + Solving C_2 system ... 11 iterations. + Solving Stokes system (GMG)... 5+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00009 + +*** Timestep 16: t=1 seconds, dt=0.0625 seconds + Solving temperature system... 12 iterations. + Solving C_1 system ... 9 iterations. + Solving C_2 system ... 11 iterations. + Solving Stokes system (GMG)... 5+0 iterations. + + Postprocessing: + Writing graphical output: output-prescribed_solution_composition_initial_condition/solution/solution-00010 + +Termination requested by criterion: end time + + +