Skip to content

Commit 23b9afd

Browse files
lhy11009lhy11009
committed
Prescribed condition from initial composition
1 parent b9a61a5 commit 23b9afd

File tree

5 files changed

+504
-0
lines changed

5 files changed

+504
-0
lines changed

‎doc/modules/changes/20260311_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 'initial composition' in the prescribed solution plugin system.
2+
<br>
3+
(Haoyuan Li, 2026/03/11)

‎include/aspect/prescribed_solution/initial_composition.h‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
/*
2+
Prescribe compositional fields using the initial composition model
3+
within a region defined by an indicator function.
4+
*/
5+
6+
#ifndef _aspect_prescribed_solution_initial_composition_h
7+
#define _aspect_prescribed_solution_initial_composition_h
8+
9+
#include <aspect/initial_composition/interface.h>
10+
#include <aspect/prescribed_solution/interface.h>
11+
#include <aspect/simulator_access.h>
12+
13+
#include <deal.II/base/parsed_function.h>
14+
15+
namespace aspect
16+
{
17+
namespace PrescribedSolution
18+
{
19+
/**
20+
* Prescribe compositional fields using the initial composition model.
21+
* The region where the constraint applies is determined by an
22+
* indicator function.
23+
*/
24+
template <int dim>
25+
class InitialComposition
26+
: public Interface<dim>,
27+
public SimulatorAccess<dim>
28+
{
29+
public:
30+
31+
InitialComposition ();
32+
33+
/**
34+
* Store shared pointer to initial composition manager.
35+
*/
36+
void initialize () override;
37+
38+
/**
39+
* Update function time.
40+
*/
41+
void update () override;
42+
43+
/**
44+
* Declare parameters.
45+
*/
46+
static void declare_parameters (ParameterHandler &prm);
47+
48+
/**
49+
* Parse parameters.
50+
*/
51+
void parse_parameters (ParameterHandler &prm) override;
52+
53+
/**
54+
* Constrain compositional solution.
55+
*/
56+
void constrain_solution (const typename DoFHandler<dim>::active_cell_iterator &cell,
57+
const std::vector<Point<dim>> &positions,
58+
const std::vector<unsigned int> &component_indices,
59+
std::vector<bool> &should_be_constrained,
60+
std::vector<double> &solution) override;
61+
62+
private:
63+
64+
/**
65+
* Indicator function defining region where composition
66+
* is prescribed.
67+
*/
68+
Functions::ParsedFunction<dim> indicator_function;
69+
70+
/**
71+
* Coordinate system used for evaluating indicator.
72+
*/
73+
Utilities::Coordinates::CoordinateSystem coordinate_system;
74+
75+
/**
76+
* Pointer to initial composition manager.
77+
*/
78+
std::shared_ptr<const aspect::InitialComposition::Manager<dim>> initial_composition_manager;
79+
};
80+
}
81+
}
82+
83+
#endif

‎source/prescribed_solution/initial_composition.cc‎

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
#include <aspect/initial_composition/interface.h>
2+
#include <aspect/prescribed_solution/initial_composition.h>
3+
4+
namespace aspect
5+
{
6+
namespace PrescribedSolution
7+
{
8+
9+
template <int dim>
10+
InitialComposition<dim>::InitialComposition ()
11+
:
12+
indicator_function(1)
13+
{}
14+
15+
16+
17+
template <int dim>
18+
void
19+
InitialComposition<dim>::initialize ()
20+
{
21+
initial_composition_manager =
22+
this->get_initial_composition_manager_pointer();
23+
}
24+
25+
26+
27+
template <int dim>
28+
void
29+
InitialComposition<dim>::update ()
30+
{
31+
if (this->convert_output_to_years())
32+
indicator_function.set_time(this->get_time() / year_in_seconds);
33+
else
34+
indicator_function.set_time(this->get_time());
35+
}
36+
37+
38+
39+
template <int dim>
40+
void
41+
InitialComposition<dim>::constrain_solution (const typename DoFHandler<dim>::active_cell_iterator &,
42+
const std::vector<Point<dim>> &positions,
43+
const std::vector<unsigned int> &component_indices,
44+
std::vector<bool> &should_be_constrained,
45+
std::vector<double> &solution)
46+
{
47+
// Determine the component range corresponding to compositional fields.
48+
// These component indices are originally mapped from local DoFs and
49+
// determine which DoFs in the system belong to compositional fields.
50+
const unsigned int first_comp =
51+
this->introspection().component_indices.compositional_fields[0];
52+
53+
const unsigned int last_comp =
54+
this->introspection().component_indices.compositional_fields[this->introspection().n_compositional_fields-1];
55+
56+
const auto &geometry_model = this->get_geometry_model();
57+
58+
for (unsigned int q=0; q<positions.size(); ++q)
59+
{
60+
const unsigned int component = component_indices[q];
61+
62+
if (component < first_comp || component > last_comp)
63+
continue;
64+
65+
const auto point = geometry_model.cartesian_to_other_coordinates(
66+
positions[q], coordinate_system);
67+
68+
const double indicator = indicator_function.value(
69+
Utilities::convert_array_to_point<dim>(point.get_coordinates()));
70+
71+
if (indicator > 0.5)
72+
{
73+
const unsigned int field_index = component - first_comp;
74+
75+
solution[q] =
76+
initial_composition_manager->initial_composition(
77+
positions[q], field_index);
78+
79+
should_be_constrained[q] = true;
80+
}
81+
}
82+
}
83+
84+
85+
86+
template <int dim>
87+
void
88+
InitialComposition<dim>::declare_parameters (ParameterHandler &prm)
89+
{
90+
prm.enter_subsection("Prescribed solution");
91+
{
92+
prm.enter_subsection("Initial composition");
93+
{
94+
95+
prm.declare_entry("Coordinate system",
96+
"cartesian",
97+
Patterns::Selection("cartesian|spherical|depth"),
98+
"Coordinate system used for evaluating the indicator function.");
99+
100+
prm.enter_subsection("Indicator function");
101+
{
102+
Functions::ParsedFunction<dim>::declare_parameters(prm,1);
103+
}
104+
prm.leave_subsection();
105+
106+
}
107+
prm.leave_subsection();
108+
}
109+
prm.leave_subsection();
110+
}
111+
112+
113+
114+
template <int dim>
115+
void
116+
InitialComposition<dim>::parse_parameters (ParameterHandler &prm)
117+
{
118+
prm.enter_subsection("Prescribed solution");
119+
{
120+
prm.enter_subsection("Initial composition");
121+
{
122+
123+
coordinate_system =
124+
Utilities::Coordinates::string_to_coordinate_system(
125+
prm.get("Coordinate system"));
126+
127+
prm.enter_subsection("Indicator function");
128+
{
129+
indicator_function.parse_parameters(prm);
130+
}
131+
prm.leave_subsection();
132+
133+
}
134+
prm.leave_subsection();
135+
}
136+
prm.leave_subsection();
137+
}
138+
}
139+
}
140+
141+
142+
143+
namespace aspect
144+
{
145+
namespace PrescribedSolution
146+
{
147+
148+
ASPECT_REGISTER_PRESCRIBED_SOLUTION(InitialComposition,
149+
"initial composition",
150+
"Prescribe compositional fields from the initial composition model "
151+
"within a region defined by an indicator function.")
152+
153+
}
154+
}

