Skip to content

Commit d976c6d

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

File tree

5 files changed

+460
-0
lines changed

5 files changed

+460
-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 initial' in the prescribed solution plugin system.
2+
<br>
3+
(Haoyuan Li, 2025/12/30)

‎include/aspect/prescribed_solution/initial_temperature.h‎

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
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 InitialTemperature
45+
: public Interface<dim>,
46+
public SimulatorAccess<dim>
47+
{
48+
public:
49+
/**
50+
* Constructor.
51+
*/
52+
InitialTemperature ();
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+
* Shared pointer to the initial temperature manager. We keep this
100+
* alive because the simulator may release its own pointer after
101+
* initialization.
102+
*/
103+
std::shared_ptr<const aspect::InitialTemperature::Manager<dim>> initial_temperature_manager;
104+
};
105+
}
106+
}
107+
108+
#endif

‎source/prescribed_solution/initial_temperature.cc‎

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
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/initial_temperature.h>
23+
#include <aspect/geometry_model/interface.h>
24+
25+
namespace aspect
26+
{
27+
namespace PrescribedSolution
28+
{
29+
template <int dim>
30+
InitialTemperature<dim>::InitialTemperature ()
31+
:
32+
prescribed_temperature_indicator_function(1)
33+
{}
34+
35+
36+
37+
template <int dim>
38+
void
39+
InitialTemperature<dim>::initialize_simulator (const Simulator<dim> &simulator_object)
40+
{
41+
SimulatorAccess<dim>::initialize_simulator(simulator_object);
42+
43+
initial_temperature_manager =
44+
this->get_initial_temperature_manager_pointer();
45+
}
46+
47+
48+
49+
template <int dim>
50+
void
51+
InitialTemperature<dim>::update ()
52+
{
53+
if (this->convert_output_to_years())
54+
prescribed_temperature_indicator_function.set_time(this->get_time() / year_in_seconds);
55+
else
56+
prescribed_temperature_indicator_function.set_time(this->get_time());
57+
}
58+
59+
60+
61+
template <int dim>
62+
void
63+
InitialTemperature<dim>::declare_parameters (ParameterHandler &prm)
64+
{
65+
prm.enter_subsection("Prescribed solution");
66+
{
67+
prm.enter_subsection("Initial temperature");
68+
{
69+
prm.declare_entry("Coordinate system",
70+
"cartesian",
71+
Patterns::Selection("cartesian|spherical|depth"),
72+
"A selection that determines the assumed coordinate "
73+
"system for the indicator function variables. "
74+
"Allowed values are `cartesian', `spherical', and `depth'.");
75+
76+
prm.enter_subsection("Indicator function");
77+
{
78+
Functions::ParsedFunction<dim>::declare_parameters(prm, 1);
79+
}
80+
prm.leave_subsection();
81+
}
82+
prm.leave_subsection();
83+
}
84+
prm.leave_subsection();
85+
}
86+
87+
88+
89+
template <int dim>
90+
void
91+
InitialTemperature<dim>::parse_parameters (ParameterHandler &prm)
92+
{
93+
prm.enter_subsection("Prescribed solution");
94+
{
95+
prm.enter_subsection("Initial temperature");
96+
{
97+
coordinate_system =
98+
Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system"));
99+
100+
prm.enter_subsection("Indicator function");
101+
{
102+
try
103+
{
104+
prescribed_temperature_indicator_function.parse_parameters(prm);
105+
}
106+
catch (...)
107+
{
108+
std::cerr << "ERROR: FunctionParser failed to parse\n"
109+
<< "\t'Prescribed solution.Initial temperature.Indicator function'\n"
110+
<< "with expression\n"
111+
<< "\t'" << prm.get("Function expression") << "'\n";
112+
throw;
113+
}
114+
}
115+
prm.leave_subsection();
116+
}
117+
prm.leave_subsection();
118+
}
119+
prm.leave_subsection();
120+
}
121+
122+
123+
124+
template <int dim>
125+
void
126+
InitialTemperature<dim>::constrain_solution (
127+
const typename DoFHandler<dim>::active_cell_iterator &/*cell*/,
128+
const std::vector<Point<dim>> &positions,
129+
const std::vector<unsigned int> &component_indices,
130+
std::vector<bool> &should_be_constrained,
131+
std::vector<double> &solution)
132+
{
133+
134+
const unsigned int temperature_index =
135+
this->introspection().component_indices.temperature;
136+
137+
for (unsigned int q=0; q<positions.size(); ++q)
138+
{
139+
if (component_indices[q] != temperature_index)
140+
continue;
141+
142+
const auto point =
143+
this->get_geometry_model().cartesian_to_other_coordinates(positions[q], coordinate_system);
144+
145+
const double indicator =
146+
prescribed_temperature_indicator_function.value(Utilities::convert_array_to_point<dim>(point.get_coordinates()), 0);
147+
148+
if (indicator > 0.5)
149+
{
150+
should_be_constrained[q] = true;
151+
solution[q] = initial_temperature_manager->initial_temperature(positions[q]);
152+
}
153+
}
154+
}
155+
}
156+
}
157+
158+
159+
160+
// explicit instantiations
161+
namespace aspect
162+
{
163+
namespace PrescribedSolution
164+
{
165+
ASPECT_REGISTER_PRESCRIBED_SOLUTION(
166+
InitialTemperature,
167+
"initial temperature",
168+
"Prescribe the temperature in a selected region using the active "
169+
"initial temperature model. The selected region is defined through "
170+
"an indicator function. At locations where the indicator value is greater than 0.5, "
171+
"the temperature is constrained to the initial temperature evaluated "
172+
"at that position.")
173+
}
174+
}

0 commit comments

Comments
 (0)