gas (the so-called Euler equations of gas dynamics) are
a basic system of conservation laws. In spatial dimension $d$ they read
@f[
-\partial_t \mathbf{w} + \nabla \cdot \mathbf{F}(\mathbf{w}) = \mathbf{0},
+\partial_t \mathbf{w} + \nabla \cdot \mathbf{F}(\mathbf{w}) =
+\mathbf{G}(\mathbf w),
@f]
with the solution $\mathbf{w}=(\rho v_1,\ldots,\rho v_d,\rho,
E)^{\top}$ consisting of $\rho$ the fluid density, ${\mathbf v}=(v_1,\ldots v_d)^T$ the
flow velocity (and thus $\rho\mathbf v$ being the linear momentum
density), and
$E$ the energy density of the gas. We interpret the equations above as
-$\partial_t \mathbf{w}_i + \nabla \cdot \mathbf{F}_i(\mathbf{w}) = 0$, $i=1,\ldots,dim+2$.
+$\partial_t \mathbf{w}_i + \nabla \cdot \mathbf{F}_i(\mathbf{w}) = \mathbf
+G_i(\mathbf w)$, $i=1,\ldots,dim+2$.
For the Euler equations, the flux matrix $\mathbf F$ (or system of flux functions)
is defined as (shown here for the case $d=3$)
@f{eqnarray*}
- \mathbf F
+ \mathbf F(\mathbf w)
=
\left(
\begin{array}{ccc}
\end{array}
\right),
@f}
-such that the entire system of equations reads:
+and we will choose as particular right hand side forcing only the effects of
+gravity, described by
+@f{eqnarray*}
+ \mathbf G(\mathbf w)
+ =
+ \left(
+ \begin{array}{c}
+ g_1\rho \\
+ g_2\rho \\
+ g_3\rho \\
+ 0 \\
+ \rho \mathbf g \cdot \mathbf v
+ \end{array}
+ \right),
+@f}
+where $\mathbf g=(g_1,g_2,g_3)^T$ denotes the gravity vector.
+With this, the entire system of equations reads:
@f{eqnarray*}
- \partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0, \\
\partial_t (\rho v_i) + \sum_{s=1}^d \frac{\partial(\rho v_i v_s +
- \delta_{is} p)}{\partial x_s} &=& 0, \qquad i=1,\dots,d, \\
- \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=& 0.
+ \delta_{is} p)}{\partial x_s} &=& g_i \rho, \qquad i=1,\dots,d, \\
+ \partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0, \\
+ \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=&
+ \rho \mathbf g \cdot \mathbf v.
@f}
These equations describe, respectively, the conservation of momentum,
mass, and energy.
\mathbf{z}\right)
- \int_{\Omega} \left(\mathbf{F}(\tilde{\mathbf{w}}),
\nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z})
+\\
+&& \qquad
+
\int_{\partial \Omega} \left(\mathbf{H}(\tilde{\mathbf{w}}^+),
\mathbf{w}^-(\tilde{\mathbf{w}}^+), \mathbf{n}), \mathbf{z}\right)
+-
+\int_{\partial \Omega} \left(\mathbf{G}(\mathbf{w}),
+\mathbf{z}\right)
\\
& = & 0
@f}
// well:
template <typename InputVector, typename number>
static
- void flux_matrix (const InputVector &W,
- number (&flux)[n_components][dim])
+ void compute_flux_matrix (const InputVector &W,
+ number (&flux)[n_components][dim])
{
// First compute the pressure that
// appears in the flux matrix, and
Sacado::Fad::DFad<double> iflux[n_components][dim];
Sacado::Fad::DFad<double> oflux[n_components][dim];
- flux_matrix (Wplus, iflux);
- flux_matrix (Wminus, oflux);
+ compute_flux_matrix (Wplus, iflux);
+ compute_flux_matrix (Wminus, oflux);
for (unsigned int di=0; di<n_components; ++di)
{
}
}
+
+ template <typename InputVector, typename number>
+ static
+ void compute_forcing_vector (const InputVector &W,
+ number (&forcing)[n_components])
+ {
+ const double gravity = -1.0;
+
+ for (unsigned int c=0; c<n_components; ++c)
+ switch (c)
+ {
+ case first_momentum_component+dim-1:
+ forcing[c] = gravity * W[density_component];
+ break;
+ case energy_component:
+ forcing[c] = gravity *
+ W[density_component] *
+ W[first_momentum_component+dim-1];
+ break;
+ default:
+ forcing[c] = 0;
+ }
+ }
+
+
// Finally, we declare a class that
// implements a postprocessing of data
// \left(\mathbf{F}(\tilde{\mathbf{w}}),
// \nabla\mathbf{z}_i\right)_K +
// h^{\eta}(\nabla \mathbf{w} , \nabla
+ // \mathbf{z}_i)_K -
+ // (\mathbf{G}(\tilde{\mathbf w}),
// \mathbf{z}_i)_K$ where $\tilde{\mathbf w}$
// is represented by the variable
// <code>W_theta</code>, $\mathbf{z}_i$ is
// Next, in order to compute the cell
// contributions, we need to evaluate
- // $F(\tilde{\mathbf w})$ at all quadrature
+ // $F(\tilde{\mathbf w})$ and
+ // $G(\tilde{\mathbf w})$ at all quadrature
// points. To store these, we also need to
// allocate a bit of memory. Note that we
- // compute the flux matrices in terms of
- // autodifferentiation variables, so that
- // the Jacobian contributions can later
- // easily be computed from it:
+ // compute the flux matrices and right hand
+ // sides in terms of autodifferentiation
+ // variables, so that the Jacobian
+ // contributions can later easily be
+ // computed from it:
typedef Sacado::Fad::DFad<double> FluxMatrix[EulerEquations<dim>::n_components][dim];
FluxMatrix *flux = new FluxMatrix[n_q_points];
+
+ typedef Sacado::Fad::DFad<double> ForcingVector[EulerEquations<dim>::n_components];
+ ForcingVector *forcing = new ForcingVector[n_q_points];
for (unsigned int q=0; q<n_q_points; ++q)
- EulerEquations<dim>::flux_matrix (W_theta[q], flux[q]);
+ {
+ EulerEquations<dim>::compute_flux_matrix (W_theta[q], flux[q]);
+ EulerEquations<dim>::compute_forcing_vector (W_theta[q], forcing[q]);
+ }
// We now have all of the pieces in place,
// \mathbf{w}_{\text{component\_i}}}{\partial
// x_d} , \frac{\partial
// (\mathbf{z}_i)_{\text{component\_i}}}{\partial
- // x_d} \right)_K$, where integrals are
- // understood to be evaluated through
- // summation over quadrature points.
+ // x_d} \right)_K$
+ // $-(\mathbf{G}(\tilde{\mathbf{w}}
+ // )_{\text{component\_i}},
+ // (\mathbf{z}_i)_{\text{component\_i}})_K$,
+ // where integrals are understood to be
+ // evaluated through summation over
+ // quadrature points.
//
// We initialy sum all contributions of the
// residual in the positive sense, so that
fe_v.shape_grad_component(i, point, component_i)[d] *
fe_v.JxW(point);
- for (unsigned int d = 0; d < dim; d++)
+ for (unsigned int d=0; d<dim; d++)
F_i += 1.0*std::pow(fe_v.get_cell()->diameter(),
parameters.diffusion_power) *
grad_W[point][component_i][d] *
fe_v.shape_grad_component(i, point, component_i)[d] *
fe_v.JxW(point);
-
- // The gravity component only
- // enters into the energy equation
- // and into the vertical component
- // of the velocity.
- if (component_i == dim - 1)
- F_i += parameters.gravity *
- W_theta[point][EulerEquations<dim>::density_component] *
- fe_v.shape_value_component(i,point, component_i) *
- fe_v.JxW(point);
- else if (component_i == EulerEquations<dim>::energy_component)
- F_i += parameters.gravity *
- W_theta[point][EulerEquations<dim>::density_component] *
- W_theta[point][dim-1] *
- fe_v.shape_value_component(i,point, component_i) *
- fe_v.JxW(point);
+
+ F_i -= forcing[point][component_i] *
+ fe_v.shape_value_component(i, point, component_i) *
+ fe_v.JxW(point);
}
// At the end of the loop, we have to
projects / dealii-svn.git / commitdiff