off centre supernova 0.5 radius

#include <iostream>
#include "source/tessellation/VoronoiMesh.hpp"
#include "source/newtonian/two_dimensional/hdsim2d.hpp"
#include "source/newtonian/two_dimensional/interpolations/pcm2d.hpp"
#include "source/newtonian/two_dimensional/spatial_distributions/uniform2d.hpp"
#include "source/newtonian/common/ideal_gas.hpp"
#include "source/newtonian/two_dimensional/geometric_outer_boundaries/SquareBox.hpp"
#include "source/newtonian/two_dimensional/hydro_boundary_conditions/RigidWallHydro.hpp"
#include "source/newtonian/common/hllc.hpp"
#include "source/newtonian/two_dimensional/source_terms/zero_force.hpp"
#include "source/newtonian/two_dimensional/point_motions/eulerian.hpp"
#include "source/newtonian/two_dimensional/point_motions/lagrangian.hpp"
#include "source/newtonian/two_dimensional/point_motions/round_cells.hpp"
#include "source/newtonian/two_dimensional/diagnostics.hpp"
#include "source/newtonian/two_dimensional/source_terms/cylindrical_complementary.hpp"
#include "source/misc/simple_io.hpp"
#include "source/newtonian/test_2d/main_loop_2d.hpp"
#include "source/newtonian/two_dimensional/hdf5_diagnostics.hpp"
#include "source/misc/mesh_generator.hpp"
#include "source/tessellation/shape_2d.hpp"
#include "source/newtonian/test_2d/piecewise.hpp"
#include "source/newtonian/two_dimensional/simple_flux_calculator.hpp"
#include "source/newtonian/two_dimensional/simple_cell_updater.hpp"
#include "source/newtonian/two_dimensional/simple_extensive_updater.hpp"
#include "source/newtonian/test_2d/consecutive_snapshots.hpp"
#include "source/newtonian/test_2d/multiple_diagnostics.hpp"

using namespace std;
using namespace simulation2d;

namespace {

  class PhysicalConstants
  {
  public:

    // Metric prefixes
    const double kilo;
    const double centi;

    // Length
    const double solar_radius;
    const double meter;

    // Mass
    const double solar_mass;
    const double gram;

    // Time
    const double second;

    // Energy
    const double erg;

    PhysicalConstants(void):
      kilo(1e3),
      centi(1e-2),
      solar_radius(1),
      meter(1.438e-9),
      solar_mass(1),
      gram(5.029e-34),
      second(1),
      erg(gram*pow(centi*meter/second,2)) {}
  };

  vector<ComputationalCell> calc_init_cond(const Tessellation& tess,
                       const PhysicalConstants& c)
  {
    vector<ComputationalCell> res(static_cast<size_t>(tess.GetPointNo()));
    const double star_mass = 10*c.solar_mass;
    const double star_radius = 50*c.solar_radius;
    const double average_density = star_mass/pow(star_radius,3);
    const double hot_spot_radius = 1*c.solar_radius;
    const double hot_spot_volume = (4.*M_PI/3.)*pow(hot_spot_radius,3);
    const double hot_spot_energy = 1e51*c.erg;
    const double hot_spot_pressure = hot_spot_energy/hot_spot_volume;
    const Circle hot_spot(Vector2D(0,0.5*star_radius), hot_spot_radius);
    for(size_t i=0;i<res.size();++i){
      res[i].density = 1e-9*average_density;
      res[i].pressure = 1e-9*hot_spot_pressure;
      res[i].velocity = Vector2D(0,0);
      const Vector2D r = tess.GetMeshPoint(static_cast<int>(i));
      if(abs(r)<star_radius)
    res[i].density = average_density*pow(1-abs(r)/star_radius,1.5);
      if(hot_spot(r))
    res[i].pressure = hot_spot_pressure;
    }
    return res;
  }

  double calc_tracer_flux(const Edge& edge,
              const Tessellation& tess,
              const vector<ComputationalCell>& cells,
              const string& name,
              const Conserved& hf)
  {
    if(hf.Mass>0 &&
       edge.neighbors.first>0 &&
       edge.neighbors.first<tess.GetPointNo())
      return hf.Mass*
    cells[static_cast<size_t>(edge.neighbors.first)].tracers.find(name)->second;
    if(hf.Mass<0 &&
       edge.neighbors.second>0 &&
       edge.neighbors.second<tess.GetPointNo())
      return hf.Mass*
    cells[static_cast<size_t>(edge.neighbors.second)].tracers.find(name)->second;
    return 0;
  }

  class CustomFluxCalculator: public FluxCalculator
  {
  public:

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

    vector<Extensive> operator()
    (const Tessellation& tess,
     const vector<Vector2D>& point_velocities,
     const vector<ComputationalCell>& cells,
     const EquationOfState& eos,
     const double /*time*/,
     const double /*dt*/) const
    {
      vector<Extensive> res(tess.getAllEdges().size());
      for(size_t i=0;i<tess.getAllEdges().size();++i){
    const Conserved hydro_flux =
      calcHydroFlux(tess, point_velocities,
            cells, eos, i);
    res[i].mass = hydro_flux.Mass;
    res[i].momentum = hydro_flux.Momentum;
    res[i].energy = hydro_flux.Energy;
    for(std::map<std::string,double>::const_iterator it =
          cells.front().tracers.begin();
        it!=cells.front().tracers.end();++it)
      res[i].tracers[it->first] =
        calc_tracer_flux(tess.getAllEdges()[i],
                 tess,cells,it->first,hydro_flux);

      }
      return res;
    }

