From 711a5c4cb9e710b7eed00375c6a30747ca97aa28 Mon Sep 17 00:00:00 2001 From: Jiaqi Zhang Date: Thu, 11 Aug 2022 00:30:31 -0400 Subject: [PATCH] stokes test --- .../numerics/{no_flux_16.cc => no_flux_17.cc} | 0 .../{no_flux_16.output => no_flux_17.output} | 0 tests/numerics/no_flux_18.cc | 1519 +++++++++++++++++ ...apack=true.with_p4est=true.mpirun=4.output | 19 + 4 files changed, 1538 insertions(+) rename tests/numerics/{no_flux_16.cc => no_flux_17.cc} (100%) rename tests/numerics/{no_flux_16.output => no_flux_17.output} (100%) create mode 100644 tests/numerics/no_flux_18.cc create mode 100644 tests/numerics/no_flux_18.with_lapack=true.with_p4est=true.mpirun=4.output diff --git a/tests/numerics/no_flux_16.cc b/tests/numerics/no_flux_17.cc similarity index 100% rename from tests/numerics/no_flux_16.cc rename to tests/numerics/no_flux_17.cc diff --git a/tests/numerics/no_flux_16.output b/tests/numerics/no_flux_17.output similarity index 100% rename from tests/numerics/no_flux_16.output rename to tests/numerics/no_flux_17.output diff --git a/tests/numerics/no_flux_18.cc b/tests/numerics/no_flux_18.cc new file mode 100644 index 0000000000..6813fdec13 --- /dev/null +++ b/tests/numerics/no_flux_18.cc @@ -0,0 +1,1519 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2012 - 2022 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. +// +// --------------------------------------------------------------------- + + +// copied from matrix_free/stokes_computation.cc, but use no normal flux +// constraints on all boundaries instead of zero velocity. + + +#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 "../tests.h" + + +unsigned int minlevel = 0; +const unsigned int velocity_degree = 2; + +double pressure_scaling = 1.0; + +namespace StokesClass +{ + class QuietException + {}; + + namespace StokesSolver + { + /** + * Implement the block Schur preconditioner for the Stokes system. + */ + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + /** + * brief Constructor + * + * S The entire Stokes matrix + * Spre The matrix whose blocks are used in the definition of + * the preconditioning of the Stokes matrix, i.e. containing + *approximations of the A and S blocks. Mppreconditioner Preconditioner + *object for the Schur complement, typically chosen as the mass matrix. + * Apreconditioner Preconditioner object for the matrix A. + * do_solve_A A flag indicating whether we should actually solve with + * the matrix $A$, or only apply one preconditioner step with it. + * A_block_tolerance The tolerance for the CG solver which computes + * the inverse of the A block. + * S_block_tolerance The tolerance for the CG solver which computes + * the inverse of the S block (Schur complement matrix). + **/ + BlockSchurPreconditioner(const StokesMatrixType &S, + const MassMatrixType & Mass, + const PreconditionerMp &Mppreconditioner, + const PreconditionerA & Apreconditioner, + const bool do_solve_A, + const double A_block_tolerance, + const double S_block_tolerance); + + /** + * Matrix vector product with this preconditioner object. + */ + void + vmult(LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src) const; + + unsigned int + n_iterations_A() const; + unsigned int + n_iterations_S() const; + + private: + /** + * References to the various matrix object this preconditioner works on. + */ + const StokesMatrixType &stokes_matrix; + const MassMatrixType & mass_matrix; + const PreconditionerMp &mp_preconditioner; + const PreconditionerA & a_preconditioner; + + /** + * Whether to actually invert the $\tilde A$ part of the preconditioner + *matrix or to just apply a single preconditioner step with it. + **/ + const bool do_solve_A; + mutable unsigned int n_iterations_A_; + mutable unsigned int n_iterations_S_; + const double A_block_tolerance; + const double S_block_tolerance; + }; + + + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner(const StokesMatrixType &S, + const MassMatrixType & Mass, + const PreconditionerMp &Mppreconditioner, + const PreconditionerA & Apreconditioner, + const bool do_solve_A, + const double A_block_tolerance, + const double S_block_tolerance) + : stokes_matrix(S) + , mass_matrix(Mass) + , mp_preconditioner(Mppreconditioner) + , a_preconditioner(Apreconditioner) + , do_solve_A(do_solve_A) + , n_iterations_A_(0) + , n_iterations_S_(0) + , A_block_tolerance(A_block_tolerance) + , S_block_tolerance(S_block_tolerance) + {} + + template + unsigned int + BlockSchurPreconditioner::n_iterations_A() const + { + return n_iterations_A_; + } + + template + unsigned int + BlockSchurPreconditioner::n_iterations_S() const + { + return n_iterations_S_; + } + + template + void + BlockSchurPreconditioner:: + vmult(LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src) const + { + LinearAlgebra::distributed::BlockVector utmp(src); + + // first solve with the bottom left block, which we have built + // as a mass matrix with the inverse of the viscosity + { + SolverControl solver_control(1000, + src.block(1).l2_norm() * S_block_tolerance, + false, + false); + + SolverCG> solver( + solver_control); + try + { + dst.block(1) = 0.0; + solver.solve(mass_matrix, + dst.block(1), + src.block(1), + mp_preconditioner); + n_iterations_S_ += solver_control.last_step(); + } + // if the solver fails, report the error from processor 0 with some + // additional information about its location, and throw a quiet + // exception on all other processors + catch (const std::exception &exc) + { + if (Utilities::MPI::this_mpi_process( + src.block(0).get_mpi_communicator()) == 0) + AssertThrow( + false, + ExcMessage( + std::string( + "The iterative (bottom right) solver in BlockSchurPreconditioner::vmult " + "did not converge to a tolerance of " + + Utilities::to_string(solver_control.tolerance()) + + ". It reported the following error:\n\n") + + exc.what())) else throw QuietException(); + } + dst.block(1) *= -1.0; + } + + // apply the top right block + { + LinearAlgebra::distributed::BlockVector dst_tmp(dst); + dst_tmp.block(0) *= 0.0; + stokes_matrix.vmult(utmp, dst_tmp); // B^T + utmp.block(0) *= -1.0; + utmp.block(0) += src.block(0); + } + + // now either solve with the top left block (if do_solve_A==true) + // or just apply one preconditioner sweep (for the first few + // iterations of our two-stage outer GMRES iteration) + if (do_solve_A == true) + { + Assert(false, ExcNotImplemented()); + } + else + { + a_preconditioner.vmult(dst.block(0), utmp.block(0)); + n_iterations_A_ += 1; + } + } + } // namespace StokesSolver + + + // Parameters for Sinker example + double beta = 10.0; + double delta = 200.0; + double omega = 0.1; + + template + struct Sinker + { + unsigned int problem_dim; + unsigned int n_sinkers; + std::vector> centers; + double DR_mu; + double mu_min; + double mu_max; + }; + + template + class Viscosity + { + public: + Viscosity(const Sinker &sink); + virtual double + value(const Point &p, const unsigned int component = 0) const; + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const; + + Sinker sinker; + }; + template + Viscosity::Viscosity(const Sinker &sink) + { + sinker = sink; + } + template + double + Viscosity::value(const Point &p, + const unsigned int /*component*/) const + { + double Chi = 1.0; + for (unsigned int s = 0; s < sinker.n_sinkers; ++s) + { + double dist = p.distance(sinker.centers[s]); + double temp = + 1 - std::exp(-delta * std::pow(std::max(0.0, dist - omega / 2.0), 2)); + Chi *= temp; + } + return (sinker.mu_max - sinker.mu_min) * (1 - Chi) + sinker.mu_min; + } + template + void + Viscosity::value_list(const std::vector> &points, + std::vector & values, + const unsigned int component) const + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + Assert(component == 0, ExcIndexRange(component, 0, 1)); + const unsigned int n_points = points.size(); + for (unsigned int i = 0; i < n_points; ++i) + values[i] = value(points[i], component); + } + + template + class RightHandSide + { + public: + RightHandSide(const Sinker &sink); + Sinker sinker; + virtual void + vector_value(const Point &p, Vector &value) const; + }; + template + RightHandSide::RightHandSide(const Sinker &sink) + { + sinker = sink; + } + template + void + RightHandSide::vector_value(const Point &p, + Vector & values) const + { + double Chi = 1.0; + for (unsigned int s = 0; s < sinker.n_sinkers; ++s) + { + double dist = p.distance(sinker.centers[s]); + double temp = + 1 - std::exp(-delta * std::pow(std::max(0.0, dist - omega / 2.0), 2)); + Chi *= temp; + } + + if (sinker.problem_dim == 2) + { + values[0] = 0; + values[1] = beta * (Chi - 1.0); + values[2] = 0; + } + else if (sinker.problem_dim == 3) + { + values[0] = 0; + values[1] = 0; + values[2] = beta * (Chi - 1.0); + values[3] = 0; + } + return; + } + + template + class ExactSolution_BoundaryValues : public Function + { + public: + ExactSolution_BoundaryValues() + : Function(dim + 1) + {} + virtual void + vector_value(const Point &p, Vector &value) const; + }; + template + void + ExactSolution_BoundaryValues::vector_value(const Point &p, + Vector &values) const + { + (void)p; + for (unsigned int i = 0; i < values.size(); ++i) + values(i) = 0.0; + return; + } + + template + class ExactSolution_BoundaryValues_u : public Function + { + public: + ExactSolution_BoundaryValues_u() + : Function(dim) + {} + virtual void + vector_value(const Point &p, Vector &value) const; + }; + template + void + ExactSolution_BoundaryValues_u::vector_value( + const Point &p, + Vector & values) const + { + (void)p; + for (unsigned int i = 0; i < values.size(); ++i) + values(i) = 0.0; + return; + } + + + + template + class StokesOperator + : public MatrixFreeOperators:: + Base> + { + public: + StokesOperator() + : MatrixFreeOperators:: + Base>() + {} + void + clear(); + void + evaluate_2_x_viscosity(const Viscosity &viscosity_function); + virtual void + compute_diagonal(); + + private: + virtual void + apply_add(LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src) const; + + void + local_apply(const dealii::MatrixFree & data, + LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src, + const std::pair &cell_range) const; + + Table<2, VectorizedArray> viscosity_x_2; + }; + template + void + StokesOperator::clear() + { + viscosity_x_2.reinit(0, 0); + MatrixFreeOperators:: + Base>::clear(); + } + template + void + StokesOperator::evaluate_2_x_viscosity( + const Viscosity &viscosity_function) + { + const unsigned int n_cells = this->data->n_cell_batches(); + FEEvaluation velocity(*this->data, + 0); + viscosity_x_2.reinit(n_cells, velocity.n_q_points); + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + velocity.reinit(cell); + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + VectorizedArray return_value = + make_vectorized_array(1.); + for (unsigned int i = 0; i < VectorizedArray::size(); ++i) + { + Point p; + for (unsigned int d = 0; d < dim; ++d) + { + p(d) = velocity.quadrature_point(q)(d)[i]; + } + return_value[i] = 2.0 * viscosity_function.value(p); + } + viscosity_x_2(cell, q) = return_value; + } + } + } + template + void + StokesOperator::local_apply( + const dealii::MatrixFree & data, + LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src, + const std::pair & cell_range) const + { + using vector_t = VectorizedArray; + FEEvaluation velocity(data, 0); + FEEvaluation pressure(data, 1); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + velocity.reinit(cell); + velocity.read_dof_values(src.block(0)); + velocity.evaluate(EvaluationFlags::gradients); + pressure.reinit(cell); + pressure.read_dof_values(src.block(1)); + pressure.evaluate(EvaluationFlags::values); + + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + SymmetricTensor<2, dim, vector_t> sym_grad_u = + velocity.get_symmetric_gradient(q); + vector_t pres = pressure.get_value(q); + vector_t div = -trace(sym_grad_u); + pressure.submit_value(div, q); + + sym_grad_u *= viscosity_x_2(cell, q); + // subtract p * I + for (unsigned int d = 0; d < dim; ++d) + sym_grad_u[d][d] -= pres; + + velocity.submit_symmetric_gradient(sym_grad_u, q); + } + + velocity.integrate(EvaluationFlags::gradients); + velocity.distribute_local_to_global(dst.block(0)); + pressure.integrate(EvaluationFlags::values); + pressure.distribute_local_to_global(dst.block(1)); + } + } + template + void + StokesOperator::apply_add( + LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src) const + { + MatrixFreeOperators:: + Base>::data + ->cell_loop(&StokesOperator::local_apply, this, dst, src); + } + template + void + StokesOperator::compute_diagonal() + { + Assert(false, ExcNotImplemented()); + } + + + template + class MassMatrixOperator + : public MatrixFreeOperators:: + Base> + { + public: + MassMatrixOperator() + : MatrixFreeOperators::Base>() + {} + void + clear(); + void + evaluate_1_over_viscosity(const Viscosity &viscosity_function); + virtual void + compute_diagonal(); + + private: + virtual void + apply_add(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const; + + void + local_apply(const dealii::MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const; + + void + local_compute_diagonal( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const unsigned int & dummy, + const std::pair &cell_range) const; + + Table<2, VectorizedArray> one_over_viscosity; + }; + template + void + MassMatrixOperator::clear() + { + one_over_viscosity.reinit(0, 0); + MatrixFreeOperators::Base>:: + clear(); + } + template + void + MassMatrixOperator::evaluate_1_over_viscosity( + const Viscosity &viscosity_function) + { + const unsigned int n_cells = this->data->n_cell_batches(); + FEEvaluation pressure(*this->data, + 0); + one_over_viscosity.reinit(n_cells, pressure.n_q_points); + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + pressure.reinit(cell); + for (unsigned int q = 0; q < pressure.n_q_points; ++q) + { + VectorizedArray return_value = + make_vectorized_array(1.); + for (unsigned int i = 0; i < VectorizedArray::size(); ++i) + { + Point p; + for (unsigned int d = 0; d < dim; ++d) + p(d) = pressure.quadrature_point(q)(d)[i]; + return_value[i] = 1.0 / viscosity_function.value(p); + } + one_over_viscosity(cell, q) = return_value; + } + } + } + template + void + MassMatrixOperator::local_apply( + const dealii::MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const + { + FEEvaluation pressure(data); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + AssertDimension(one_over_viscosity.size(0), data.n_cell_batches()); + AssertDimension(one_over_viscosity.size(1), pressure.n_q_points); + + pressure.reinit(cell); + pressure.read_dof_values(src); + pressure.evaluate(EvaluationFlags::values); + for (unsigned int q = 0; q < pressure.n_q_points; ++q) + pressure.submit_value(one_over_viscosity(cell, q) * + pressure.get_value(q), + q); + pressure.integrate(EvaluationFlags::values); + pressure.distribute_local_to_global(dst); + } + } + template + void + MassMatrixOperator::apply_add( + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + MatrixFreeOperators::Base>::data + ->cell_loop(&MassMatrixOperator::local_apply, this, dst, src); + } + template + void + MassMatrixOperator::compute_diagonal() + { + this->inverse_diagonal_entries.reset( + new DiagonalMatrix>()); + this->diagonal_entries.reset( + new DiagonalMatrix>()); + + LinearAlgebra::distributed::Vector &inverse_diagonal = + this->inverse_diagonal_entries->get_vector(); + LinearAlgebra::distributed::Vector &diagonal = + this->diagonal_entries->get_vector(); + + unsigned int dummy = 0; + this->data->initialize_dof_vector(inverse_diagonal); + this->data->initialize_dof_vector(diagonal); + + this->data->cell_loop(&MassMatrixOperator::local_compute_diagonal, + this, + diagonal, + dummy); + + this->set_constrained_entries_to_one(diagonal); + inverse_diagonal = diagonal; + const unsigned int local_size = inverse_diagonal.local_size(); + for (unsigned int i = 0; i < local_size; ++i) + { + Assert(inverse_diagonal.local_element(i) > 0., + ExcMessage("No diagonal entry in a positive definite operator " + "should be zero")); + inverse_diagonal.local_element(i) = + 1. / inverse_diagonal.local_element(i); + } + } + template + void + MassMatrixOperator::local_compute_diagonal( + const MatrixFree & data, + LinearAlgebra::distributed::Vector &dst, + const unsigned int &, + const std::pair &cell_range) const + { + FEEvaluation pressure(data, 0); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + pressure.reinit(cell); + AlignedVector> diagonal(pressure.dofs_per_cell); + for (unsigned int i = 0; i < pressure.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < pressure.dofs_per_cell; ++j) + pressure.begin_dof_values()[j] = VectorizedArray(); + pressure.begin_dof_values()[i] = make_vectorized_array(1.); + + pressure.evaluate(EvaluationFlags::values); + for (unsigned int q = 0; q < pressure.n_q_points; ++q) + pressure.submit_value(one_over_viscosity(cell, q) * + pressure.get_value(q), + q); + pressure.integrate(EvaluationFlags::values); + + diagonal[i] = pressure.begin_dof_values()[i]; + } + + for (unsigned int i = 0; i < pressure.dofs_per_cell; ++i) + pressure.begin_dof_values()[i] = diagonal[i]; + pressure.distribute_local_to_global(dst); + } + } + + + template + class ABlockOperator : public MatrixFreeOperators:: + Base> + { + public: + ABlockOperator() + : MatrixFreeOperators::Base>() + {} + void + clear(); + void + evaluate_2_x_viscosity(const Viscosity &viscosity_function); + virtual void + compute_diagonal(); + + private: + virtual void + apply_add(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const; + + void + local_apply(const dealii::MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const; + + void + local_compute_diagonal( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const unsigned int & dummy, + const std::pair &cell_range) const; + + Table<2, VectorizedArray> viscosity_x_2; + }; + template + void + ABlockOperator::clear() + { + viscosity_x_2.reinit(0, 0); + MatrixFreeOperators::Base>:: + clear(); + } + template + void + ABlockOperator::evaluate_2_x_viscosity( + const Viscosity &viscosity_function) + { + const unsigned int n_cells = this->data->n_cell_batches(); + FEEvaluation velocity(*this->data, + 0); + viscosity_x_2.reinit(n_cells, velocity.n_q_points); + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + velocity.reinit(cell); + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + VectorizedArray return_value = + make_vectorized_array(1.); + for (unsigned int i = 0; i < VectorizedArray::size(); ++i) + { + Point p; + for (unsigned int d = 0; d < dim; ++d) + { + p(d) = velocity.quadrature_point(q)(d)[i]; + } + return_value[i] = 2.0 * viscosity_function.value(p); + } + viscosity_x_2(cell, q) = return_value; + } + } + } + template + void + ABlockOperator::local_apply( + const dealii::MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const + { + FEEvaluation velocity(data); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + AssertDimension(viscosity_x_2.size(0), data.n_cell_batches()); + AssertDimension(viscosity_x_2.size(1), velocity.n_q_points); + + velocity.reinit(cell); + velocity.read_dof_values(src); + velocity.evaluate(EvaluationFlags::gradients); + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + velocity.submit_symmetric_gradient( + viscosity_x_2(cell, q) * velocity.get_symmetric_gradient(q), q); + } + velocity.integrate(EvaluationFlags::gradients); + velocity.distribute_local_to_global(dst); + } + } + template + void + ABlockOperator::apply_add( + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + MatrixFreeOperators::Base>::data + ->cell_loop(&ABlockOperator::local_apply, this, dst, src); + } + template + void + ABlockOperator::compute_diagonal() + { + this->inverse_diagonal_entries.reset( + new DiagonalMatrix>()); + LinearAlgebra::distributed::Vector &inverse_diagonal = + this->inverse_diagonal_entries->get_vector(); + this->data->initialize_dof_vector(inverse_diagonal); + unsigned int dummy = 0; + this->data->cell_loop(&ABlockOperator::local_compute_diagonal, + this, + inverse_diagonal, + dummy); + + this->set_constrained_entries_to_one(inverse_diagonal); + + for (unsigned int i = 0; i < inverse_diagonal.local_size(); ++i) + { + Assert(inverse_diagonal.local_element(i) > 0., + ExcMessage("No diagonal entry in a positive definite operator " + "should be zero")); + inverse_diagonal.local_element(i) = + 1. / inverse_diagonal.local_element(i); + } + } + template + void + ABlockOperator::local_compute_diagonal( + const MatrixFree & data, + LinearAlgebra::distributed::Vector &dst, + const unsigned int &, + const std::pair &cell_range) const + { + FEEvaluation velocity(data, 0); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + velocity.reinit(cell); + AlignedVector> diagonal(velocity.dofs_per_cell); + for (unsigned int i = 0; i < velocity.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < velocity.dofs_per_cell; ++j) + velocity.begin_dof_values()[j] = VectorizedArray(); + velocity.begin_dof_values()[i] = make_vectorized_array(1.); + + velocity.evaluate(EvaluationFlags::gradients); + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + velocity.submit_symmetric_gradient( + viscosity_x_2(cell, q) * velocity.get_symmetric_gradient(q), + q); + } + velocity.integrate(EvaluationFlags::gradients); + + diagonal[i] = velocity.begin_dof_values()[i]; + } + + for (unsigned int i = 0; i < velocity.dofs_per_cell; ++i) + velocity.begin_dof_values()[i] = diagonal[i]; + velocity.distribute_local_to_global(dst); + } + } + + + + template + class StokesProblem + { + public: + StokesProblem(); + + void + run(); + + private: + void + make_grid(const unsigned int ref = 4); + void + create_sinker(const unsigned int n_sinkers, const double visc_jump); + void + setup_system(); + void + assemble_rhs(); + void + solve(); + + using vector_t = LinearAlgebra::distributed::Vector; + using block_vector_t = LinearAlgebra::distributed::BlockVector; + + using StokesMatrixType = StokesOperator; + using MassMatrixType = MassMatrixOperator; + using LevelMatrixType = ABlockOperator; + + unsigned int degree_u; + + FESystem fe_u; + FE_Q fe_p; + + parallel::distributed::Triangulation triangulation; + DoFHandler dof_handler_u; + DoFHandler dof_handler_p; + + std::vector owned_partitioning; + std::vector relevant_partitioning; + + AffineConstraints constraints_u; + AffineConstraints constraints_p; + + block_vector_t solution; + block_vector_t system_rhs; + + StokesMatrixType stokes_matrix; + MassMatrixType mass_matrix; + + MGLevelObject mg_matrices; + MGConstrainedDoFs mg_constrained_dofs; + + Sinker sinker; + }; + + + + template + StokesProblem::StokesProblem() + : degree_u(velocity_degree) + , fe_u(FE_Q(degree_u), dim) + , fe_p(FE_Q(degree_u - 1)) + , triangulation(MPI_COMM_WORLD, + typename Triangulation::MeshSmoothing( + Triangulation::limit_level_difference_at_vertices | + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening), + parallel::distributed::Triangulation< + dim>::construct_multigrid_hierarchy) + , dof_handler_u(triangulation) + , dof_handler_p(triangulation) + {} + + template + void + StokesProblem::create_sinker(const unsigned int n_sinkers, + const double visc_jump) + { + sinker.problem_dim = dim; + sinker.n_sinkers = n_sinkers; + std::srand(171); + for (unsigned int s = 0; s < sinker.n_sinkers; ++s) + { + std::vector coords(dim); + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + for (unsigned int i = 0; i < dim; ++i) + coords[i] = std::rand() / (double)RAND_MAX; + + MPI_Bcast(&(coords[0]), dim, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + Tensor<1, dim, double> coords_tens; + for (unsigned int i = 0; i < dim; ++i) + coords_tens[i] = coords[i]; + sinker.centers.push_back(Point(coords_tens)); + } + + sinker.DR_mu = visc_jump; + sinker.mu_max = std::sqrt(sinker.DR_mu); + sinker.mu_min = 1.0 / std::sqrt(sinker.DR_mu); + } + + + template + void + StokesProblem::make_grid(const unsigned int ref) + { + GridGenerator::hyper_cube(triangulation, 0, 1); + triangulation.refine_global(ref); + } + + + template + void + StokesProblem::setup_system() + { + dof_handler_u.clear(); + dof_handler_u.distribute_dofs(fe_u); + dof_handler_u.distribute_mg_dofs(); + + dof_handler_p.clear(); + dof_handler_p.distribute_dofs(fe_p); + + IndexSet locally_relevant_dofs_u; + DoFTools::extract_locally_relevant_dofs(dof_handler_u, + locally_relevant_dofs_u); + constraints_u.reinit(locally_relevant_dofs_u); + DoFTools::make_hanging_node_constraints(dof_handler_u, constraints_u); + + // VectorTools::interpolate_boundary_values( + // dof_handler_u, 0, ExactSolution_BoundaryValues_u(), + // constraints_u); + + const std::set dirichlet_boundary = {0}; + VectorTools::compute_no_normal_flux_constraints(dof_handler_u, + 0, + dirichlet_boundary, + constraints_u); + constraints_u.close(); + + IndexSet locally_relevant_dofs_p; + DoFTools::extract_locally_relevant_dofs(dof_handler_p, + locally_relevant_dofs_p); + constraints_p.reinit(locally_relevant_dofs_p); + DoFTools::make_hanging_node_constraints(dof_handler_p, constraints_p); + constraints_p.close(); + + + // Stokes matrix stuff... + typename MatrixFree::AdditionalData additional_data_stokes; + additional_data_stokes.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data_stokes.mapping_update_flags = + (update_values | update_gradients | update_JxW_values | + update_quadrature_points); + + std::vector *> stokes_dofs; + stokes_dofs.push_back(&dof_handler_u); + stokes_dofs.push_back(&dof_handler_p); + std::vector *> stokes_constraints; + stokes_constraints.push_back(&constraints_u); + stokes_constraints.push_back(&constraints_p); + + std::shared_ptr> stokes_mf_storage( + new MatrixFree()); + stokes_mf_storage->reinit(MappingQ1{}, + stokes_dofs, + stokes_constraints, + QGauss<1>(degree_u + 1), + additional_data_stokes); + + stokes_matrix.initialize(stokes_mf_storage); + stokes_matrix.evaluate_2_x_viscosity(Viscosity(sinker)); + + // Mass matrix stuff... + typename MatrixFree::AdditionalData additional_data_mass; + additional_data_mass.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data_mass.mapping_update_flags = + (update_values | update_JxW_values | update_quadrature_points); + std::shared_ptr> mass_mf_storage( + new MatrixFree()); + mass_mf_storage->reinit(MappingQ1{}, + dof_handler_p, + constraints_p, + QGauss<1>(degree_u + 1), + additional_data_mass); + + mass_matrix.initialize(mass_mf_storage); + mass_matrix.evaluate_1_over_viscosity(Viscosity(sinker)); + mass_matrix.compute_diagonal(); + + // GMG stuff... + const unsigned int n_levels = triangulation.n_global_levels(); + mg_matrices.clear_elements(); + + mg_matrices.resize(0, n_levels - 1); + + mg_constrained_dofs.clear(); + // const std::set dirichlet_boundary = {0}; + mg_constrained_dofs.initialize(dof_handler_u); + // mg_constrained_dofs.make_zero_boundary_constraints(dof_handler_u, + // dirichlet_boundary); + + + for (unsigned int level = 0; level < n_levels; ++level) + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler_u, + level, + relevant_dofs); + AffineConstraints level_constraints; + level_constraints.reinit(relevant_dofs); + // level_constraints.add_lines( + // mg_constrained_dofs.get_boundary_indices(level)); + + MappingQ mapping(1); + const auto & refinement_edge_indices = + mg_constrained_dofs.get_refinement_edge_indices(level); + + VectorTools::compute_no_normal_flux_constraints(dof_handler_u, + 0, + dirichlet_boundary, + level_constraints, + mapping, + refinement_edge_indices, + level); + level_constraints.close(); + + mg_constrained_dofs.add_user_constraints(level, level_constraints); + + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + additional_data.mg_level = level; + std::shared_ptr> mg_mf_storage_level( + new MatrixFree()); + mg_mf_storage_level->reinit(MappingQ1{}, + dof_handler_u, + level_constraints, + QGauss<1>(degree_u + 1), + additional_data); + + mg_matrices[level].initialize(mg_mf_storage_level, + mg_constrained_dofs, + level); + mg_matrices[level].evaluate_2_x_viscosity(Viscosity(sinker)); + mg_matrices[level].compute_diagonal(); + } + + solution.reinit(2); + system_rhs.reinit(2); + + stokes_matrix.initialize_dof_vector(solution); + stokes_matrix.initialize_dof_vector(system_rhs); + + solution.update_ghost_values(); + solution.collect_sizes(); + + system_rhs.update_ghost_values(); + system_rhs.collect_sizes(); + } + + + + template + void + StokesProblem::assemble_rhs() + { + system_rhs = 0.0; + + // Create operator with no Dirchlet info + StokesMatrixType operator_homogeneous; + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = (update_values | update_gradients | + update_JxW_values | update_quadrature_points); + + // Create constraints with no Dirchlet info + AffineConstraints constraints_u_no_dirchlet; + IndexSet locally_relevant_dofs_u; + DoFTools::extract_locally_relevant_dofs(dof_handler_u, + locally_relevant_dofs_u); + constraints_u_no_dirchlet.reinit(locally_relevant_dofs_u); + DoFTools::make_hanging_node_constraints(dof_handler_u, + constraints_u_no_dirchlet); + constraints_u_no_dirchlet.close(); + + std::vector *> constraints_no_dirchlet; + constraints_no_dirchlet.push_back(&constraints_u_no_dirchlet); + constraints_no_dirchlet.push_back(&constraints_p); + std::vector *> dofs; + dofs.push_back(&dof_handler_u); + dofs.push_back(&dof_handler_p); + + std::shared_ptr> matrix_free_homogeneous( + new MatrixFree()); + matrix_free_homogeneous->reinit(MappingQ1{}, + dofs, + constraints_no_dirchlet, + QGauss<1>(degree_u + 1), + data); + + operator_homogeneous.initialize(matrix_free_homogeneous); + operator_homogeneous.evaluate_2_x_viscosity(Viscosity(sinker)); + LinearAlgebra::distributed::BlockVector inhomogeneity(2); + operator_homogeneous.initialize_dof_vector(inhomogeneity); + constraints_u.distribute(inhomogeneity.block(0)); + operator_homogeneous.vmult(system_rhs, inhomogeneity); + system_rhs *= -1.; + + // Normal apply boundary + RightHandSide right_hand_side(sinker); + + FEEvaluation + velocity(*stokes_matrix.get_matrix_free(), 0); + FEEvaluation + pressure(*stokes_matrix.get_matrix_free(), 1); + + for (unsigned int cell = 0; + cell < stokes_matrix.get_matrix_free()->n_cell_batches(); + ++cell) + { + velocity.reinit(cell); + pressure.reinit(cell); + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + Tensor<1, dim, VectorizedArray> rhs_u; + for (unsigned int d = 0; d < dim; ++d) + rhs_u[d] = make_vectorized_array(1.0); + VectorizedArray rhs_p = make_vectorized_array(1.0); + for (unsigned int i = 0; i < VectorizedArray::size(); ++i) + { + Point p; + for (unsigned int d = 0; d < dim; ++d) + p(d) = velocity.quadrature_point(q)(d)[i]; + + Vector rhs_temp(dim + 1); + right_hand_side.vector_value(p, rhs_temp); + + for (unsigned int d = 0; d < dim; ++d) + rhs_u[d][i] = rhs_temp(d); + rhs_p[i] = rhs_temp(dim); + } + velocity.submit_value(rhs_u, q); + pressure.submit_value(rhs_p, q); + } + velocity.integrate(EvaluationFlags::values); + velocity.distribute_local_to_global(system_rhs.block(0)); + pressure.integrate(EvaluationFlags::values); + pressure.distribute_local_to_global(system_rhs.block(1)); + } + system_rhs.compress(VectorOperation::add); + } + + template + void + StokesProblem::solve() + { + const double solver_tolerance = 1e-6 * system_rhs.l2_norm(); + const unsigned int n_cheap_stokes_solver_steps = 1000; + const double linear_solver_A_block_tolerance = 1e-2; + const double linear_solver_S_block_tolerance = 1e-6; + + // extract Stokes parts of solution vector, without any ghost elements + block_vector_t distributed_stokes_solution(solution); + + const unsigned int block_vel = 0; + const unsigned int block_p = 1; + + // extract Stokes parts of rhs vector + block_vector_t distributed_stokes_rhs(system_rhs); + + PrimitiveVectorMemory mem; + + SolverControl solver_control_cheap(n_cheap_stokes_solver_steps, + solver_tolerance, + false, + false); + + using Transfer = MGTransferMatrixFree; + + Transfer mg_transfer(mg_constrained_dofs); + mg_transfer.initialize_constraints(mg_constrained_dofs); + mg_transfer.build(dof_handler_u); + + LevelMatrixType & coarse_matrix = mg_matrices[0]; + SolverControl coarse_solver_control(1000, 1e-12, false, false); + SolverCG coarse_solver(coarse_solver_control); + + PreconditionIdentity coarse_prec_identity; + MGCoarseGridIterativeSolver, + LevelMatrixType, + PreconditionIdentity> + mg_coarse; + mg_coarse.initialize(coarse_solver, coarse_matrix, coarse_prec_identity); + + using SmootherType = PreconditionChebyshev; + mg::SmootherRelaxation mg_smoother; + MGLevelObject smoother_data; + smoother_data.resize(0, triangulation.n_global_levels() - 1); + for (unsigned int level = 0; level < triangulation.n_global_levels(); + ++level) + { + if (level > 0) + { + smoother_data[level].smoothing_range = 15.; + smoother_data[level].degree = 4; + smoother_data[level].eig_cg_n_iterations = 10; + } + else + { + smoother_data[0].smoothing_range = 1e-3; + smoother_data[0].degree = numbers::invalid_unsigned_int; + smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m(); + } + smoother_data[level].preconditioner = + mg_matrices[level].get_matrix_diagonal_inverse(); + } + mg_smoother.initialize(mg_matrices, smoother_data); + + mg::Matrix mg_matrix(mg_matrices); + + MGLevelObject> + mg_interface_matrices; + mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1); + for (unsigned int level = 0; level < triangulation.n_global_levels(); + ++level) + mg_interface_matrices[level].initialize(mg_matrices[level]); + mg::Matrix mg_interface(mg_interface_matrices); + + Multigrid mg( + mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother, 0); + mg.set_edge_matrices(mg_interface, mg_interface); + + + PreconditionMG prec_A(dof_handler_u, + mg, + mg_transfer); + + using MassPrec = PreconditionChebyshev; + MassPrec prec_S; + typename MassPrec::AdditionalData prec_S_data; + prec_S_data.smoothing_range = 1e-3; + prec_S_data.degree = numbers::invalid_unsigned_int; + prec_S_data.eig_cg_n_iterations = mass_matrix.m(); + prec_S_data.preconditioner = mass_matrix.get_matrix_diagonal_inverse(); + prec_S.initialize(mass_matrix, prec_S_data); + + using A_prec_type = PreconditionMG; + + // create a cheap preconditioner that consists of only a single V-cycle + const StokesSolver::BlockSchurPreconditioner + preconditioner_cheap(stokes_matrix, + mass_matrix, + prec_S, + prec_A, + false, + linear_solver_A_block_tolerance, + linear_solver_S_block_tolerance); + try + { + SolverFGMRES solver( + solver_control_cheap, + mem, + SolverFGMRES::AdditionalData(50)); + + solver.solve(stokes_matrix, + distributed_stokes_solution, + distributed_stokes_rhs, + preconditioner_cheap); + } + catch (const SolverControl::NoConvergence &) + { + deallog + << "********************************************************************" + << std::endl + << "SOLVER DID NOT CONVERGE AFTER " << n_cheap_stokes_solver_steps + << " ITERATIONS. res=" << solver_control_cheap.last_value() + << std::endl + << "********************************************************************" + << std::endl; + } + + constraints_u.distribute(distributed_stokes_solution.block(0)); + + distributed_stokes_solution.block(block_p) *= pressure_scaling; + + solution.block(block_vel) = distributed_stokes_solution.block(block_vel); + solution.block(block_p) = distributed_stokes_solution.block(block_p); + + deallog << "Solved-in " + << (solver_control_cheap.last_step() != + numbers::invalid_unsigned_int ? + solver_control_cheap.last_step() : + 0) + << " iterations, final residual: " + << solver_control_cheap.last_value() << std::endl; + } + + + template + void + StokesProblem::run() + { + deallog << "Sinker problem in " << dim << "D." << std::endl; + + create_sinker(4, 1000); + deallog << "n_sinker: " << sinker.n_sinkers + << " max/min viscosity ratio: " << sinker.DR_mu << std::endl + << std::endl; + + unsigned int initial_ref; + if (dim == 2) + { + initial_ref = 5; + } + else if (dim == 3) + { + initial_ref = 3; + } + + unsigned int n_cycles = 1; + if (dim == 2) + n_cycles = 2; + for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) + { + if (cycle == 0) + make_grid(initial_ref); + else + triangulation.refine_global(); + + setup_system(); + deallog << "Number of active cells: " + << triangulation.n_global_active_cells() << " (on " + << triangulation.n_global_levels() << " levels)" << std::endl; + deallog << "Number of degrees of freedom: " + << dof_handler_u.n_dofs() + dof_handler_p.n_dofs() << " (" + << dof_handler_u.n_dofs() << '+' << dof_handler_p.n_dofs() + << ')' << std::endl; + + assemble_rhs(); + solve(); + + deallog << std::endl; + } + } +} // namespace StokesClass + + +int +main(int argc, char *argv[]) +{ + dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + mpi_initlog(); + try + { + { + deallog.push("2d"); + StokesClass::StokesProblem<2> problem; + problem.run(); + deallog.pop(); + } + { + deallog.push("3d"); + StokesClass::StokesProblem<3> problem; + problem.run(); + deallog.pop(); + } + } + 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; + throw; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + throw; + } + + return 0; +} diff --git a/tests/numerics/no_flux_18.with_lapack=true.with_p4est=true.mpirun=4.output b/tests/numerics/no_flux_18.with_lapack=true.with_p4est=true.mpirun=4.output new file mode 100644 index 0000000000..dfc261478d --- /dev/null +++ b/tests/numerics/no_flux_18.with_lapack=true.with_p4est=true.mpirun=4.output @@ -0,0 +1,19 @@ + +DEAL:2d::Sinker problem in 2D. +DEAL:2d::n_sinker: 4 max/min viscosity ratio: 1000.00 +DEAL:2d:: +DEAL:2d::Number of active cells: 1024 (on 6 levels) +DEAL:2d::Number of degrees of freedom: 9539 (8450+1089) +DEAL:2d::Solved-in 48 iterations, final residual: 4.42329e-08 +DEAL:2d:: +DEAL:2d::Number of active cells: 4096 (on 7 levels) +DEAL:2d::Number of degrees of freedom: 37507 (33282+4225) +DEAL:2d::Solved-in 48 iterations, final residual: 2.31457e-08 +DEAL:2d:: +DEAL:3d::Sinker problem in 3D. +DEAL:3d::n_sinker: 4 max/min viscosity ratio: 1000.00 +DEAL:3d:: +DEAL:3d::Number of active cells: 512 (on 4 levels) +DEAL:3d::Number of degrees of freedom: 15468 (14739+729) +DEAL:3d::Solved-in 29 iterations, final residual: 1.16604e-08 +DEAL:3d:: -- 2.39.5