‎tests/prescribed_solution_composition_initial_condition.prm‎

Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
#########################################################
2+
# This is modified from the composition_passive.prm
3+
# parameter file by prescribing the compositional solution
4+
# from initial compositional fields.
5+
6+
set Dimension = 2
7+
set Start time = 0
8+
set End time = 1
9+
set Use years instead of seconds = false
10+
11+
subsection Geometry model
12+
set Model name = box
13+
14+
subsection Box
15+
set X extent = 2
16+
set Y extent = 1
17+
end
18+
end
19+
20+
subsection Boundary temperature model
21+
set Fixed temperature boundary indicators = 2, 3
22+
set List of model names = box
23+
24+
subsection Box
25+
set Bottom temperature = 1
26+
set Top temperature = 0
27+
end
28+
end
29+
30+
subsection Boundary velocity model
31+
set Tangential velocity boundary indicators = 0, 1, 2
32+
set Prescribed velocity boundary indicators = 3: function
33+
34+
subsection Function
35+
set Variable names = x,z,t
36+
set Function constants = pi=3.1415926
37+
set Function expression = if(x>1+sin(0.5*pi*t), 1, -1); 0
38+
end
39+
end
40+
41+
subsection Gravity model
42+
set Model name = vertical
43+
end
44+
45+
subsection Initial temperature model
46+
set Model name = function
47+
48+
subsection Function
49+
set Variable names = x,z
50+
set Function expression = (1-z)
51+
end
52+
end
53+
54+
subsection Material model
55+
set Model name = simple
56+
57+
subsection Simple model
58+
set Thermal conductivity = 1e-6
59+
set Thermal expansion coefficient = 1e-4
60+
set Viscosity = 1
61+
end
62+
end
63+
64+
subsection Mesh refinement
65+
set Initial adaptive refinement = 0
66+
set Initial global refinement = 3
67+
set Time steps between mesh refinement = 0
68+
end
69+
70+
# The next section prescribes the composition inside the domain
71+
# using a spatially and temporally varying function. An indicator
72+
# function selects the region where the composition is constrained.
73+
subsection Prescribed solution
74+
set List of model names = initial composition
75+
subsection Initial composition
76+
subsection Indicator function
77+
set Variable names = x, y
78+
set Function expression = (x<1) ? 1:0
79+
end
80+
end
81+
end
82+
83+
subsection Postprocess
84+
set List of postprocessors = visualization
85+
86+
subsection Visualization
87+
set Interpolate output = false
88+
set Time between graphical output = 0.1
89+
set Output format = gnuplot
90+
end
91+
end
92+
93+
# There will be two compositional fields that will be
94+
# advected along. Their initial conditions are given by
95+
# a function that is one for the lowermost 0.2 height
96+
# units of the domain and zero otherwise in the first case,
97+
# and one in the top most 0.2 height units in the latter.
98+
subsection Compositional fields
99+
set Number of fields = 2
100+
end
101+
102+
subsection Initial composition model
103+
set Model name = function
104+
105+
subsection Function
106+
set Variable names = x,y
107+
set Function expression = if(y<0.2, 1, 0) ; if(y>0.8, 1, 0)
108+
end
109+
end

0 commit comments

Comments
 (0)