Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

step-69: Remove boundary c_ij fix #10926

Merged
merged 2 commits into from
Sep 23, 2020
Merged

step-69: Remove boundary c_ij fix #10926

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bffc0c7
Remove boundary fix
tamiko Sep 16, 2020
8c4eb3e
Modified dox and cc. Now they are consistent.
Sep 23, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
81 changes: 50 additions & 31 deletions examples/step-69/doc/intro.dox
View file
Open in desktop
Expand Up @@ -163,9 +163,9 @@ $\mathbf{u}(\mathbf{x},t)$ remains in $\mathcal{B}$.

<h4>Variational versus collocation-type discretizations</h4>

Following Step-9, Step-12, and Step-33, at this point it might look tempting
to base a discretization of Euler's equations on a (semi-discrete) variational
formulation:
Following Step-9, Step-12, Step-33, and Step-67, at this point it might look
tempting to base a discretization of Euler's equations on a (semi-discrete)
variational formulation:
@f{align*}
(\partial_t\mathbf{u}_{h},\textbf{v}_h)_{L^2(\Omega)}
- ( \mathbb{f}(\mathbf{u}_{h}) ,\text{grad} \, \textbf{v}_{h})_{L^2(\Omega)}
Expand Down Expand Up @@ -193,7 +193,10 @@ techniques (see the early work of @cite Brooks1982 and @cite Johnson1986). They
have proven to be some of the best approaches for simulations in the subsonic
shockless regime and similarly benign situations.

However, in the transonic and supersonic regime, and shock-hydrodynamics
<!-- In particular, tutorial Step-67 focuses on Euler's equation of gas
dynamics in the subsonic regime using dG techniques. -->

However, in the transonic and supersonic regimes, and shock-hydrodynamics
applications the use of variational schemes might be questionable. In fact,
at the time of this writing, most shock-hydrodynamics codes are still
firmly grounded on finite volume methods. The main reason for failure of
Expand Down Expand Up @@ -244,7 +247,7 @@ vector-valued DoFHandler (i.e. initialized with an FESystem) it is possible
to compute the block-matrices (required in order to advance the solution)
with relative ease. However, the interactions that have to be computed in
the context of time-explicit collocation-type schemes (such as finite
differences and/or the scheme presented in this tutorial) can be are
differences and/or the scheme presented in this tutorial) can be
better described as <i>interactions between nodes</i> (not between DOFs).
In addition, in our case we do not solve a linear equation in order to
advance the solution. This leaves very little reason to use vector-valued
Expand Down Expand Up @@ -282,7 +285,7 @@ where
(\mathbf{U}_i^{n},\mathbf{U}_j^{n}, \textbf{n}_{ij}),
\lambda_{\text{max}} (\mathbf{U}_j^{n}, \mathbf{U}_i^{n},
\textbf{n}_{ji}) \} \|\mathbf{c}_{ij}\|$ if $i \not = j$ is the so
called <i>graph viscosity</i>. The graph viscosity serves as a
called <i>graph-viscosity</i>. The graph-viscosity serves as a
stabilization term, it is somewhat the discrete counterpart of
$\epsilon \Delta \mathbf{u}$ that appears in the notion of viscosity
solution described above. We will base our construction of $d_{ij}$ on
Expand Down Expand Up @@ -363,15 +366,15 @@ compute all $d_{ij}$ in a separate step prior to actually performing above
update. The core principle remains unchanged, though: we do not loop over
cells but rather over all edges of the sparsity graph.