  private:
    const RiemannSolver& rs_;

    const Conserved calcHydroFlux
    (const Tessellation& tess,
     const vector<Vector2D>& point_velocities,
     const vector<ComputationalCell>& cells,
     const EquationOfState& eos,
     const size_t i) const
    {
      const Edge& edge = tess.GetEdge(static_cast<int>(i));
      const std::pair<bool,bool> flags
    (edge.neighbors.first>=0 && edge.neighbors.first<tess.GetPointNo(),
     edge.neighbors.second>=0 && edge.neighbors.second<tess.GetPointNo());
      assert(flags.first || flags.second);
      if(!flags.first){
    const size_t right_index =
      static_cast<size_t>(edge.neighbors.second);
    const ComputationalCell& right_cell = cells[right_index];
    const Vector2D p = Parallel(edge);
    const Primitive right = convert_to_primitive(right_cell,eos);
    const Primitive left = right;
    const Vector2D n = remove_parallel_component
      (tess.GetMeshPoint(edge.neighbors.second) -
       edge.vertices.first, p);
    return rotate_solve_rotate_back
      (rs_, left, right, 0, n, p);
      }
      if(!flags.second){
    const size_t left_index =
      static_cast<size_t>(edge.neighbors.first);
    const ComputationalCell& left_cell = cells[left_index];
    const Primitive left = convert_to_primitive(left_cell, eos);
    const Vector2D p = Parallel(edge);
    const Primitive right = left;
    const Vector2D n = remove_parallel_component
      (edge.vertices.second -
       tess.GetMeshPoint(edge.neighbors.first), p);
    return rotate_solve_rotate_back
      (rs_, left, right, 0, n, p);
      }
      const size_t left_index =
    static_cast<size_t>(edge.neighbors.first);
      const size_t right_index =
    static_cast<size_t>(edge.neighbors.second);
      const ComputationalCell& left_cell = cells[left_index];
      const ComputationalCell& right_cell = cells[right_index];
      const Vector2D p = Parallel(edge);
      const Vector2D n =
    tess.GetMeshPoint(edge.neighbors.second) -
    tess.GetMeshPoint(edge.neighbors.first);
      const double velocity = Projection
    (tess.CalcFaceVelocity
     (point_velocities[left_index],
      point_velocities[right_index],
      tess.GetCellCM(edge.neighbors.first),
      tess.GetCellCM(edge.neighbors.second),
      calc_centroid(edge)),n);
      const Primitive left =
    convert_to_primitive(left_cell, eos);
      const Primitive right =
    convert_to_primitive(right_cell, eos);
      return rotate_solve_rotate_back
    (rs_,left,right,velocity,n,p);
    }
  };

  class SimData
  {
  public:
    SimData(const PhysicalConstants& c):
      pg_(Vector2D(0,0), Vector2D(0,1)),
      outer_(Vector2D(0,-100*c.solar_radius),
         Vector2D(200*c.solar_radius,
              300*c.solar_radius)),
      mesh_(cartesian_mesh(2*200,2*300,
               outer_.getBoundary().first,
               outer_.getBoundary().second)),
      tess_(mesh_,outer_),
      interpm_(),
      eos_(5./3.),
      rs_(),
      raw_point_motion_(),
      point_motion_(raw_point_motion_,eos_),
      alt_point_motion_(),
      force_(pg_.getAxis()),
      tsf_(0.3),
      fc_(rs_),
      eu_(),
      cu_(),
      sim_(tess_,
       outer_,
       pg_,
       calc_init_cond(tess_,c),
       eos_,
       alt_point_motion_,
       force_,
       tsf_,
       fc_,
       eu_,
       cu_) {}

    hdsim& getSim(void)
    {
      return sim_;
    }

  private:
    const CylindricalSymmetry pg_;
    const SquareBox outer_;
    const vector<Vector2D> mesh_;
    VoronoiMesh tess_;
    PCM2D interpm_;
    const IdealGas eos_;
    const Hllc rs_;
    //    const RigidWallHydro hbc_;
    Lagrangian raw_point_motion_;
    RoundCells point_motion_;
    Eulerian alt_point_motion_;
    CylindricalComplementary force_;
    const SimpleCFL tsf_;
    //    const SimpleFluxCalculator fc_;
    const CustomFluxCalculator fc_;
    const SimpleExtensiveUpdater eu_;
    const SimpleCellUpdater cu_;
    hdsim sim_;
  };

  void my_main_loop(hdsim& sim,
            const PhysicalConstants& c)
  {
    const double tf = 20000*c.second;
    SafeTimeTermination term_cond(tf,1e6);
    WriteTime diag1("time.txt");
    ConsecutiveSnapshots diag2(new ConstantTimeInterval(tf/1000),
                   new Rubric("snapshot_",".h5"));
    MultipleDiagnostics diag;
    diag.diag_list.push_back(&diag1);
    diag.diag_list.push_back(&diag2);
    main_loop(sim,
          term_cond,
          &hdsim::TimeAdvance,
          &diag);
  }
}

int main(void)
{
  const PhysicalConstants c;
  SimData sim_data(c);
  hdsim& sim = sim_data.getSim();

  write_snapshot_to_hdf5(sim, "initial.h5");

  my_main_loop(sim,c);

  write_snapshot_to_hdf5(sim, "final.h5");

  return 0;
}