From 644fa451bab71ccd4b4ccf1cf0e3805a2a8212a9 Mon Sep 17 00:00:00 2001 From: tcclevenger Date: Mon, 9 Apr 2018 14:21:36 -0400 Subject: [PATCH] add test stokes computation --- tests/matrix_free/stokes_computation.cc | 1371 +++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=4.output | 19 + 2 files changed, 1390 insertions(+) create mode 100644 tests/matrix_free/stokes_computation.cc create mode 100644 tests/matrix_free/stokes_computation.with_mpi=true.with_p4est=true.mpirun=4.output diff --git a/tests/matrix_free/stokes_computation.cc b/tests/matrix_free/stokes_computation.cc new file mode 100644 index 0000000000..75054db18d --- /dev/null +++ b/tests/matrix_free/stokes_computation.cc @@ -0,0 +1,1371 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2012 - 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test for correctness of matrix free implementation for multigrid stokes + + +#include "../tests.h" + +#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 + + +unsigned int minlevel = 0; +const unsigned int velocity_degree = 2; + +double pressure_scaling = 1.0; + +namespace StokesClass +{ + using namespace dealii; + 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; + } + } + } + + +// 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 + 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 + 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 + 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 + 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 + 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_macro_cells(); + FEEvaluation velocity (*this->data, 0); + viscosity_x_2.reinit (n_cells, velocity.n_q_points); + for (unsigned int cell=0; cell return_value = make_vectorized_array(1.); + for (unsigned int i=0; i::n_array_elements; ++i) + { + Point p; + for (unsigned int d=0; d + void + StokesOperator + ::local_apply (const dealii::MatrixFree &data, + LinearAlgebra::distributed::BlockVector &dst, + const LinearAlgebra::distributed::BlockVector &src, + const std::pair &cell_range) const + { + typedef VectorizedArray vector_t; + FEEvaluation velocity (data, 0); + FEEvaluation pressure (data, 1); + + for (unsigned int cell=cell_range.first; cell 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 + 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_macro_cells(); + FEEvaluation pressure (*this->data, 0); + one_over_viscosity.reinit (n_cells, pressure.n_q_points); + for (unsigned int cell=0; cell return_value = make_vectorized_array(1.); + for (unsigned int i=0; i::n_array_elements; ++i) + { + Point p; + for (unsigned int d=0; d + 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 + 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 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 > diagonal(pressure.dofs_per_cell); + for (unsigned int i=0; i(); + pressure.begin_dof_values()[i] = make_vectorized_array (1.); + + pressure.evaluate (true,false,false); + for (unsigned int q=0; q + 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_macro_cells(); + FEEvaluation velocity (*this->data, 0); + viscosity_x_2.reinit (n_cells, velocity.n_q_points); + for (unsigned int cell=0; cell return_value = make_vectorized_array(1.); + for (unsigned int i=0; i::n_array_elements; ++i) + { + Point p; + for (unsigned int d=0; d + 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 + 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 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 > diagonal(velocity.dofs_per_cell); + for (unsigned int i=0; i(); + velocity.begin_dof_values()[i] = make_vectorized_array (1.); + + velocity.evaluate (false,true,false); + for (unsigned int q=0; q + 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 (); + + typedef LinearAlgebra::distributed::Vector vector_t; + typedef LinearAlgebra::distributed::BlockVector block_vector_t; + + typedef StokesOperator StokesMatrixType; + typedef MassMatrixOperator MassMatrixType; + typedef ABlockOperator LevelMatrixType; + + 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; + + ConstraintMatrix constraints_u; + ConstraintMatrix 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::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 coords(dim); + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + for (unsigned int i=0; i coords_tens; + for (unsigned int i=0; i(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); + 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(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(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(); + std::set dirichlet_boundary; + dirichlet_boundary.insert(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::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.level_mg_handler = level; + std::shared_ptr > + mg_mf_storage_level(new MatrixFree()); + mg_mf_storage_level->reinit(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 + ConstraintMatrix 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(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; celln_macro_cells(); ++cell) + { + velocity.reinit (cell); + pressure.reinit (cell); + for (unsigned int q=0; q > rhs_u; + for (unsigned int d=0; d (1.0); + VectorizedArray rhs_p = make_vectorized_array (1.0); + for (unsigned int i=0; i::n_array_elements; ++i) + { + Point p; + for (unsigned int d=0; d rhs_temp(dim+1); + right_hand_side.vector_value(p,rhs_temp); + + for (unsigned int d=0; d + 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< block_vector_t > mem; + + SolverControl solver_control_cheap (n_cheap_stokes_solver_steps, + solver_tolerance, false, false); + + typedef MGTransferMatrixFree Transfer; + + 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); + + typedef PreconditionChebyshev SmootherType; + mg::SmootherRelaxation mg_smoother; + MGLevelObject smoother_data; + smoother_data.resize(0, triangulation.n_global_levels()-1); + for (unsigned int level = 0; 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 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); + + typedef PreconditionChebyshev MassPrec; + 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); + + typedef PreconditionMG A_prec_type; + + // 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, true)); + + solver.solve (stokes_matrix, + distributed_stokes_solution, + distributed_stokes_rhs, + preconditioner_cheap); + } + catch (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 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/matrix_free/stokes_computation.with_mpi=true.with_p4est=true.mpirun=4.output b/tests/matrix_free/stokes_computation.with_mpi=true.with_p4est=true.mpirun=4.output new file mode 100644 index 0000000000..c33433cf25 --- /dev/null +++ b/tests/matrix_free/stokes_computation.with_mpi=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 62 iterations, final residual: 5.37736e-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 61 iterations, final residual: 2.94930e-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 54 iterations, final residual: 2.12652e-08 +DEAL:3d:: -- 2.39.5