Skip to content

Commit e10b32f

Browse files
lhy11009lhy11009
committed
Prescribed Solution for initial temperature condition
1 parent b9a61a5 commit e10b32f

File tree

5 files changed

+470
-0
lines changed

5 files changed

+470
-0
lines changed

‎doc/modules/changes/20260307_lhy11009‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
Added: a new instance of 'temperature from intial' in the prescribed solution plugin system.
2+
<br>
3+
(Haoyuan Li, 2025/12/30)

‎include/aspect/prescribed_solution/temperature_from_initial.h‎

Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
/*
2+
Copyright (C) 2011 - 2026 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#ifndef _aspect_prescribed_solution_temperature_from_initial_h
22+
#define _aspect_prescribed_solution_temperature_from_initial_h
23+
24+
#include <aspect/initial_temperature/interface.h>
25+
#include <aspect/prescribed_solution/interface.h>
26+
#include <aspect/simulator_access.h>
27+
#include <deal.II/base/parsed_function.h>
28+
29+
namespace aspect
30+
{
31+
namespace PrescribedSolution
32+
{
33+
/**
34+
* A class that prescribes temperature in a selected region using
35+
* the value returned by the active initial temperature model.
36+
*
37+
* The region is selected through an indicator function. Where the
38+
* indicator is greater than 0.5, the temperature is constrained
39+
* to the initial temperature value evaluated at that position.
40+
*
41+
* @ingroup PrescribedSolution
42+
*/
43+
template <int dim>
44+
class TemperatureFromInitial
45+
: public Interface<dim>,
46+
public SimulatorAccess<dim>
47+
{
48+
public:
49+
/**
50+
* Constructor.
51+
*/
52+
TemperatureFromInitial ();
53+
54+
/**
55+
* Store a shared pointer to the initial temperature manager so the
56+
* plugin can safely access it after initialization.
57+
*/
58+
void initialize_simulator (const Simulator<dim> &simulator_object) override;
59+
60+
/**
61+
* Update the current time in the indicator function.
62+
*/
63+
void update () override;
64+
65+
/**
66+
* Declare the parameters this class takes through input files.
67+
*/
68+
static void declare_parameters (ParameterHandler &prm);
69+
70+
/**
71+
* Read the parameters this class declares from the parameter file.
72+
*/
73+
void parse_parameters (ParameterHandler &prm) override;
74+
75+
/**
76+
* Decide and assign cell-wise constraints for temperature DoFs.
77+
*/
78+
void constrain_solution (
79+
const typename DoFHandler<dim>::active_cell_iterator &cell,
80+
const std::vector<Point<dim>> &positions,
81+
const std::vector<unsigned int> &component_indices,
82+
std::vector<bool> &should_be_constrained,
83+
std::vector<double> &solution) override;
84+
85+
private:
86+
/**
87+
* Indicator function for selecting where the temperature should
88+
* be prescribed.
89+
*/
90+
Functions::ParsedFunction<dim> prescribed_temperature_indicator_function;
91+
92+
/**
93+
* The coordinate representation used to evaluate the indicator
94+
* function. Possible choices are cartesian, spherical, and depth.
95+
*/
96+
Utilities::Coordinates::CoordinateSystem coordinate_system;
97+
98+
99+
/**
100+
* Shared pointer to the initial temperature manager. We keep this
101+
* alive because the simulator may release its own pointer after
102+
* initialization.
103+
*/
104+
std::shared_ptr<const InitialTemperature::Manager<dim>> initial_temperature_manager;
105+
};
106+
}
107+
}
108+
109+
#endif

‎source/prescribed_solution/temperature_from_initial.cc‎