@note It is not uncommon to encounter such fully algebraic schemes (i.e.
@note It is not uncommon to encounter such fully-algebraic schemes (i.e.
no bilinear forms, no cell loops, and no quadrature) outside of the finite
element community in the wider CFD community. There is a rich history of
application of this kind of schemes, also called <i>edge-based</i> or
<i>graph-based</i> finite element schemes (see for instance
@cite Rainald2008 for a historical overview). However, it is important to
highlight that the algebraic structure of the scheme (presented in this
tutorial) and the node-loops are not just a performance gimmick. Actually, the
structure of this scheme was borne out of theoretical necessity: the proof of
structure of this scheme was born out of theoretical necessity: the proof of
pointwise stability of the scheme hinges on the specific algebraic structure of
the scheme. In addition, it is not possible to compute the algebraic
viscosities $d_{ij}$ using cell-loops since they depend nonlinearly on
Expand Down Expand Up @@ -413,12 +416,17 @@ For the case of the reflecting boundary conditions we will proceed as follows:
@f{align*}
\mathbf{m}_i \dealcoloneq \mathbf{m}_i - (\widehat{\boldsymbol{\nu}}_i
\cdot \mathbf{m}_i) \widehat{\boldsymbol{\nu}}_i \ \
\text{for all }\mathbf{x}_i \in \partial\Omega^r
\text{where} \ \
\widehat{\boldsymbol{\nu}}_i \dealcoloneq
\frac{\int_{\partial\Omega} \phi_i \widehat{\boldsymbol{\nu}} \,
\, \mathrm{d}\mathbf{s}}{\big|\int_{\partial\Omega} \phi_i
\widehat{\boldsymbol{\nu}} \, \mathrm{d}\mathbf{s}\big|}
\ \ \text{for all }\mathbf{x}_i \in \partial\Omega^r
\ \ \ \ \boldsymbol{(1)}
@f}
that removes the normal component of $\mathbf{m}$. Here the definition of
nodal normal $\widehat{\boldsymbol{\nu}}_i$ is very much arbitrary (there is
no unique definition) but it should be consistent upon refinement with the
underlying geometry.
that removes the normal component of $\mathbf{m}$. This is a somewhat
naive idea that preserves a few fundamental properties of the PDE as we
explain below.

