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

Changes from 1 commit
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
116 changes: 0 additions & 116 deletions examples/step-69/step-69.cc
View file
Open in desktop
Expand Up @@ -1295,122 +1295,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