The Euler Equations and Air Flow Over a Symmetric Airfoil: Solution

Table of Contents

• Python Implementation
• C++ Implementation

Python Implementation

 
import peano4, exahype2
 
domain_size = [120.0, 120.0]
domain_offset = [-10.0, -60.0]
 
min_level = 5  # Increase for finer mesh
 
unknowns = {"p": 1, "v": 2, "E": 1}
auxiliary_variables = {}
max_h = 1.1 * min(domain_size) / (3.0**min_level)
min_h = max_h
 
my_solver = exahype2.solvers.fv.godunov.GlobalAdaptiveTimeStep(
    name="EulerSolver",
    patch_size=22,
    min_volume_h=min_h,
    max_volume_h=max_h,
    time_step_relaxation=0.5,
    unknowns=unknowns,
    auxiliary_variables=auxiliary_variables,
)
 
my_solver.set_implementation(
    initial_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
    boundary_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
    flux=exahype2.solvers.PDETerms.User_Defined_Implementation,
)
 
exahype2_project = exahype2.Project(
    namespace=["tutorials", "exahype2", "euler"],
    directory=".",
    project_name="Airfoil",
    executable="Airfoil",
)
 
exahype2_project.add_solver(my_solver)
exahype2_project.set_output_path("solutions")
 
exahype2_project.set_global_simulation_parameters(
    dimensions=2,
    size=domain_size,
    offset=domain_offset,
    min_end_time=10.0,
    max_end_time=10.0,
    first_plot_time_stamp=0.0,
    time_in_between_plots=0.5,
    periodic_BC=[False, False],
)
 
exahype2_project.set_load_balancer("new ::exahype2::LoadBalancingConfiguration")
exahype2_project.set_Peano4_installation("../../../", peano4.output.CompileMode.Release)
peano4_project = exahype2_project.generate_Peano4_project()
peano4_project.generate()
peano4_project.build(make=True, make_clean_first=False, throw_away_data_after_build=False)

C++ Implementation

//
// A sample implementation file for this problem:
//
#include "EulerSolver.h"
 
tarch::logging::Log tutorials::exahype2::euler::EulerSolver::_log("tutorials::exahype2::euler::EulerSolver");
 
// c chord length, t maximum relative thickness.
double airfoilSymmetric(double x, double c = 100, double t = 0.12) {
  return 5.0 * t * c
         * (0.2969 * sqrt(x / c) + ((((-0.1015) * (x / c) + 0.2843) * (x / c) - 0.3516) * (x / c) - 0.1260) * (x / c));
}
 
void tutorials::exahype2::euler::EulerSolver::initialCondition(
  [[maybe_unused]] double* const NOALIAS                        Q, // Q[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const bool                                   gridIsConstructed
) {
  const double y     = airfoilSymmetric(x[0]);
  const double alpha = (y > x[1] && -y < x[1]) ? 0.0 : 1.0;
 
  Q[0] = (y > x[1] && -y < x[1]) ? 1000.0 : 1.0;
  Q[1] = alpha * 1.0;
  Q[2] = 0.0;
  Q[3] = alpha * 2.5;
}
 
void tutorials::exahype2::euler::EulerSolver::boundaryConditions(
  [[maybe_unused]] const double* const NOALIAS                  Qinside,  // Qinside[4+0]
  [[maybe_unused]] double* const NOALIAS                        Qoutside, // Qoutside[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const double                                 t,
  [[maybe_unused]] const int                                    normal
) {
  if (normal == 0 && x[0] == DomainOffset[0]) {
    Qoutside[0] = 1.0;
    Qoutside[1] = 1.0;
    Qoutside[2] = 0.0;
    Qoutside[3] = 2.5;
  } else {
    Qoutside[0] = Qinside[0];
    Qoutside[1] = Qinside[1];
    Qoutside[2] = Qinside[2];
    Qoutside[3] = Qinside[3];
  }
}
 
double tutorials::exahype2::euler::EulerSolver::maxEigenvalue(
  [[maybe_unused]] const double* const NOALIAS                  Q, // Q[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const double                                 t,
  [[maybe_unused]] const double                                 dt,
  [[maybe_unused]] const int                                    normal
) {
  const double     irho  = 1.0 / Q[0];
  constexpr double gamma = 1.4;
  const double     p     = (gamma - 1) * (Q[3] - 0.5 * irho * (Q[1] * Q[1] + Q[2] * Q[2]));
 
  const double c = std::sqrt(gamma * p * irho);
  const double u = Q[normal + 1] * irho;
 
  return std::max(std::abs(u - c), std::abs(u + c));
}
 
void tutorials::exahype2::euler::EulerSolver::flux(
  [[maybe_unused]] const double* const NOALIAS                  Q, // Q[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const double                                 t,
  [[maybe_unused]] const double                                 dt,
  [[maybe_unused]] const int                                    normal,
  [[maybe_unused]] double* const NOALIAS                        F // F[4]
) {
  const double     irho  = 1.0 / Q[0];
  constexpr double gamma = 1.4;
  const double     p     = (gamma - 1) * (Q[3] - 0.5 * irho * (Q[1] * Q[1] + Q[2] * Q[2]));
 
  F[0] = Q[normal + 1];
  F[1] = Q[normal + 1] * Q[1] * irho;
  F[2] = Q[normal + 1] * Q[2] * irho;
  F[3] = Q[normal + 1] * irho * (Q[3] + p);
 
  F[normal + 1] += p;
}

These should yield following solution for the final time: