From: Martin Kronbichler Date: Tue, 1 Nov 2022 16:33:23 +0000 (+0100) Subject: Add performance test for compressible Navier-Stokes equations X-Git-Tag: v9.5.0-rc1~860^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4caf49f7dc022916360090df751c62b461141010;p=dealii.git Add performance test for compressible Navier-Stokes equations --- diff --git a/tests/performance/timing_navier_stokes.cc b/tests/performance/timing_navier_stokes.cc new file mode 100644 index 0000000000..e6293f4bc7 --- /dev/null +++ b/tests/performance/timing_navier_stokes.cc @@ -0,0 +1,2298 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// +// Description: +// +// A performance benchmark based on step 67 and step 76, extended to the +// compressible Navier-Stokes equations. We measure timings for grid creation, +// setup of unknowns, explicit Runge-Kutta time stepping with face-centric +// loop through MatrixFree::loop() (step-67 style) as well as cell-centric +// loop through MatrixFree::loop_cell_centric() (step-76 style). We use a +// problem with periodic boundary conditions to avoid +// +// Status: experimental +// + +#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 + + +#define ENABLE_MPI + +#include "performance_test_driver.h" + +namespace NavierStokes_DG +{ + using namespace dealii; + + constexpr unsigned int testcase = 2; + constexpr unsigned int dimension = 3; + constexpr unsigned int fe_degree = 4; + constexpr unsigned int n_q_points_1d = fe_degree + 2; + + constexpr unsigned int group_size = numbers::invalid_unsigned_int; + + using Number = double; + + using VectorizedArrayType = VectorizedArray; + + constexpr double gamma = 1.4; + constexpr double R = 287.; + constexpr double c_v = R / (gamma - 1.); + constexpr double c_p = gamma / c_v; + constexpr double viscosity = 1. / 1600; + constexpr double lambda = viscosity * c_p / 0.71; + constexpr double Ma = 0.1; + + const double courant_number = 0.07 / std::pow(fe_degree, 1.5); + + enum LowStorageRungeKuttaScheme + { + stage_3_order_3, + stage_5_order_4, + stage_7_order_4, + stage_9_order_5, + }; + constexpr LowStorageRungeKuttaScheme lsrk_scheme = stage_3_order_3; + + + + class LowStorageRungeKuttaIntegrator + { + public: + LowStorageRungeKuttaIntegrator(const LowStorageRungeKuttaScheme scheme) + { + TimeStepping::runge_kutta_method lsrk; + switch (scheme) + { + case stage_3_order_3: + lsrk = TimeStepping::LOW_STORAGE_RK_STAGE3_ORDER3; + break; + case stage_5_order_4: + lsrk = TimeStepping::LOW_STORAGE_RK_STAGE5_ORDER4; + break; + case stage_7_order_4: + lsrk = TimeStepping::LOW_STORAGE_RK_STAGE7_ORDER4; + break; + case stage_9_order_5: + lsrk = TimeStepping::LOW_STORAGE_RK_STAGE9_ORDER5; + break; + + default: + AssertThrow(false, ExcNotImplemented()); + } + TimeStepping::LowStorageRungeKutta< + LinearAlgebra::distributed::Vector> + rk_integrator(lsrk); + std::vector ci; // not used + rk_integrator.get_coefficients(ai, bi, ci); + } + + unsigned int + n_stages() const + { + return bi.size(); + } + + template + void + perform_time_step(const Operator &pde_operator, + const double current_time, + const double time_step, + VectorType & solution, + VectorType & vec_ri, + VectorType & vec_ki) const + { + AssertDimension(ai.size() + 1, bi.size()); + + vec_ki.swap(solution); + + double sum_previous_bi = 0; + for (unsigned int stage = 0; stage < bi.size(); ++stage) + { + const double c_i = stage == 0 ? 0 : sum_previous_bi + ai[stage - 1]; + + pde_operator.perform_stage(stage, + current_time + c_i * time_step, + bi[stage] * time_step, + (stage == bi.size() - 1 ? + 0 : + ai[stage] * time_step), + (stage % 2 == 0 ? vec_ki : vec_ri), + (stage % 2 == 0 ? vec_ri : vec_ki), + solution); + + if (stage > 0) + sum_previous_bi += bi[stage - 1]; + } + } + + private: + std::vector bi; + std::vector ai; + }; + + + enum EulerNumericalFlux + { + lax_friedrichs_modified, + harten_lax_vanleer, + }; + constexpr EulerNumericalFlux numerical_flux_type = harten_lax_vanleer; + + + + template + class ExactSolution : public Function + { + public: + ExactSolution(const double time) + : Function(dim + 2, time) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const override; + }; + + + + template + double + ExactSolution::value(const Point & x, + const unsigned int component) const + { + const double t = this->get_time(); + + switch (testcase) + { + case 0: + { + Assert(dim == 2, ExcNotImplemented()); + const double beta = 5; + + Point x0; + x0[0] = 5.; + const double radius_sqr = + (x - x0).norm_square() - 2. * (x[0] - x0[0]) * t + t * t; + const double factor = + beta / (numbers::PI * 2) * std::exp(1. - radius_sqr); + const double density_log = std::log2( + std::abs(1. - (gamma - 1.) / gamma * 0.25 * factor * factor)); + const double density = std::exp2(density_log * (1. / (gamma - 1.))); + const double u = 1. - factor * (x[1] - x0[1]); + const double v = factor * (x[0] - t - x0[0]); + + if (component == 0) + return density; + else if (component == 1) + return density * u; + else if (component == 2) + return density * v; + else + { + const double pressure = + std::exp2(density_log * (gamma / (gamma - 1.))); + return pressure / (gamma - 1.) + + 0.5 * (density * u * u + density * v * v); + } + } + + case 1: + { + if (component == 0) + return 1.; + else if (component == 1) + return 0.4; + else if (component == dim + 1) + return 3.097857142857143; + else + return 0.; + } + + case 2: + { + AssertThrow(dim == 3, ExcNotImplemented()); + const double c0 = 1. / Ma; + const double T0 = c0 * c0 / (gamma * R); + if (component == 0) + return 1 + 1. / (R * T0) * 1. / 16. * + (std::cos(2 * x[0]) + std::cos(2 * x[1])) * + (std::cos(2 * x[2]) + 2.); + else if (component == 1) + return std::sin(x[0]) * std::cos(x[1]) * std::cos(x[2]); + else if (component == 2) + return -std::cos(x[0]) * std::sin(x[1]) * std::cos(x[2]); + else if (component == 3) + return 0.; + else + return c_v * T0 + + 0.5 * + (Utilities::fixed_power<2>( + std::sin(x[0]) * std::cos(x[1]) * std::cos(x[2])) + + Utilities::fixed_power<2>( + std::cos(x[0]) * std::sin(x[1]) * std::cos(x[2]))); + } + + default: + Assert(false, ExcNotImplemented()); + return 0.; + } + } + + + + template + inline DEAL_II_ALWAYS_INLINE // + Tensor<1, dim, Number> + fluid_velocity(const Tensor<1, dim + 2, Number> &conserved_variables) + { + const Number inverse_density = Number(1.) / conserved_variables[0]; + + Tensor<1, dim, Number> velocity; + for (unsigned int d = 0; d < dim; ++d) + velocity[d] = conserved_variables[1 + d] * inverse_density; + + return velocity; + } + + template + inline DEAL_II_ALWAYS_INLINE // + Number + fluid_pressure(const Tensor<1, dim + 2, Number> &conserved_variables) + { + const Tensor<1, dim, Number> velocity = + fluid_velocity(conserved_variables); + + Number rho_u_dot_u = conserved_variables[1] * velocity[0]; + for (unsigned int d = 1; d < dim; ++d) + rho_u_dot_u += conserved_variables[1 + d] * velocity[d]; + + return (gamma - 1.) * (conserved_variables[dim + 1] - 0.5 * rho_u_dot_u); + } + + template + inline DEAL_II_ALWAYS_INLINE // + Tensor<1, dim + 2, Tensor<1, dim, Number>> + euler_flux(const Tensor<1, dim + 2, Number> &conserved_variables) + { + const Tensor<1, dim, Number> velocity = + fluid_velocity(conserved_variables); + const Number pressure = fluid_pressure(conserved_variables); + + Tensor<1, dim + 2, Tensor<1, dim, Number>> flux; + for (unsigned int d = 0; d < dim; ++d) + { + flux[0][d] = conserved_variables[1 + d]; + for (unsigned int e = 0; e < dim; ++e) + flux[e + 1][d] = conserved_variables[e + 1] * velocity[d]; + flux[d + 1][d] += pressure; + flux[dim + 1][d] = + velocity[d] * (conserved_variables[dim + 1] + pressure); + } + + return flux; + } + + template + inline DEAL_II_ALWAYS_INLINE // + Tensor<1, n_components, Number> + operator*(const Tensor<1, n_components, Tensor<1, dim, Number>> &matrix, + const Tensor<1, dim, Number> & vector) + { + Tensor<1, n_components, Number> result; + for (unsigned int d = 0; d < n_components; ++d) + result[d] = matrix[d] * vector; + return result; + } + + template + inline DEAL_II_ALWAYS_INLINE // + Tensor<1, dim + 2, Number> + euler_numerical_flux(const Tensor<1, dim + 2, Number> &u_m, + const Tensor<1, dim + 2, Number> &u_p, + const Tensor<1, dim, Number> & normal) + { + const auto velocity_m = fluid_velocity(u_m); + const auto velocity_p = fluid_velocity(u_p); + + const auto pressure_m = fluid_pressure(u_m); + const auto pressure_p = fluid_pressure(u_p); + + const auto flux_m = euler_flux(u_m); + const auto flux_p = euler_flux(u_p); + + switch (numerical_flux_type) + { + case lax_friedrichs_modified: + { + const auto Lambda = + 0.5 * std::sqrt(std::max(velocity_p.norm_square() + + gamma * pressure_p * (1. / u_p[0]), + velocity_m.norm_square() + + gamma * pressure_m * (1. / u_m[0]))); + + return 0.5 * (flux_m * normal + flux_p * normal) + + 0.5 * Lambda * (u_m - u_p); + } + + case harten_lax_vanleer: + { + const auto avg_velocity_normal = + 0.5 * ((velocity_m + velocity_p) * normal); + const auto avg_c = std::sqrt(std::abs( + 0.5 * gamma * + (pressure_p * (1. / u_p[0]) + pressure_m * (1. / u_m[0])))); + const Number s_pos = + std::max(Number(), avg_velocity_normal + avg_c); + const Number s_neg = + std::min(Number(), avg_velocity_normal - avg_c); + const Number inverse_s = Number(1.) / (s_pos - s_neg); + + return inverse_s * + ((s_pos * (flux_m * normal) - s_neg * (flux_p * normal)) - + s_pos * s_neg * (u_m - u_p)); + } + + default: + { + Assert(false, ExcNotImplemented()); + return {}; + } + } + } + + template + inline DEAL_II_ALWAYS_INLINE // + Tensor<2, dim, Number> + fluid_velocity_gradient( + const Tensor<1, dim + 2, Number> & conserved_variables, + const Tensor<1, dim + 2, Tensor<1, dim, Number>> &gradients) + { + const Number inverse_density = Number(1.) / conserved_variables[0]; + const Tensor<1, dim, Number> velocity = + fluid_velocity(conserved_variables); + + Tensor<2, dim, Number> gradient; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + gradient[d][e] = inverse_density * + (gradients[d + 1][e] - velocity[d] * gradients[0][e]); + + return gradient; + } + + template + inline DEAL_II_ALWAYS_INLINE // + Number + fluid_temperature(const Tensor<1, dim + 2, Number> &conserved_variables) + { + const Number inverse_density = Number(1.) / conserved_variables[0]; + const Number inverse_R = 1. / R; + return fluid_pressure(conserved_variables) * inverse_density * inverse_R; + } + + template + inline DEAL_II_ALWAYS_INLINE // + Tensor<1, dim, Number> + fluid_temperature_gradient( + const Tensor<1, dim + 2, Number> & conserved_variables, + const Tensor<1, dim + 2, Tensor<1, dim, Number>> &gradients) + { + const Number inverse_R = 1. / R; + return (gamma - 1.) * inverse_R * + (gradients[dim + 1] - + fluid_velocity(conserved_variables) * + fluid_velocity_gradient(conserved_variables, gradients)); + } + + template + inline DEAL_II_ALWAYS_INLINE // + Tensor<1, dim + 2, Tensor<1, dim, Number>> + viscous_flux(const Tensor<1, dim + 2, Number> &conserved_variables, + const Tensor<1, dim + 2, Tensor<1, dim, Number>> &gradients) + { + const Tensor<1, dim, Number> velocity = + fluid_velocity(conserved_variables); + const Tensor<2, dim, Number> grad_u = + fluid_velocity_gradient(conserved_variables, gradients); + const Number scaled_div_u = viscosity * (2. / 3.) * trace(grad_u); + + Tensor<1, dim + 2, Tensor<1, dim, Number>> result; + for (unsigned int d = 0; d < dim; ++d) + { + for (unsigned int e = d; e < dim; ++e) + { + result[d + 1][e] = viscosity * (grad_u[d][e] + grad_u[e][d]); + result[e + 1][d] = result[d + 1][e]; + } + result[d + 1][d] -= scaled_div_u; + } + + result[dim + 1] = + lambda * fluid_temperature_gradient(conserved_variables, gradients); + for (unsigned int d = 0; d < dim; ++d) + result[dim + 1][d] += result[d + 1] * velocity; + + return result; + } + + + + template + VectorizedArrayType + evaluate_function(const Function & function, + const Point &p_vectorized, + const unsigned int component) + { + VectorizedArrayType result; + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + { + Point p; + for (unsigned int d = 0; d < dim; ++d) + p[d] = p_vectorized[d][v]; + result[v] = function.value(p, component); + } + return result; + } + + + + template + Tensor<1, n_components, VectorizedArrayType> + evaluate_function(const Function & function, + const Point &p_vectorized) + { + AssertDimension(function.n_components, n_components); + Tensor<1, n_components, VectorizedArrayType> result; + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + { + Point p; + for (unsigned int d = 0; d < dim; ++d) + p[d] = p_vectorized[d][v]; + for (unsigned int d = 0; d < n_components; ++d) + result[d][v] = function.value(p, d); + } + return result; + } + + + + template + class NavierStokesOperator + { + public: + static constexpr unsigned int n_quadrature_points_1d = n_points_1d; + + NavierStokesOperator(); + + ~NavierStokesOperator(); + + void + reinit(const Mapping &mapping, const DoFHandler &dof_handler); + + void + set_inflow_boundary(const types::boundary_id boundary_id, + std::unique_ptr> inflow_function); + + void + set_subsonic_outflow_boundary( + const types::boundary_id boundary_id, + std::unique_ptr> outflow_energy); + + void + set_wall_boundary(const types::boundary_id boundary_id); + + void + set_body_force(std::unique_ptr> body_force); + + void + perform_stage(const unsigned int stage, + const Number cur_time, + const Number bi, + const Number ai, + const LinearAlgebra::distributed::Vector ¤t_ri, + LinearAlgebra::distributed::Vector & vec_ki, + LinearAlgebra::distributed::Vector &solution) const; + + void + perform_stage_face( + const unsigned int stage, + const Number cur_time, + const Number bi, + const Number ai, + const LinearAlgebra::distributed::Vector ¤t_ri, + LinearAlgebra::distributed::Vector & vec_ki, + LinearAlgebra::distributed::Vector & solution) const; + + void + project(const Function & function, + LinearAlgebra::distributed::Vector &solution) const; + + std::array + compute_errors( + const Function & function, + const LinearAlgebra::distributed::Vector &solution) const; + + std::array + compute_kinetic_energy( + const LinearAlgebra::distributed::Vector &solution) const; + + double + compute_cell_transport_speed( + const LinearAlgebra::distributed::Vector &solution) const; + + void + initialize_vector(LinearAlgebra::distributed::Vector &vector) const; + + mutable double time_loop; + mutable double time_rk_update; + + private: + MPI_Comm subcommunicator; + + MatrixFree data; + + std::map>> + inflow_boundaries; + std::map>> + subsonic_outflow_boundaries; + std::set wall_boundaries; + std::unique_ptr> body_force; + + void + operation_on_cell(const MatrixFree &mf, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector & src, + const std::pair &range) const; + + void + operation_cell(const MatrixFree &mf, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector & src, + const std::pair &range) const; + + void + operation_face(const MatrixFree &mf, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector & src, + const std::pair &range) const; + + void + operation_boundary( + const MatrixFree &mf, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector & src, + const std::pair & range) const; + + void + local_apply_inverse_mass_matrix( + const MatrixFree &, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const; + + mutable double ai; + mutable double bi; + mutable LinearAlgebra::distributed::Vector *solution; + mutable unsigned int stage; + }; + + + + template + NavierStokesOperator::NavierStokesOperator() + { +#ifdef DEAL_II_WITH_MPI + if (group_size == 1) + { + this->subcommunicator = MPI_COMM_SELF; + } + else if (group_size == numbers::invalid_unsigned_int) + { + const auto rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + MPI_Comm_split_type(MPI_COMM_WORLD, + MPI_COMM_TYPE_SHARED, + rank, + MPI_INFO_NULL, + &subcommunicator); + } + else + { + Assert(false, ExcNotImplemented()); + } +#else + (void)subcommunicator; + (void)group_size; + this->subcommunicator = MPI_COMM_SELF; +#endif + + time_loop = 0.; + time_rk_update = 0.; + } + + + template + NavierStokesOperator::~NavierStokesOperator() + { +#ifdef DEAL_II_WITH_MPI + if (this->subcommunicator != MPI_COMM_SELF) + MPI_Comm_free(&subcommunicator); +#endif + } + + + template + void + NavierStokesOperator::reinit( + const Mapping & mapping, + const DoFHandler &dof_handler) + { + const std::vector *> dof_handlers = {&dof_handler}; + const AffineConstraints dummy; + const std::vector *> constraints = {&dummy}; + const std::vector> quadratures = {QGauss<1>(n_q_points_1d), + QGauss<1>(fe_degree + 1)}; + + typename MatrixFree::AdditionalData + additional_data; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points | + update_values); + additional_data.mapping_update_flags_inner_faces = + (update_JxW_values | update_quadrature_points | update_normal_vectors | + update_values); + additional_data.mapping_update_flags_boundary_faces = + (update_JxW_values | update_quadrature_points | update_normal_vectors | + update_values); + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + MatrixFreeTools::categorize_by_boundary_ids(dof_handler.get_triangulation(), + additional_data); + + additional_data.communicator_sm = subcommunicator; + + data.reinit( + mapping, dof_handlers, constraints, quadratures, additional_data); + } + + + + template + void + NavierStokesOperator::perform_stage( + const unsigned int stage, + const Number current_time, + const Number bi, + const Number ai, + const LinearAlgebra::distributed::Vector ¤t_ri, + LinearAlgebra::distributed::Vector & vec_ki, + LinearAlgebra::distributed::Vector & solution) const + { + for (auto &i : inflow_boundaries) + i.second->set_time(current_time); + for (auto &i : subsonic_outflow_boundaries) + i.second->set_time(current_time); + + this->ai = ai; + this->bi = bi; + this->solution = &solution; + this->stage = stage; + + data.loop_cell_centric( + &NavierStokesOperator::operation_on_cell, + this, + vec_ki, + current_ri, + true, + MatrixFree::DataAccessOnFaces::values); + } + + + + template + void + NavierStokesOperator::operation_on_cell( + const MatrixFree &data, + LinearAlgebra::distributed::Vector & vec_ki, + const LinearAlgebra::distributed::Vector & current_ri, + const std::pair & cell_range) const + { + using FECellIntegral = FEEvaluation; + using FEFaceIntegral = FEFaceEvaluation; + + FECellIntegral phi(data); + FECellIntegral phi_temp(data); + FEFaceIntegral phi_m(data, true); + FEFaceIntegral phi_p(data, false); + + Tensor<1, dim, VectorizedArrayType> constant_body_force; + const Functions::ConstantFunction *constant_function = + dynamic_cast *>(body_force.get()); + + if (constant_function) + constant_body_force = evaluate_function( + *constant_function, Point()); + + const dealii::internal::EvaluatorTensorProduct< + dealii::internal::EvaluatorVariant::evaluate_evenodd, + dim, + n_points_1d, + n_points_1d, + VectorizedArrayType> + eval(AlignedVector(), + data.get_shape_info().data[0].shape_gradients_collocation_eo, + AlignedVector()); + + internal::EvaluatorTensorProduct< + internal::EvaluatorVariant::evaluate_evenodd, + dim - 1, + n_q_points_1d, + n_q_points_1d, + VectorizedArrayType> + eval_face({}, + data.get_shape_info().data[0].shape_gradients_collocation_eo, + {}); + + AlignedVector buffer; + buffer.resize_fast(phi.static_n_q_points * phi.n_components); + AlignedVector buffer_face; + buffer_face.resize_fast(phi_m.static_n_q_points * 2); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + + if (ai != Number()) + phi_temp.reinit(cell); + + if (ai != Number() && stage == 0) + { + phi.read_dof_values(current_ri); + + for (unsigned int i = 0; + i < phi.static_dofs_per_component * (dim + 2); + ++i) + phi_temp.begin_dof_values()[i] = phi.begin_dof_values()[i]; + + phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + } + else + { + phi.gather_evaluate(current_ri, + EvaluationFlags::values | + EvaluationFlags::gradients); + } + + for (unsigned int i = 0; i < phi.static_n_q_points * (dim + 2); ++i) + buffer[i] = phi.begin_values()[i]; + + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + const auto w_q = phi.get_value(q); + const auto grad_w_q = phi.get_gradient(q); + auto flux = euler_flux(w_q); + const auto viscous = viscous_flux(w_q, grad_w_q); + for (unsigned int d = 0; d < dim + 2; ++d) + flux[d] = flux[d] - viscous[d]; + phi.submit_gradient(flux, q); + if (body_force.get() != nullptr) + { + const Tensor<1, dim, VectorizedArrayType> force = + constant_function ? + constant_body_force : + evaluate_function( + *body_force, phi.quadrature_point(q)); + + Tensor<1, dim + 2, VectorizedArrayType> forcing; + for (unsigned int d = 0; d < dim; ++d) + forcing[d + 1] = w_q[0] * force[d]; + for (unsigned int d = 0; d < dim; ++d) + forcing[dim + 1] += force[d] * w_q[d + 1]; + + phi.submit_value(forcing, q); + } + } + + { + auto *values_ptr = phi.begin_values(); + auto *gradient_ptr = phi.begin_gradients(); + + for (unsigned int c = 0; c < dim + 2; ++c) + { + if (dim >= 1 && body_force.get() == nullptr) + eval.template gradients<0, false, false>( + gradient_ptr + phi.static_n_q_points * 0, values_ptr); + else if (dim >= 1) + eval.template gradients<0, false, true>( + gradient_ptr + phi.static_n_q_points * 0, values_ptr); + if (dim >= 2) + eval.template gradients<1, false, true>( + gradient_ptr + phi.static_n_q_points * 1, values_ptr); + if (dim >= 3) + eval.template gradients<2, false, true>( + gradient_ptr + phi.static_n_q_points * 2, values_ptr); + + values_ptr += phi.static_n_q_points; + gradient_ptr += phi.static_n_q_points * dim; + } + } + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + { + const auto boundary_ids = + data.get_faces_by_cells_boundary_id(cell, face); + + Assert(std::equal(boundary_ids.begin(), + boundary_ids.begin() + + data.n_active_entries_per_cell_batch(cell), + boundary_ids.begin()), + ExcMessage("Boundary IDs of lanes differ.")); + + const auto boundary_id = boundary_ids[0]; + + phi_m.reinit(cell, face); + + internal::EvaluatorTensorProduct< + internal::EvaluatorVariant::evaluate_general, + dim, + n_points_1d, + 0, + VectorizedArrayType> + evalf(data.get_shape_info() + .data.front() + .quadrature_data_on_face[face % 2], + {}, + {}, + n_points_1d, + 0); + + for (unsigned int d = 0; d < dim + 2; ++d) + { + const unsigned int n_q_points_face = phi_m.static_n_q_points; + if (face / 2 == 0) + evalf.template apply_face<0, true, false, 1>( + buffer.data() + d * phi.static_n_q_points, + buffer_face.data()); + else if (face / 2 == 1) + evalf.template apply_face<1, true, false, 1>( + buffer.data() + d * phi.static_n_q_points, + buffer_face.data()); + else if (face / 2 == 2) + evalf.template apply_face<2, true, false, 1>( + buffer.data() + d * phi.static_n_q_points, + buffer_face.data()); + + if (dim > 1) + eval_face.template gradients<0, true, false>( + buffer_face.data(), + phi_m.begin_gradients() + (d * dim) * n_q_points_face); + if (dim > 2) + eval_face.template gradients<1, true, false>( + buffer_face.data(), + phi_m.begin_gradients() + (d * dim + 1) * n_q_points_face); + + for (unsigned int i = 0; i < n_q_points_face; ++i) + { + phi_m.begin_values()[d * n_q_points_face + i] = + buffer_face[i]; + + phi_m + .begin_gradients()[(d * dim + dim - 1) * n_q_points_face + + i] = buffer_face[n_q_points_face + i]; + } + } + + if (boundary_id == numbers::internal_face_boundary_id) + { + phi_p.reinit(cell, face); + phi_p.gather_evaluate(current_ri, + EvaluationFlags::values | + EvaluationFlags::gradients); + + const auto tau_ip = + (std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs((phi_p.get_normal_vector(0) * + phi_p.inverse_jacobian(0))[dim - 1])) * + Number(viscosity * (degree + 1) * (degree + 1)); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + const auto w_m = phi_m.get_value(q); + const auto w_p = phi_p.get_value(q); + const auto normal = phi_m.get_normal_vector(q); + auto numerical_flux = + -euler_numerical_flux(w_m, w_p, normal); + const auto grad_w_m = phi_m.get_gradient(q); + const auto grad_w_p = phi_p.get_gradient(q); + + const auto flux_q1 = viscous_flux(w_m, grad_w_m); + for (unsigned int d = 0; d < dim + 2; ++d) + numerical_flux[d] += 0.5 * (flux_q1[d] * normal); + const auto flux_q2 = viscous_flux(w_p, grad_w_p); + for (unsigned int d = 0; d < dim + 2; ++d) + numerical_flux[d] += 0.5 * (flux_q2[d] * normal); + numerical_flux -= tau_ip * (w_m - w_p); + phi_m.submit_value(numerical_flux, q); + + Tensor<1, dim + 2, Tensor<1, dim, VectorizedArrayType>> + w_jump; + for (unsigned int d = 0; d < dim + 2; ++d) + for (unsigned int e = 0; e < dim; ++e) + w_jump[d][e] = + (w_m[d] - w_p[d]) * (Number(0.5) * normal[e]); + phi_m.submit_gradient(viscous_flux(w_m, w_jump), q); + } + } + else + { + const auto tau_ip = + std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) * + Number(2. * viscosity * (degree + 1) * (degree + 1)); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + const auto w_m = phi_m.get_value(q); + const auto normal = phi_m.get_normal_vector(q); + const auto grad_w_m = phi_m.get_gradient(q); + const auto grad_w_p = grad_w_m; + + auto rho_u_dot_n = w_m[1] * normal[0]; + for (unsigned int d = 1; d < dim; ++d) + rho_u_dot_n += w_m[1 + d] * normal[d]; + + bool at_outflow = false; + + Tensor<1, dim + 2, VectorizedArrayType> w_p; + + if (wall_boundaries.find(boundary_id) != + wall_boundaries.end()) + { + w_p[0] = w_m[0]; + for (unsigned int d = 0; d < dim; ++d) + w_p[d + 1] = + w_m[d + 1] - 2. * rho_u_dot_n * normal[d]; + w_p[dim + 1] = w_m[dim + 1]; + } + else if (inflow_boundaries.find(boundary_id) != + inflow_boundaries.end()) + w_p = evaluate_function( + *inflow_boundaries.find(boundary_id)->second, + phi_m.quadrature_point(q)); + else if (subsonic_outflow_boundaries.find(boundary_id) != + subsonic_outflow_boundaries.end()) + { + w_p = w_m; + w_p[dim + 1] = + evaluate_function(*subsonic_outflow_boundaries + .find(boundary_id) + ->second, + phi_m.quadrature_point(q), + dim + 1); + at_outflow = true; + } + else + AssertThrow(false, + ExcMessage( + "Unknown boundary id, did " + "you set a boundary condition for " + "this part of the domain boundary?")); + + auto flux = -euler_numerical_flux(w_m, w_p, normal); + + if (at_outflow) + for (unsigned int v = 0; v < VectorizedArrayType::size(); + ++v) + { + if (rho_u_dot_n[v] < -1e-12) + for (unsigned int d = 0; d < dim; ++d) + flux[d + 1][v] = 0.; + } + + const auto flux_q1 = viscous_flux(w_m, grad_w_m); + for (unsigned int d = 0; d < dim + 2; ++d) + flux[d] += 0.5 * (flux_q1[d] * normal); + const auto flux_q2 = viscous_flux(w_p, grad_w_p); + for (unsigned int d = 0; d < dim + 2; ++d) + flux[d] += 0.5 * (flux_q2[d] * normal); + flux -= tau_ip * (w_m - w_p); + phi_m.submit_value(flux, q); + + Tensor<1, dim + 2, Tensor<1, dim, VectorizedArrayType>> + w_jump; + for (unsigned int d = 0; d < dim + 2; ++d) + for (unsigned int e = 0; e < dim; ++e) + w_jump[d][e] = + (w_m[d] - w_p[d]) * (Number(0.5) * normal[e]); + phi_m.submit_gradient(viscous_flux(w_m, w_jump), q); + } + } + + for (unsigned int d = 0; d < dim + 2; ++d) + { + const unsigned int n_q_points_face = phi_m.static_n_q_points; + for (unsigned int i = 0; i < n_q_points_face; ++i) + { + buffer_face[i] = + phi_m.begin_values()[d * n_q_points_face + i]; + buffer_face[n_q_points_face + i] = + phi_m.begin_gradients()[(d * dim + dim - 1) * + n_q_points_face + + i]; + } + + if (dim > 2) + eval_face.template gradients<1, false, true>( + phi_m.begin_gradients() + (d * dim + 1) * n_q_points_face, + buffer_face.data()); + if (dim > 1) + eval_face.template gradients<0, false, true>( + phi_m.begin_gradients() + (d * dim) * n_q_points_face, + buffer_face.data()); + + if (face / 2 == 0) + evalf.template apply_face<0, false, true, 1>( + buffer_face.data(), + phi.begin_values() + d * phi.static_n_q_points); + else if (face / 2 == 1) + evalf.template apply_face<1, false, true, 1>( + buffer_face.data(), + phi.begin_values() + d * phi.static_n_q_points); + else if (face / 2 == 2) + evalf.template apply_face<2, false, true, 1>( + buffer_face.data(), + phi.begin_values() + d * phi.static_n_q_points); + } + } + + for (unsigned int q = 0; q < phi.static_n_q_points; ++q) + { + const auto factor = VectorizedArrayType(1.0) / phi.JxW(q); + for (unsigned int c = 0; c < dim + 2; ++c) + phi.begin_values()[c * phi.static_n_q_points + q] = + phi.begin_values()[c * phi.static_n_q_points + q] * factor; + } + + internal::FEEvaluationImplBasisChange< + dealii::internal::EvaluatorVariant::evaluate_evenodd, + internal::EvaluatorQuantity::hessian, + dim, + degree + 1, + n_points_1d, + VectorizedArrayType, + VectorizedArrayType>::do_backward(dim + 2, + data.get_shape_info() + .data[0] + .inverse_shape_values_eo, + false, + phi.begin_values(), + phi.begin_dof_values()); + + if (ai == Number()) + { + for (unsigned int q = 0; q < phi.static_dofs_per_cell; ++q) + phi.begin_dof_values()[q] = bi * phi.begin_dof_values()[q]; + phi.distribute_local_to_global(*solution); + } + else + { + if (stage != 0) + phi_temp.read_dof_values(*solution); + + for (unsigned int q = 0; q < phi.static_dofs_per_cell; ++q) + { + const auto K_i = phi.begin_dof_values()[q]; + + phi.begin_dof_values()[q] = + phi_temp.begin_dof_values()[q] + (ai * K_i); + + phi_temp.begin_dof_values()[q] += bi * K_i; + } + phi.set_dof_values(vec_ki); + phi_temp.set_dof_values(*solution); + } + } + } + + + + template + void + NavierStokesOperator::perform_stage_face( + const unsigned int stage, + const Number current_time, + const Number bi, + const Number ai, + const LinearAlgebra::distributed::Vector ¤t_ri, + LinearAlgebra::distributed::Vector & vec_ki, + LinearAlgebra::distributed::Vector & solution) const + { + for (auto &i : inflow_boundaries) + i.second->set_time(current_time); + for (auto &i : subsonic_outflow_boundaries) + i.second->set_time(current_time); + + { + Timer timer; + data.loop(&NavierStokesOperator::operation_cell, + &NavierStokesOperator::operation_face, + &NavierStokesOperator::operation_boundary, + this, + vec_ki, + current_ri, + true, + MatrixFree:: + DataAccessOnFaces::gradients, + MatrixFree:: + DataAccessOnFaces::gradients); + time_loop += timer.wall_time(); + } + + { + Timer timer; + data.cell_loop( + &NavierStokesOperator::local_apply_inverse_mass_matrix, + this, + vec_ki, + vec_ki, + std::function(), + [&](const unsigned int start_range, const unsigned int end_range) { + if (ai == Number()) + { + /* DEAL_II_OPENMP_SIMD_PRAGMA */ + for (unsigned int i = start_range; i < end_range; ++i) + { + const Number k_i = vec_ki.local_element(i); + const Number sol_i = solution.local_element(i); + solution.local_element(i) = sol_i + bi * k_i; + } + } + else + { + /* DEAL_II_OPENMP_SIMD_PRAGMA */ + if (stage == 0) + for (unsigned int i = start_range; i < end_range; ++i) + { + const Number k_i = vec_ki.local_element(i); + const Number sol_i = current_ri.local_element(i); + solution.local_element(i) = sol_i + bi * k_i; + vec_ki.local_element(i) = sol_i + ai * k_i; + } + else + for (unsigned int i = start_range; i < end_range; ++i) + { + const Number k_i = vec_ki.local_element(i); + const Number sol_i = solution.local_element(i); + solution.local_element(i) = sol_i + bi * k_i; + vec_ki.local_element(i) = sol_i + ai * k_i; + } + } + }); + time_rk_update += timer.wall_time(); + } + } + + + + template + void + NavierStokesOperator::operation_cell( + const MatrixFree &data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector & src, + const std::pair & cell_range) const + { + using FECellIntegral = FEEvaluation; + + FECellIntegral phi(data); + + Tensor<1, dim, VectorizedArrayType> constant_body_force; + const Functions::ConstantFunction *constant_function = + dynamic_cast *>(body_force.get()); + + if (constant_function) + constant_body_force = evaluate_function( + *constant_function, Point()); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + const auto w_q = phi.get_value(q); + const auto grad_w_q = phi.get_gradient(q); + auto flux = euler_flux(w_q); + const auto viscous = viscous_flux(w_q, grad_w_q); + for (unsigned int d = 0; d < dim + 2; ++d) + flux[d] = flux[d] - viscous[d]; + phi.submit_gradient(flux, q); + if (body_force.get() != nullptr) + { + const Tensor<1, dim, VectorizedArrayType> force = + constant_function ? + constant_body_force : + evaluate_function( + *body_force, phi.quadrature_point(q)); + + Tensor<1, dim + 2, VectorizedArrayType> forcing; + for (unsigned int d = 0; d < dim; ++d) + forcing[d + 1] = w_q[0] * force[d]; + for (unsigned int d = 0; d < dim; ++d) + forcing[dim + 1] += force[d] * w_q[d + 1]; + + phi.submit_value(forcing, q); + } + } + + phi.integrate_scatter(((body_force.get() != nullptr) ? + EvaluationFlags::values : + EvaluationFlags::nothing) | + EvaluationFlags::gradients, + dst); + } + } + + + + template + void + NavierStokesOperator::operation_face( + const MatrixFree &data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector & src, + const std::pair & face_range) const + { + using FEFaceIntegral = FEFaceEvaluation; + FEFaceIntegral phi_m(data, true); + FEFaceIntegral phi_p(data, false); + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + phi_p.reinit(face); + phi_p.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + + phi_m.reinit(face); + phi_m.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + + const auto tau_ip = (std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs((phi_p.get_normal_vector(0) * + phi_p.inverse_jacobian(0))[dim - 1])) * + Number(viscosity * (degree + 1) * (degree + 1)); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + const auto w_m = phi_m.get_value(q); + const auto w_p = phi_p.get_value(q); + const auto normal = phi_m.get_normal_vector(q); + auto numerical_flux = -euler_numerical_flux(w_m, w_p, normal); + const auto grad_w_m = phi_m.get_gradient(q); + const auto grad_w_p = phi_p.get_gradient(q); + + const auto flux_q1 = viscous_flux(w_m, grad_w_m); + for (unsigned int d = 0; d < dim + 2; ++d) + numerical_flux[d] += 0.5 * (flux_q1[d] * normal); + const auto flux_q2 = viscous_flux(w_p, grad_w_p); + for (unsigned int d = 0; d < dim + 2; ++d) + numerical_flux[d] += 0.5 * (flux_q2[d] * normal); + numerical_flux -= tau_ip * (w_m - w_p); + phi_m.submit_value(numerical_flux, q); + phi_p.submit_value(-numerical_flux, q); + + Tensor<1, dim + 2, Tensor<1, dim, VectorizedArrayType>> w_jump; + for (unsigned int d = 0; d < dim + 2; ++d) + for (unsigned int e = 0; e < dim; ++e) + w_jump[d][e] = (w_m[d] - w_p[d]) * (Number(0.5) * normal[e]); + phi_m.submit_gradient(viscous_flux(w_m, w_jump), q); + phi_p.submit_gradient(viscous_flux(w_p, w_jump), q); + } + + phi_m.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + phi_p.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + } + + + + template + void + NavierStokesOperator::operation_boundary( + const MatrixFree &data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector & src, + const std::pair & face_range) const + { + AssertThrow(false, ExcNotImplemented()); + FEFaceEvaluation + phi_m(data, true); + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + phi_m.reinit(face); + phi_m.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + + const auto tau_ip = + std::abs( + (phi_m.get_normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) * + Number(2. * viscosity * (degree + 1) * (degree + 1)); + + const auto boundary_id = data.get_boundary_id(face); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + const auto w_m = phi_m.get_value(q); + const auto normal = phi_m.get_normal_vector(q); + const auto grad_w_m = phi_m.get_gradient(q); + const auto grad_w_p = grad_w_m; + + auto rho_u_dot_n = w_m[1] * normal[0]; + for (unsigned int d = 1; d < dim; ++d) + rho_u_dot_n += w_m[1 + d] * normal[d]; + + bool at_outflow = false; + + Tensor<1, dim + 2, VectorizedArrayType> w_p; + + if (wall_boundaries.find(boundary_id) != wall_boundaries.end()) + { + w_p[0] = w_m[0]; + for (unsigned int d = 0; d < dim; ++d) + w_p[d + 1] = w_m[d + 1] - 2. * rho_u_dot_n * normal[d]; + w_p[dim + 1] = w_m[dim + 1]; + } + else if (inflow_boundaries.find(boundary_id) != + inflow_boundaries.end()) + w_p = + evaluate_function(*inflow_boundaries.find(boundary_id)->second, + phi_m.quadrature_point(q)); + else if (subsonic_outflow_boundaries.find(boundary_id) != + subsonic_outflow_boundaries.end()) + { + w_p = w_m; + w_p[dim + 1] = evaluate_function( + *subsonic_outflow_boundaries.find(boundary_id)->second, + phi_m.quadrature_point(q), + dim + 1); + at_outflow = true; + } + else + AssertThrow(false, + ExcMessage("Unknown boundary id, did " + "you set a boundary condition for " + "this part of the domain boundary?")); + + auto flux = -euler_numerical_flux(w_m, w_p, normal); + + if (at_outflow) + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + { + if (rho_u_dot_n[v] < -1e-12) + for (unsigned int d = 0; d < dim; ++d) + flux[d + 1][v] = 0.; + } + + const auto flux_q1 = viscous_flux(w_m, grad_w_m); + for (unsigned int d = 0; d < dim + 2; ++d) + flux[d] += 0.5 * (flux_q1[d] * normal); + const auto flux_q2 = viscous_flux(w_p, grad_w_p); + for (unsigned int d = 0; d < dim + 2; ++d) + flux[d] += 0.5 * (flux_q2[d] * normal); + flux -= tau_ip * (w_m - w_p); + phi_m.submit_value(flux, q); + + Tensor<1, dim + 2, Tensor<1, dim, VectorizedArrayType>> w_jump; + for (unsigned int d = 0; d < dim + 2; ++d) + for (unsigned int e = 0; e < dim; ++e) + w_jump[d][e] = (w_m[d] - w_p[d]) * (Number(0.5) * normal[e]); + + phi_m.submit_gradient(viscous_flux(w_m, w_jump), q); + } + + phi_m.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + } + + + + template + void + NavierStokesOperator:: + local_apply_inverse_mass_matrix( + const MatrixFree &, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const + { + FEEvaluation phi(data, 0, 1); + MatrixFreeOperators::CellwiseInverseMassMatrix + inverse(phi); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + + inverse.apply(phi.begin_dof_values(), phi.begin_dof_values()); + + phi.set_dof_values(dst); + } + } + + + + template + void + NavierStokesOperator::initialize_vector( + LinearAlgebra::distributed::Vector &vector) const + { + data.initialize_dof_vector(vector); + } + + + + template + void + NavierStokesOperator::set_inflow_boundary( + const types::boundary_id boundary_id, + std::unique_ptr> inflow_function) + { + AssertThrow(subsonic_outflow_boundaries.find(boundary_id) == + subsonic_outflow_boundaries.end() && + wall_boundaries.find(boundary_id) == wall_boundaries.end(), + ExcMessage("You already set the boundary with id " + + std::to_string(static_cast(boundary_id)) + + " to another type of boundary before now setting " + + "it as inflow")); + AssertThrow(inflow_function->n_components == dim + 2, + ExcMessage("Expected function with dim+2 components")); + + inflow_boundaries[boundary_id] = std::move(inflow_function); + } + + + + template + void + NavierStokesOperator::set_subsonic_outflow_boundary( + const types::boundary_id boundary_id, + std::unique_ptr> outflow_function) + { + AssertThrow(inflow_boundaries.find(boundary_id) == + inflow_boundaries.end() && + wall_boundaries.find(boundary_id) == wall_boundaries.end(), + ExcMessage("You already set the boundary with id " + + std::to_string(static_cast(boundary_id)) + + " to another type of boundary before now setting " + + "it as subsonic outflow")); + AssertThrow(outflow_function->n_components == dim + 2, + ExcMessage("Expected function with dim+2 components")); + + subsonic_outflow_boundaries[boundary_id] = std::move(outflow_function); + } + + + + template + void + NavierStokesOperator::set_wall_boundary( + const types::boundary_id boundary_id) + { + AssertThrow(inflow_boundaries.find(boundary_id) == + inflow_boundaries.end() && + subsonic_outflow_boundaries.find(boundary_id) == + subsonic_outflow_boundaries.end(), + ExcMessage("You already set the boundary with id " + + std::to_string(static_cast(boundary_id)) + + " to another type of boundary before now setting " + + "it as wall boundary")); + + wall_boundaries.insert(boundary_id); + } + + + + template + void + NavierStokesOperator::set_body_force( + std::unique_ptr> body_force) + { + AssertDimension(body_force->n_components, dim); + + this->body_force = std::move(body_force); + } + + + + template + void + NavierStokesOperator::project( + const Function & function, + LinearAlgebra::distributed::Vector &solution) const + { + FEEvaluation + phi(data, 0, 1); + MatrixFreeOperators::CellwiseInverseMassMatrix + inverse(phi); + solution.zero_out_ghost_values(); + for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell) + { + phi.reinit(cell); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_dof_value(evaluate_function(function, + phi.quadrature_point(q)), + q); + inverse.transform_from_q_points_to_basis(dim + 2, + phi.begin_dof_values(), + phi.begin_dof_values()); + phi.set_dof_values(solution); + } + } + + + + template + std::array + NavierStokesOperator::compute_errors( + const Function & function, + const LinearAlgebra::distributed::Vector &solution) const + { + double errors_squared[3] = {}; + FEEvaluation + phi(data, 0, 0); + + for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(solution, EvaluationFlags::values); + VectorizedArrayType local_errors_squared[3] = {}; + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + const auto error = + evaluate_function(function, phi.quadrature_point(q)) - + phi.get_value(q); + const auto JxW = phi.JxW(q); + + local_errors_squared[0] += error[0] * error[0] * JxW; + for (unsigned int d = 0; d < dim; ++d) + local_errors_squared[1] += (error[d + 1] * error[d + 1]) * JxW; + local_errors_squared[2] += (error[dim + 1] * error[dim + 1]) * JxW; + } + for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell); + ++v) + for (unsigned int d = 0; d < 3; ++d) + errors_squared[d] += local_errors_squared[d][v]; + } + + Utilities::MPI::sum(errors_squared, MPI_COMM_WORLD, errors_squared); + + std::array errors; + for (unsigned int d = 0; d < 3; ++d) + errors[d] = std::sqrt(errors_squared[d]); + + return errors; + } + + + + template + std::array + NavierStokesOperator::compute_kinetic_energy( + const LinearAlgebra::distributed::Vector &solution) const + { + double squared[2] = {}; + FEEvaluation + phi(data, 0, 0); + + for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(solution, + EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArrayType local_squared[2] = {}; + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + const auto JxW = phi.JxW(q); + const auto w_q = phi.get_value(q); + const auto velocity = fluid_velocity(w_q); + const auto velocity_grad = + fluid_velocity_gradient(w_q, phi.get_gradient(q)); + local_squared[0] += velocity.norm_square() * JxW; + local_squared[1] += + scalar_product(velocity_grad, velocity_grad) * JxW; + } + for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell); + ++v) + for (unsigned int d = 0; d < 2; ++d) + squared[d] += local_squared[d][v]; + } + + Utilities::MPI::sum(squared, MPI_COMM_WORLD, squared); + + std::array result{ + {0.5 * squared[0] / Utilities::fixed_power(2. * numbers::PI), + viscosity * squared[1] / Utilities::fixed_power(2. * numbers::PI)}}; + + return result; + } + + + + template + double + NavierStokesOperator::compute_cell_transport_speed( + const LinearAlgebra::distributed::Vector &solution) const + { + Number max_transport = 0; + FEEvaluation + phi(data, 0, 1); + + for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(solution, EvaluationFlags::values); + VectorizedArrayType local_max = 0.; + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + const auto solution = phi.get_value(q); + const auto velocity = fluid_velocity(solution); + const auto pressure = fluid_pressure(solution); + + const auto inverse_jacobian = phi.inverse_jacobian(q); + const auto convective_speed = inverse_jacobian * velocity; + VectorizedArrayType convective_limit = 0.; + for (unsigned int d = 0; d < dim; ++d) + convective_limit = + std::max(convective_limit, std::abs(convective_speed[d])); + + const auto speed_of_sound = + std::sqrt(gamma * pressure * (1. / solution[0])); + + Tensor<1, dim, VectorizedArrayType> eigenvector; + for (unsigned int d = 0; d < dim; ++d) + eigenvector[d] = 1.; + for (unsigned int i = 0; i < 5; ++i) + { + eigenvector = transpose(inverse_jacobian) * + (inverse_jacobian * eigenvector); + VectorizedArrayType eigenvector_norm = 0.; + for (unsigned int d = 0; d < dim; ++d) + eigenvector_norm = + std::max(eigenvector_norm, std::abs(eigenvector[d])); + eigenvector /= eigenvector_norm; + } + const auto jac_times_ev = inverse_jacobian * eigenvector; + const auto max_eigenvalue = std::sqrt( + (jac_times_ev * jac_times_ev) / (eigenvector * eigenvector)); + local_max = + std::max(local_max, + max_eigenvalue * speed_of_sound + convective_limit); + } + + for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell); + ++v) + for (unsigned int d = 0; d < 3; ++d) + max_transport = std::max(max_transport, local_max[v]); + } + + max_transport = Utilities::MPI::max(max_transport, MPI_COMM_WORLD); + + return max_transport; + } + + + + template + class NavierStokesOperatorFaceCentric + { + public: + NavierStokesOperatorFaceCentric( + const NavierStokesOperator &ns_operator) + : ns_operator(ns_operator) + {} + + void + perform_stage(const unsigned int stage, + const Number current_time, + const Number bi, + const Number ai, + const LinearAlgebra::distributed::Vector ¤t_ri, + LinearAlgebra::distributed::Vector & vec_ki, + LinearAlgebra::distributed::Vector &solution) const + { + ns_operator.perform_stage_face( + stage, current_time, bi, ai, current_ri, vec_ki, solution); + } + + private: + const NavierStokesOperator &ns_operator; + }; + + + + template + class FlowProblem + { + public: + FlowProblem(); + + Measurement + run(); + + private: + void + make_grid(); + + void + output_results(const unsigned int result_number); + + LinearAlgebra::distributed::Vector solution; + +#ifdef DEAL_II_WITH_P4EST + parallel::distributed::Triangulation triangulation; +#else + Triangulation triangulation; +#endif + + FESystem fe; + MappingQ mapping; + DoFHandler dof_handler; + + NavierStokesOperator flow_operator; + + double time, time_step; + + class Postprocessor : public DataPostprocessor + { + public: + Postprocessor(); + + virtual void + evaluate_vector_field( + const DataPostprocessorInputs::Vector &inputs, + std::vector> &computed_quantities) const override; + + virtual std::vector + get_names() const override; + + virtual std::vector< + DataComponentInterpretation::DataComponentInterpretation> + get_data_component_interpretation() const override; + + virtual UpdateFlags + get_needed_update_flags() const override; + + private: + const bool do_schlieren_plot; + }; + }; + + + + template + FlowProblem::Postprocessor::Postprocessor() + : do_schlieren_plot(dim == 2) + {} + + + + template + void + FlowProblem::Postprocessor::evaluate_vector_field( + const DataPostprocessorInputs::Vector &inputs, + std::vector> & computed_quantities) const + { + const unsigned int n_evaluation_points = inputs.solution_values.size(); + + if (do_schlieren_plot == true) + Assert(inputs.solution_gradients.size() == n_evaluation_points, + ExcInternalError()); + + Assert(computed_quantities.size() == n_evaluation_points, + ExcInternalError()); + Assert(inputs.solution_values[0].size() == dim + 2, ExcInternalError()); + Assert(computed_quantities[0].size() == + dim + 2 + (do_schlieren_plot == true ? 1 : 0), + ExcInternalError()); + + for (unsigned int q = 0; q < n_evaluation_points; ++q) + { + Tensor<1, dim + 2> solution; + for (unsigned int d = 0; d < dim + 2; ++d) + solution[d] = inputs.solution_values[q](d); + + const double density = solution[0]; + const Tensor<1, dim> velocity = fluid_velocity(solution); + const double pressure = fluid_pressure(solution); + + for (unsigned int d = 0; d < dim; ++d) + computed_quantities[q](d) = velocity[d]; + computed_quantities[q](dim) = pressure; + computed_quantities[q](dim + 1) = std::sqrt(gamma * pressure / density); + + if (do_schlieren_plot == true) + computed_quantities[q](dim + 2) = + inputs.solution_gradients[q][0] * inputs.solution_gradients[q][0]; + } + } + + + + template + std::vector + FlowProblem::Postprocessor::get_names() const + { + std::vector names; + for (unsigned int d = 0; d < dim; ++d) + names.emplace_back("velocity"); + names.emplace_back("pressure"); + names.emplace_back("speed_of_sound"); + + if (do_schlieren_plot == true) + names.emplace_back("schlieren_plot"); + + return names; + } + + + + template + std::vector + FlowProblem::Postprocessor::get_data_component_interpretation() const + { + std::vector + interpretation; + for (unsigned int d = 0; d < dim; ++d) + interpretation.push_back( + DataComponentInterpretation::component_is_part_of_vector); + interpretation.push_back(DataComponentInterpretation::component_is_scalar); + interpretation.push_back(DataComponentInterpretation::component_is_scalar); + + if (do_schlieren_plot == true) + interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + + return interpretation; + } + + + + template + UpdateFlags + FlowProblem::Postprocessor::get_needed_update_flags() const + { + if (do_schlieren_plot == true) + return update_values | update_gradients; + else + return update_values; + } + + + + template + FlowProblem::FlowProblem() + : +#ifdef DEAL_II_WITH_P4EST + triangulation(MPI_COMM_WORLD) +#endif + , fe(FE_DGQHermite(fe_degree), dim + 2) + , mapping(fe_degree) + , dof_handler(triangulation) + , flow_operator() + , time(0) + , time_step(0) + {} + + + + template + void + FlowProblem::make_grid() + { + switch (testcase) + { + case 0: + { + Point lower_left; + for (unsigned int d = 1; d < dim; ++d) + lower_left[d] = -5; + + Point upper_right; + upper_right[0] = 10; + for (unsigned int d = 1; d < dim; ++d) + upper_right[d] = 5; + + GridGenerator::hyper_rectangle(triangulation, + lower_left, + upper_right); + triangulation.refine_global(2); + + flow_operator.set_inflow_boundary( + 0, std::make_unique>(0)); + + break; + } + + case 1: + { + GridGenerator::channel_with_cylinder( + triangulation, 0.03, 1, 0, true); + + flow_operator.set_inflow_boundary( + 0, std::make_unique>(0)); + flow_operator.set_subsonic_outflow_boundary( + 1, std::make_unique>(0)); + + flow_operator.set_wall_boundary(2); + flow_operator.set_wall_boundary(3); + + if (dim == 3) + flow_operator.set_body_force( + std::make_unique>( + std::vector({0., 0., -0.2}))); + + break; + } + + case 2: + { + Point lower_left, upper_right; + for (unsigned int d = 0; d < dim; ++d) + lower_left[d] = -numbers::PI; + + for (unsigned int d = 0; d < dim; ++d) + upper_right[d] = numbers::PI; + + GridGenerator::hyper_rectangle(triangulation, + lower_left, + upper_right); + for (const auto &cell : triangulation.cell_iterators()) + for (unsigned int face : cell->face_indices()) + if (cell->at_boundary(face)) + cell->face(face)->set_boundary_id(face); + std::vector::cell_iterator>> + periodic_faces; + for (unsigned int d = 0; d < dim; ++d) + GridTools::collect_periodic_faces( + triangulation, 2 * d, 2 * d + 1, d, periodic_faces); + triangulation.add_periodicity(periodic_faces); + + triangulation.refine_global(2); + + break; + } + + default: + Assert(false, ExcNotImplemented()); + } + + switch (get_testing_environment()) + { + case TestingEnvironment::light: + triangulation.refine_global(1); + break; + case TestingEnvironment::medium: + triangulation.refine_global(2); + break; + case TestingEnvironment::heavy: + triangulation.refine_global(3); + break; + } + } + + + + template + void + FlowProblem::output_results(const unsigned int) + { + Postprocessor postprocessor; + DataOut data_out; + + DataOutBase::VtkFlags flags; + flags.write_higher_order_cells = true; + data_out.set_flags(flags); + + data_out.attach_dof_handler(dof_handler); + { + std::vector names; + names.emplace_back("density"); + for (unsigned int d = 0; d < dim; ++d) + names.emplace_back("momentum"); + names.emplace_back("energy"); + + std::vector + interpretation; + interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + for (unsigned int d = 0; d < dim; ++d) + interpretation.push_back( + DataComponentInterpretation::component_is_part_of_vector); + interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + + data_out.add_data_vector(dof_handler, solution, names, interpretation); + } + data_out.add_data_vector(solution, postprocessor); + + LinearAlgebra::distributed::Vector reference; + if (testcase == 0 && dim == 2) + { + reference.reinit(solution); + flow_operator.project(ExactSolution(time), reference); + reference.sadd(-1., 1, solution); + std::vector names; + names.emplace_back("error_density"); + for (unsigned int d = 0; d < dim; ++d) + names.emplace_back("error_momentum"); + names.emplace_back("error_energy"); + + std::vector + interpretation; + interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + for (unsigned int d = 0; d < dim; ++d) + interpretation.push_back( + DataComponentInterpretation::component_is_part_of_vector); + interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + + data_out.add_data_vector(dof_handler, reference, names, interpretation); + } + + Vector mpi_owner(triangulation.n_active_cells()); + mpi_owner = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + data_out.add_data_vector(mpi_owner, "owner"); + + data_out.build_patches(mapping, + fe.degree, + DataOut::curved_inner_cells); + + // Do not write a file here to be independent of file system + } + + + + template + Measurement + FlowProblem::run() + { + std::map timer; + + timer["setup_grid"].start(); + make_grid(); + timer["setup_grid"].stop(); + + timer["setup_dofs"].start(); + dof_handler.distribute_dofs(fe); + timer["setup_dofs"].stop(); + + timer["setup_matrix_free"].start(); + flow_operator.reinit(mapping, dof_handler); + flow_operator.initialize_vector(solution); + LinearAlgebra::distributed::Vector rk_register_1; + LinearAlgebra::distributed::Vector rk_register_2; + rk_register_1.reinit(solution); + rk_register_2.reinit(solution); + timer["setup_matrix_free"].stop(); + + const LowStorageRungeKuttaIntegrator integrator(lsrk_scheme); + + timer["project_initial"].start(); + flow_operator.project(ExactSolution(time), solution); + + double min_vertex_distance = std::numeric_limits::max(); + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->is_locally_owned()) + min_vertex_distance = + std::min(min_vertex_distance, cell->minimum_vertex_distance()); + min_vertex_distance = + Utilities::MPI::min(min_vertex_distance, MPI_COMM_WORLD); + + time_step = courant_number * integrator.n_stages() / + flow_operator.compute_cell_transport_speed(solution); + timer["project_initial"].stop(); + + time = 0; + + timer["write_output"].start(); + output_results(0); + timer["write_output"].stop(); + + unsigned int timestep_number = 0; + while (timestep_number < 20) + { + timer["rk_timestep_cellbased"].start(); + integrator.perform_time_step(flow_operator, + time, + time_step, + solution, + rk_register_1, + rk_register_2); + timer["rk_timestep_cellbased"].stop(); + + timer["analyze_solution"].start(); + const std::array energy = + flow_operator.compute_kinetic_energy(solution); + AssertThrow(energy[0] > 0 && energy[1] > 0, ExcInternalError()); + timer["analyze_solution"].stop(); + + time += time_step; + ++timestep_number; + } + + while (timestep_number < 40) + { + timer["rk_timestep_facebased"].start(); + NavierStokesOperatorFaceCentric + flow_operator_face(flow_operator); + integrator.perform_time_step(flow_operator_face, + time, + time_step, + solution, + rk_register_1, + rk_register_2); + timer["rk_timestep_facebased"].stop(); + + timer["analyze_solution"].start(); + const std::array energy = + flow_operator.compute_kinetic_energy(solution); + AssertThrow(energy[0] > 0 && energy[1] > 0, ExcInternalError()); + timer["analyze_solution"].stop(); + + time += time_step; + ++timestep_number; + } + + return {timer["setup_grid"].wall_time(), + timer["setup_dofs"].wall_time(), + timer["setup_matrix_free"].wall_time(), + timer["project_initial"].wall_time(), + timer["write_output"].wall_time(), + timer["analyze_solution"].wall_time(), + timer["rk_timestep_cellbased"].wall_time(), + timer["rk_timestep_facebased"].wall_time(), + flow_operator.time_loop, + flow_operator.time_rk_update}; + } + +} // namespace NavierStokes_DG + + + +std::tuple> +describe_measurements() +{ + return {Metric::timing, + 4, + {"setup_grid", + "setup_dofs", + "setup_matrix_free", + "project_initial", + "write_output", + "analyze_solution", + "rk_timestep_cellbased", + "rk_timestep_facebased", + "rk_timestep_facebased_loop", + "rk_timestep_facebased_update"}}; +} + + + +Measurement +perform_single_measurement() +{ + return NavierStokes_DG::FlowProblem().run(); +} diff --git a/tests/performance/timing_navier_stokes.threads=1.mpirun=max.exclusive.release.run_only b/tests/performance/timing_navier_stokes.threads=1.mpirun=max.exclusive.release.run_only new file mode 100644 index 0000000000..e69de29bb2