Lines changed: 168 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,168 @@
1+
/*
2+
Copyright (C) 2011 - 2026 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#include <aspect/initial_temperature/interface.h>
22+
#include <aspect/prescribed_solution/temperature_from_initial.h>
23+
24+
namespace aspect
25+
{
26+
namespace PrescribedSolution
27+
{
28+
template <int dim>
29+
TemperatureFromInitial<dim>::TemperatureFromInitial ()
30+
:
31+
prescribed_temperature_indicator_function(1)
32+
{}
33+
34+
35+
36+
template <int dim>
37+
void
38+
TemperatureFromInitial<dim>::initialize_simulator (const Simulator<dim> &simulator_object)
39+
{
40+
SimulatorAccess<dim>::initialize_simulator(simulator_object);
41+
42+
initial_temperature_manager =
43+
this->get_initial_temperature_manager_pointer();
44+
45+
AssertThrow(initial_temperature_manager,
46+
ExcMessage("TemperatureFromInitial: Failed to obtain a valid "
47+
"shared pointer to the initial temperature manager."));
48+
}
49+
50+
template <int dim>
51+
void
52+
TemperatureFromInitial<dim>::update ()
53+
{
54+
if (this->convert_output_to_years())
55+
prescribed_temperature_indicator_function.set_time(this->get_time() / year_in_seconds);
56+
else
57+
prescribed_temperature_indicator_function.set_time(this->get_time());
58+
}
59+
60+
template <int dim>
61+
void
62+
TemperatureFromInitial<dim>::declare_parameters (ParameterHandler &prm)
63+
{
64+
prm.enter_subsection("Prescribed solution");
65+
{
66+
prm.enter_subsection("Temperature from initial");
67+
{
68+
prm.declare_entry("Coordinate system",
69+
"cartesian",
70+
Patterns::Selection("cartesian|spherical|depth"),
71+
"A selection that determines the assumed coordinate "
72+
"system for the indicator function variables. "
73+
"Allowed values are `cartesian', `spherical', and `depth'.");
74+
75+
prm.enter_subsection("Indicator function");
76+
{
77+
Functions::ParsedFunction<dim>::declare_parameters(prm, 1);
78+
}
79+
prm.leave_subsection();
80+
}
81+
prm.leave_subsection();
82+
}
83+
prm.leave_subsection();
84+
}
85+
86+
template <int dim>
87+
void
88+
TemperatureFromInitial<dim>::parse_parameters (ParameterHandler &prm)
89+
{
90+
prm.enter_subsection("Prescribed solution");
91+
{
92+
prm.enter_subsection("Temperature from initial");
93+
{
94+
coordinate_system =
95+
Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system"));
96+
97+
prm.enter_subsection("Indicator function");
98+
{
99+
try
100+
{
101+
prescribed_temperature_indicator_function.parse_parameters(prm);
102+
}
103+
catch (...)
104+
{
105+
std::cerr << "ERROR: FunctionParser failed to parse\n"
106+
<< "\t'Prescribed solution.Temperature from initial.Indicator function'\n"
107+
<< "with expression\n"
108+
<< "\t'" << prm.get("Function expression") << "'\n";
109+
throw;
110+
}
111+
}
112+
prm.leave_subsection();
113+
}
114+
prm.leave_subsection();
115+
}
116+
prm.leave_subsection();
117+
}
118+
119+
template <int dim>
120+
void
121+
TemperatureFromInitial<dim>::constrain_solution (
122+
const typename DoFHandler<dim>::active_cell_iterator &/*cell*/,
123+
const std::vector<Point<dim>> &positions,
124+
const std::vector<unsigned int> &component_indices,
125+
std::vector<bool> &should_be_constrained,
126+
std::vector<double> &solution)
127+
{
128+
129+
const unsigned int temperature_index =
130+
this->introspection().component_indices.temperature;
131+
132+
for (unsigned int q=0; q<positions.size(); ++q)
133+
{
134+
if (component_indices[q] != temperature_index)
135+
continue;
136+
137+
const auto point =
138+
this->get_geometry_model().cartesian_to_other_coordinates(positions[q], coordinate_system);
139+
140+
const double indicator =
141+
prescribed_temperature_indicator_function.value(Utilities::convert_array_to_point<dim>(point.get_coordinates()), 0);
142+
143+
if (indicator > 0.5)
144+
{
145+
should_be_constrained[q] = true;
146+
solution[q] = initial_temperature_manager->initial_temperature(positions[q]);
147+
}
148+
}
149+
}
150+
}
151+
}
152+
153+
154+
// explicit instantiations
155+
namespace aspect
156+
{
157+
namespace PrescribedSolution
158+
{
159+
ASPECT_REGISTER_PRESCRIBED_SOLUTION(
160+
TemperatureFromInitial,
161+
"temperature from initial",
162+
"Prescribe the temperature in a selected region using the active "
163+
"initial temperature model. The selected region is defined through "
164+
"an indicator function. Where the indicator is greater than 0.5, "
165+
"the temperature is constrained to the initial temperature evaluated "
166+
"at that position.")
167+
}
168+
}

0 commit comments

Comments
 (0)