-
Notifications
You must be signed in to change notification settings - Fork 232
/
inner_core_assembly.cc
133 lines (115 loc) · 6.18 KB
/
inner_core_assembly.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#include <aspect/simulator_access.h>
#include <aspect/global.h>
#include <aspect/simulator.h>
#include <aspect/simulator/assemblers/interface.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/fe/fe_values.h>
#include "inner_core_convection.cc"
namespace aspect
{
/**
* A new assembler class that implements boundary conditions for the
* normal stress and the normal velocity that take into account the
* rate of phase change (melting/freezing) at the inner-outer core
* boundary. The model is based on Deguen, Alboussiere, and Cardin
* (2013), Thermal convection in Earth’s inner core with phase change
* at its boundary. GJI, 194, 1310-133.
*
* The mechanical boundary conditions for the inner core are
* tangential stress-free and continuity of the normal stress at the
* inner-outer core boundary. For the non-dimensional equations, that
* means that we can define a 'phase change number' $\mathcal{P}$ so
* that the normal stress at the boundary is $-\mathcal{P} u_r$ with
* the radial velocity $u_r$. This number characterizes the resistance
* to phase change at the boundary, with $\mathcal{P}\rightarrow\infty$
* corresponding to infinitely slow melting/freezing (free slip
* boundary), and $\mathcal{P}\rightarrow0$ corresponding to
* instantaneous melting/freezing (zero normal stress, open boundary).
*
* In the weak form, this results in boundary conditions of the form
* of a surface integral:
* $$\int_S \mathcal{P} (\mathbf u \cdot \mathbf n) (\mathbf v \cdot \mathbf n) dS,$$,
* with the normal vector $\mathbf n$.
*
* The function value of $\mathcal{P}$ is taken from the inner core
* material model.
*/
template <int dim>
class PhaseBoundaryAssembler :
public aspect::Assemblers::Interface<dim>,
public SimulatorAccess<dim>
{
public:
virtual
void
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
{
internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>& > (scratch_base);
internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>& > (data_base);
const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.face_finite_element_values.n_quadrature_points;
const typename DoFHandler<dim>::active_cell_iterator cell (&this->get_triangulation(),
scratch.face_finite_element_values.get_cell()->level(),
scratch.face_finite_element_values.get_cell()->index(),
&this->get_dof_handler());
unsigned int face_no = numbers::invalid_unsigned_int;
for (face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
if (scratch.face_finite_element_values.get_face_index() == cell->face_index(face_no))
break;
Assert(face_no != numbers::invalid_unsigned_int,ExcInternalError());
//assemble force terms for the matrix for all boundary faces
if (cell->face(face_no)->at_boundary())
{
scratch.face_finite_element_values.reinit (cell, face_no);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double P = dynamic_cast<const MaterialModel::InnerCore<dim>&>
(this->get_material_model()).resistance_to_phase_change
.value(scratch.material_model_inputs.position[q]);
for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection
.extractors.velocities].value(i, q);
++i_stokes;
}
++i;
}
const Tensor<1,dim> normal_vector = scratch.face_finite_element_values.normal_vector(q);
const double JxW = scratch.face_finite_element_values.JxW(q);
// boundary term: P*u*n*v*n*JxW(q)
for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
data.local_matrix(i,j) += P *
scratch.phi_u[i] *
normal_vector *
scratch.phi_u[j] *
normal_vector *
JxW;
}
}
}
};
template <int dim>
void set_assemblers_phase_boundary(const SimulatorAccess<dim> &simulator_access,
Assemblers::Manager<dim> &assemblers)
{
AssertThrow (dynamic_cast<const MaterialModel::InnerCore<dim>*>
(&simulator_access.get_material_model()) != 0,
ExcMessage ("The phase boundary assembler can only be used with the "
"material model 'inner core material'!"));
PhaseBoundaryAssembler<dim> *phase_boundary_assembler = new PhaseBoundaryAssembler<dim>();
assemblers.stokes_system_on_boundary_face.push_back (std_cxx11::unique_ptr<PhaseBoundaryAssembler<dim> > (phase_boundary_assembler));
}
}
template <int dim>
void signal_connector (aspect::SimulatorSignals<dim> &signals)
{
signals.set_assemblers.connect (&aspect::set_assemblers_phase_boundary<dim>);
}
ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>,
signal_connector<3>)