diff --git a/doc/modules/changes/20260307_lhy11009 b/doc/modules/changes/20260307_lhy11009 new file mode 100644 index 00000000000..8733127577e --- /dev/null +++ b/doc/modules/changes/20260307_lhy11009 @@ -0,0 +1,3 @@ +Added: a new instance of 'initial temperature' in the prescribed solution plugin system. +
+(Haoyuan Li, 2025/12/30) diff --git a/doc/sphinx/parameters/Material_20model.md b/doc/sphinx/parameters/Material_20model.md index 0ae46de72da..7a2f7638032 100644 --- a/doc/sphinx/parameters/Material_20model.md +++ b/doc/sphinx/parameters/Material_20model.md @@ -200,7 +200,7 @@ Viscous stress may also be limited by a non-linear stress limiter that has a for The value for the components of this formula and additional parameters are read from the parameter file in subsection ’Material model/Visco Plastic’. -‘viscoelastic’: An implementation of a simple linear viscoelastic rheology that only includes the deviatoric components of elasticity. Specifically, the viscoelastic rheology only takes into account the elastic shear strength (e.g., shear modulus), while the tensile and volumetric strength (e.g., Young’s and bulk modulus) are not considered. The model is incompressible and allows specifying an arbitrary number of compositional fields, where each field represents a different rock type or component of the viscoelastic stress tensor. The stress tensor in 2d and 3d, respectively, contains 3 or 6 components. The compositional fields representing these components must be named and listed in a very specific format, which is designed to minimize mislabeling stress tensor components as distinct ’compositional rock types’ (or vice versa). For 2d models, the first six compositional fields must be labeled ’stress\_xx’, ’stress\_yy’ and ’stress\_xy’, ’stress\_xx\_old’, ’stress\_yy\_old’ and ’stress\_xy\_old’, In 3d, the first twelve compositional fields must be labeled ’stress\_xx’, ’stress\_yy’, ’stress\_zz’, ’stress\_xy’, ’stress\_xz’, ’stress\_yz’, ’stress\_xx\_old’, ’stress\_yy\_old’, ’stress\_zz\_old’, ’stress\_xy\_old’, ’stress\_xz\_old’, ’stress\_yz\_old’. +‘viscoelastic’: An implementation of a simple linear viscoelastic rheology that only includes the deviatoric components of elasticity. Specifically, the viscoelastic rheology only takes into account the elastic shear strength (e.g., shear modulus), while the tensile and volumetric strength (e.g., Young’s and bulk modulus) are not considered. The model is incompressible and allows specifying an arbitrary number of compositional fields, where each field represents a different rock type or component of the viscoelastic stress tensor. The stress tensor in 2d and 3d, respectively, contains 3 or 6 components. The compositional fields representing these components must be named and listed in a very specific format, which is designed to minimize mislabeling stress tensor components as distinct ’compositional rock types’ (or vice versa). For 2d models, the first six compositional fields of type stress must be labeled ’ve\_stress\_xx’, ’ve\_stress\_yy’ and ’ve\_stress\_xy’, ’ve\_stress\_xx\_old’, ’ve\_stress\_yy\_old’ and ’ve\_stress\_xy\_old’, In 3d, the first twelve compositional fields of type stress must be labeled ’ve\_stress\_xx’, ’ve\_stress\_yy’, ’ve\_stress\_zz’, ’ve\_stress\_xy’, ’ve\_stress\_xz’, ’ve\_stress\_yz’, ’ve\_stress\_xx\_old’, ’ve\_stress\_yy\_old’, ’ve\_stress\_zz\_old’, ’ve\_stress\_xy\_old’, ’ve\_stress\_xz\_old’, ’ve\_stress\_yz\_old’. Expanding the model to include non-linear viscous flow (e.g., diffusion/dislocation creep) and plasticity would produce a constitutive relationship commonly referred to as partial elastoviscoplastic (e.g., pEVP) in the geodynamics community. While extensively discussed and applied within the geodynamics literature, notable references include: Moresi et al. (2003), J. Comp. Phys., v. 184, p. 476-497. Gerya and Yuen (2007), Phys. Earth. Planet. Inter., v. 163, p. 83-105. Gerya (2010), Introduction to Numerical Geodynamic Modeling. Kaus (2010), Tectonophysics, v. 484, p. 36-47. Choi et al. (2013), J. Geophys. Res., v. 118, p. 2429-2444. Keller et al. (2013), Geophys. J. Int., v. 195, p. 1406-1442. diff --git a/doc/sphinx/parameters/Prescribed_20solution.md b/doc/sphinx/parameters/Prescribed_20solution.md index 598b0179538..4742bd37279 100644 --- a/doc/sphinx/parameters/Prescribed_20solution.md +++ b/doc/sphinx/parameters/Prescribed_20solution.md @@ -9,17 +9,63 @@ :name: parameters:Prescribed_20solution/List_20of_20model_20names **Default value:** -**Pattern:** [MultipleSelection temperature function|velocity function ] +**Pattern:** [MultipleSelection 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 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. ‘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_20temperature)= +## **Subsection:** Prescribed solution / Initial temperature +::::{dropdown} __Parameter:__ {ref}`Coordinate system` +:name: parameters:Prescribed_20solution/Initial_20temperature/Coordinate_20system +**Default value:** cartesian + +**Pattern:** [Selection cartesian|spherical|depth ] + +**Documentation:** A selection that determines the assumed coordinate system for the indicator function variables. Allowed values are ‘cartesian’, ‘spherical’, and ‘depth’. +:::: + +(parameters:Prescribed_20solution/Initial_20temperature/Indicator_20function)= +## **Subsection:** Prescribed solution / Initial temperature / Indicator function +::::{dropdown} __Parameter:__ {ref}`Function constants` +:name: parameters:Prescribed_20solution/Initial_20temperature/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_20temperature/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_20temperature/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/Temperature_20function)= ## **Subsection:** Prescribed solution / Temperature function ::::{dropdown} __Parameter:__ {ref}`Coordinate system` diff --git a/include/aspect/prescribed_solution/initial_temperature.h b/include/aspect/prescribed_solution/initial_temperature.h new file mode 100644 index 00000000000..73f3ddb0194 --- /dev/null +++ b/include/aspect/prescribed_solution/initial_temperature.h @@ -0,0 +1,108 @@ +/* + Copyright (C) 2011 - 2026 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_temperature_h +#define _aspect_prescribed_solution_initial_temperature_h + +#include +#include +#include +#include + +namespace aspect +{ + namespace PrescribedSolution + { + /** + * A class that prescribes temperature in a selected region using + * the value returned by the active initial temperature model. + * + * The region is selected through an indicator function. Where the + * indicator is greater than 0.5, the temperature is constrained + * to the initial temperature value evaluated at that position. + * + * @ingroup PrescribedSolution + */ + template + class InitialTemperature + : public Interface, + public SimulatorAccess + { + public: + /** + * Constructor. + */ + InitialTemperature (); + + /** + * Store a shared pointer to the initial temperature manager so the + * plugin can safely access it after initialization. + */ + void initialize () override; + + /** + * Update the current time in the indicator function. + */ + void update () override; + + /** + * Declare the parameters this class takes through input files. + */ + static void declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + void parse_parameters (ParameterHandler &prm) override; + + /** + * Decide and assign cell-wise constraints for temperature DoFs. + */ + 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 for selecting where the temperature should + * be prescribed. + */ + Functions::ParsedFunction prescribed_temperature_indicator_function; + + /** + * The coordinate representation used to evaluate the indicator + * function. Possible choices are cartesian, spherical, and depth. + */ + Utilities::Coordinates::CoordinateSystem coordinate_system; + + /** + * Shared pointer to the initial temperature manager. We keep this + * alive because the simulator may release its own pointer after + * initialization. + */ + std::shared_ptr> initial_temperature_manager; + }; + } +} + +#endif diff --git a/source/prescribed_solution/initial_temperature.cc b/source/prescribed_solution/initial_temperature.cc new file mode 100644 index 00000000000..fb988e56223 --- /dev/null +++ b/source/prescribed_solution/initial_temperature.cc @@ -0,0 +1,170 @@ +/* + Copyright (C) 2011 - 2026 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 + InitialTemperature::InitialTemperature () + : + prescribed_temperature_indicator_function(1) + {} + + + + template + void + InitialTemperature::initialize () + { + initial_temperature_manager = + this->get_initial_temperature_manager_pointer(); + } + + + + template + void + InitialTemperature::update () + { + if (this->convert_output_to_years()) + prescribed_temperature_indicator_function.set_time(this->get_time() / year_in_seconds); + else + prescribed_temperature_indicator_function.set_time(this->get_time()); + } + + + + template + void + InitialTemperature::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) + { + + const unsigned int temperature_index = + this->introspection().component_indices.temperature; + + for (unsigned int q=0; qget_geometry_model().cartesian_to_other_coordinates(positions[q], coordinate_system); + + const double indicator = + prescribed_temperature_indicator_function.value(Utilities::convert_array_to_point(point.get_coordinates()), 0); + + if (indicator > 0.5) + { + should_be_constrained[q] = true; + solution[q] = initial_temperature_manager->initial_temperature(positions[q]); + } + } + } + + + + template + void + InitialTemperature::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Prescribed solution"); + { + prm.enter_subsection("Initial temperature"); + { + prm.declare_entry("Coordinate system", + "cartesian", + Patterns::Selection("cartesian|spherical|depth"), + "A selection that determines the assumed coordinate " + "system for the indicator function variables. " + "Allowed values are `cartesian', `spherical', and `depth'."); + + prm.enter_subsection("Indicator function"); + { + Functions::ParsedFunction::declare_parameters(prm, 1); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + + + template + void + InitialTemperature::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Prescribed solution"); + { + prm.enter_subsection("Initial temperature"); + { + coordinate_system = + Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system")); + + prm.enter_subsection("Indicator function"); + { + try + { + prescribed_temperature_indicator_function.parse_parameters(prm); + } + catch (...) + { + std::cerr << "ERROR: FunctionParser failed to parse\n" + << "\t'Prescribed solution.Initial temperature.Indicator function'\n" + << "with expression\n" + << "\t'" << prm.get("Function expression") << "'\n"; + throw; + } + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + } +} + + + +// explicit instantiations +namespace aspect +{ + namespace PrescribedSolution + { + ASPECT_REGISTER_PRESCRIBED_SOLUTION(InitialTemperature, + "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.") + } +} diff --git a/tests/prescribed_solution_temperature_initial_condition.prm b/tests/prescribed_solution_temperature_initial_condition.prm new file mode 100644 index 00000000000..431359517e5 --- /dev/null +++ b/tests/prescribed_solution_temperature_initial_condition.prm @@ -0,0 +1,145 @@ +# This test verifies implementation of prescribed solutions interface. +# It sets up a 2D convection problem in a rectangular box with +# prescribed temperature conditions from the initial temperature fields. +# A prescribed temperature indicator function is to set the region of +# prescribed temperature. + +# At the top, we define the number of space dimensions we would like to +# work in: +set Dimension = 2 + +# There are several global variables that have to do with what +# time system we want to work in and what the end time is. We +# also designate an output directory. +set Use years instead of seconds = false +set End time = 1 +set Output directory = output +set Resume computation = false + + +# Then there are variables that describe how the pressure should +# be normalized. Here, we choose a zero average pressure +# at the surface of the domain (for the current geometry, the +# surface is defined as the top boundary). +set Pressure normalization = surface +set Surface pressure = 0 +set Adiabatic surface temperature = 1573.0 + +subsection Solver parameters + set Temperature solver tolerance = 1e-10 +end + +# Then come a number of sections that deal with the setup +# of the problem to solve. The first one deals with the +# geometry of the domain within which we want to solve. +# The sections that follow all have the same basic setup +# where we select the name of a particular model (here, +# the box geometry) and then, in a further subsection, +# set the parameters that are specific to this particular +# model. +subsection Geometry model + set Model name = box + subsection Box + set X extent = 2e6 + set Y extent = 1.0000e+06 + set X repetitions = 2 + end +end + +subsection Mesh refinement + set Initial global refinement = 3 +end + +# The next section deals with the initial conditions for the +# temperature with an adiabatic condition +subsection Initial temperature model + set Model name = adiabatic +end + +# Then follows a section that describes the boundary conditions +# for the temperature. The model we choose is called 'box' and +# allows to set a constant temperature on each of the four sides +# of the box geometry. In our case, we choose something that is +# heated from below and cooled from above, whereas all other +# parts of the boundary are insulated (i.e., no heat flux through +# these boundaries; this is also often used to specify symmetry +# boundaries). +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top + set List of model names = box + subsection Box + set Bottom temperature = 2.0136e+03 + set Top temperature = 273 + end +end + +# The next section prescribes the temperature from the initial +# temperature field inside the domain told by an indicator function. +subsection Prescribed solution + set List of model names = initial temperature + subsection Initial temperature + subsection Indicator function + set Variable names = x, y + set Function expression = (x<2e6*0.75 && x>2e6*0.25 && y<1e6*0.4) ? 1:0 + end + end +end + +# The next parameters then describe on which parts of the +# boundary we prescribe a zero or nonzero velocity and +# on which parts the flow is allowed to be tangential. +# Here, all four sides of the box allow tangential +# unrestricted flow but with a zero normal component: +subsection Boundary velocity model + set Tangential velocity boundary indicators = left, right, bottom, top +end + +# The following two sections describe first the +# direction (vertical) and magnitude of gravity and the +# material model (i.e., density, viscosity, etc). We have +# discussed the settings used here in the introduction to +# this cookbook in the manual already. +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 10.0 + end +end + +subsection Material model + set Model name = simple + + subsection Simple model + set Reference density = 3400.0 + set Reference specific heat = 1250.0 + set Reference temperature = 1673.0 + set Thermal conductivity = 1e22 + set Thermal expansion coefficient = 3.1e-5 + set Viscosity = 5e21 + end +end + +# We also have to specify that we want to use the Boussinesq +# approximation (assuming the density in the temperature +# equation to be constant, and incompressibility). +subsection Formulation + set Formulation = Boussinesq approximation +end + +# The final part is to specify what ASPECT should do with the +# solution once computed at the end of every time step. The +# process of evaluating the solution is called `postprocessing' +# and we choose to compute velocity and temperature statistics, +# statistics about the heat flux through the boundaries of the +# domain, and to generate graphical output files for later +# visualization. +subsection Postprocess + set List of postprocessors = velocity statistics, temperature statistics, heat flux statistics, visualization, heat flux map + + subsection Visualization + set List of output variables = material properties, nonadiabatic pressure, strain rate, stress + set Time between graphical output = 1 + set Output format = gnuplot + end +end diff --git a/tests/prescribed_solution_temperature_initial_condition/screen-output b/tests/prescribed_solution_temperature_initial_condition/screen-output new file mode 100644 index 00000000000..36beea52ccc --- /dev/null +++ b/tests/prescribed_solution_temperature_initial_condition/screen-output @@ -0,0 +1,30 @@ + +Number of active cells: 128 (on 4 levels) +Number of degrees of freedom: 1,836 (1,122+153+561) + +*** Timestep 0: t=0 seconds, dt=0 seconds + Solving temperature system... 0 iterations. + Solving Stokes system (GMG)... 7+0 iterations. + + Postprocessing: + RMS, max velocity: 2.96e-12 m/s, 1.43e-11 m/s + Temperature min/avg/max: 273 K, 1555 K, 2014 K + Heat fluxes through boundary parts: 0 W, 0 W, -1.645e+26 W, 4.853e+26 W + Writing graphical output: output-prescribed_solution_temperature_initial_condition/solution/solution-00000 + Writing heat flux map output-prescribed_solution_temperature_initial_condition/heat_flux.00000 + +*** Timestep 1: t=1 seconds, dt=1 seconds + Solving temperature system... 36 iterations. + Solving Stokes system (GMG)... 8+0 iterations. + + Postprocessing: + RMS, max velocity: 1.62e-10 m/s, 2.5e-10 m/s + Temperature min/avg/max: 273 K, 1167 K, 2014 K + Heat fluxes through boundary parts: 0 W, 0 W, -1.011e+26 W, 3.882e+25 W + Writing graphical output: output-prescribed_solution_temperature_initial_condition/solution/solution-00001 + Writing heat flux map output-prescribed_solution_temperature_initial_condition/heat_flux.00001 + +Termination requested by criterion: end time + + +