initial_state_astro_jet.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2023 - 2024 by the ryujin authors
//
 
#pragma once
 
#include <initial_state_library.h>
 
namespace ryujin
{
  namespace EulerInitialStates
  {
    template <typename Description, int dim, typename Number>
    class AstroJet : public InitialState<Description, dim, Number>
    {
    public:
      using HyperbolicSystem = typename Description::HyperbolicSystem;
      using View =
          typename Description::template HyperbolicSystemView<dim, Number>;
      using state_type = typename View::state_type;
 
      AstroJet(const HyperbolicSystem &hyperbolic_system,
               const std::string subsection)
          : InitialState<Description, dim, Number>("astro jet", subsection)
          , hyperbolic_system_(hyperbolic_system)
      {
        gamma_ = 5. / 3.;
        if constexpr (!View::have_gamma) {
          this->add_parameter("gamma", gamma_, "The ratio of specific heats");
        }
 
        jet_width_ = 0.05;
        this->add_parameter("jet width",
                            jet_width_,
                            "The width of the jet coming out of boundary");
 
        jet_state_[0] = 5.0;
        jet_state_[1] = 30.0;
        jet_state_[2] = 0.4127;
        this->add_parameter(
            "primitive jet state",
            jet_state_,
            "Initial primitive state (rho, u, p) for jet state");
 
        ambient_state_[0] = 5.0;
        ambient_state_[1] = 0.0;
        ambient_state_[2] = 0.4127;
        this->add_parameter(
            "primitive ambient right",
            ambient_state_,
            "Initial primitive state (rho, u, p) for ambient state");
 
        const auto convert_states = [&]() {
          const auto view = hyperbolic_system_.template view<dim, Number>();
          state_left_ = view.from_initial_state(jet_state_);
          state_right_ = view.from_initial_state(ambient_state_);
        };
        this->parse_parameters_call_back.connect(convert_states);
        convert_states();
      }
 
      state_type compute(const dealii::Point<dim> &point, Number /*t*/) final
      {
        return (point[0] < 1.e-12 && std::abs(point[1]) <= jet_width_
                    ? state_left_
                    : state_right_);
      }
 
    private:
      const HyperbolicSystem &hyperbolic_system_;
 
      Number gamma_;
      Number jet_width_;
 
      dealii::Tensor<1, 3, Number> jet_state_;
      dealii::Tensor<1, 3, Number> ambient_state_;
 
      state_type state_left_;
      state_type state_right_;
    };
  } // namespace EulerInitialStates
} // namespace ryujin