capricorn_breakout

#include <iostream>
#include <cmath>
#include "hydro_snapshot.hpp"
#include "hdsim.hpp"
#include "ideal_gas.hpp"
#include "hll.hpp"
//#include "free_flow.hpp"
#include "vacuum.hpp"
#include "pcm.hpp"
#include "hdf5_diagnostics.hpp"
#include <fenv.h>

using namespace std;

namespace {

  double cell_width(double position)
  {
    return abs(position)*0.01+0.01;
  }

  vector<double> generate_grid(void)
  {
    const double left = -2000;
    const double right = 0;
    vector<double> res;
    for(double x=left;x<right;x+=cell_width(x))
      res.push_back(x);
    return res;
  }

  Primitive generate_cell(double position)
  {
    return Primitive(pow(-min(position,-1e-12),1.5),
             exp(-pow(position+1000,2)/10),
             0);
  }

  vector<Primitive> generate_cells(const vector<double>& grid)
  {
    vector<Primitive> cells(grid.size()-1);
    for(size_t i=0;i<grid.size()-1;++i)
      cells[i] = generate_cell(0.5*(grid[i]+grid[i+1]));
    return cells;
  }

  HydroSnapshot initial_conditions(void)
  {
    const vector<double> grid = generate_grid();
    const vector<Primitive> cells = generate_cells(grid);
    return HydroSnapshot(grid,cells);
  }

  class SimData
  {
  public:

    SimData(void):
      eos_(5./3.),
      rs_(eos_),
      bc_(rs_,1e-30),
      sr_(),
      sim_(initial_conditions(),
       eos_, rs_, bc_, sr_) {}

    HydroSimulation& getSim(void)
    {
      return sim_;
    }

  private:
    const IdealGas eos_;
    const HLL rs_;
    //    const FreeFlow bc_;
    const Vacuum bc_;
    const PCM sr_;
    HydroSimulation sim_;
  };

  /*
  bool detect_motion(const HydroSimulation& sim,
             double pos, double vel_thres)
  {
    size_t i = 0;
    while(sim.getSnapshot().grid[i]<pos)
      ++i;
    return sim.getSnapshot().cells[i].velocity<vel_thres;
  }
  */

  bool finish_line(const HydroSimulation& sim,
           const double pos)
  {
    return sim.getSnapshot().grid.back()<pos;
  }

  void main_loop(HydroSimulation& sim)
  {
    write_snapshot_to_hdf5(sim,"initial.h5");
    //    while(sim.getTime()<5e6){
    //while(detect_motion(sim,8e4,1e-6)){
    while(finish_line(sim,1e9)){
      sim.timeAdvance();
      cout << sim.getTime() << endl;
    }
    write_snapshot_to_hdf5(sim,"final.h5");
  }
}

int main(void)
{
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);

  main_loop(SimData().getSim());

  return 0;
}