planar_relativistic_breakout

#include <cassert>
#include <iostream>
#include <fstream>
#include <vector>
#include "srhd_sim.hpp"
#include "spatial_distribution.hpp"
#include "ideal_gas.hpp"
#include "hll.hpp"
#include "pcm.hpp"
#include "rigid_wall.hpp"
#include "main_loop.hpp"
#include "hdf5_snapshot.hpp"
#include "imgrs.hpp"
#include "rubric.hpp"
//#include <fenv.h>

using namespace std;

namespace {

  vector<double> generate_grid(double left, double right, double rat)
  {
    assert(left<0 && right<0);
    vector<double> res;
    for(double x=left*(1-0.5*rat);
    x<right*(1+0.5*rat);x-=rat*x)
      res.push_back(x);
    return res;
  }

  /*
  class Gaussian: public SpatialDistribution
  {
  public:

    Gaussian(double prefactor,
         double mean,
         double std,
         double minval):
      prefactor_(prefactor),
      mean_(mean), std_(std), minval_(minval) {}

    double operator()(double x) const
    {
      if(abs(x-mean_)/std_>20)
    return 0;
      return max(prefactor_*exp(-pow((x-mean_)/std_,2)),
         minval_);
    }

  private:
    const double prefactor_;
    const double mean_;
    const double std_;
    const double minval_;
  };
  */

  class PowerLaw: public SpatialDistribution
  {
  public:

    PowerLaw(double prefactor, double index):
      prefactor_(prefactor), index_(index) {}

    double operator()(double x) const
    {
      return prefactor_*pow(-x,index_);
    }

  private:
    const double prefactor_;
    const double index_;
  };

  class MockVacuum: public BoundaryCondition
  {
  public:

    MockVacuum(const RiemannSolver& rs):
      rs_(rs) {}

    RiemannSolution CalcRS
    (size_t idx, const vector<Primitive>& cells) const
    {
      assert(0==idx || cells.size()==idx);
      if(0==idx){
    const Primitive right = cells.front();
    const Primitive left(right.Density*0.99,
                 right.Pressure*0.99,
                 right.Celerity);
    return rs_(left,right);
      }
      const Primitive left = cells.back();
      const Primitive right(left.Density/100,
                left.Pressure/100,
                left.Celerity);
      return rs_(left, right);
    }

  private:
    const RiemannSolver& rs_;
  };

  /*
  class IdealGasVacuum: public BoundaryCondition
  {
  public:

    IdealGasVacuum(double g): g_(g) {}

    RiemannSolution CalcRS
    (size_t idx, const vector<Primitive>& cells) const
    {
      assert(0==idx || cells.size()==idx);
      return RiemannSolution
    (0, calc_vacuum_velocity(idx==cells.size() ?
                 cells.back() : cells.front(),
                 idx==cells.size()));
    }

  private:

    double calc_vacuum_velocity(const Primitive& cell,
                bool dir) const
    {
      const double ba = IdealGas(g_).dp2ba(cell.Density, cell.Pressure);
      if(dir)
    return tanh(atanh(celerity2velocity(cell.Celerity))+
            (2./sqrt(g_-1))*
            atanh(ba/sqrt(g_-1)));
      else
    return tanh(atanh(celerity2velocity(cell.Celerity))-
            (2./sqrt(g_-1))*
            atanh(ba/sqrt(g_-1)));
    }

    const double g_;
  };
  */

  /*
  class ConstGhost: public BoundaryCondition
  {
  public:

    ConstGhost(const Primitive& ghost,
           const RiemannSolver& rs):
      ghost_(ghost), rs_(rs) {}

    RiemannSolution CalcRS
    (size_t idx, const vector<Primitive>& cells) const
    {
      assert(0==idx || cells.size()==idx);
      return idx==0 ? rs_(ghost_, cells.front()) :
    rs_(cells.back(), ghost_);
    }

  private:
    const Primitive ghost_;
    const RiemannSolver& rs_;
  };
  */

  class SimData
  {
  public:

    SimData(void):
      eos_(4./3.),
      rs_(eos_),
      //rs_(eos_.getAdiabaticIndex()),
      //bc_(Primitive(1e-14,1e-14,0), rs_),
      //      bc_(eos_.getAdiabaticIndex()),
      bc_left_(rs_),
      bc_right_(rs_),
      sr_(),
      geometry_(),
      sim_(generate_grid(-1,-2e-3,0.01),
       PowerLaw(1,1.5),
       //Gaussian(1e2,-1,0.05,1e-14),
       //      TwoSteps(1e-30,-1.01,1e3,-0.99,1e-30),
       Step(1e2,1e-30,-0.99),
       Uniform(0),
       //      bc_, bc_,
       bc_left_,
       bc_right_,
       eos_,
       rs_,
       sr_,
       geometry_) {}

    SRHDSimulation& getSim(void)
    {
      return sim_;
    }


  private:
    IdealGas eos_;
    HLL rs_;
    //    IdealGasRiemannSolver rs_;
    RigidWall bc_left_;
    MockVacuum bc_right_;
    //    ConstGhost bc_;
    //    IdealGasVacuum bc_;
    PCM sr_;
    const Planar geometry_;
    SRHDSimulation sim_;
  };

  /*
  class ConditionalSnapshot: public DiagnosticFunction
  {
  public:

    ConditionalSnapshot(void):
      rubric_("snapshot_",".h5") {}

    void operator()(const SRHDSimulation& sim) const
    {
      if(sim.getHydroSnapshot().edges.back()<0)
    write_hdf5_snapshot(sim, rubric_(sim.GetCycle()));
    }

  private:
    mutable Rubric rubric_;
  };
  */

  class MyTimeTermination: public TerminationCondition
  {
  public:

    MyTimeTermination(void) {}

    bool operator()(const SRHDSimulation& sim) const
    {
      //      return sim.getHydroSnapshot().edges.back() > -1e-9;
      //      const size_t n = sim.getHydroSnapshot().cells.size();
      return sim.getHydroSnapshot().cells.back().Celerity >= 2;
      //      return sim.GetTime()>0.1;
    }
  };

  void my_main_loop(SRHDSimulation& sim)
  {
    //    WriteTime diag("time.txt");
    MultipleDiagnostics diag;
    diag.diag_list.push_back(new WriteTime("time.txt"));
    //    diag.diag_list.push_back(new ConditionalSnapshot());
    write_hdf5_snapshot(sim,"initial.h5");
    main_loop(sim, MyTimeTermination(),
          &SRHDSimulation::TimeAdvance, diag);
    write_hdf5_snapshot(sim,"mid.h5");
    for(int tf=10;tf<=1e6;tf*=10){
      SafeTimeTermination stt(static_cast<double>(tf),1e9);
      main_loop(sim, stt,
        &SRHDSimulation::TimeAdvance, diag);
      write_hdf5_snapshot(sim,"milestone_"+int2str(tf)+".h5");
    }
  }
}

int main()
{
  /*
  feenableexcept(FE_INVALID   |
         FE_DIVBYZERO |
         FE_OVERFLOW  |
         FE_UNDERFLOW);
  */

  my_main_loop(SimData().getSim());

  return 0;
}