acoustic_planar_waves.py

Go to the documentation of this file.

# This file is part of the ExaHyPE2 project. For conditions of distribution and
# use, please see the copyright notice at www.peano-framework.org
from .scenario import Scenario
 
import os
import sys
from math import sqrt
 
sys.path.insert(0, os.path.abspath("../equations"))
from equations import Acoustic
 
 
class AcousticPlanarWaves(Scenario):
    """
    Scenario reproduced from Dumbser & Käser, https://doi.org/10.1111/j.1365-246X.2006.03120.x (Chapter 5, p. 328)
    """
 
    _plot_dt = 0.0
    _offset = -1.0
    _domain_size = 2.0
    _periodic_bc = True
 
    def __init__(self, dimensions, iterations=2):
        self._dimensions = dimensions
        self._end_time = iterations * sqrt(dimensions)
        self._K0 = 4.0
        self._rho = 1.0
        self._equation = Acoustic(dimensions, rho=self._rho, K0=self._K0)
 
    def initial_conditions(self):
        return (
            """
  // simple translation in positive diagonal direction
  const double val = cos( - std::numbers::pi*(
      x[0] + x[1]
#if DIMENSIONS==3
      + x[2]
#endif
  ));
 
  Q[1] = val;
  Q[2] = val;
#if DIMENSIONS==3
  Q[3] = val;
#endif
 
  constexpr double K0  = """
            + str(self._K0)
            + """;
  constexpr double rho = """
            + str(self._rho)
            + """;
 
  // These are defined by the eigenvector of the plane wave operator
#if DIMENSIONS==3
  const double kr = K0*std::sqrt(rho/(3.0*K0)) + 2*std::sqrt(K0*rho/3.0);
#else
  const double kr = std::sqrt(2*K0*rho);
#endif
 
  Q[0] = kr*val;
"""
        )
 
    def analytical_solution(self):
        return (
            """
  const double w = 2*std::sqrt(DIMENSIONS)*M_PI;
 
  const double val = cos( w*t - std::numbers::pi*(
      x[0] + x[1]
#if DIMENSIONS==3
    + x[2]
#endif
  ));
 
  solution[1] = val;
  solution[2] = val;
 
  constexpr double K0  = """
            + str(self._K0)
            + """;
  constexpr double rho = """
            + str(self._rho)
            + """;
 
#if DIMENSIONS==3
  solution[3] = val;
  const double kr = K0*std::sqrt(rho/(3.0*K0)) + 2*std::sqrt(K0*rho/3.0);
#else
  const double kr = std::sqrt(2*K0*rho);
#endif
 
  solution[0] = kr*val;
"""
        )