//========================================================================================
// Athena++ astrophysical MHD code
// Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
// Licensed under the 3-clause BSD License, see LICENSE file for details
//========================================================================================
//! \file two_ibw.c
//  \brief Problem generator for two interacting blast waves test.
//
// REFERENCE: P. Woodward & P. Colella, "The numerical simulation of two-dimensional
//   fluid flow with strong shocks", JCP, 54, 115, sect. IVa
//========================================================================================

// C/C++ headers
#include <iostream>   // endl
#include <sstream>    // stringstream
#include <stdexcept>  // runtime_error
#include <string>     // c_str()

// Athena++ headers
#include "../athena.hpp"
#include "../athena_arrays.hpp"
#include "../parameter_input.hpp"
#include "../coordinates/coordinates.hpp"
#include "../eos/eos.hpp"
#include "../hydro/hydro.hpp"
#include "../mesh/mesh.hpp"

#if MAGNETIC_FIELDS_ENABLED
#error "This problem generator does not support magnetic fields"
#endif

//========================================================================================
//! \fn void MeshBlock::ProblemGenerator(ParameterInput *pin)
//  \brief Problem Generator for the two interacting blast waves
//========================================================================================

void MeshBlock::ProblemGenerator(ParameterInput *pin)
{
  std::stringstream msg;

  // parse shock direction: {1,2,3} -> {x1,x2,x3}
  int shk_dir = pin->GetInteger("problem","shock_dir");
  if (shk_dir < 1 || shk_dir > 3) {
    msg << "### FATAL ERROR in Problem Generator" << std::endl
        << "shk_dir = " << shk_dir << " must be either 1,2, or 3" << std::endl;
    throw std::runtime_error(msg.str().c_str());
  }

  for (int k=ks; k<=ke; ++k) {
  for (int j=js; j<=je; ++j) {
#pragma simd
    for (int i=is; i<=ie; ++i) {
      phydro->u(IDN,k,j,i) = 1.0;
      phydro->u(IM1,k,j,i) = 0.0;
      phydro->u(IM2,k,j,i) = 0.0;
      phydro->u(IM3,k,j,i) = 0.0;
      if ((shk_dir==1 && pcoord->x1v(i) < 0.1) ||
          (shk_dir==2 && pcoord->x2v(j) < 0.1) ||
          (shk_dir==3 && pcoord->x3v(k) < 0.1)) {
        phydro->u(IEN,k,j,i)= 1.0e3/(peos->GetGamma() - 1.0);
      }
      else if ((shk_dir==1 && pcoord->x1v(i) > 0.9) ||
               (shk_dir==2 && pcoord->x2v(j) > 0.9) ||
               (shk_dir==3 && pcoord->x3v(k) > 0.9)) {
        phydro->u(IEN,k,j,i)= 1.0e2/(peos->GetGamma() - 1.0);
      }
      else {
        phydro->u(IEN,k,j,i)= 0.01/(peos->GetGamma() - 1.0);
      }
    }
  }}

  return;
}
