From: Peter Munch Date: Mon, 11 Jan 2021 13:28:04 +0000 (+0100) Subject: Convert step-31 to a simplex test X-Git-Tag: v9.3.0-rc1~382^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b54da6e810f6e9a323e0b7d25e430f793784c820;p=dealii.git Convert step-31 to a simplex test --- diff --git a/tests/simplex/step-31.cc b/tests/simplex/step-31.cc new file mode 100644 index 0000000000..5fa28f8ad1 --- /dev/null +++ b/tests/simplex/step-31.cc @@ -0,0 +1,1128 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 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. +// +// --------------------------------------------------------------------- + + +// Step-31 on a simplex mesh. Following incompatible modifications had to be +// made: +// - number of maximum quadrature points had to be limited +// - no local refinement + +#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" + +namespace Step31 +{ + using namespace dealii; + namespace EquationData + { + constexpr double eta = 1; + constexpr double kappa = 1e-6; + constexpr double beta = 10; + constexpr double density = 1; + template + class TemperatureInitialValues : public Function + { + public: + TemperatureInitialValues() + : Function(1) + {} + virtual double + value(const Point & /*p*/, + const unsigned int /*component*/ = 0) const override + { + return 0; + } + virtual void + vector_value(const Point &p, Vector &value) const override + { + for (unsigned int c = 0; c < this->n_components; ++c) + value(c) = TemperatureInitialValues::value(p, c); + } + }; + template + class TemperatureRightHandSide : public Function + { + public: + TemperatureRightHandSide() + : Function(1) + {} + virtual double + value(const Point & p, + const unsigned int component = 0) const override + { + (void)component; + Assert(component == 0, + ExcMessage("Invalid operation for a scalar function.")); + Assert((dim == 2) || (dim == 3), ExcNotImplemented()); + static const Point source_centers[3] = { + (dim == 2 ? Point(.3, .1) : Point(.3, .5, .1)), + (dim == 2 ? Point(.45, .1) : Point(.45, .5, .1)), + (dim == 2 ? Point(.75, .1) : Point(.75, .5, .1))}; + static const double source_radius = (dim == 2 ? 1. / 32 : 1. / 8); + return ((source_centers[0].distance(p) < source_radius) || + (source_centers[1].distance(p) < source_radius) || + (source_centers[2].distance(p) < source_radius) ? + 1 : + 0); + } + virtual void + vector_value(const Point &p, Vector &value) const override + { + for (unsigned int c = 0; c < this->n_components; ++c) + value(c) = TemperatureRightHandSide::value(p, c); + } + }; + } // namespace EquationData + namespace LinearSolvers + { + template + class InverseMatrix : public Subscriptor + { + public: + InverseMatrix(const MatrixType & m, + const PreconditionerType &preconditioner); + template + void + vmult(VectorType &dst, const VectorType &src) const; + + private: + const SmartPointer matrix; + const PreconditionerType & preconditioner; + }; + template + InverseMatrix::InverseMatrix( + const MatrixType & m, + const PreconditionerType &preconditioner) + : matrix(&m) + , preconditioner(preconditioner) + {} + template + template + void + InverseMatrix::vmult( + VectorType & dst, + const VectorType &src) const + { + SolverControl solver_control(src.size(), 1e-7 * src.l2_norm()); + SolverCG cg(solver_control); + dst = 0; + try + { + cg.solve(*matrix, dst, src, preconditioner); + } + catch (std::exception &e) + { + Assert(false, ExcMessage(e.what())); + } + } + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner( + const TrilinosWrappers::BlockSparseMatrix &S, + const InverseMatrix &Mpinv, + const PreconditionerTypeA & Apreconditioner); + void + vmult(TrilinosWrappers::MPI::BlockVector & dst, + const TrilinosWrappers::MPI::BlockVector &src) const; + + private: + const SmartPointer + stokes_matrix; + const SmartPointer> + m_inverse; + const PreconditionerTypeA & a_preconditioner; + mutable TrilinosWrappers::MPI::Vector tmp; + }; + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner( + const TrilinosWrappers::BlockSparseMatrix &S, + const InverseMatrix &Mpinv, + const PreconditionerTypeA & Apreconditioner) + : stokes_matrix(&S) + , m_inverse(&Mpinv) + , a_preconditioner(Apreconditioner) + , tmp(complete_index_set(stokes_matrix->block(1, 1).m())) + {} + template + void + BlockSchurPreconditioner::vmult( + TrilinosWrappers::MPI::BlockVector & dst, + const TrilinosWrappers::MPI::BlockVector &src) const + { + a_preconditioner.vmult(dst.block(0), src.block(0)); + stokes_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1)); + tmp *= -1; + m_inverse->vmult(dst.block(1), tmp); + } + } // namespace LinearSolvers + template + class BoussinesqFlowProblem + { + public: + BoussinesqFlowProblem(); + void + run(); + + private: + void + setup_dofs(); + void + assemble_stokes_preconditioner(); + void + build_stokes_preconditioner(); + void + assemble_stokes_system(); + void + assemble_temperature_system(const double maximal_velocity); + void + assemble_temperature_matrix(); + double + get_maximal_velocity() const; + std::pair + get_extrapolated_temperature_range() const; + void + solve(); + void + output_results() const; + void + refine_mesh(const unsigned int max_grid_level); + double + compute_viscosity( + const std::vector & old_temperature, + const std::vector & old_old_temperature, + const std::vector> &old_temperature_grads, + const std::vector> &old_old_temperature_grads, + const std::vector & old_temperature_laplacians, + const std::vector & old_old_temperature_laplacians, + const std::vector> &old_velocity_values, + const std::vector> &old_old_velocity_values, + const std::vector & gamma_values, + const double global_u_infty, + const double global_T_variation, + const double cell_diameter) const; + Triangulation triangulation; + double global_Omega_diameter; + MappingFE mapping; + const unsigned int stokes_degree; + FESystem stokes_fe; + DoFHandler stokes_dof_handler; + AffineConstraints stokes_constraints; + std::vector stokes_partitioning; + TrilinosWrappers::BlockSparseMatrix stokes_matrix; + TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; + TrilinosWrappers::MPI::BlockVector stokes_solution; + TrilinosWrappers::MPI::BlockVector old_stokes_solution; + TrilinosWrappers::MPI::BlockVector stokes_rhs; + const unsigned int temperature_degree; + FE_SimplexP temperature_fe; + DoFHandler temperature_dof_handler; + + QGaussSimplex quadrature_temperature; + QGaussSimplex quadrature_stokes; + + AffineConstraints temperature_constraints; + TrilinosWrappers::SparseMatrix temperature_mass_matrix; + TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; + TrilinosWrappers::SparseMatrix temperature_matrix; + TrilinosWrappers::MPI::Vector temperature_solution; + TrilinosWrappers::MPI::Vector old_temperature_solution; + TrilinosWrappers::MPI::Vector old_old_temperature_solution; + TrilinosWrappers::MPI::Vector temperature_rhs; + double time_step; + double old_time_step; + unsigned int timestep_number; + std::shared_ptr Amg_preconditioner; + std::shared_ptr Mp_preconditioner; + bool rebuild_stokes_matrix; + bool rebuild_temperature_matrices; + bool rebuild_stokes_preconditioner; + }; + template + BoussinesqFlowProblem::BoussinesqFlowProblem() + : triangulation(/*Triangulation::maximum_smoothing*/) + , global_Omega_diameter(std::numeric_limits::quiet_NaN()) + , mapping(FE_SimplexP(1)) + , stokes_degree(1) + , stokes_fe(FE_SimplexP(stokes_degree + 1), + dim, + FE_SimplexP(stokes_degree), + 1) + , stokes_dof_handler(triangulation) + , temperature_degree(2) + , temperature_fe(temperature_degree) + , temperature_dof_handler(triangulation) + , quadrature_temperature(temperature_degree + 1) + , quadrature_stokes(stokes_degree + 2) + , time_step(0) + , old_time_step(0) + , timestep_number(0) + , rebuild_stokes_matrix(true) + , rebuild_temperature_matrices(true) + , rebuild_stokes_preconditioner(true) + {} + template + double + BoussinesqFlowProblem::get_maximal_velocity() const + { + const auto & quadrature_formula = quadrature_stokes; + const unsigned int n_q_points = quadrature_formula.size(); + FEValues fe_values(mapping, + stokes_fe, + quadrature_formula, + update_values); + std::vector> velocity_values(n_q_points); + double max_velocity = 0; + const FEValuesExtractors::Vector velocities(0); + for (const auto &cell : stokes_dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + fe_values[velocities].get_function_values(stokes_solution, + velocity_values); + for (unsigned int q = 0; q < n_q_points; ++q) + max_velocity = std::max(max_velocity, velocity_values[q].norm()); + } + return max_velocity; + } + template + std::pair + BoussinesqFlowProblem::get_extrapolated_temperature_range() const + { + const auto & quadrature_formula = quadrature_temperature; + const unsigned int n_q_points = quadrature_formula.size(); + FEValues fe_values(mapping, + temperature_fe, + quadrature_formula, + update_values); + std::vector old_temperature_values(n_q_points); + std::vector old_old_temperature_values(n_q_points); + if (timestep_number != 0) + { + double min_temperature = std::numeric_limits::max(), + max_temperature = -std::numeric_limits::max(); + for (const auto &cell : temperature_dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + fe_values.get_function_values(old_temperature_solution, + old_temperature_values); + fe_values.get_function_values(old_old_temperature_solution, + old_old_temperature_values); + for (unsigned int q = 0; q < n_q_points; ++q) + { + const double temperature = + (1. + time_step / old_time_step) * old_temperature_values[q] - + time_step / old_time_step * old_old_temperature_values[q]; + min_temperature = std::min(min_temperature, temperature); + max_temperature = std::max(max_temperature, temperature); + } + } + return std::make_pair(min_temperature, max_temperature); + } + else + { + double min_temperature = std::numeric_limits::max(), + max_temperature = -std::numeric_limits::max(); + for (const auto &cell : temperature_dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + fe_values.get_function_values(old_temperature_solution, + old_temperature_values); + for (unsigned int q = 0; q < n_q_points; ++q) + { + const double temperature = old_temperature_values[q]; + min_temperature = std::min(min_temperature, temperature); + max_temperature = std::max(max_temperature, temperature); + } + } + return std::make_pair(min_temperature, max_temperature); + } + } + template + double + BoussinesqFlowProblem::compute_viscosity( + const std::vector & old_temperature, + const std::vector & old_old_temperature, + const std::vector> &old_temperature_grads, + const std::vector> &old_old_temperature_grads, + const std::vector & old_temperature_laplacians, + const std::vector & old_old_temperature_laplacians, + const std::vector> &old_velocity_values, + const std::vector> &old_old_velocity_values, + const std::vector & gamma_values, + const double global_u_infty, + const double global_T_variation, + const double cell_diameter) const + { + constexpr double beta = 0.017 * dim; + constexpr double alpha = 1.0; + if (global_u_infty == 0) + return 5e-3 * cell_diameter; + const unsigned int n_q_points = old_temperature.size(); + double max_residual = 0; + double max_velocity = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + const Tensor<1, dim> u = + (old_velocity_values[q] + old_old_velocity_values[q]) / 2; + const double dT_dt = + (old_temperature[q] - old_old_temperature[q]) / old_time_step; + const double u_grad_T = + u * (old_temperature_grads[q] + old_old_temperature_grads[q]) / 2; + const double kappa_Delta_T = + EquationData::kappa * + (old_temperature_laplacians[q] + old_old_temperature_laplacians[q]) / + 2; + const double residual = + std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) * + std::pow((old_temperature[q] + old_old_temperature[q]) / 2, + alpha - 1.)); + max_residual = std::max(residual, max_residual); + max_velocity = std::max(std::sqrt(u * u), max_velocity); + } + const double c_R = std::pow(2., (4. - 2 * alpha) / dim); + const double global_scaling = c_R * global_u_infty * global_T_variation * + std::pow(global_Omega_diameter, alpha - 2.); + return ( + beta * max_velocity * + std::min(cell_diameter, + std::pow(cell_diameter, alpha) * max_residual / global_scaling)); + } + template + void + BoussinesqFlowProblem::setup_dofs() + { + std::vector stokes_sub_blocks(dim + 1, 0); + stokes_sub_blocks[dim] = 1; + { + stokes_dof_handler.distribute_dofs(stokes_fe); + DoFRenumbering::component_wise(stokes_dof_handler, stokes_sub_blocks); + stokes_constraints.clear(); + DoFTools::make_hanging_node_constraints(stokes_dof_handler, + stokes_constraints); + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert(0); + VectorTools::compute_no_normal_flux_constraints(stokes_dof_handler, + 0, + no_normal_flux_boundaries, + stokes_constraints, + mapping); + stokes_constraints.close(); + } + { + temperature_dof_handler.distribute_dofs(temperature_fe); + temperature_constraints.clear(); + DoFTools::make_hanging_node_constraints(temperature_dof_handler, + temperature_constraints); + temperature_constraints.close(); + } + const std::vector stokes_dofs_per_block = + DoFTools::count_dofs_per_fe_block(stokes_dof_handler, stokes_sub_blocks); + const unsigned int n_u = stokes_dofs_per_block[0], + n_p = stokes_dofs_per_block[1], + n_T = temperature_dof_handler.n_dofs(); + deallog << "Number of active cells: " << triangulation.n_active_cells() + << " (on " << triangulation.n_levels() << " levels)" << std::endl + << "Number of degrees of freedom: " << n_u + n_p + n_T << " (" + << n_u << '+' << n_p << '+' << n_T << ')' << std::endl + << std::endl; + stokes_partitioning.resize(2); + stokes_partitioning[0] = complete_index_set(n_u); + stokes_partitioning[1] = complete_index_set(n_p); + { + stokes_matrix.clear(); + BlockDynamicSparsityPattern dsp(2, 2); + dsp.block(0, 0).reinit(n_u, n_u); + dsp.block(0, 1).reinit(n_u, n_p); + dsp.block(1, 0).reinit(n_p, n_u); + dsp.block(1, 1).reinit(n_p, n_p); + dsp.collect_sizes(); + Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1); + for (unsigned int c = 0; c < dim + 1; ++c) + for (unsigned int d = 0; d < dim + 1; ++d) + if (!((c == dim) && (d == dim))) + coupling[c][d] = DoFTools::always; + else + coupling[c][d] = DoFTools::none; + DoFTools::make_sparsity_pattern( + stokes_dof_handler, coupling, dsp, stokes_constraints, false); + stokes_matrix.reinit(dsp); + } + { + Amg_preconditioner.reset(); + Mp_preconditioner.reset(); + stokes_preconditioner_matrix.clear(); + BlockDynamicSparsityPattern dsp(2, 2); + dsp.block(0, 0).reinit(n_u, n_u); + dsp.block(0, 1).reinit(n_u, n_p); + dsp.block(1, 0).reinit(n_p, n_u); + dsp.block(1, 1).reinit(n_p, n_p); + dsp.collect_sizes(); + Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1); + for (unsigned int c = 0; c < dim + 1; ++c) + for (unsigned int d = 0; d < dim + 1; ++d) + if (c == d) + coupling[c][d] = DoFTools::always; + else + coupling[c][d] = DoFTools::none; + DoFTools::make_sparsity_pattern( + stokes_dof_handler, coupling, dsp, stokes_constraints, false); + stokes_preconditioner_matrix.reinit(dsp); + } + { + temperature_mass_matrix.clear(); + temperature_stiffness_matrix.clear(); + temperature_matrix.clear(); + DynamicSparsityPattern dsp(n_T, n_T); + DoFTools::make_sparsity_pattern(temperature_dof_handler, + dsp, + temperature_constraints, + false); + temperature_matrix.reinit(dsp); + temperature_mass_matrix.reinit(temperature_matrix); + temperature_stiffness_matrix.reinit(temperature_matrix); + } + IndexSet temperature_partitioning = complete_index_set(n_T); + stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD); + old_stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD); + stokes_rhs.reinit(stokes_partitioning, MPI_COMM_WORLD); + temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD); + old_temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD); + old_old_temperature_solution.reinit(temperature_partitioning, + MPI_COMM_WORLD); + temperature_rhs.reinit(temperature_partitioning, MPI_COMM_WORLD); + } + template + void + BoussinesqFlowProblem::assemble_stokes_preconditioner() + { + stokes_preconditioner_matrix = 0; + const auto & quadrature_formula = quadrature_stokes; + FEValues stokes_fe_values(mapping, + stokes_fe, + quadrature_formula, + update_JxW_values | update_values | + update_gradients); + const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell(); + const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + std::vector> grad_phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + for (const auto &cell : stokes_dof_handler.active_cell_iterators()) + { + stokes_fe_values.reinit(cell); + local_matrix = 0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + grad_phi_u[k] = stokes_fe_values[velocities].gradient(k, q); + phi_p[k] = stokes_fe_values[pressure].value(k, q); + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + local_matrix(i, j) += + (EquationData::eta * + scalar_product(grad_phi_u[i], grad_phi_u[j]) + + (1. / EquationData::eta) * phi_p[i] * phi_p[j]) * + stokes_fe_values.JxW(q); + } + cell->get_dof_indices(local_dof_indices); + stokes_constraints.distribute_local_to_global( + local_matrix, local_dof_indices, stokes_preconditioner_matrix); + } + } + template + void + BoussinesqFlowProblem::build_stokes_preconditioner() + { + if (rebuild_stokes_preconditioner == false) + return; + deallog << " Rebuilding Stokes preconditioner..." << std::flush; + assemble_stokes_preconditioner(); + Amg_preconditioner = std::make_shared(); + std::vector> constant_modes; + FEValuesExtractors::Vector velocity_components(0); + DoFTools::extract_constant_modes(stokes_dof_handler, + stokes_fe.component_mask( + velocity_components), + constant_modes); + TrilinosWrappers::PreconditionAMG::AdditionalData amg_data; + amg_data.constant_modes = constant_modes; + amg_data.elliptic = true; + amg_data.higher_order_elements = true; + amg_data.smoother_sweeps = 2; + amg_data.aggregation_threshold = 0.02; + Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0, 0), + amg_data); + Mp_preconditioner = std::make_shared(); + Mp_preconditioner->initialize(stokes_preconditioner_matrix.block(1, 1)); + deallog << std::endl; + rebuild_stokes_preconditioner = false; + } + template + void + BoussinesqFlowProblem::assemble_stokes_system() + { + deallog << " Assembling..." << std::flush; + if (rebuild_stokes_matrix == true) + stokes_matrix = 0; + stokes_rhs = 0; + const auto & quadrature_formula = quadrature_stokes; + FEValues stokes_fe_values( + mapping, + stokes_fe, + quadrature_formula, + update_values | update_quadrature_points | update_JxW_values | + (rebuild_stokes_matrix == true ? update_gradients : UpdateFlags(0))); + FEValues temperature_fe_values(mapping, + temperature_fe, + quadrature_formula, + update_values); + const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell(); + const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + Vector local_rhs(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + std::vector old_temperature_values(n_q_points); + std::vector> phi_u(dofs_per_cell); + std::vector> grads_phi_u(dofs_per_cell); + std::vector div_phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + auto cell = stokes_dof_handler.begin_active(); + const auto endc = stokes_dof_handler.end(); + auto temperature_cell = temperature_dof_handler.begin_active(); + for (; cell != endc; ++cell, ++temperature_cell) + { + stokes_fe_values.reinit(cell); + temperature_fe_values.reinit(temperature_cell); + local_matrix = 0; + local_rhs = 0; + temperature_fe_values.get_function_values(old_temperature_solution, + old_temperature_values); + for (unsigned int q = 0; q < n_q_points; ++q) + { + const double old_temperature = old_temperature_values[q]; + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + phi_u[k] = stokes_fe_values[velocities].value(k, q); + if (rebuild_stokes_matrix) + { + grads_phi_u[k] = + stokes_fe_values[velocities].symmetric_gradient(k, q); + div_phi_u[k] = + stokes_fe_values[velocities].divergence(k, q); + phi_p[k] = stokes_fe_values[pressure].value(k, q); + } + } + if (rebuild_stokes_matrix) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + local_matrix(i, j) += + (EquationData::eta * 2 * (grads_phi_u[i] * grads_phi_u[j]) - + div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) * + stokes_fe_values.JxW(q); + const Point gravity = + -((dim == 2) ? (Point(0, 1)) : (Point(0, 0, 1))); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + local_rhs(i) += (-EquationData::density * EquationData::beta * + gravity * phi_u[i] * old_temperature) * + stokes_fe_values.JxW(q); + } + cell->get_dof_indices(local_dof_indices); + if (rebuild_stokes_matrix == true) + stokes_constraints.distribute_local_to_global(local_matrix, + local_rhs, + local_dof_indices, + stokes_matrix, + stokes_rhs); + else + stokes_constraints.distribute_local_to_global(local_rhs, + local_dof_indices, + stokes_rhs); + } + rebuild_stokes_matrix = false; + deallog << std::endl; + } + template + void + BoussinesqFlowProblem::assemble_temperature_matrix() + { + if (rebuild_temperature_matrices == false) + return; + temperature_mass_matrix = 0; + temperature_stiffness_matrix = 0; + const auto & quadrature_formula = quadrature_temperature; + FEValues temperature_fe_values(mapping, + temperature_fe, + quadrature_formula, + update_values | update_gradients | + update_JxW_values); + const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell(); + const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix local_mass_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix local_stiffness_matrix(dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + std::vector phi_T(dofs_per_cell); + std::vector> grad_phi_T(dofs_per_cell); + for (const auto &cell : temperature_dof_handler.active_cell_iterators()) + { + local_mass_matrix = 0; + local_stiffness_matrix = 0; + temperature_fe_values.reinit(cell); + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + grad_phi_T[k] = temperature_fe_values.shape_grad(k, q); + phi_T[k] = temperature_fe_values.shape_value(k, q); + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + local_mass_matrix(i, j) += + (phi_T[i] * phi_T[j] * temperature_fe_values.JxW(q)); + local_stiffness_matrix(i, j) += + (EquationData::kappa * grad_phi_T[i] * grad_phi_T[j] * + temperature_fe_values.JxW(q)); + } + } + cell->get_dof_indices(local_dof_indices); + temperature_constraints.distribute_local_to_global( + local_mass_matrix, local_dof_indices, temperature_mass_matrix); + temperature_constraints.distribute_local_to_global( + local_stiffness_matrix, + local_dof_indices, + temperature_stiffness_matrix); + } + rebuild_temperature_matrices = false; + } + template + void + BoussinesqFlowProblem::assemble_temperature_system( + const double maximal_velocity) + { + const bool use_bdf2_scheme = (timestep_number != 0); + if (use_bdf2_scheme == true) + { + temperature_matrix.copy_from(temperature_mass_matrix); + temperature_matrix *= + (2 * time_step + old_time_step) / (time_step + old_time_step); + temperature_matrix.add(time_step, temperature_stiffness_matrix); + } + else + { + temperature_matrix.copy_from(temperature_mass_matrix); + temperature_matrix.add(time_step, temperature_stiffness_matrix); + } + temperature_rhs = 0; + const auto & quadrature_formula = quadrature_temperature; + FEValues temperature_fe_values(mapping, + temperature_fe, + quadrature_formula, + update_values | update_gradients | + update_hessians | + update_quadrature_points | + update_JxW_values); + FEValues stokes_fe_values(stokes_fe, + quadrature_formula, + update_values); + const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell(); + const unsigned int n_q_points = quadrature_formula.size(); + Vector local_rhs(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + std::vector> old_velocity_values(n_q_points); + std::vector> old_old_velocity_values(n_q_points); + std::vector old_temperature_values(n_q_points); + std::vector old_old_temperature_values(n_q_points); + std::vector> old_temperature_grads(n_q_points); + std::vector> old_old_temperature_grads(n_q_points); + std::vector old_temperature_laplacians(n_q_points); + std::vector old_old_temperature_laplacians(n_q_points); + EquationData::TemperatureRightHandSide temperature_right_hand_side; + std::vector gamma_values(n_q_points); + std::vector phi_T(dofs_per_cell); + std::vector> grad_phi_T(dofs_per_cell); + const std::pair global_T_range = + get_extrapolated_temperature_range(); + const FEValuesExtractors::Vector velocities(0); + auto cell = temperature_dof_handler.begin_active(); + const auto endc = temperature_dof_handler.end(); + auto stokes_cell = stokes_dof_handler.begin_active(); + for (; cell != endc; ++cell, ++stokes_cell) + { + local_rhs = 0; + temperature_fe_values.reinit(cell); + stokes_fe_values.reinit(stokes_cell); + temperature_fe_values.get_function_values(old_temperature_solution, + old_temperature_values); + temperature_fe_values.get_function_values(old_old_temperature_solution, + old_old_temperature_values); + temperature_fe_values.get_function_gradients(old_temperature_solution, + old_temperature_grads); + temperature_fe_values.get_function_gradients( + old_old_temperature_solution, old_old_temperature_grads); + temperature_fe_values.get_function_laplacians( + old_temperature_solution, old_temperature_laplacians); + temperature_fe_values.get_function_laplacians( + old_old_temperature_solution, old_old_temperature_laplacians); + temperature_right_hand_side.value_list( + temperature_fe_values.get_quadrature_points(), gamma_values); + stokes_fe_values[velocities].get_function_values(stokes_solution, + old_velocity_values); + stokes_fe_values[velocities].get_function_values( + old_stokes_solution, old_old_velocity_values); + const double nu = + compute_viscosity(old_temperature_values, + old_old_temperature_values, + old_temperature_grads, + old_old_temperature_grads, + old_temperature_laplacians, + old_old_temperature_laplacians, + old_velocity_values, + old_old_velocity_values, + gamma_values, + maximal_velocity, + global_T_range.second - global_T_range.first, + cell->diameter()); + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + grad_phi_T[k] = temperature_fe_values.shape_grad(k, q); + phi_T[k] = temperature_fe_values.shape_value(k, q); + } + const double T_term_for_rhs = + (use_bdf2_scheme ? + (old_temperature_values[q] * (1 + time_step / old_time_step) - + old_old_temperature_values[q] * (time_step * time_step) / + (old_time_step * (time_step + old_time_step))) : + old_temperature_values[q]); + const Tensor<1, dim> ext_grad_T = + (use_bdf2_scheme ? + (old_temperature_grads[q] * (1 + time_step / old_time_step) - + old_old_temperature_grads[q] * time_step / old_time_step) : + old_temperature_grads[q]); + const Tensor<1, dim> extrapolated_u = + (use_bdf2_scheme ? + (old_velocity_values[q] * (1 + time_step / old_time_step) - + old_old_velocity_values[q] * time_step / old_time_step) : + old_velocity_values[q]); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + local_rhs(i) += + (T_term_for_rhs * phi_T[i] - + time_step * extrapolated_u * ext_grad_T * phi_T[i] - + time_step * nu * ext_grad_T * grad_phi_T[i] + + time_step * gamma_values[q] * phi_T[i]) * + temperature_fe_values.JxW(q); + } + cell->get_dof_indices(local_dof_indices); + temperature_constraints.distribute_local_to_global(local_rhs, + local_dof_indices, + temperature_rhs); + } + } + template + void + BoussinesqFlowProblem::solve() + { + deallog << " Solving..." << std::endl; + { + const LinearSolvers::InverseMatrix + mp_inverse(stokes_preconditioner_matrix.block(1, 1), + *Mp_preconditioner); + const LinearSolvers::BlockSchurPreconditioner< + TrilinosWrappers::PreconditionAMG, + TrilinosWrappers::PreconditionIC> + preconditioner(stokes_matrix, mp_inverse, *Amg_preconditioner); + SolverControl solver_control(stokes_matrix.m(), + 1e-6 * stokes_rhs.l2_norm()); + SolverGMRES gmres( + solver_control, + SolverGMRES::AdditionalData(100)); + for (unsigned int i = 0; i < stokes_solution.size(); ++i) + if (stokes_constraints.is_constrained(i)) + stokes_solution(i) = 0; + gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner); + stokes_constraints.distribute(stokes_solution); + deallog << " " << solver_control.last_step() + << " GMRES iterations for Stokes subsystem." << std::endl; + } + old_time_step = time_step; + const double maximal_velocity = get_maximal_velocity(); + if (maximal_velocity >= 0.01) + time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree * + GridTools::minimal_cell_diameter(triangulation) / + maximal_velocity; + else + time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree * + GridTools::minimal_cell_diameter(triangulation) / .01; + deallog << " " + << "Time step: " << time_step << std::endl; + temperature_solution = old_temperature_solution; + assemble_temperature_system(maximal_velocity); + { + SolverControl solver_control(temperature_matrix.m(), + 1e-8 * temperature_rhs.l2_norm()); + SolverCG cg(solver_control); + TrilinosWrappers::PreconditionIC preconditioner; + preconditioner.initialize(temperature_matrix); + cg.solve(temperature_matrix, + temperature_solution, + temperature_rhs, + preconditioner); + temperature_constraints.distribute(temperature_solution); + deallog << " " << solver_control.last_step() + << " CG iterations for temperature." << std::endl; + double min_temperature = temperature_solution(0), + max_temperature = temperature_solution(0); + for (unsigned int i = 0; i < temperature_solution.size(); ++i) + { + min_temperature = + std::min(min_temperature, temperature_solution(i)); + max_temperature = + std::max(max_temperature, temperature_solution(i)); + } + deallog << " Temperature range: " << min_temperature << ' ' + << max_temperature << std::endl; + } + } + template + void + BoussinesqFlowProblem::output_results() const + { + if (timestep_number % 10 != 0) + return; + std::vector stokes_names(dim, "velocity"); + stokes_names.emplace_back("p"); + std::vector + stokes_component_interpretation( + dim + 1, DataComponentInterpretation::component_is_scalar); + for (unsigned int i = 0; i < dim; ++i) + stokes_component_interpretation[i] = + DataComponentInterpretation::component_is_part_of_vector; + DataOut data_out; + data_out.add_data_vector(stokes_dof_handler, + stokes_solution, + stokes_names, + stokes_component_interpretation); + data_out.add_data_vector(temperature_dof_handler, + temperature_solution, + "T"); + data_out.build_patches(mapping, + std::max(stokes_degree, temperature_degree)); + std::ofstream output("solution-" + + Utilities::int_to_string(timestep_number, 4) + ".vtk"); + data_out.write_vtk(output); + } + template + void + BoussinesqFlowProblem::refine_mesh(const unsigned int max_grid_level) + { + return; + + Vector estimated_error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator::estimate(temperature_dof_handler, + QGauss(temperature_degree + 1), + {}, + temperature_solution, + estimated_error_per_cell); + GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, + estimated_error_per_cell, + 0.8, + 0.1); + if (triangulation.n_levels() > max_grid_level) + for (auto &cell : + triangulation.active_cell_iterators_on_level(max_grid_level)) + cell->clear_refine_flag(); + std::vector x_temperature(2); + x_temperature[0] = temperature_solution; + x_temperature[1] = old_temperature_solution; + TrilinosWrappers::MPI::BlockVector x_stokes = stokes_solution; + SolutionTransfer temperature_trans( + temperature_dof_handler); + SolutionTransfer stokes_trans( + stokes_dof_handler); + triangulation.prepare_coarsening_and_refinement(); + temperature_trans.prepare_for_coarsening_and_refinement(x_temperature); + stokes_trans.prepare_for_coarsening_and_refinement(x_stokes); + triangulation.execute_coarsening_and_refinement(); + setup_dofs(); + std::vector tmp(2); + tmp[0].reinit(temperature_solution); + tmp[1].reinit(temperature_solution); + temperature_trans.interpolate(x_temperature, tmp); + temperature_solution = tmp[0]; + old_temperature_solution = tmp[1]; + temperature_constraints.distribute(temperature_solution); + temperature_constraints.distribute(old_temperature_solution); + stokes_trans.interpolate(x_stokes, stokes_solution); + stokes_constraints.distribute(stokes_solution); + rebuild_stokes_matrix = true; + rebuild_temperature_matrices = true; + rebuild_stokes_preconditioner = true; + } + template + void + BoussinesqFlowProblem::run() + { + const unsigned int initial_refinement = (dim == 2 ? 4 : 2); + const unsigned int n_pre_refinement_steps = (dim == 2 ? 0 : 0); + + Triangulation tria_temp; + GridGenerator::hyper_cube(tria_temp); + + GridGenerator::convert_hypercube_to_simplex_mesh(tria_temp, triangulation); + + global_Omega_diameter = GridTools::diameter(triangulation); + triangulation.refine_global(initial_refinement); + setup_dofs(); + unsigned int pre_refinement_step = 0; + start_time_iteration: + VectorTools::project(mapping, + temperature_dof_handler, + temperature_constraints, + quadrature_temperature, + EquationData::TemperatureInitialValues(), + old_temperature_solution); + timestep_number = 0; + time_step = old_time_step = 0; + double time = 0; + do + { + deallog << "Timestep " << timestep_number << ": t=" << time + << std::endl; + assemble_stokes_system(); + build_stokes_preconditioner(); + assemble_temperature_matrix(); + solve(); + output_results(); + deallog << std::endl; + if ((timestep_number == 0) && + (pre_refinement_step < n_pre_refinement_steps)) + { + refine_mesh(initial_refinement + n_pre_refinement_steps); + ++pre_refinement_step; + goto start_time_iteration; + } + else if ((timestep_number > 0) && (timestep_number % 5 == 0)) + refine_mesh(initial_refinement + n_pre_refinement_steps); + time += time_step; + ++timestep_number; + old_stokes_solution = stokes_solution; + old_old_temperature_solution = old_temperature_solution; + old_temperature_solution = temperature_solution; + } + while (time <= 2); + } +} // namespace Step31 +int +main(int argc, char *argv[]) +{ + try + { + using namespace dealii; + using namespace Step31; + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + initlog(); + deallog.depth_file(1); + + AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1, + ExcMessage( + "This program can only be run in serial, use ./step-31")); + BoussinesqFlowProblem<2> flow_problem; + flow_problem.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + return 0; +} diff --git a/tests/simplex/step-31.with_simplex_support=on.with_trilinos=true.output b/tests/simplex/step-31.with_simplex_support=on.with_trilinos=true.output new file mode 100644 index 0000000000..07eb3ff8f2 --- /dev/null +++ b/tests/simplex/step-31.with_simplex_support=on.with_trilinos=true.output @@ -0,0 +1,45 @@ + +DEAL::Number of active cells: 2048 (on 5 levels) +DEAL::Number of degrees of freedom: 13764 (8450+1089+4225) +DEAL:: +DEAL::Timestep 0: t=0.00000 +DEAL:: Assembling... +DEAL:: Rebuilding Stokes preconditioner... +DEAL:: Solving... +DEAL:: 0 GMRES iterations for Stokes subsystem. +DEAL:: Time step: 0.459559 +DEAL:: 9 CG iterations for temperature. +DEAL:: Temperature range: -0.141015 0.586975 +DEAL:: +DEAL::Timestep 1: t=0.459559 +DEAL:: Assembling... +DEAL:: Solving... +DEAL:: 24 GMRES iterations for Stokes subsystem. +DEAL:: Time step: 0.459559 +DEAL:: 9 CG iterations for temperature. +DEAL:: Temperature range: -0.301941 1.17152 +DEAL:: +DEAL::Timestep 2: t=0.919118 +DEAL:: Assembling... +DEAL:: Solving... +DEAL:: 23 GMRES iterations for Stokes subsystem. +DEAL:: Time step: 0.459559 +DEAL:: 9 CG iterations for temperature. +DEAL:: Temperature range: -0.467381 1.71946 +DEAL:: +DEAL::Timestep 3: t=1.37868 +DEAL:: Assembling... +DEAL:: Solving... +DEAL:: 23 GMRES iterations for Stokes subsystem. +DEAL:: Time step: 0.459559 +DEAL:: 9 CG iterations for temperature. +DEAL:: Temperature range: -0.612772 2.18114 +DEAL:: +DEAL::Timestep 4: t=1.83824 +DEAL:: Assembling... +DEAL:: Solving... +DEAL:: 22 GMRES iterations for Stokes subsystem. +DEAL:: Time step: 0.459559 +DEAL:: 9 CG iterations for temperature. +DEAL:: Temperature range: -0.705842 2.76592 +DEAL::