From: Luca Heltai Date: Thu, 21 Nov 2019 19:28:22 +0000 (+0100) Subject: Step-70 X-Git-Tag: v9.2.0-rc2~3^2~35 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ba54847b47ab35c3d2084848a3c49b3b51de60ae;p=dealii.git Step-70 Work in progress for step-70. Co-authored-by: Bruno Blais --- diff --git a/examples/step-70/CMakeLists.txt b/examples/step-70/CMakeLists.txt new file mode 100644 index 0000000000..760967fdff --- /dev/null +++ b/examples/step-70/CMakeLists.txt @@ -0,0 +1,62 @@ +## +# CMake script +## + +# Set the name of the project and target: +SET(TARGET "step-70") + +# 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_PETSC AND NOT DEAL_II_PETSC_WITH_COMPLEX) OR DEAL_II_WITH_TRILINOS) OR NOT DEAL_II_WITH_P4EST) # 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_PETSC = ON + DEAL_II_PETSC_WITH_COMPLEX = OFF + DEAL_II_WITH_P4EST = ON +or + DEAL_II_WITH_TRILINOS = ON + DEAL_II_WITH_P4EST = ON +However, the deal.II library found at ${DEAL_II_PATH} was configured with these options + DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC} + DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX} + DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST} + DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS} +which conflict with the requirements. +One or both of the aforementioned combinations of prerequisites are not met by your installation, but at least one is required for this tutorial step." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +SET(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.vtk *.ucd *.d2 *.vtu *.pvtu) +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-70/step-70.cc b/examples/step-70/step-70.cc new file mode 100644 index 0000000000..fa0d265fb4 --- /dev/null +++ b/examples/step-70/step-70.cc @@ -0,0 +1,1137 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2020 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * + * Authors: Luca Heltai, Bruno Blais, 2019 + */ + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#define FORCE_USE_OF_TRILINOS + +namespace LA +{ +#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \ + !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS)) + using namespace dealii::LinearAlgebraPETSc; +# define USE_PETSC_LA +#elif defined(DEAL_II_WITH_TRILINOS) + using namespace dealii::LinearAlgebraTrilinos; +#else +# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required +#endif +} // namespace LA + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +namespace Step70 +{ + using namespace dealii; + + template + void interpolate_field_on_particles( + const DoFHandler & field_dh, + const Particles::ParticleHandler &particle_handler, + const InputVectorType & field_vector, + OutputVectorType & interpolated_field, + const ComponentMask & field_comps) + { + if (particle_handler.n_locally_owned_particles() == 0) + { + interpolated_field.compress(VectorOperation::add); + return; // nothing else to do here + } + + const auto &tria = field_dh.get_triangulation(); + const auto &fe = field_dh.get_fe(); + auto particle = particle_handler.begin(); + + // Take care of components + const ComponentMask comps = + (field_comps.size() == 0 ? ComponentMask(fe.n_components(), true) : + field_comps); + AssertDimension(comps.size(), fe.n_components()); + const auto n_comps = comps.n_selected_components(); + + AssertDimension(field_vector.size(), field_dh.n_dofs()); + AssertDimension(interpolated_field.size(), + particle_handler.get_next_free_particle_index() * n_comps); + // Add check on locally owned indices + + // Global to local indices + std::vector space_gtl(fe.n_components(), + numbers::invalid_unsigned_int); + for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i) + if (comps[i]) + space_gtl[i] = j++; + + std::vector dof_indices(fe.dofs_per_cell); + + while (particle != particle_handler.end()) + { + const auto &cell = particle->get_surrounding_cell(tria); + const auto &dh_cell = + typename DoFHandler::cell_iterator(*cell, &field_dh); + dh_cell->get_dof_indices(dof_indices); + const auto pic = particle_handler.particles_in_cell(cell); + Assert(pic.begin() == particle, ExcInternalError()); + for (unsigned int i = 0; particle != pic.end(); ++particle, ++i) + { + const auto &reference_location = particle->get_reference_location(); + + const auto id = particle->get_id(); + + for (unsigned int j = 0; j < fe.dofs_per_cell; ++j) + { + const auto comp_j = + space_gtl[fe.system_to_component_index(j).first]; + if (comp_j != numbers::invalid_unsigned_int) + interpolated_field[id * n_comps + comp_j] += + fe.shape_value(j, reference_location) * + field_vector(dof_indices[j]); + } + } + } + interpolated_field.compress(VectorOperation::add); + } + + template + class StokesImmersedProblemParameters : public ParameterAcceptor + { + public: + StokesImmersedProblemParameters() + : ParameterAcceptor("Stokes Immersed Problem/") + , rhs("Right hand side", spacedim + 1) + , angular_velocity("Angular velocity", spacedim == 3 ? spacedim : 1) + { + add_parameter("Velocity degree", + velocity_degree, + "", + this->prm, + Patterns::Integer(1)); + + add_parameter("Number of time steps", number_of_time_steps); + + add_parameter("Final time", final_time); + + add_parameter("Viscosity", viscosity); + + add_parameter("Nitsche penalty term", penalty_term); + + add_parameter("Initial fluid refinement", + initial_fluid_refinement, + "Initial mesh refinement used for the fluid domain Omega"); + + add_parameter("Initial solid refinement", + initial_solid_refinement, + "Initial mesh refinement used for the solid domain Gamma"); + + add_parameter( + "Particle insertion refinement", + particle_insertion_refinement, + "Refinement of the volumetric mesh used to insert the particles"); + + add_parameter( + "Homogeneous Dirichlet boundary ids", + homogeneous_dirichlet_ids, + "Boundary Ids over which homogeneous Dirichlet boundary conditions are applied"); + + enter_my_subsection(this->prm); + this->prm.enter_subsection("Grid generation"); + this->prm.add_parameter("Grid one generator", name_of_grid1); + this->prm.add_parameter("Grid one generator arguments", + arguments_for_grid1); + + this->prm.add_parameter("Grid two generator", name_of_grid2); + this->prm.add_parameter("Grid two generator arguments", + arguments_for_grid2); + + this->prm.add_parameter("Particle grid generator", name_of_particle_grid); + this->prm.add_parameter("Particle grid generator arguments", + arguments_for_particle_grid); + this->prm.leave_subsection(); + + leave_my_subsection(this->prm); + + // correct the default dimension for the functions + rhs.declare_parameters_call_back.connect([&]() { + Functions::ParsedFunction::declare_parameters(this->prm, + spacedim + 1); + }); + angular_velocity.declare_parameters_call_back.connect([&]() { + Functions::ParsedFunction::declare_parameters( + this->prm, spacedim == 3 ? spacedim : 1); + }); + } + + void set_time(const double &time) const + { + rhs.set_time(time); + angular_velocity.set_time(time); + } + + unsigned int velocity_degree = 2; + unsigned int number_of_time_steps = 1; + double viscosity = 1.0; + double final_time = 1.0; + unsigned int initial_fluid_refinement = 3; + unsigned int initial_solid_refinement = 3; + unsigned int particle_insertion_refinement = 1; + double penalty_term = 1e3; + std::list homogeneous_dirichlet_ids{0, 1, 2, 3}; + std::string name_of_grid1 = "hyper_cube"; + std::string arguments_for_grid1 = "-1: 1: false"; + std::string name_of_grid2 = "hyper_rectangle"; + std::string arguments_for_grid2 = + dim == 2 ? "-.5, -.1: .5, .1: false" : "-.5, -.1, -.1: .5, .1, .1: false"; + std::string name_of_particle_grid = "hyper_ball"; + std::string arguments_for_particle_grid = + dim == 2 ? "0.3, 0.3: 0.1: false" : "0.3, 0.3, 0.3 : 0.1: false"; + + mutable ParameterAcceptorProxy> rhs; + mutable ParameterAcceptorProxy> + angular_velocity; + }; // namespace Step70 + + + template + class SolidVelocity : public Function + { + public: + SolidVelocity(const Functions::ParsedFunction &angular_velocity) + : angular_velocity(angular_velocity) + {} + + virtual double value(const Point &p, + unsigned int component = 0) const + { + Tensor<1, spacedim> velocity; + if (spacedim == 3) + { + Tensor<1, spacedim> omega; + for (unsigned int i = 0; i < spacedim; ++i) + omega[i] = angular_velocity.value(p, i); + + velocity = cross_product_3d(p, omega); + } + + if (spacedim == 2) + { + double omega = angular_velocity.value(p, 0); + + velocity[0] = -omega * p[1]; + velocity[1] = omega * p[0]; + } + + return velocity[component]; + } + + private: + const Functions::ParsedFunction &angular_velocity; + }; + + + template + class SolidDisplacement : public Function + { + public: + SolidDisplacement( + const Functions::ParsedFunction &angular_velocity, + const double time_step) + : Function(spacedim) + , angular_velocity(angular_velocity) + , time_step(time_step) + {} + + virtual double value(const Point &p, + unsigned int component = 0) const + { + Tensor<1, spacedim> displacement; + + double dtheta = angular_velocity.value(p, 0) * time_step; + + displacement[0] = std::cos(dtheta) * p[0] - std::sin(dtheta) * p[1]; + displacement[1] = std::sin(dtheta) * p[0] + std::cos(dtheta) * p[1]; + + return displacement[component]; + } + + void set_time_step(const double new_time_step) + { + time_step = new_time_step; + } + + private: + const Functions::ParsedFunction &angular_velocity; + double time_step; + }; + + template + class StokesImmersedProblem + { + public: + StokesImmersedProblem( + const StokesImmersedProblemParameters &par); + + void run(); + + private: + void make_grid(); + void setup_tracer_particles(); + void setup_solid_particles(); + void setup_system(); + void assemble_stokes_system(); + void assemble_nitche_restriction(); + void solve(); + void refine_grid(); + void output_results(const unsigned int cycle) const; + + void + output_particles(const Particles::ParticleHandler &particles, + std::string fprefix, + const unsigned int iter) const; + + const StokesImmersedProblemParameters ∥ + + MPI_Comm mpi_communicator; + + std::unique_ptr> fe1; + std::unique_ptr> fe2; + + parallel::distributed::Triangulation tria1; + parallel::distributed::Triangulation tria2; + + DoFHandler dh1; + DoFHandler dh2; + + std::unique_ptr> mapping2; + + std::vector owned1; + std::vector owned2; + + std::vector relevant1; + std::vector relevant2; + + IndexSet owned_tracer_particles; + IndexSet relevant_tracer_particles; + + AffineConstraints constraints; + + LA::MPI::BlockSparseMatrix system_matrix; + LA::MPI::BlockSparseMatrix coupling_matrix; + + LA::MPI::BlockSparseMatrix preconditioner_matrix; + LA::MPI::BlockVector solution; + LA::MPI::BlockVector locally_relevant_solution; + LA::MPI::BlockVector system_rhs; + + LA::MPI::Vector tracer_particle_velocities; + LA::MPI::Vector relevant_tracer_particle_displacements; + + std::unique_ptr> quadrature_formula; + + // std::unique_ptr> + // dof_coupling; + + Particles::ParticleHandler tracer_particle_handler; + Particles::ParticleHandler solid_particle_handler; + + ConditionalOStream pcout; + mutable TimerOutput computing_timer; + }; + + + + template + StokesImmersedProblem::StokesImmersedProblem( + const StokesImmersedProblemParameters &par) + : par(par) + , mpi_communicator(MPI_COMM_WORLD) + , tria1(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , tria2(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , dh1(tria1) + , dh2(tria2) + , pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + , computing_timer(mpi_communicator, + pcout, + TimerOutput::summary, + TimerOutput::wall_times) + {} + + + template + void StokesImmersedProblem::make_grid() + { + GridGenerator::generate_from_name_and_arguments(tria1, + par.name_of_grid1, + par.arguments_for_grid1); + tria1.refine_global(par.initial_fluid_refinement); + + GridGenerator::generate_from_name_and_arguments(tria2, + par.name_of_grid2, + par.arguments_for_grid2); + tria2.refine_global(par.initial_solid_refinement); + } + + template + void StokesImmersedProblem::setup_tracer_particles() + { + // Generate a triangulation that will be used to decide the position + // of the particles to insert. In this case we choose an hyper_ball, a + // circle (spacedim==2) or a sphere (spacedim==3) filled with particles are + // the position of the support points of the triangulation + parallel::distributed::Triangulation particle_insert_tria( + mpi_communicator); + GridGenerator::generate_from_name_and_arguments( + particle_insert_tria, + par.name_of_particle_grid, + par.arguments_for_particle_grid); + particle_insert_tria.refine_global(par.particle_insertion_refinement); + + // Generate the support point on the triangulation that will be used as + // particle insertion point + DoFHandler particles_dof_handler(particle_insert_tria); + FE_Q particles_fe(1); + particles_dof_handler.distribute_dofs(particles_fe); + + // Create the particle handler associated with the fluid triangulation + tracer_particle_handler.initialize(tria1, + StaticMappingQ1::mapping); + + + // Generate the necessary local and global bounding boxes for the generator. + // The generation of the global bounding boxes requires an all-to-all + // communication + auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box( + tria1, IteratorFilters::LocallyOwnedCell()); + auto global_bounding_boxes = + Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box); + + + // Finally generate the particles from the support point of the + // particle_insert_tria triangulation + Particles::Generators::dof_support_points(particles_dof_handler, + global_bounding_boxes, + tracer_particle_handler); + + owned_tracer_particles = + tracer_particle_handler.locally_relevant_ids().tensor_product( + complete_index_set(spacedim)); + + relevant_tracer_particles = owned_tracer_particles; + } + + template + void StokesImmersedProblem::setup_solid_particles() + { + QGauss quadrature(fe1->degree + 1); + // In codimension one case, we store also the normal, else only the + // quadrature weight. + const unsigned int n_properties = (dim == spacedim) ? 1 : spacedim + 1; + + solid_particle_handler.initialize(tria1, + StaticMappingQ1::mapping, + n_properties); + + std::vector> quadrature_points_vec( + quadrature.size() * tria2.n_locally_owned_active_cells()); + + std::vector> properties( + quadrature.size() * tria2.n_locally_owned_active_cells(), + std::vector(n_properties)); + + UpdateFlags flags = update_JxW_values | update_quadrature_points; + if (spacedim > dim) + flags |= update_normal_vectors; + FEValues fe_v(*fe2, quadrature, flags); + + unsigned int cell_index = 0; + for (const auto &cell : dh2.active_cell_iterators()) + if (cell->is_locally_owned()) + { + fe_v.reinit(cell); + const auto &points = fe_v.get_quadrature_points(); + const auto &JxW = fe_v.get_JxW_values(); + + for (unsigned int q = 0; q < points.size(); ++q) + { + const auto i = cell_index * points.size() + q; + quadrature_points_vec[i] = points[q]; + properties[i][0] = JxW[q]; + if (dim < spacedim) + for (unsigned int d = 0; d < spacedim; ++d) + { + properties[i][d + 1] = fe_v.normal_vector(q)[d]; + } + } + ++cell_index; + } + + // Distribute the local points to the processor that owns + // them on the triangulation + auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box( + tria1, IteratorFilters::LocallyOwnedCell()); + + auto global_bounding_boxes = + Utilities::MPI::all_gather(mpi_communicator, my_bounding_box); + + auto cpu_to_index = + solid_particle_handler.insert_global_particles(quadrature_points_vec, + global_bounding_boxes, + properties); + } + + template + void StokesImmersedProblem::setup_system() + { + TimerOutput::Scope t(computing_timer, "setup"); + + fe1 = + std::make_unique>(FE_Q(par.velocity_degree), + spacedim, + FE_Q(par.velocity_degree - + 1), + 1); + + fe2 = std::make_unique>( + FE_Q(par.velocity_degree), spacedim); + + quadrature_formula = std::make_unique>(par.velocity_degree + 1); + + dh1.distribute_dofs(*fe1); + dh2.distribute_dofs(*fe2); + + std::vector stokes_sub_blocks(dim + 1, 0); + stokes_sub_blocks[dim] = 1; + DoFRenumbering::component_wise(dh1, stokes_sub_blocks); + + auto dofs_per_block = + DoFTools::count_dofs_per_fe_block(dh1, stokes_sub_blocks); + + const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1]; + + pcout << " Number of degrees of freedom: " << dh1.n_dofs() << " (" << n_u + << '+' << n_p << ')' << std::endl; + + owned1.resize(2); + owned1[0] = dh1.locally_owned_dofs().get_view(0, n_u); + owned1[1] = dh1.locally_owned_dofs().get_view(n_u, n_u + n_p); + + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dh1, locally_relevant_dofs); + relevant1.resize(2); + relevant1[0] = locally_relevant_dofs.get_view(0, n_u); + relevant1[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p); + + { + constraints.reinit(locally_relevant_dofs); + + FEValuesExtractors::Vector velocities(0); + DoFTools::make_hanging_node_constraints(dh1, constraints); + VectorTools::interpolate_boundary_values(dh1, + 0, + ZeroFunction(spacedim + + 1), + constraints, + fe1->component_mask(velocities)); + constraints.close(); + } + + { + system_matrix.clear(); + + 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::none; + else if (c == dim || d == dim || c == d) + coupling[c][d] = DoFTools::always; + else + coupling[c][d] = DoFTools::none; + + BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block); + + DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false); + + SparsityTools::distribute_sparsity_pattern( + dsp, + dh1.compute_locally_owned_dofs_per_processor(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(owned1, dsp, mpi_communicator); + } + + { + preconditioner_matrix.clear(); + + 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; + + BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block); + + DoFTools::make_sparsity_pattern(dh1, coupling, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern( + dsp, + dh1.compute_locally_owned_dofs_per_processor(), + mpi_communicator, + locally_relevant_dofs); + preconditioner_matrix.reinit(owned1, dsp, mpi_communicator); + } + + locally_relevant_solution.reinit(owned1, relevant1, mpi_communicator); + system_rhs.reinit(owned1, mpi_communicator); + solution.reinit(owned1, mpi_communicator); + } + + + + template + void StokesImmersedProblem::assemble_stokes_system() + { + system_matrix = 0; + preconditioner_matrix = 0; + system_rhs = 0; + + TimerOutput::Scope t(computing_timer, "Stokes_assembly"); + + + FEValues fe_values(*fe1, + *quadrature_formula, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe1->dofs_per_cell; + const unsigned int n_q_points = quadrature_formula->size(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix cell_matrix2(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector> rhs_values(n_q_points, + Vector(spacedim + 1)); + + std::vector> grad_phi_u(dofs_per_cell); + std::vector div_phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(spacedim); + + for (const auto &cell : dh1.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell_matrix = 0; + cell_matrix2 = 0; + cell_rhs = 0; + + fe_values.reinit(cell); + par.rhs.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) + { + grad_phi_u[k] = fe_values[velocities].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 < dofs_per_cell; ++j) + { + cell_matrix(i, j) += + (par.viscosity * + scalar_product(grad_phi_u[i], grad_phi_u[j]) - + div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) * + fe_values.JxW(q); + + cell_matrix2(i, j) += 1.0 / par.viscosity * phi_p[i] * + phi_p[j] * fe_values.JxW(q); + } + + const unsigned int component_i = + fe1->system_to_component_index(i).first; + cell_rhs(i) += fe_values.shape_value(i, q) * + rhs_values[q](component_i) * fe_values.JxW(q); + } + } + + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + + constraints.distribute_local_to_global(cell_matrix2, + local_dof_indices, + preconditioner_matrix); + } + + system_matrix.compress(VectorOperation::add); + preconditioner_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + + + template + void StokesImmersedProblem::assemble_nitche_restriction() + { + TimerOutput::Scope t(computing_timer, "Nitsche_assembly"); + + SolidVelocity solid_velocity(par.angular_velocity); + + std::vector dof_indices1(fe1->dofs_per_cell); + + FullMatrix local_matrix(fe1->dofs_per_cell, fe1->dofs_per_cell); + dealii::Vector local_rhs(fe1->dofs_per_cell); + + auto particle = solid_particle_handler.begin(); + while (particle != solid_particle_handler.end()) + { + local_matrix = 0; + local_rhs = 0; + const auto &cell = particle->get_surrounding_cell(tria1); + const auto &dh_cell = + typename DoFHandler::cell_iterator(*cell, &dh1); + dh_cell->get_dof_indices(dof_indices1); + + const auto pic = solid_particle_handler.particles_in_cell(cell); + Assert(pic.begin() == particle, ExcInternalError()); + for (const auto &p : pic) + { + const auto ref_q = p.get_reference_location(); + const auto real_q = p.get_location(); + const auto properties = p.get_properties(); + const auto &JxW = properties[0]; + for (unsigned int i = 0; i < fe1->dofs_per_cell; ++i) + { + const auto comp_i = fe1->system_to_component_index(i).first; + if (comp_i < spacedim) + { + for (unsigned int j = 0; j < fe1->dofs_per_cell; ++j) + { + const auto comp_j = + fe1->system_to_component_index(j).first; + if (comp_i == comp_j) + local_matrix(i, j) += + par.penalty_term * fe1->shape_value(i, ref_q) * + fe1->shape_value(j, ref_q) * JxW; + } + local_rhs(i) += par.penalty_term * + solid_velocity.value(real_q, comp_i) * + fe1->shape_value(i, ref_q) * JxW; + } + } + } + constraints.distribute_local_to_global( + local_matrix, local_rhs, dof_indices1, system_matrix, system_rhs); + particle = pic.end(); + } + } + + + + template + void StokesImmersedProblem::solve() + { + TimerOutput::Scope t(computing_timer, "solve"); + + LA::MPI::PreconditionAMG prec_A; + { + LA::MPI::PreconditionAMG::AdditionalData data; + +#ifdef USE_PETSC_LA + data.symmetric_operator = true; +#endif + prec_A.initialize(system_matrix.block(0, 0), data); + } + + LA::MPI::PreconditionAMG prec_S; + { + LA::MPI::PreconditionAMG::AdditionalData data; + +#ifdef USE_PETSC_LA + data.symmetric_operator = true; +#endif + prec_S.initialize(preconditioner_matrix.block(1, 1), data); + } + + const auto A = linear_operator(system_matrix.block(0, 0)); + const auto amgA = linear_operator(A, prec_A); + + const auto S = + linear_operator(preconditioner_matrix.block(1, 1)); + const auto amgS = linear_operator(S, prec_S); + + ReductionControl inner_solver_control(10, + 1e-8 * system_rhs.l2_norm(), + 1.e-2); + SolverCG cg(inner_solver_control); + + const auto invS = inverse_operator(S, cg, amgS); + + const auto P = + block_diagonal_operator<2, LA::MPI::BlockVector>({amgA, amgS}); + + SolverControl solver_control(system_matrix.m(), + 1e-10 * system_rhs.l2_norm()); + + SolverMinRes solver(solver_control); + + constraints.set_zero(solution); + + solver.solve(system_matrix, solution, system_rhs, P); + + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(solution); + + locally_relevant_solution = solution; + const double mean_pressure = + VectorTools::compute_mean_value(dh1, + QGauss(par.velocity_degree + 2), + locally_relevant_solution, + spacedim); + solution.block(1).add(-mean_pressure); + locally_relevant_solution.block(1) = solution.block(1); + } + + + + template + void StokesImmersedProblem::refine_grid() + { + TimerOutput::Scope t(computing_timer, "refine"); + + tria1.refine_global(); + } + + + + template + void StokesImmersedProblem::output_results( + const unsigned int cycle) const + { + TimerOutput::Scope t(computing_timer, "Output fluid"); + std::vector solution_names(spacedim, "velocity"); + solution_names.emplace_back("pressure"); + std::vector + data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.attach_dof_handler(dh1); + data_out.add_data_vector(locally_relevant_solution, + solution_names, + DataOut::type_dof_data, + data_component_interpretation); + + LA::MPI::BlockVector interpolated; + interpolated.reinit(owned1, MPI_COMM_WORLD); + VectorTools::interpolate(dh1, + ConstantFunction(1.0, spacedim + 1), + interpolated); + + LA::MPI::BlockVector interpolated_relevant(owned1, + relevant1, + MPI_COMM_WORLD); + interpolated_relevant = interpolated; + { + std::vector solution_names(dim, "ref_u"); + solution_names.emplace_back("ref_p"); + data_out.add_data_vector(interpolated_relevant, + solution_names, + DataOut::type_dof_data, + data_component_interpretation); + } + + + Vector subdomain(tria1.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = tria1.locally_owned_subdomain(); + data_out.add_data_vector(subdomain, "subdomain"); + + data_out.build_patches(); + + const std::string filename = + ("solution-" + Utilities::int_to_string(cycle, 2) + "." + + Utilities::int_to_string(tria1.locally_owned_subdomain(), 4)); + std::ofstream output((filename + ".vtu")); + data_out.write_vtu(output); + + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + std::vector filenames; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + filenames.push_back("solution-" + Utilities::int_to_string(cycle, 2) + + "." + Utilities::int_to_string(i, 4) + ".vtu"); + + std::ofstream master_output( + "solution-" + Utilities::int_to_string(cycle, 2) + ".pvtu"); + data_out.write_pvtu_record(master_output, filenames); + } + } + + template + void StokesImmersedProblem::output_particles( + const Particles::ParticleHandler &particles, + std::string fprefix, + const unsigned int iter) const + { + Particles::DataOut particles_out; + particles_out.build_patches(particles); + const std::string filename = + (fprefix + "-" + Utilities::int_to_string(iter, 2) + "." + + Utilities::int_to_string(tria1.locally_owned_subdomain(), 4)); + std::ofstream output((filename + ".vtu")); + particles_out.write_vtu(output); + + if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + std::vector filenames; + for (unsigned int i = 0; + i < Utilities::MPI::n_mpi_processes(mpi_communicator); + ++i) + filenames.push_back(fprefix + "-" + + Utilities::int_to_string(iter, 2) + "." + + Utilities::int_to_string(i, 4) + ".vtu"); + + std::ofstream master_output( + fprefix + "-" + Utilities::int_to_string(iter, 2) + ".pvtu"); + particles_out.write_pvtu_record(master_output, filenames); + } + } + + + template + void StokesImmersedProblem::run() + { +#ifdef USE_PETSC_LA + pcout << "Running using PETSc." << std::endl; +#else + pcout << "Running using Trilinos." << std::endl; +#endif + + ComponentMask velocity_mask(spacedim + 1, true); + velocity_mask.set(spacedim, false); + + const double time_step = par.final_time / (par.number_of_time_steps - 1); + double time = 0; + for (unsigned int cycle = 0; cycle < par.number_of_time_steps; + ++cycle, time += time_step) + { + par.set_time(time); + pcout << "Cycle " << cycle << ':' << std::endl + << "Time : " << time << ", time step: " << time_step << std::endl; + + if (cycle == 0) + { + make_grid(); + setup_system(); + setup_tracer_particles(); + setup_solid_particles(); + tracer_particle_velocities.reinit(owned_tracer_particles, + mpi_communicator); + } + else + { + TimerOutput::Scope t(computing_timer, + "Set solid particle position"); + + SolidDisplacement solid_displacement(par.angular_velocity, + time_step); + solid_particle_handler.set_particle_positions(solid_displacement, + false); + } + { + TimerOutput::Scope t(computing_timer, "Set tracer particle motion"); + interpolate_field_on_particles(dh1, + tracer_particle_handler, + locally_relevant_solution, + tracer_particle_velocities, + velocity_mask); + + tracer_particle_velocities *= time_step; + + relevant_tracer_particles = + tracer_particle_handler.locally_relevant_ids().tensor_product( + complete_index_set(spacedim)); + + relevant_tracer_particle_displacements.reinit( + owned_tracer_particles, + relevant_tracer_particles, + mpi_communicator); + + relevant_tracer_particle_displacements = tracer_particle_velocities; + + tracer_particle_handler.set_particle_positions( + relevant_tracer_particle_displacements); + } + assemble_stokes_system(); + assemble_nitche_restriction(); + solve(); + + if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32) + { + output_results(cycle); + { + TimerOutput::Scope t(computing_timer, "Output tracer particles"); + output_particles(tracer_particle_handler, "tracer", cycle); + } + { + TimerOutput::Scope t(computing_timer, "Output solid particles"); + output_particles(solid_particle_handler, "solid", cycle); + } + } + } + } +} // namespace Step70 + + + +int main(int argc, char *argv[]) +{ + using namespace Step70; + using namespace dealii; + deallog.depth_console(1); + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + StokesImmersedProblemParameters<2> par; + par.declare_all_parameters(); + std::ofstream out("default.prm"); + par.prm.print_parameters(out, ParameterHandler::ShortText); + ParameterAcceptor::initialize("parameters.prm", "used_parameters.prm"); + + StokesImmersedProblem<2> problem(par); + 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; +}