Hi Wolfgang, Bruno and all community members,
My apologies for a delayed reply as I was pondering on a good way to
replicate the problem without bogging you down with the unnecessary
dependencies in my program. So what I come up with is to adapt the step-22
by adding a particle handler object and mimic what is done in step-19. I
was able to reproduce the error for refinement. I am working on a way to
reproduce the error for coarsening. The overall idea is:
- spawn one particle in one cell
- refine that cell
- An exception will be thrown in sort_particles_into_domains_and_cells():
--------------------------------------------------------
An error occurred in line <495> of file
</home/wuqihang/.local/src/dealii-9.2.0/include/deal.II/grid/tria_iterator.templates.h>
in function
dealii::TriaActiveIterator<Accessor>&
dealii::TriaActiveIterator<Accessor>::operator=(const
dealii::TriaIterator<Accessor>&) [with Accessor = dealii::CellAccessor<2,
2>]
The violated condition was:
this->accessor.has_children() == false
Additional information:
(none)
Stacktrace:
-----------
#0 /home/wuqihang/.local/lib/libdeal_II.g.so.9.2.0:
dealii::TriaActiveIterator<dealii::CellAccessor<2, 2>
>::operator=(dealii::TriaIterator<dealii::CellAccessor<2, 2> > const&)
#1 /home/wuqihang/.local/lib/libdeal_II.g.so.9.2.0:
dealii::Particles::ParticleHandler<2, 2>::update_cached_numbers()
#2 /home/wuqihang/.local/lib/libdeal_II.g.so.9.2.0:
dealii::Particles::ParticleHandler<2,
2>::sort_particles_into_subdomains_and_cells()
#3
/home/wuqihang/Dropbox/Coding/DEAL.II/Step22ParticlesAMR/build/step-22_PtclAMR:
Step22::StokesParticles<2, double>::move_particles()
#4
/home/wuqihang/Dropbox/Coding/DEAL.II/Step22ParticlesAMR/build/step-22_PtclAMR:
Step22::StokesProblem<2>::second_step()
#5
/home/wuqihang/Dropbox/Coding/DEAL.II/Step22ParticlesAMR/build/step-22_PtclAMR:
Step22::StokesProblem<2>::run()
#6
/home/wuqihang/Dropbox/Coding/DEAL.II/Step22ParticlesAMR/build/step-22_PtclAMR:
main
--------------------------------------------------------
My limited understanding is that after the grid refinement the particle is
not correctly associated with the finer child cell. I have not setup the
callback function in dealii::Triangulation for this test because the
dynamic_cast in the callback functions will throw an error anyways. What I
will try now is to follow Bruno's advice to remove dynamic_cast in the
callback functions and see what happens. I will report back if I shall make
any progress.
I very much appreciate your time bearing with me and I look forward to
helping out to fix the problem.
Thank you,
Qihang.
On Wednesday, December 16, 2020 at 7:37:01 p.m. UTC-5 [email protected]
wrote:
> Dear Wu,
> Have you tried without the dynamic cast? Can you bind the callback of a
> regular triangulation?
> I must admit I have never tried to use particles with regular
> triangulation (and I think that outside of step-19, not much people have
> done it ).
>
>
> On Tuesday, December 15, 2020 at 4:49:39 p.m. UTC-5 WU, Qihang wrote:
>
>> Hi deal.II Community,
>>
>> I am new to the library and I'm trying to implement a prototype Stokes
>> flow solver with tracer particles using shared memory computing. The
>> overall procedure is something like:
>>
>> 1) Solve velocities and pressure analogous to Step-22.
>> 2) Execute coarsening and refinement
>> 3) Interpolate velocities to particles and move particles.
>>
>> What I found is after grid coarsening, the particles which originally
>> belong to the finer grid loses its host cell and an exception is thrown
>> from ParticleAccessor::get_surrounding_cell().
>>
>> I read through step-19 on Github. From my understanding the
>> implementation there does not involve AMR and therefore such problems do
>> not occur. I also consulted step-68 and 70. Those two implementations use
>> parallel::distributed::Triangulation and attach pre- and post-refinement
>> callback functions to save and reload particles. However in these functions
>> the dynamic_cast from Triangulation to p::d::Triangulation is invalid and
>> returns a nullptr.
>>
>> My question is am I correct to assume that the library does not yet fully
>> support particles on dealii::Triangulation with AMR? If yes can I ask for
>> some guidance on how to implement it myself? One naive way out seems to be
>> to store and reload procedure as is done in the library for
>> p::d::Triangulation, but is it too costly?
>>
>> I really appreciate your time and help.
>>
>> Cheers,
>> Qihang.
>>
>>
>>
>>
>>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/15564e56-acb0-4a33-b658-2ce295638667n%40googlegroups.com.
##
# CMake script for the step-22 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "step-22_PtclAMR")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
FIND_PACKAGE(deal.II 9.2.0 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
#
# Are all dependencies fulfilled?
#
IF(NOT DEAL_II_WITH_UMFPACK) # keep in one line
MESSAGE(FATAL_ERROR "
Error! This tutorial requires a deal.II library that was configured with the
following options:
DEAL_II_WITH_UMFPACK = ON
However, the deal.II library found at ${DEAL_II_PATH} was configured with these
options
DEAL_II_WITH_UMFPACK = ${DEAL_II_WITH_UMFPACK}
which conflict with the requirements."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_ilu.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <fstream>
#include <iostream>
#include <memory>
#include "stokes_particles.h"
namespace Step22
{
using namespace dealii;
template <int dim>
struct InnerPreconditioner;
template <>
struct InnerPreconditioner<2>
{
using type = SparseDirectUMFPACK;
};
template <>
struct InnerPreconditioner<3>
{
using type = SparseILU<double>;
};
template <int dim>
class StokesProblem
{
friend class StokesParticles<dim>;
public:
StokesProblem(const unsigned int degree);
void run();
private:
void setup_dofs();
void assemble_system();
void solve();
void output_results(const unsigned int refinement_cycle) const;
using ActiveCellIterator = typename DoFHandler<dim>::active_cell_iterator;
void first_step();
void second_step();
const unsigned int degree;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
MappingQGeneric<dim> mapping;
AffineConstraints<double> constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockSparsityPattern preconditioner_sparsity_pattern;
BlockSparseMatrix<double> preconditioner_matrix;
BlockVector<double> solution;
BlockVector<double> system_rhs;
std::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
std::unique_ptr<StokesParticles<dim>> ptr_particles;
constexpr static double floating_point_tolerance = 1.0e-10;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues()
: Function<dim>(dim + 1)
{}
virtual double value(const Point<dim>& p,
const unsigned int component = 0) const override;
virtual void vector_value(const Point<dim>& p,
Vector<double>& value) const override;
};
template <int dim>
double BoundaryValues<dim>::value(const Point<dim>& p,
const unsigned int component) const
{
Assert(component < this->n_components,
ExcIndexRange(component, 0, this->n_components));
if (component == 0)
return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0));
return 0;
}
template <int dim>
void BoundaryValues<dim>::vector_value(const Point<dim>& p,
Vector<double>& values) const
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = BoundaryValues<dim>::value(p, c);
}
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide()
: Function<dim>(dim + 1)
{}
virtual double value(const Point<dim>& p,
const unsigned int component = 0) const override;
virtual void vector_value(const Point<dim>& p,
Vector<double>& value) const override;
};
template <int dim>
double RightHandSide<dim>::value(const Point<dim>& /*p*/,
const unsigned int /*component*/) const
{
return 0;
}
template <int dim>
void RightHandSide<dim>::vector_value(const Point<dim>& p,
Vector<double>& values) const
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = RightHandSide<dim>::value(p, c);
}
template <class MatrixType, class PreconditionerType>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix(const MatrixType& m,
const PreconditionerType& preconditioner);
void vmult(Vector<double>& dst, const Vector<double>& src) const;
private:
const SmartPointer<const MatrixType> matrix;
const SmartPointer<const PreconditionerType> preconditioner;
};
template <class MatrixType, class PreconditionerType>
InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
const MatrixType& m,
const PreconditionerType& preconditioner)
: matrix(&m)
, preconditioner(&preconditioner)
{}
template <class MatrixType, class PreconditionerType>
void InverseMatrix<MatrixType, PreconditionerType>::vmult(
Vector<double>& dst,
const Vector<double>& src) const
{
SolverControl solver_control(src.size(), 1e-6 * src.l2_norm());
SolverCG<Vector<double>> cg(solver_control);
dst = 0;
cg.solve(*matrix, dst, src, *preconditioner);
}
template <class PreconditionerType>
class SchurComplement : public Subscriptor
{
public:
SchurComplement(
const BlockSparseMatrix<double>& system_matrix,
const InverseMatrix<SparseMatrix<double>, PreconditionerType>& A_inverse);
void vmult(Vector<double>& dst, const Vector<double>& src) const;
private:
const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
const SmartPointer<
const InverseMatrix<SparseMatrix<double>, PreconditionerType>>
A_inverse;
mutable Vector<double> tmp1, tmp2;
};
template <class PreconditionerType>
SchurComplement<PreconditionerType>::SchurComplement(
const BlockSparseMatrix<double>& system_matrix,
const InverseMatrix<SparseMatrix<double>, PreconditionerType>& A_inverse)
: system_matrix(&system_matrix)
, A_inverse(&A_inverse)
, tmp1(system_matrix.block(0, 0).m())
, tmp2(system_matrix.block(0, 0).m())
{}
template <class PreconditionerType>
void
SchurComplement<PreconditionerType>::vmult(Vector<double>& dst,
const Vector<double>& src) const
{
system_matrix->block(0, 1).vmult(tmp1, src);
A_inverse->vmult(tmp2, tmp1);
system_matrix->block(1, 0).vmult(dst, tmp2);
}
template <int dim>
StokesProblem<dim>::StokesProblem(const unsigned int degree)
: degree(degree)
, triangulation()
, fe(FE_Q<dim>(degree + 1), dim, FE_Q<dim>(degree), 1)
, dof_handler(triangulation)
, mapping(degree)
, ptr_particles(std::make_unique<StokesParticles<dim>>(*this))
{}
template <int dim>
void StokesProblem<dim>::setup_dofs()
{
A_preconditioner.reset();
system_matrix.clear();
preconditioner_matrix.clear();
dof_handler.distribute_dofs(fe);
DoFRenumbering::Cuthill_McKee(dof_handler);
std::vector<unsigned int> block_component(dim + 1, 0);
block_component[dim] = 1;
DoFRenumbering::component_wise(dof_handler, block_component);
{
constraints.clear();
FEValuesExtractors::Vector velocities(0);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
1,
BoundaryValues<dim>(),
constraints,
fe.component_mask(velocities));
}
constraints.close();
const std::vector<types::global_dof_index> dofs_per_block =
DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
const unsigned int n_u = dofs_per_block[0];
const unsigned int n_p = dofs_per_block[1];
std::cout << " Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< " (" << n_u << '+' << n_p << ')' << std::endl;
{
BlockDynamicSparsityPattern dsp(2, 2);
dsp.block(0, 0).reinit(n_u, n_u);
dsp.block(1, 0).reinit(n_p, n_u);
dsp.block(0, 1).reinit(n_u, n_p);
dsp.block(1, 1).reinit(n_p, n_p);
dsp.collect_sizes();
Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
for (unsigned int c = 0; c < dim + 1; ++c)
for (unsigned int d = 0; d < dim + 1; ++d)
if (!((c == dim) && (d == dim)))
coupling[c][d] = DoFTools::always;
else
coupling[c][d] = DoFTools::none;
DoFTools::make_sparsity_pattern(
dof_handler, coupling, dsp, constraints, false);
sparsity_pattern.copy_from(dsp);
}
{
BlockDynamicSparsityPattern preconditioner_dsp(2, 2);
preconditioner_dsp.block(0, 0).reinit(n_u, n_u);
preconditioner_dsp.block(1, 0).reinit(n_p, n_u);
preconditioner_dsp.block(0, 1).reinit(n_u, n_p);
preconditioner_dsp.block(1, 1).reinit(n_p, n_p);
preconditioner_dsp.collect_sizes();
Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1);
for (unsigned int c = 0; c < dim + 1; ++c)
for (unsigned int d = 0; d < dim + 1; ++d)
if (((c == dim) && (d == dim)))
preconditioner_coupling[c][d] = DoFTools::always;
else
preconditioner_coupling[c][d] = DoFTools::none;
DoFTools::make_sparsity_pattern(dof_handler,
preconditioner_coupling,
preconditioner_dsp,
constraints,
false);
preconditioner_sparsity_pattern.copy_from(preconditioner_dsp);
}
system_matrix.reinit(sparsity_pattern);
preconditioner_matrix.reinit(preconditioner_sparsity_pattern);
solution.reinit(2);
solution.block(0).reinit(n_u);
solution.block(1).reinit(n_p);
solution.collect_sizes();
system_rhs.reinit(2);
system_rhs.block(0).reinit(n_u);
system_rhs.block(1).reinit(n_p);
system_rhs.collect_sizes();
}
template <int dim>
void StokesProblem<dim>::assemble_system()
{
system_matrix = 0;
system_rhs = 0;
preconditioner_matrix = 0;
QGauss<dim> quadrature_formula(degree + 2);
FEValues<dim> fe_values(mapping,
fe,
quadrature_formula,
update_values | update_quadrature_points |
update_JxW_values | update_gradients);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> local_preconditioner_matrix(dofs_per_cell,
dofs_per_cell);
Vector<double> local_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const RightHandSide<dim> right_hand_side;
std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim + 1));
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(dim);
std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
std::vector<double> div_phi_u(dofs_per_cell);
std::vector<double> phi_p(dofs_per_cell);
for (const auto& cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
local_matrix = 0;
local_preconditioner_matrix = 0;
local_rhs = 0;
right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
rhs_values);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
symgrad_phi_u[k] =
fe_values[velocities].symmetric_gradient(k, q);
div_phi_u[k] = fe_values[velocities].divergence(k, q);
phi_p[k] = fe_values[pressure].value(k, q);
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j <= i; ++j)
{
local_matrix(i, j) +=
(2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) // (1)
- div_phi_u[i] * phi_p[j] // (2)
- phi_p[i] * div_phi_u[j]) // (3)
* fe_values.JxW(q); // * dx
local_preconditioner_matrix(i, j) +=
(phi_p[i] * phi_p[j]) // (4)
* fe_values.JxW(q); // * dx
}
const unsigned int component_i =
fe.system_to_component_index(i).first;
local_rhs(i) += (fe_values.shape_value(i, q) // (phi_u_i(x_q)
* rhs_values[q](component_i)) // * f(x_q))
* fe_values.JxW(q); // * dx
}
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = i + 1; j < dofs_per_cell; ++j)
{
local_matrix(i, j) = local_matrix(j, i);
local_preconditioner_matrix(i, j) =
local_preconditioner_matrix(j, i);
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(local_matrix,
local_rhs,
local_dof_indices,
system_matrix,
system_rhs);
constraints.distribute_local_to_global(local_preconditioner_matrix,
local_dof_indices,
preconditioner_matrix);
}
std::cout << " Computing preconditioner..." << std::endl << std::flush;
A_preconditioner =
std::make_shared<typename InnerPreconditioner<dim>::type>();
A_preconditioner->initialize(
system_matrix.block(0, 0),
typename InnerPreconditioner<dim>::type::AdditionalData());
}
template <int dim>
void StokesProblem<dim>::solve()
{
const InverseMatrix<SparseMatrix<double>,
typename InnerPreconditioner<dim>::type>
A_inverse(system_matrix.block(0, 0), *A_preconditioner);
Vector<double> tmp(solution.block(0).size());
{
Vector<double> schur_rhs(solution.block(1).size());
A_inverse.vmult(tmp, system_rhs.block(0));
system_matrix.block(1, 0).vmult(schur_rhs, tmp);
schur_rhs -= system_rhs.block(1);
SchurComplement<typename InnerPreconditioner<dim>::type> schur_complement(
system_matrix, A_inverse);
SolverControl solver_control(solution.block(1).size(),
1e-6 * schur_rhs.l2_norm());
SolverCG<Vector<double>> cg(solver_control);
SparseILU<double> preconditioner;
preconditioner.initialize(preconditioner_matrix.block(1, 1),
SparseILU<double>::AdditionalData());
InverseMatrix<SparseMatrix<double>, SparseILU<double>> m_inverse(
preconditioner_matrix.block(1, 1), preconditioner);
cg.solve(schur_complement, solution.block(1), schur_rhs, m_inverse);
constraints.distribute(solution);
std::cout << " " << solver_control.last_step()
<< " outer CG Schur complement iterations for pressure"
<< std::endl;
}
{
system_matrix.block(0, 1).vmult(tmp, solution.block(1));
tmp *= -1;
tmp += system_rhs.block(0);
A_inverse.vmult(solution.block(0), tmp);
constraints.distribute(solution);
}
}
template <int dim>
void
StokesProblem<dim>::output_results(const unsigned int refinement_cycle) const
{
std::vector<std::string> solution_names(dim, "velocity");
solution_names.emplace_back("pressure");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(
dim, DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation.push_back(
DataComponentInterpretation::component_is_scalar);
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution,
solution_names,
DataOut<dim>::type_dof_data,
data_component_interpretation);
data_out.build_patches();
std::ofstream output(
"solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtk");
data_out.write_vtk(output);
// add output for particles
std::ofstream ptcl_output(
"particles-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtk");
ptr_particles->output_vtk(ptcl_output);
}
template <int dim>
void StokesProblem<dim>::first_step()
{
// define a lambda function to
// find one cell that has a vertex overlapping with origin
auto find_origin_cell = [&]() {
for (auto cell = dof_handler.begin_active(); cell != dof_handler.end();
++cell)
for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
if (cell->vertex(i).norm() < floating_point_tolerance)
return cell;
return dof_handler.begin_active();
};
setup_dofs();
// We will spawn one and only one particle
// in an origin-overlapping cell
ActiveCellIterator cell = find_origin_cell();
ptr_particles->create_particles(cell);
// Let Stokes solver do its work
assemble_system();
solve();
output_results(0);
// advect particles
ptr_particles->move_particles();
// only refine the cell that contains the particle
cell->set_refine_flag();
triangulation.execute_coarsening_and_refinement();
}
template <int dim>
void StokesProblem<dim>::second_step()
{
setup_dofs();
assemble_system();
solve();
output_results(1);
ptr_particles->move_particles();
}
template <int dim>
void StokesProblem<dim>::run()
{
{
std::vector<unsigned int> subdivisions(dim, 1);
subdivisions[0] = 4;
const Point<dim> bottom_left = (dim == 2 ? //
Point<dim>(-2, -1) : // 2d case
Point<dim>(-2, 0, -1)); // 3d case
const Point<dim> top_right = (dim == 2 ? //
Point<dim>(2, 0) : // 2d case
Point<dim>(2, 1, 0)); // 3d case
GridGenerator::subdivided_hyper_rectangle(triangulation,
subdivisions,
bottom_left,
top_right);
}
for (const auto& cell : triangulation.active_cell_iterators())
for (const auto& face : cell->face_iterators())
if (face->center()[dim - 1] == 0)
face->set_all_boundary_ids(1);
triangulation.refine_global(4 - dim);
first_step();
second_step();
}
} // namespace Step22
int main()
{
try
{
using namespace Step22;
StokesProblem<2> flow_problem(1);
flow_problem.run();
}
catch (std::exception& exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
#ifndef _STEP_22_PTCLAMR_STOKES_PARTICLES_H_
#define _STEP_22_PTCLAMR_STOKES_PARTICLES_H_
#include <deal.II/base/array_view.h>
#include <deal.II/base/types.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/particles/data_out.h>
#include <deal.II/particles/particle_accessor.h>
#include <deal.II/particles/particle_handler.h>
namespace Step22
{
using namespace dealii;
/* Forward declarations */
template <int dim>
class StokesProblem;
/**
* Tracer particles for Stokes solver.
*/
template <int dim, typename NumberType = double>
class StokesParticles
{
public:
using value_type = NumberType;
/**
* @brief Construct Particles and set CFL.
* We choose a small CFL so that the particles will not
* leave the cell after one step of advection.
*/
StokesParticles(const StokesProblem<dim>& stokes_solver,
value_type CFL = 0.1);
/**
* @brief Create particles in the given cell.
* Spawn 1 particle at the center at the mid-point integration point
*/
void create_particles(
const typename DoFHandler<dim>::active_cell_iterator& cell);
/**
* @brief Create particles in each cell.
* Will call the function above at every active cell.
*/
void create_particles();
/**
* @brief Advect particles.
* The particle velocities will be interpolated from
* the solution of \c StokesProblem flow solver.
* @return the time step that is advanced by this advection step.
*/
value_type move_particles();
/**
* @brief Write the particles to file.
*/
void output_vtk(std::ostream& os) const;
private:
const StokesProblem<dim>* ptr_stokes_solver;
Particles::ParticleHandler<dim> particle_handler;
value_type cfl_number;
types::particle_index next_unused_id;
};
/****************************************************/
/* IMPLEMENTATIONS */
/****************************************************/
template <int dim, typename NumberType>
StokesParticles<dim, NumberType>::StokesParticles(
const StokesProblem<dim>& stokes_solver,
value_type CFL)
: ptr_stokes_solver(&stokes_solver)
, particle_handler(stokes_solver.triangulation, stokes_solver.mapping, dim)
, cfl_number(CFL)
, next_unused_id(0)
{
Assert(CFL > 0.0, ExcMessage("CFL number must be positive."));
}
template <int dim, typename NumberType>
void StokesParticles<dim, NumberType>::create_particles(
const typename DoFHandler<dim>::active_cell_iterator& cell)
{
using namespace dealii::Particles;
FEValues<dim> feval(ptr_stokes_solver->mapping,
ptr_stokes_solver->fe,
QMidpoint<dim>(),
update_quadrature_points);
feval.reinit(cell);
// spawn new particles at the quadrature point(s)
for (const auto iqpt : feval.quadrature_point_indices())
{
Particle<dim> new_ptcl;
new_ptcl.set_location(feval.quadrature_point(iqpt));
new_ptcl.set_id(next_unused_id++);
particle_handler.insert_particle(new_ptcl, cell);
}
particle_handler.update_cached_numbers();
}
template <int dim, typename NumberType>
void StokesParticles<dim, NumberType>::create_particles()
{
using namespace dealii::Particles;
for (const auto& cell :
ptr_stokes_solver->dof_handler.active_cell_iterators())
create_particles(cell);
}
template <int dim, typename NumberType>
auto StokesParticles<dim, NumberType>::move_particles() -> value_type
{
using namespace dealii::Particles;
FEValuesExtractors::Vector v(0);
value_type max_velocity_over_diameter = 0.0;
// **************************************************
// If there is a refinement in the previous time step,
// An expcetion will be thrown in the next line.
// **************************************************
particle_handler.sort_particles_into_subdomains_and_cells();
// run CFL estimation
for (const auto& cell :
ptr_stokes_solver->dof_handler.active_cell_iterators())
{
if (particle_handler.n_particles_in_cell(cell) > 0)
{
auto particles_range = particle_handler.particles_in_cell(cell);
std::vector<Point<dim>> particle_positions;
std::transform(particles_range.begin(),
particles_range.end(),
std::back_inserter(particle_positions),
[](const ParticleAccessor<dim>& p) {
return p.get_location();
});
// Use the particle location to construct a quadrature rule
Quadrature<dim> quadrature_on_particles(particle_positions);
// Now use the FEValus object to
// interpolate velocity solution (from Stokes solver)
// to particle(s)
FEValues<dim> fevals(ptr_stokes_solver->mapping,
ptr_stokes_solver->fe,
quadrature_on_particles,
update_values);
fevals.reinit(cell);
std::vector<dealii::Tensor<1, dim, value_type>> particle_velocities(
particle_handler.n_particles_in_cell(cell));
fevals[v].get_function_values(ptr_stokes_solver->solution,
particle_velocities);
// copy velocities into PropertyPool
size_t iptcl = 0;
for (auto& ptcl : particles_range)
{
std::vector<value_type> ptcl_velo_vector(dim);
std::copy(particle_velocities[iptcl].begin_raw(),
particle_velocities[iptcl].end_raw(),
ptcl_velo_vector.begin());
ptcl.set_properties(ptcl_velo_vector);
Assert(cell->diameter() > 0.0, ExcInternalError());
max_velocity_over_diameter =
std::max(particle_velocities[iptcl].norm() / cell->diameter(),
max_velocity_over_diameter);
++iptcl;
}
}
}
// Estimate the size of next CFL-complying time step
const value_type dt = cfl_number / max_velocity_over_diameter;
// Execute the particle advection
std::for_each(particle_handler.begin(),
particle_handler.end(),
[&](ParticleAccessor<dim>& p) {
dealii::Point<dim, value_type> new_position =
p.get_location();
for (int idim = 0; idim < dim; ++idim)
new_position[idim] += dt * p.get_properties()[idim];
p.set_location(new_position);
});
particle_handler.sort_particles_into_subdomains_and_cells();
return dt;
}
template <int dim, typename NumberType>
void StokesParticles<dim, NumberType>::output_vtk(std::ostream& os) const
{
Particles::DataOut<dim, dim> ptcl_out;
ptcl_out.build_patches(particle_handler);
ptcl_out.write_vtk(os);
}
} // namespace Step22
#endif // _STEP_22_PTCLAMR_STOKES_PARTICLES_H_ //
