Dear Michal, This is expected: the matrix-free operator evaluation cannot apply non-homogeneous boundary conditions while solving (at least I have never figured out how to do that). To solve such a problem, you need to bring the non-homogeneous part on the right hand side first. I often solve this by first creating a vector that is filled with the Dirichlet values and apply the operator on that (without setting the zero constraints) and then solve for an increment that has zero boundary conditions. I attach an example which is an extension of step-37 (and uses coefficients described here, if you want some more background: http://www.it.uu.se/research/publications/reports/2017-006/2017-006-nc.pdf
In the program, you will find that LaplaceOperator not only contains a vmult() function, but also a compute_residual() function that reads the data in from an FEEvaluation object without Dirichlet conditions, due the operation, and assembles the right hand side into an FEEvaluation object with Dirichlet conditions (-> local_compute_residual()). Please let me know in case you have questions. Best, Martin On 26.10.2017 22:07, Michał Wichrowski wrote: > Dear deal.II developers, > I think I found problem with non-homogeneus boundary conditions. > By changing: > VectorTools::interpolate_boundary_values (dof_handler, > 0, > Functions::ZeroFunction<dim>(), > constraints); > to > VectorTools::interpolate_boundary_values (dof_handler, > 0, > Functions::ConstantFunction<dim>(5), > constraints); > > > The results shown on included picture were outputted. I've tried to > figure out what is going on, and it looks like matrix-free framework > is solving problem with homogeneous constrains and then Dirichlet > nodes are overwritten. > -- > The deal.II project is located at http://www.dealii.org/ > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > --- > You received this message because you are subscribed to the Google > Groups "deal.II User Group" group. > To unsubscribe from this group and stop receiving emails from it, send > an email to dealii+unsubscr...@googlegroups.com > <mailto:dealii+unsubscr...@googlegroups.com>. > For more options, visit https://groups.google.com/d/optout. -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
/* --------------------------------------------------------------------- * * Copyright (C) 2009 - 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. * * --------------------------------------------------------------------- * * Authors: Katharina Kormann, Martin Kronbichler, Uppsala University, * 2009-2012, updated to MPI version with parallel vectors in 2016 */ #include <deal.II/base/quadrature_lib.h> #include <deal.II/base/function.h> #include <deal.II/base/function_lib.h> #include <deal.II/base/logstream.h> #include <deal.II/base/timer.h> #include <deal.II/base/convergence_table.h> #include <deal.II/lac/full_matrix.h> #include <deal.II/lac/solver_cg.h> #include <deal.II/lac/la_parallel_vector.h> #include <deal.II/lac/precondition.h> #include <deal.II/fe/fe_q.h> #include <deal.II/fe/fe_values.h> #include <deal.II/grid/tria.h> #include <deal.II/grid/tria_accessor.h> #include <deal.II/grid/tria_iterator.h> #include <deal.II/grid/grid_generator.h> #include <deal.II/grid/manifold_lib.h> #include <deal.II/multigrid/multigrid.h> #include <deal.II/multigrid/mg_transfer_matrix_free.h> #include <deal.II/multigrid/mg_tools.h> #include <deal.II/multigrid/mg_coarse.h> #include <deal.II/multigrid/mg_smoother.h> #include <deal.II/multigrid/mg_matrix.h> #include <deal.II/numerics/data_out.h> #include <deal.II/numerics/vector_tools.h> #include <deal.II/numerics/error_estimator.h> #include <deal.II/matrix_free/matrix_free.h> #include <deal.II/matrix_free/operators.h> #include <deal.II/matrix_free/fe_evaluation.h> #include <deal.II/distributed/solution_transfer.h> #include <deal.II/distributed/grid_refinement.h> #include <deal.II/distributed/tria.h> #include <iostream> #include <fstream> #include <sstream> namespace Step37 { using namespace dealii; const unsigned int degree_finite_element = 3; const unsigned int dimension = 3; typedef double number; typedef float level_number; // ========================================================================== // Functions // ========================================================================== template <int dim> class Solution : public Function<dim> { public: Solution () : Function<dim>() {} virtual double value (const Point<dim> &p, const unsigned int component = 0) const; virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component = 0) const; virtual Tensor<1,dim,VectorizedArray<double>> gradient (const Point<dim,VectorizedArray<double>> &p, const unsigned int component = 0) const; virtual double laplacian (const Point<dim> &p, const unsigned int component = 0) const; virtual VectorizedArray<double> laplacian (const Point<dim,VectorizedArray<double>> &p, const unsigned int component = 0) const; }; template <int dim> double Solution<dim>::value (const Point<dim> &p, const unsigned int) const { return std::sin(2*numbers::PI*(p[0] + p[1])); } template <int dim> Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &p, const unsigned int) const { Tensor<1,dim> grad_tensor; grad_tensor[0] = 1; grad_tensor[1] = 1; return 2*numbers::PI*std::cos(2*numbers::PI*(p[0] + p[1]))*grad_tensor; } template <int dim> Tensor<1,dim,VectorizedArray<double>> Solution<dim>::gradient (const Point<dim,VectorizedArray<double>> &p, const unsigned int ) const { Tensor<1,dim,VectorizedArray<double>> direction; for(int i = 0; i < 2; ++i) direction[i] = make_vectorized_array(1.0); return 2*numbers::PI*std::cos(2*numbers::PI*(p[0] + p[1]))*direction; } template <int dim> double Solution<dim>::laplacian (const Point<dim> &p, const unsigned int) const { return -2*2*numbers::PI*2*numbers::PI*std::sin(2*numbers::PI*(p[0] + p[1])); } template <int dim> VectorizedArray<double> Solution<dim>::laplacian (const Point<dim,VectorizedArray<double>> &p, const unsigned int) const { return -2*2*numbers::PI*2*numbers::PI*std::sin(2*numbers::PI*(p[0] + p[1])); } // function computing the coefficient template <int dim, typename Number = double> class Coefficient : Function<dim,Number> { public: Coefficient () : Function<dim,Number>() {} virtual Number value (const Point<dim,Number> &p, const unsigned int component = 0) const; virtual VectorizedArray<Number> value (const Point<dim,VectorizedArray<Number>> &p, const unsigned int component = 0) const; virtual Tensor<1,dim,Number> gradient (const Point<dim,Number> &p, const unsigned int component = 0) const; virtual Tensor<1,dim,VectorizedArray<Number>> gradient (const Point<dim,VectorizedArray<Number>> &p, const unsigned int component = 0) const; }; template <int dim, typename Number> Number Coefficient<dim,Number>::value (const Point<dim,Number> &p, const unsigned int) const { Number prod = 1.0; for(int e = 0; e < dim; ++e) { Number c = std::cos(2*numbers::PI*p[e]+0.1*e); prod *= c*c; } return 1.0 + 1.0e6*prod; } template <int dim, typename Number> VectorizedArray<Number> Coefficient<dim,Number>::value (const Point<dim,VectorizedArray<Number>> &p, const unsigned int) const { VectorizedArray<Number> prod = make_vectorized_array<Number>(1.0); for(int e = 0; e < dim; ++e) { VectorizedArray<Number> c = std::cos(2*numbers::PI*p[e]+0.1*e); prod *= c*c; } return 1.0 + 1.0e6*prod; } template <int dim, typename Number> Tensor<1,dim,Number> Coefficient<dim,Number>::gradient (const Point<dim,Number> &p, const unsigned int) const { Tensor<1,dim,Number> prod; for(int d = 0; d < dim; ++d) { prod[d] = 1.0; for(int e = 0; e < dim; ++e) { Number c = std::cos(2*numbers::PI*p[e]+0.1*e); if(e==d) { prod[d] *= -4*numbers::PI* c * std::sin(2*numbers::PI*p[e]+0.1*e); } else { prod[d] *= c*c; } } } return 1.0e6*prod; } template <int dim, typename Number> Tensor<1,dim,VectorizedArray<Number>> Coefficient<dim,Number>::gradient (const Point<dim,VectorizedArray<Number>> &p, const unsigned int ) const { Tensor<1,dim,VectorizedArray<Number>> prod; for(int d = 0; d < dim; ++d) { prod[d] = make_vectorized_array<Number>(1.0); for(int e = 0; e < dim; ++e) { VectorizedArray<Number> c = std::cos(2*numbers::PI*p[e]+0.1*e); if(e==d) { prod[d] *= -4*numbers::PI* c * std::sin(2*numbers::PI*p[e]+0.1*e); } else { prod[d] *= c*c; } } } return 1.0e6*prod; } // function computing the right-hand side template <int dim> class RightHandSide : public Function<dim> { private: Solution<dim> solution; Coefficient<dim> coefficient; public: RightHandSide () : Function<dim>() {} virtual double value (const Point<dim> &p, const unsigned int component = 0) const; virtual VectorizedArray<double> value (const Point<dim,VectorizedArray<double>> &p, const unsigned int component = 0) const; }; template <int dim> double RightHandSide<dim>::value (const Point<dim> &p, const unsigned int) const { return -(solution.laplacian(p)*coefficient.value(p) + coefficient.gradient(p)*solution.gradient(p)); } template <int dim> VectorizedArray<double> RightHandSide<dim>::value (const Point<dim,VectorizedArray<double>> &p, const unsigned int ) const { return -(solution.laplacian(p)*coefficient.value(p) + coefficient.gradient(p)*solution.gradient(p)); } // ========================================================================== // Operator // ========================================================================== template <int dim, int fe_degree, typename number> class LaplaceOperator : public MatrixFreeOperators::Base<dim,LinearAlgebra::distributed::Vector<number> > { public: typedef number value_type; LaplaceOperator (); void compute_residual (LinearAlgebra::distributed::Vector<number> &dst, const LinearAlgebra::distributed::Vector<number> &src) const; void evaluate_coefficient (); void clear(); virtual void compute_diagonal(); private: virtual void apply_add(LinearAlgebra::distributed::Vector<number> &dst, const LinearAlgebra::distributed::Vector<number> &src) const; void local_apply (const MatrixFree<dim,number> &data, LinearAlgebra::distributed::Vector<number> &dst, const LinearAlgebra::distributed::Vector<number> &src, const std::pair<unsigned int,unsigned int> &cell_range) const; void local_compute_residual (const MatrixFree<dim,number> &data, LinearAlgebra::distributed::Vector<number> &dst, const LinearAlgebra::distributed::Vector<number> &src, const std::pair<unsigned int,unsigned int> &cell_range) const; void local_compute_diagonal (const MatrixFree<dim,number> &data, LinearAlgebra::distributed::Vector<number> &dst, const unsigned int &dummy, const std::pair<unsigned int,unsigned int> &cell_range) const; void local_compute_coefficient (const MatrixFree<dim,number> &data, unsigned int &dummy, const unsigned int &dummy2, const std::pair<unsigned int,unsigned int> &cell_range); Table<2, VectorizedArray<number> > coefficient; }; template <int dim, int fe_degree, typename number> LaplaceOperator<dim,fe_degree,number>::LaplaceOperator () : MatrixFreeOperators::Base<dim, LinearAlgebra::distributed::Vector<number> >() {} template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>::clear () { MatrixFreeOperators::Base<dim,LinearAlgebra::distributed::Vector<number> >::clear(); coefficient.reinit(0, 0); } template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>:: evaluate_coefficient () { coefficient.reinit(this->data->n_macro_cells(), this->data->get_n_q_points()); unsigned int dummy; this->data->cell_loop (&LaplaceOperator::local_compute_coefficient, this, dummy, dummy); } template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>:: local_compute_coefficient (const MatrixFree<dim,number> &, unsigned int &, const unsigned int &, const std::pair<unsigned int,unsigned int> &cell_range) { Coefficient<dim,number> coefficient_function; FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (*this->data, 1); for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell) { phi.reinit (cell); for (unsigned int q=0; q<phi.n_q_points; ++q) { coefficient(cell,q) = coefficient_function.value(phi.quadrature_point(q)); } } } template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number> ::local_apply (const MatrixFree<dim,number> &data, LinearAlgebra::distributed::Vector<number> &dst, const LinearAlgebra::distributed::Vector<number> &src, const std::pair<unsigned int,unsigned int> &cell_range) const { FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (data); for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell) { phi.reinit (cell); phi.read_dof_values(src); phi.evaluate (false, true); for (unsigned int q=0; q<phi.n_q_points; ++q) phi.submit_gradient (coefficient(cell,q) * phi.get_gradient(q), q); phi.integrate (false, true); phi.distribute_local_to_global (dst); } } template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number> ::apply_add (LinearAlgebra::distributed::Vector<number> &dst, const LinearAlgebra::distributed::Vector<number> &src) const { this->data->cell_loop (&LaplaceOperator::local_apply, this, dst, src); } template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>:: compute_residual (LinearAlgebra::distributed::Vector<number> &dst, const LinearAlgebra::distributed::Vector<number> &src) const { dst = 0; this->data->cell_loop (&LaplaceOperator::local_compute_residual, this, dst, src); } template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number>:: local_compute_residual (const MatrixFree<dim,number> &, LinearAlgebra::distributed::Vector<number> &dst, const LinearAlgebra::distributed::Vector<number> &solution, const std::pair<unsigned int,unsigned int> &cell_range) const { FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (*this->data, 0); FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi_nodirichlet (*this->data, 1); for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell) { phi.reinit (cell); phi_nodirichlet.reinit (cell); phi_nodirichlet.read_dof_values(solution); phi_nodirichlet.evaluate(false, true); for (unsigned int q=0; q<phi.n_q_points; ++q) phi.submit_gradient (-coefficient(cell,q) * phi_nodirichlet.get_gradient(q), q); phi.integrate(false, true); phi.distribute_local_to_global(dst); } } template <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number> ::compute_diagonal () { this->inverse_diagonal_entries. reset(new DiagonalMatrix<LinearAlgebra::distributed::Vector<number> >()); LinearAlgebra::distributed::Vector<number> &inverse_diagonal = this->inverse_diagonal_entries->get_vector(); this->data->initialize_dof_vector(inverse_diagonal); unsigned int dummy = 0; this->data->cell_loop (&LaplaceOperator::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 <int dim, int fe_degree, typename number> void LaplaceOperator<dim,fe_degree,number> ::local_compute_diagonal (const MatrixFree<dim,number> &data, LinearAlgebra::distributed::Vector<number> &dst, const unsigned int &, const std::pair<unsigned int,unsigned int> &cell_range) const { FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (data); AlignedVector<VectorizedArray<number> > diagonal(phi.dofs_per_cell); for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell) { phi.reinit (cell); for (unsigned int i=0; i<phi.dofs_per_cell; ++i) { for (unsigned int j=0; j<phi.dofs_per_cell; ++j) phi.submit_dof_value(VectorizedArray<number>(), j); phi.submit_dof_value(make_vectorized_array<number>(1.), i); phi.evaluate (false, true); for (unsigned int q=0; q<phi.n_q_points; ++q) phi.submit_gradient (coefficient(cell,q) * phi.get_gradient(q), q); phi.integrate (false, true); diagonal[i] = phi.get_dof_value(i); } for (unsigned int i=0; i<phi.dofs_per_cell; ++i) phi.submit_dof_value(diagonal[i], i); phi.distribute_local_to_global (dst); } } template <int dim> class LaplaceProblem { public: LaplaceProblem (); void run (); private: void setup_system (); void assemble_rhs(); std::pair<unsigned int,std::pair<double,double> > solve (); void output_results (const unsigned int cycle) const; void interpolate_boundary_values(); #ifdef DEAL_II_WITH_P4EST parallel::distributed::Triangulation<dim> triangulation; #else Triangulation<dim> triangulation; #endif MappingQGeneric<dim> mapping; FE_Q<dim> fe; DoFHandler<dim> dof_handler; ConstraintMatrix constraints_without_dirichlet; ConstraintMatrix constraints; std::shared_ptr<MatrixFree<dim,number> > system_matrix_free; typedef LaplaceOperator<dim,degree_finite_element,number> SystemMatrixType; SystemMatrixType system_matrix; MGConstrainedDoFs mg_constrained_dofs; MGTransferMatrixFree<dim,level_number> mg_transfer; typedef LaplaceOperator<dim,degree_finite_element,level_number> LevelMatrixType; MGLevelObject<LevelMatrixType> mg_matrices; LinearAlgebra::distributed::Vector<number> solution; LinearAlgebra::distributed::Vector<number> solution_update; LinearAlgebra::distributed::Vector<number> system_rhs; double setup_time; ConditionalOStream pcout; ConditionalOStream time_details; }; template <int dim> LaplaceProblem<dim>::LaplaceProblem () : #ifdef DEAL_II_WITH_P4EST triangulation(MPI_COMM_WORLD, Triangulation<dim>::limit_level_difference_at_vertices, parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy), #else triangulation (Triangulation<dim>::limit_level_difference_at_vertices), #endif mapping (degree_finite_element), fe (degree_finite_element), dof_handler (triangulation), pcout (std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0), time_details (std::cout, false && Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {} template <int dim> void LaplaceProblem<dim>::setup_system () { Timer time; time.start (); setup_time = 0; system_matrix_free.reset(new MatrixFree<dim,number>()); system_matrix.clear(); mg_matrices.clear_elements(); dof_handler.distribute_dofs (fe); dof_handler.distribute_mg_dofs (); pcout << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; IndexSet locally_relevant_dofs; DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs); constraints.clear(); constraints.reinit(locally_relevant_dofs); DoFTools::make_hanging_node_constraints(dof_handler, constraints); VectorTools::interpolate_boundary_values (dof_handler, 0, ZeroFunction<dim>(), constraints); constraints.close(); constraints_without_dirichlet.clear(); constraints_without_dirichlet.reinit(locally_relevant_dofs); DoFTools::make_hanging_node_constraints(dof_handler, constraints_without_dirichlet); constraints_without_dirichlet.close(); setup_time += time.wall_time(); time_details << "Distribute DoFs & B.C. (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << "s" << std::endl; time.restart(); std::vector<const DoFHandler<dim> *> dof_handlers(2, &dof_handler); { std::vector<const ConstraintMatrix *> constraint(2); constraint[0] = &constraints; constraint[1] = &constraints_without_dirichlet; typename MatrixFree<dim,number>::AdditionalData additional_data; additional_data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none; additional_data.mapping_update_flags = (update_gradients | update_JxW_values | update_quadrature_points); system_matrix_free->reinit (mapping, dof_handlers, constraint, QGauss<1>(fe.degree+1), additional_data); std::vector<unsigned int> selected_blocks(1, 0); system_matrix.initialize (system_matrix_free, selected_blocks); system_matrix.evaluate_coefficient(); } system_matrix_free->initialize_dof_vector(solution, 1); system_matrix_free->initialize_dof_vector(solution_update, 0); system_matrix_free->initialize_dof_vector(system_rhs, 0); setup_time += time.wall_time(); time_details << "Setup matrix-free system (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << "s" << std::endl; time.restart(); const unsigned int nlevels = triangulation.n_global_levels(); mg_matrices.resize(0, nlevels-1); std::set<types::boundary_id> dirichlet_boundary; dirichlet_boundary.insert(0); mg_constrained_dofs.initialize(dof_handler); mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, dirichlet_boundary); for (unsigned int level=0; level<nlevels; ++level) { IndexSet relevant_dofs; DoFTools::extract_locally_relevant_level_dofs(dof_handler, level, relevant_dofs); ConstraintMatrix level_constraints; level_constraints.reinit(relevant_dofs); level_constraints.add_lines(mg_constrained_dofs.get_boundary_indices(level)); level_constraints.close(); std::vector<const ConstraintMatrix *> constraint(2); ConstraintMatrix dummy; dummy.close(); constraint[0] = &level_constraints; constraint[1] = &dummy; typename MatrixFree<dim,level_number>::AdditionalData additional_data; additional_data.tasks_parallel_scheme = MatrixFree<dim,level_number>::AdditionalData::none; additional_data.mapping_update_flags = (update_gradients | update_JxW_values | update_quadrature_points); additional_data.level_mg_handler = level; std::shared_ptr<MatrixFree<dim,level_number> > mg_matrix_free_level(new MatrixFree<dim,level_number>()); mg_matrix_free_level->reinit(mapping, dof_handlers, constraint, QGauss<1>(fe.degree+1), additional_data); std::vector<unsigned int> selected_blocks(1, 0); mg_matrices[level].initialize(mg_matrix_free_level, mg_constrained_dofs, level, selected_blocks); mg_matrices[level].evaluate_coefficient(); } setup_time += time.wall_time(); time_details << "Setup matrix-free levels (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << "s" << std::endl; time.restart(); mg_transfer.clear(); mg_transfer.initialize_constraints(mg_constrained_dofs); mg_transfer.build(dof_handler); setup_time += time.wall_time(); time_details << "MG build transfer time (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << "s\n"; pcout << "Total setup time (wall) " << setup_time << "s\n"; } template <int dim> void LaplaceProblem<dim>::assemble_rhs() { Timer time; RightHandSide<dim> right_hand_side; system_rhs = 0; FEEvaluation<dim,degree_finite_element> phi(*system_matrix.get_matrix_free()); for (unsigned int cell=0; cell<system_matrix.get_matrix_free()->n_macro_cells(); ++cell) { phi.reinit(cell); for (unsigned int q=0; q<phi.n_q_points; ++q) phi.submit_value(right_hand_side.value(phi.quadrature_point(q)), q); phi.integrate(true, false); phi.distribute_local_to_global(system_rhs); } system_rhs.compress(VectorOperation::add); setup_time += time.wall_time(); time_details << "Assemble right hand side (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << "s" << std::endl; } template <int dim> void LaplaceProblem<dim>::interpolate_boundary_values() { std::map<types::global_dof_index, double> boundary_values; VectorTools::interpolate_boundary_values(mapping, dof_handler, 0, Solution<dim>(), boundary_values); for (typename std::map<types::global_dof_index, double>::iterator it = boundary_values.begin(); it != boundary_values.end(); ++it) if (solution.locally_owned_elements().is_element(it->first)) solution(it->first) = it->second; constraints_without_dirichlet.distribute(solution); } template <int dim> std::pair<unsigned int,std::pair<double,double> > LaplaceProblem<dim>::solve () { Timer time; typedef PreconditionChebyshev<LevelMatrixType,LinearAlgebra::distributed::Vector<level_number> > SmootherType; mg::SmootherRelaxation<SmootherType, LinearAlgebra::distributed::Vector<level_number> > mg_smoother; MGLevelObject<typename SmootherType::AdditionalData> 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 = 2; smoother_data[level].eig_cg_n_iterations = 15; } 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(); } mg_matrices[level].compute_diagonal(); smoother_data[level].preconditioner = mg_matrices[level].get_matrix_diagonal_inverse(); } mg_smoother.initialize(mg_matrices, smoother_data); MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<level_number> > mg_coarse; mg_coarse.initialize(mg_smoother); mg::Matrix<LinearAlgebra::distributed::Vector<level_number> > mg_matrix(mg_matrices); MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<LevelMatrixType> > 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<LinearAlgebra::distributed::Vector<level_number> > mg_interface(mg_interface_matrices); Multigrid<LinearAlgebra::distributed::Vector<level_number> > mg(mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); mg.set_edge_matrices(mg_interface, mg_interface); PreconditionMG<dim, LinearAlgebra::distributed::Vector<level_number>, MGTransferMatrixFree<dim,level_number> > preconditioner(dof_handler, mg, mg_transfer); ReductionControl solver_control (system_matrix.m(), 1e-13, 1e-9); SolverCG<LinearAlgebra::distributed::Vector<number> > cg (solver_control); pcout << "MG build smoother time (CPU/wall) " << time.cpu_time() << "s/" << time.wall_time() << "s\n"; LinearAlgebra::distributed::Vector<number> residual; system_matrix.initialize_dof_vector(residual); system_matrix.compute_residual(residual, solution); system_rhs += residual; time.reset(); time.start(); solution_update = 0; cg.solve (system_matrix, solution_update, system_rhs, preconditioner); std::pair<unsigned int,std::pair<double,double> > stats(solver_control.last_step(), std::pair<double,double>(time.wall_time(),0)); time.reset(); time.start(); solution_update = 0; cg.solve (system_matrix, solution_update, system_rhs, preconditioner); stats.second.second = time.wall_time(); pcout << "Time solve (" << stats.first << " iterations) (CPU/wall) " << time.cpu_time() << "s/" << stats.second.second << "s\n"; constraints.distribute(solution_update); solution += solution_update; return stats; } template <int dim> void LaplaceProblem<dim>::output_results (const unsigned int cycle) const { if (triangulation.n_global_active_cells() > 1000000) return; DataOut<dim> data_out; solution.update_ghost_values(); data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); data_out.build_patches (mapping, 1, DataOut<dim>::curved_inner_cells); std::ostringstream filename; filename << "solution-" << cycle << "." << Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) << ".vtu"; std::ofstream output (filename.str().c_str()); data_out.write_vtu (output); if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) { std::vector<std::string> filenames; for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i) { std::ostringstream filename; filename << "solution-" << cycle << "." << i << ".vtu"; filenames.push_back(filename.str().c_str()); } std::string master_name = "solution-" + Utilities::to_string(cycle) + ".pvtu"; std::ofstream master_output (master_name.c_str()); data_out.write_pvtu_record (master_output, filenames); } } template <int dim> void LaplaceProblem<dim>::run () { ConvergenceTable convergence_table; for (unsigned int cycle=0; cycle<8-dim; ++cycle) { pcout << "Cycle " << cycle << std::endl; if (cycle == 0) { GridGenerator::hyper_shell (triangulation,Point<dim>(),0.5,1.0); static const SphericalManifold<dim> manifold; triangulation.set_all_manifold_ids(0); triangulation.set_manifold (0, manifold); triangulation.refine_global (8-2*dim); setup_system(); } else { triangulation.refine_global(); setup_system (); } assemble_rhs(); interpolate_boundary_values(); auto stats = solve(); output_results(cycle); solution.update_ghost_values(); Vector<float> error_per_cell(triangulation.n_active_cells()); VectorTools::integrate_difference (mapping, dof_handler, solution, Solution<dim>(), error_per_cell, QGauss<dim>(fe.degree+2), VectorTools::L2_norm); const double L2_error = std::sqrt(Utilities::MPI::sum(error_per_cell.norm_sqr(), MPI_COMM_WORLD)); //VectorTools::integrate_difference (mapping, // dof_handler, // solution, // Solution<dim>(), // error_per_cell, // QGauss<dim>(fe.degree+1), // VectorTools::H1_seminorm); //const double grad_error = // std::sqrt(Utilities::MPI::sum(error_per_cell.norm_sqr(), MPI_COMM_WORLD)); convergence_table.add_value("cells", triangulation.n_global_active_cells()); convergence_table.add_value("dofs", dof_handler.n_dofs()); convergence_table.add_value("val L2", L2_error); //convergence_table.add_value("grad L2", grad_error); convergence_table.add_value("solver_its", stats.first); convergence_table.add_value("solver_time1", stats.second.first); convergence_table.add_value("solver_time2", stats.second.second); pcout << std::endl; }; convergence_table.set_precision("val L2", 3); convergence_table.set_scientific("val L2", true); //convergence_table.set_precision("grad L2", 3); //convergence_table.set_scientific("grad L2", true); convergence_table.set_precision("solver_time1", 3); convergence_table.set_scientific("solver_time1", true); convergence_table.set_precision("solver_time2", 3); convergence_table.set_scientific("solver_time2", true); if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) { convergence_table.write_text(std::cout); convergence_table.write_tex(std::cout); } } } int main (int argc, char *argv[]) { try { using namespace Step37; Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); LaplaceProblem<dimension> laplace_problem; laplace_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; }