This is approach is usually called "explicit treatment of boundary conditions".
The well seasoned finite element person might find this approach questionable.
Expand Down Expand Up @@ -452,26 +460,37 @@ integrate in space and time $\int_{\Omega}\int_{t_1}^{t_2}$ we would obtain
@f{align*}
\int_{\Omega} \rho(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
\int_{\Omega} \rho(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ , \ \
\int_{\Omega} E(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
\int_{\Omega} E(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ , \ \
\int_{\Omega} \mathbf{m}(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
\int_{\Omega} \mathbf{m}(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x}
\int_{\Omega} \mathbf{m}(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x}
+ \int_{t_1}^{t_2} \! \int_{\partial\Omega} p \boldsymbol{\nu} \,
\mathrm{d}\mathbf{s} \mathrm{d}t\, .
\mathrm{d}\mathbf{s} \mathrm{d}t =
\int_{\Omega} \mathbf{m}(\mathbf{x},t_1) \,
\mathrm{d}\mathbf{x} \ , \ \
\int_{\Omega} E(\mathbf{x},t_2) \, \mathrm{d}\mathbf{x} =
\int_{\Omega} E(\mathbf{x},t_1) \, \mathrm{d}\mathbf{x} \ \ \ \
\boldsymbol{(2)}
@f}
Note that momentum is not a conserved quantity (interaction with walls leads to
momentum gain/loss). Even though we will not use reflecting boundary conditions
in the entirety of the domain, we would like to know that our implementation of
Note that momentum is NOT a conserved quantity (interaction with walls leads to
momentum gain/loss): however $\mathbf{m}$ has to satisfy a momentum balance.
Even though we will not use reflecting boundary conditions in the entirety of
the domain, we would like to know that our implementation of reflecting
boundary conditions is consistent with the conservation properties mentioned
above. In order to guarantee such conservation property it is necessary to
modify the values of the vectors $\mathbf{c}_{ij}$ as follows
above. In particular, if we use the projection $\boldsymbol{(1)}$ in the
entirety of the domain the following discrete mass-balance can be guaranteed:
@f{align*}
\mathbf{c}_{ij} \dealcoloneq \mathbf{c}_{ij} + \int_{\partial \Omega^r}
(\boldsymbol{\nu}_j - \boldsymbol{\nu}(\mathbf{s})) \phi_i \phi_j \,
\mathrm{d}\mathbf{s}
\ \ \text{whenever both }\mathbf{x}_i\text{ and }\mathbf{x}_j \text{ lie on the
boundary}
\sum_{i \in \mathcal{V}} m_i \rho_i^{n+1} =
\sum_{i \in \mathcal{V}} m_i \rho_i^{n} \ , \ \
\sum_{i \in \mathcal{V}} m_i \mathbf{m}_i^{n+1}
+ \tau_n \int_{\partial\Omega} \Big(\sum_{i \in \mathcal{V}} p_i^{n} \phi_i\Big)
\widehat{\boldsymbol{\nu}} \mathrm{d}\mathbf{s} =
\sum_{i \in \mathcal{V}} m_i \mathbf{m}_i^{n} \ , \ \
\sum_{i \in \mathcal{V}} m_i E_i^{n+1} = \sum_{i \in \mathcal{V}} m_i
E_i^{n} \ \ \ \
\boldsymbol{(3)}
@f}
This consistent modification of the $\mathbf{c}_{ij}$ is a direct consequence
of simple integration by parts arguments, see page 12 of @cite GuermondEtAl2018
for more details.
where $p_i$ is the pressure at the nodes that lie at the boundary. Clearly
$\boldsymbol{(3)}$ is the discrete counterpart of $\boldsymbol{(2)}$. The
proof of identity $\boldsymbol{(3)}$ is omitted, but we briefly mention that
it hinges on the definition of the <i>nodal normal</i>
$\widehat{\boldsymbol{\nu}}_i$ provided in $\boldsymbol{(1)}$. We also note that
this enforcement of reflecting boundary conditions is different from the one
originally advanced in @cite GuermondEtAl2018.
134 changes: 6 additions & 128 deletions examples/step-69/step-69.cc
View file
Open in desktop
Expand Up @@ -992,20 +992,14 @@ namespace Step69
//
// @f{align*}
// \widehat{\boldsymbol{\nu}}_i \dealcoloneq
// \frac{\boldsymbol{\nu}_i}{|\boldsymbol{\nu}_i|} \ \text{ where }
// \boldsymbol{\nu}_i \dealcoloneq \sum_{T \subset \text{supp}(\phi_i)}
// \sum_{F \subset \partial T \cap \partial \Omega}
// \sum_{\mathbf{x}_{q,F}} \nu(\mathbf{x}_{q,F})
// \phi_i(\mathbf{x}_{q,F}).
// \frac{\int_{\partial\Omega} \phi_i \widehat{\boldsymbol{\nu}} \,
// \, \mathrm{d}\mathbf{s}}{\big|\int_{\partial\Omega} \phi_i
// \widehat{\boldsymbol{\nu}} \, \mathrm{d}\mathbf{s}\big|}
// @f}
//
// Here $T$ denotes elements,
// $\text{supp}(\phi_i)$ the support of the shape function $\phi_i$,
// $F$ are faces of the element $T$, and $\mathbf{x}_{q,F}$
// are quadrature points on such face. Note that this formula for
// $\widehat{\boldsymbol{\nu}}_i$ is nothing else than some form of
// weighted averaging. Other more sophisticated definitions for $\nu_i$
// are possible but none of them have much influence in theory or practice.
// We will compute the numerator of this expression first and store it in
// <code>OfflineData<dim>::BoundaryNormalMap</code>. We will normalize these
// vectors in a posterior loop.

template <int dim>
void OfflineData<dim>::assemble()
Expand Down Expand Up @@ -1295,122 +1289,6 @@ namespace Step69
normal /= (normal.norm() + std::numeric_limits<double>::epsilon());
}
}

// As commented in the introduction (see section titled "Stable boundary
// conditions and conservation properties") we use three different types of
// boundary conditions: essential-like boundary conditions (we prescribe a
// state in the left portion of our domain), outflow boundary conditions
// (also called "do-nothing" boundary conditions) in the right portion of
// the domain, and "reflecting" boundary conditions (also called "free
// slip" boundary conditions). With these boundary conditions we should
// not expect any form of conservation to hold.
//
// However, if we were to use reflecting boundary conditions
// $\mathbf{m} \cdot \boldsymbol{\nu}_i =0$ on the entirety of the
// boundary we should preserve the density and total (mechanical)
// energy. This requires us to modify the $\mathbf{c}_{ij}$ vectors at
// the boundary as follows @cite GuermondEtAl2018 :
//
// @f{align*}
// \mathbf{c}_{ij} \, +\!\!= \int_{\partial \Omega}
// (\boldsymbol{\nu}_j - \boldsymbol{\nu}(\mathbf{s})) \phi_i \phi_j \,
// \mathrm{d}\mathbf{s} \ \text{whenever} \ \mathbf{x}_i \text{ and }
// \mathbf{x}_j \text{ lie in the boundary.}
// @f}
//
// The ideas repeat themselves: we use WorkStream in order to compute
// this correction, most of the following code is about the definition
// of the worker <code>local_assemble_system()</code>.
{
TimerOutput::Scope scope(computing_timer,
"offline_data - fix slip boundary c_ij");

const auto local_assemble_system = //
[&](const typename DoFHandler<dim>::cell_iterator &cell,
MeshWorker::ScratchData<dim> & scratch,
CopyData<dim> & copy) {
copy.is_artificial = cell->is_artificial();

if (copy.is_artificial)
return;

for (auto &matrix : copy.cell_cij_matrix)
matrix.reinit(dofs_per_cell, dofs_per_cell);

copy.local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(copy.local_dof_indices);
std::transform(copy.local_dof_indices.begin(),
copy.local_dof_indices.end(),
copy.local_dof_indices.begin(),
[&](types::global_dof_index index) {
return partitioner->global_to_local(index);
});

for (auto &matrix : copy.cell_cij_matrix)
matrix = 0.;

for (const auto f : cell->face_indices())
{
const auto face = cell->face(f);
const auto id = face->boundary_id();

if (!face->at_boundary())
continue;

if (id != Boundaries::free_slip)
continue;

const auto &fe_face_values = scratch.reinit(cell, f);

const unsigned int n_face_q_points =
fe_face_values.get_quadrature().size();

for (unsigned int q = 0; q < n_face_q_points; ++q)
{
const auto JxW = fe_face_values.JxW(q);
const auto normal_q = fe_face_values.normal_vector(q);

for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
if (!discretization->finite_element.has_support_on_face(
j, f))
continue;

const auto &normal_j = std::get<0>(
boundary_normal_map[copy.local_dof_indices[j]]);

const auto value_JxW =
fe_face_values.shape_value(j, q) * JxW;

for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const auto value = fe_face_values.shape_value(i, q);

/* This is the correction of the boundary c_ij */
for (unsigned int d = 0; d < dim; ++d)
copy.cell_cij_matrix[d](i, j) +=
(normal_j[d] - normal_q[d]) * (value * value_JxW);
} /* i */
} /* j */
} /* q */
} /* f */
};

const auto copy_local_to_global = [&](const CopyData<dim> &copy) {
if (copy.is_artificial)
return;

for (int k = 0; k < dim; ++k)
cij_matrix[k].add(copy.local_dof_indices, copy.cell_cij_matrix[k]);
};

WorkStream::run(dof_handler.begin_active(),
dof_handler.end(),
local_assemble_system,
copy_local_to_global,
scratch_data,
CopyData<dim>());
}
}

// At this point we are very much done with anything related to offline data.
Expand Down