API tools faq
paste
Advertisement
bolverk

offset_supernova_gravity_sink

bolverk
Mar 4th, 2015
434
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 23.03 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <cmath>
  3. #include <limits>
  4. #include "source/tessellation/VoronoiMesh.hpp"
  5. #include "source/newtonian/two_dimensional/hdsim2d.hpp"
  6. #include "source/newtonian/two_dimensional/interpolations/pcm2d.hpp"
  7. #include "source/newtonian/two_dimensional/spatial_distributions/uniform2d.hpp"
  8. #include "source/newtonian/common/ideal_gas.hpp"
  9. #include "source/newtonian/two_dimensional/geometric_outer_boundaries/SquareBox.hpp"
  10. #include "source/newtonian/two_dimensional/hydro_boundary_conditions/RigidWallHydro.hpp"
  11. #include "source/newtonian/common/hllc.hpp"
  12. #include "source/newtonian/two_dimensional/source_terms/zero_force.hpp"
  13. #include "source/newtonian/two_dimensional/point_motions/eulerian.hpp"
  14. #include "source/newtonian/two_dimensional/point_motions/lagrangian.hpp"
  15. #include "source/newtonian/two_dimensional/point_motions/eulerian.hpp"
  16. #include "source/newtonian/two_dimensional/point_motions/round_cells.hpp"
  17. #include "source/newtonian/two_dimensional/diagnostics.hpp"
  18. #include "source/newtonian/two_dimensional/source_terms/cylindrical_complementary.hpp"
  19. #include "source/newtonian/two_dimensional/source_terms/CenterGravity.hpp"
  20. #include "source/newtonian/two_dimensional/source_terms/SeveralSources.hpp"
  21. #include "source/misc/simple_io.hpp"
  22. #include "source/newtonian/test_2d/main_loop_2d.hpp"
  23. #include "source/newtonian/test_2d/kill_switch.hpp"
  24. #include "source/newtonian/two_dimensional/hdf5_diagnostics.hpp"
  25. #include "source/misc/mesh_generator.hpp"
  26. #include "source/tessellation/RoundGrid.hpp"
  27. #include "source/misc/int2str.hpp"
  28. #include "source/newtonian/two_dimensional/interpolations/linear_gauss_consistent.hpp"
  29. #include "source/newtonian/test_2d/consecutive_snapshots.hpp"
  30. #include <boost/foreach.hpp>
  31. #include "source/newtonian/two_dimensional/hydro_boundary_conditions/FreeFlow.hpp"
  32. #include "source/misc/utils.hpp"
  33. #include "source/tessellation/shape_2d.hpp"
  34. #include "source/newtonian/test_2d/piecewise.hpp"
  35. #include "source/mpi/MeshPointsMPI.hpp"
  36. #include "source/mpi/mpi_macro.hpp"
  37. #include "source/mpi/ConstNumberPerProc.hpp"
  38. #include "source/misc/hdf5_utils.hpp"
  39. #include "source/newtonian/test_2d/contour.hpp"
  40. #include "source/newtonian/two_dimensional/custom_evolutions/Ratchet.cpp"
  41. #include "source/newtonian/two_dimensional/diagnostics.hpp"
  42. #include <fenv.h>
  43.  
  44. using namespace std;
  45. using namespace simulation2d;
  46.  
  47. namespace {
  48.  
  49.   class WriteConserved: public DiagnosticFunction
  50.   {
  51.   public:
  52.  
  53.     WriteConserved(string const& fname):
  54.       tracer_(), cons_(), fname_(fname) {}
  55.  
  56.     void operator()(hdsim const& sim)
  57.     {
  58.       tracer_.push_back(total_tracer(sim,0));
  59.       cons_.push_back(total_conserved(sim));
  60.       time_.push_back(sim.GetTime());
  61.     }
  62.  
  63.     ~WriteConserved(void)
  64.     {
  65.     ofstream f(fname_.c_str());
  66.     for(size_t i=0;i<cons_.size();++i)
  67.       f << time_[i] << " "
  68.         << cons_[i].Mass << " "
  69.         << cons_[i].Momentum.x << " "
  70.         << cons_[i].Momentum.y << " "
  71.         << cons_[i].Energy << " "
  72.         << tracer_[i] << "\n";
  73.     f.close();
  74.     }
  75.  
  76.   private:
  77.     mutable vector<double> tracer_;
  78.     mutable vector<Conserved> cons_;
  79.     mutable vector<double> time_;
  80.     const string fname_;
  81.   };
  82.  
  83.   /*
  84.     class Veil: public Acceleration
  85.     {
  86.     public:
  87.  
  88.     Veil(Acceleration& full):
  89.     full_(full) {}
  90.  
  91.     Vector2D Calculate
  92.     (Tessellation const &  tess,
  93.     vector< Primitive > const &  cells,
  94.     int  point,
  95.     vector< Conserved > const &  fluxes,
  96.     vector< Vector2D > const &  point_velocity,
  97.     HydroBoundaryConditions const &  hbc,
  98.     vector< vector< double > > const &  tracers,
  99.     double  time,
  100.     double  dt)
  101.     {
  102.     const Primitive& cell = cells[static_cast<size_t>(point)];
  103.     if(abs(cell.Velocity.x)<std::numeric_limits<double>::epsilon()&&
  104.     abs(cell.Velocity.y)<std::numeric_limits<double>::epsilon())
  105.     return Vector2D(0,0);
  106.     else
  107.     return full_.Calculate(tess,cells,point,fluxes,point_velocity,hbc,tracers,
  108.     time,dt);
  109.     }
  110.  
  111.     private:
  112.     Acceleration& full_;
  113.     };*/
  114.  
  115.   double calc_mach_number(const Primitive& p)
  116.   {
  117.     return abs(p.Velocity)/p.SoundSpeed;
  118.   }
  119.  
  120.   Vector2D find_pos_on_line(const Vector2D& p1,
  121.                 double v1,
  122.                 const Vector2D& p2,
  123.                 double v2,
  124.                 double vt)
  125.   {
  126.     return p1+((vt-v1)/(v2-v1))*(p2-p1);
  127.   }
  128.  
  129.   /*
  130.     vector<Vector2D> process_positions(const Index2Member<Vector2D>& gg)
  131.     {
  132.     vector<Vector2D> res(get_mpi_size());
  133.     for(size_t i=0;i<res.size();++i)
  134.     res[i] = gg(i*gg.getLength()/get_mpi_size());
  135.  
  136.     return res;
  137.     }
  138.   */
  139.  
  140.   class Constants
  141.   {
  142.   public:
  143.  
  144.     // Metric prefix
  145.     const double kilo;
  146.     const double centi;
  147.  
  148.     // Length / distance
  149.     const double parsec;
  150.     const double meter;
  151.  
  152.     // Time
  153.     const double year;
  154.     const double second;
  155.  
  156.     // Mass
  157.     const double solar_mass;
  158.     const double gram;
  159.  
  160.     // Energy
  161.     const double erg;
  162.  
  163.     // Force
  164.     const double newton;
  165.  
  166.     // Parameters
  167.     const double adiabatic_index;
  168.     const double gravitation_constant;
  169.     const double rho_0;
  170.     const double R_0;
  171.     const double omega_in;
  172.     const double inner_density_prefactor;
  173.     const double R_b;
  174.     const double omega_out;
  175.     const double outer_density_prefactor;
  176.     const double offset;
  177.     const double supernova_energy;
  178.     const double supernova_radius;
  179.     const double supernova_volume;
  180.     const double supernova_mass;
  181.     const double supernova_density;
  182.     const double supernova_pressure;
  183.     const double black_hole_mass;
  184.     const Vector2D lower_left;
  185.     const Vector2D upper_right;
  186.  
  187.     Constants(void):
  188.       kilo(1e3),
  189.       centi(1e-2),
  190.       parsec(1),
  191.       meter(3.24e-17*parsec),
  192.       year(1),
  193.       second(year/3.16e7),
  194.       solar_mass(1),
  195.       gram(5.03e-34*solar_mass),
  196.       erg(gram*pow(centi*meter/second,2)),
  197.       newton(kilo*gram*meter/pow(second,2)),
  198.       adiabatic_index(5./3.),
  199.       gravitation_constant(6.673e-11*newton*pow(meter/(kilo*gram),2)),
  200.       rho_0(2.2e-22*gram/pow(centi*meter,3)),
  201.       R_0(0.04*parsec),
  202.       omega_in(1),
  203.       inner_density_prefactor(rho_0*pow(R_0,omega_in)),
  204.       R_b(0.4*parsec),
  205.       omega_out(3),
  206.       outer_density_prefactor
  207.       (inner_density_prefactor*pow(R_b,omega_out-omega_in)),
  208.       offset(1*parsec),
  209.       supernova_energy(1e51*erg),
  210.       supernova_radius(0.2*parsec),
  211.       supernova_volume((4.*M_PI/3)*pow(supernova_radius,3)),
  212.       supernova_mass(5*solar_mass),
  213.       supernova_density(supernova_mass/supernova_volume),
  214.       supernova_pressure((adiabatic_index-1)*supernova_energy/supernova_volume),
  215.       black_hole_mass(1e9*solar_mass),
  216.       lower_left(parsec*Vector2D(0,-6)),
  217.       upper_right(parsec*Vector2D(6,3)) {}
  218.   };
  219. }
  220.  
  221. namespace units
  222. {
  223.   // Time
  224.   const double year = 1;
  225.   //  const double second = year/3.16e7;
  226.  
  227.   // parameters
  228.   //  const double gravitation_constant = 6.673e-11*newton*pow(meter/(kilo*gram),2);
  229.   //  const double black_hole_mass = 1e7*solar_mass;
  230.   const double adiabatic_index = 5./3.;
  231.   const double omega_in = 1;
  232.   const double omega_out = 3;
  233.   //  const double R_0 = 0.04*parsec;
  234.   //  const double rho_0 = 2.2e-22*gram/pow(centi*meter,3);
  235.   //  const double R_b = 0.4*parsec;
  236.   //  const double inner_density_prefactor = rho_0*pow(R_0,omega_in);
  237.   /*
  238.     const double outer_density_prefactor =
  239.     inner_density_prefactor*pow(R_b,omega_out-omega_in);
  240.   */
  241.   //  const double supernova_radius = 0.2*parsec;
  242.   //  const double supernova_energy = 1e51*erg;
  243.   //  const double supernova_mass = 5*solar_mass;
  244.   //  const double supernova_volume = (4.*M_PI/3.)*pow(supernova_radius,3);
  245.   //const double supernova_density = supernova_mass / supernova_volume;
  246.   //const double supernova_pressure =
  247.   //  (adiabatic_index-1)*(supernova_energy/supernova_volume);
  248.   //const double ambient_pressure = supernova_pressure*1e-10;
  249.   //  const double ambient_pressure = 1e-30;
  250.   //  const double offset = 1*parsec;
  251.   /*
  252.     const Vector2D lower_left = parsec*Vector2D(1e-3,-6);
  253.     const Vector2D upper_right = parsec*Vector2D(6,6);
  254.   */
  255.   //  const Vector2D lower_left = parsec*Vector2D(0,-6);
  256.   //  const Vector2D upper_right = parsec*Vector2D(6,2);
  257. }
  258.  
  259. using namespace units;
  260.  
  261. namespace {
  262.  
  263.   class HydrostaticPressure: public SpatialDistribution
  264.   {
  265.   public:
  266.  
  267.     HydrostaticPressure(double gm,
  268.             double r0,
  269.             double rho_0,
  270.             double omega_1,
  271.             double R_b,
  272.             double omega_2):
  273.       gm_(gm), r0_(r0), rho_0_(rho_0),
  274.       omega_1_(omega_1), R_b_(R_b), omega_2_(omega_2) {}
  275.  
  276.     double operator()(const Vector2D& r) const
  277.     {
  278.       if(abs(r)<R_b_)
  279.     return (gm_/r0_)*(rho_0_/(omega_1_+1))*pow(abs(r)/r0_,-omega_1_-1)-
  280.       (gm_/r0_)*(rho_0_/(omega_1_+1))*pow(R_b_/r0_,-omega_1_-1)+
  281.       (gm_/r0_)*(rho_0_/(omega_2_+1))*pow(R_b_/r0_,-omega_2_-1);
  282.       else
  283.     return (gm_/r0_)*(rho_0_/(omega_2_+1))*
  284.       pow(R_b_/r0_,-omega_1_-1)*
  285.       pow(abs(r)/R_b_,-omega_2_-1);
  286.     }
  287.  
  288.   private:
  289.     const double gm_;
  290.     const double r0_;
  291.     const double rho_0_;
  292.     const double omega_1_;
  293.     const double R_b_;
  294.     const double omega_2_;
  295.   };
  296.  
  297.   class BoundedRadialPowerLaw: public SpatialDistribution
  298.   {
  299.   public:
  300.  
  301.     BoundedRadialPowerLaw(const Vector2D& center,
  302.               double index,
  303.               double prefactor,
  304.               double lower_bound,
  305.               double upper_bound):
  306.       center_(center),
  307.       index_(index),
  308.       prefactor_(prefactor),
  309.       lower_bound_(lower_bound),
  310.       upper_bound_(upper_bound) {}
  311.  
  312.     double operator()(const Vector2D& r) const
  313.     {
  314.       const double radius = max(min(abs(r-center_),
  315.                     upper_bound_),
  316.                 lower_bound_);
  317.       return prefactor_*pow(radius,index_);
  318.     }
  319.  
  320.   private:
  321.     const Vector2D center_;
  322.     const double index_;
  323.     const double prefactor_;
  324.     const double lower_bound_;
  325.     const double upper_bound_;
  326.   };
  327.  
  328.   bool is_inside_rectangle(const Vector2D& point,
  329.                const Vector2D& lower_left,
  330.                const Vector2D& upper_right)
  331.   {
  332.     return ((point.x > lower_left.x) &&
  333.         (point.x < upper_right.x) &&
  334.         (point.y > lower_left.y) &&
  335.         (point.y < upper_right.y));
  336.   }
  337.  
  338.   vector<Vector2D> rectangle_clip(const vector<Vector2D>& grid,
  339.                   const Vector2D& lower_left,
  340.                   const Vector2D& upper_right)
  341.   {
  342.     vector<Vector2D> res;
  343.     for(size_t i=0, endp=grid.size();i<endp;++i){
  344.       const Vector2D& point = grid[i];
  345.       if(is_inside_rectangle(point,lower_left,upper_right))
  346.     res.push_back(point);    
  347.     }
  348.     return res;
  349.   }
  350.  
  351.   /*
  352.     vector<double> calc_radius_list(void)
  353.     {
  354.     const double rmin = 1e-4;
  355.     const double q = 1.01;
  356.     const size_t n = 200;
  357.     vector<double> res(n);
  358.     for(size_t i=0;i<n;++i)
  359.     res[i] = rmin*pow(q,i);
  360.     write_number(0,"finished_calc_radius_list.txt");
  361.     return res;
  362.     }
  363.   */
  364.  
  365.   /*
  366.     vector<Vector2D> create_grid(Vector2D const& lower_left,
  367.     Vector2D const& upper_right,
  368.     double dx2x)
  369.     {
  370.     vector<Vector2D> res;
  371.     for(double x = lower_left.x*(1+0.5*dx2x);
  372.     x<upper_right.x; x*=(1+dx2x)){
  373.     const double dx = x*dx2x;
  374.     for(double y=lower_left.y+dx/2;
  375.     y<upper_right.y; y += dx)
  376.     res.push_back(Vector2D(x,y));
  377.     }
  378.     return res;
  379.     }
  380.   */
  381.  
  382.   /*
  383.     vector<Vector2D> centered_hexagonal_grid(double r_min,
  384.     double r_max)
  385.     {
  386.     const vector<double> r_list = arange(0,r_max,r_min);
  387.     vector<Vector2D> res;
  388.     for(size_t i=0;i<r_list.size();++i){
  389.     const size_t angle_num = max<size_t>(6*i,1);
  390.     vector<double> angle_list(angle_num,0);
  391.     for(size_t j=0;j<angle_num;++j)
  392.     angle_list.at(j) = 2*M_PI*(double)j/(double)(angle_num);
  393.     for(size_t j=0;j<angle_num;++j)
  394.     res.push_back(r_list.at(i)*Vector2D(cos(angle_list.at(j)),
  395.     sin(angle_list.at(j))));
  396.     }
  397.     return res;
  398.     }
  399.   */
  400.  
  401.   vector<Vector2D> centered_logarithmic_spiral(double r_min,
  402.                            double r_max,
  403.                            double alpha,
  404.                            const Vector2D& center)
  405.   {
  406.     const double theta_max = log(r_max/r_min)/alpha;
  407.     const vector<double> theta_list =
  408.       arange(0,theta_max,2*M_PI*alpha);
  409.  
  410.     vector<double> r_list(theta_list.size(),0);
  411.     for(size_t i=0;i<r_list.size();++i)
  412.       r_list.at(i) = r_min*exp(alpha*theta_list.at(i));
  413.  
  414.     vector<Vector2D> res(r_list.size());
  415.     for(size_t i=0;i<res.size();++i)
  416.       res[i] = center+r_list[i]*Vector2D(cos(theta_list.at(i)),
  417.                      sin(theta_list.at(i)));
  418.     return res;
  419.   }
  420.  
  421.   /*
  422.     vector<Vector2D> complete_grid(double r_inner,
  423.     double r_outer,
  424.     double alpha)
  425.     {
  426.     const vector<Vector2D> inner =
  427.     centered_hexagonal_grid(r_inner*alpha*2*M_PI,
  428.     r_inner);
  429.     const vector<Vector2D> outer =
  430.     centered_logarithmic_spiral(r_inner,
  431.     r_outer,
  432.     alpha);
  433.     return join(inner, outer);
  434.     }
  435.   */
  436.  
  437.   template<class T> class VectorInitializer
  438.   {
  439.   public:
  440.  
  441.     VectorInitializer(const T& t):
  442.       buf_(1,t) {}
  443.  
  444.     VectorInitializer& operator()(const T& t)
  445.     {
  446.       buf_.push_back(t);
  447.       return *this;
  448.     }
  449.  
  450.     vector<T> operator()(void)
  451.     {
  452.       return buf_;
  453.     }
  454.  
  455.   private:
  456.     vector<T> buf_;
  457.   };
  458. }
  459.  
  460. class SimData
  461. {
  462.  
  463. public:
  464.  
  465.   SimData(const Constants& c):
  466.     pg_(Vector2D(0,0), Vector2D(0,1)),
  467.     outer_(c.lower_left, c.upper_right),
  468.     //    init_points_(cartesian_mesh(300,450,c.lower_left,c.upper_right)),
  469.     //cgg_(11+2*150,7+2*300,lower_left, upper_right),
  470.     init_points_(rectangle_clip
  471.          (centered_logarithmic_spiral(0.01,
  472.                           abs(c.upper_right-c.lower_left),
  473.                           0.002,
  474.                           Vector2D(-0.2,-1)),
  475.           c.lower_left, c.upper_right)),
  476.     /*
  477.       init_points_(create_grid(lower_left,
  478.       upper_right,
  479.       0.1)),
  480.     */
  481.     tess_(),
  482.     //proc_tess_(process_positions(Echo<Vector2D>(init_points_)),outer_),
  483.     eos_(adiabatic_index),
  484.     rs_(),
  485.     hbc_(rs_),
  486.     //interpm_(eos_,outer_,hbc_,true,false),
  487.     interpm_(),
  488.     raw_point_motion_(),
  489.     point_motion_(raw_point_motion_,hbc_,
  490.           0.15, 0.02, false,0,&outer_),
  491.     alt_point_motion_(),
  492.     gravity_acc_(c.gravitation_constant*c.black_hole_mass,
  493.          0.1*c.parsec,
  494.          c.parsec*Vector2D(-0.05,0)),
  495.   //    veil_(gravity_acc_),
  496.     gravity_force_(gravity_acc_),
  497.     geom_force_(pg_.getAxis()),
  498.     force_(VectorInitializer<SourceTerm*>(&gravity_force_)
  499.        (&geom_force_)()),
  500.     ce_(Ratchet::in),
  501.     sim_(init_points_, //distribute_grid(proc_tess_,Echo<Vector2D>(init_points_)),
  502.      tess_,
  503.      //  proc_tess_,
  504.      interpm_,
  505.      Piecewise
  506.      (Circle(Vector2D(0,-c.offset),c.supernova_radius),
  507.       Uniform2D(c.supernova_density),
  508.       Piecewise
  509.       (Circle(Vector2D(0,0),c.R_b),
  510.        BoundedRadialPowerLaw(Vector2D(0,0),
  511.                  -omega_in,
  512.                  c.inner_density_prefactor,
  513.                  1e-6, 10),
  514.        BoundedRadialPowerLaw(Vector2D(0,0),
  515.                  -omega_out,
  516.                  c.outer_density_prefactor,
  517.                  1e-6, 10))),
  518.      Piecewise(Circle(Vector2D(0,-c.offset),c.supernova_radius),
  519.            Uniform2D(c.supernova_pressure),
  520.            HydrostaticPressure(c.gravitation_constant*c.black_hole_mass,
  521.                        c.R_0,
  522.                        c.rho_0,
  523.                        c.omega_in,
  524.                        c.R_b,
  525.                        c.omega_out)),
  526.      Uniform2D(0),
  527.      Uniform2D(0),
  528.      eos_,
  529.      rs_,
  530.      alt_point_motion_,
  531.      //  point_motion_,
  532.      force_,
  533.      outer_,
  534.      hbc_)
  535.   {
  536.     sim_.SetColdFlows(0.01,0.01);
  537.     sim_.SetDensityFloor(1e-6,1e-23);
  538.     sim_.changePhysicalGeometry(&pg_);
  539.     sim_.custom_evolution_manager.addCustomEvolution(&ce_);
  540.     for(int i=0;i<sim_.GetCellNo();++i){
  541.       if(abs(sim_.GetMeshPoint(i))<0.02*c.parsec)
  542.     sim_.custom_evolution_indices[static_cast<size_t>(i)] = 1;
  543.     }
  544.   }
  545.  
  546.   hdsim& getSim(void)
  547.   {
  548.     return sim_;
  549.   }
  550.  
  551. private:
  552.   const CylindricalSymmetry pg_;
  553.   const SquareBox outer_;
  554.   //  const CartesianGridGenerator cgg_;
  555.   const vector<Vector2D> init_points_;
  556.   VoronoiMesh tess_;
  557.   VoronoiMesh proc_tess_;
  558.   const IdealGas eos_;
  559.   const Hllc rs_;
  560.   const RigidWallHydro hbc_;
  561.   //const FreeFlow hbc_;
  562.   //  LinearGaussConsistent interpm_;
  563.   PCM2D interpm_;
  564.   Lagrangian raw_point_motion_;
  565.   RoundCells point_motion_;
  566.   Eulerian alt_point_motion_;
  567.   CenterGravity gravity_acc_;
  568.   //  Veil veil_;
  569.   ConservativeForce gravity_force_;
  570.   CylindricalComplementary geom_force_;
  571.   SeveralSources force_;
  572.   Ratchet ce_;
  573.   hdsim sim_;
  574. };
  575.  
  576. namespace {
  577.   /*
  578.     class HighestShockedPoint: public TerminationCondition
  579.     {
  580.     public:
  581.  
  582.     HighestShockedPoint(double position):
  583.     position_(position) {}
  584.  
  585.     bool operator()(hdsim const& sim)
  586.     {
  587.     const double mn_thres = 0.01;
  588.     double ys = -10;
  589.     for(int i=0;i<sim.GetCellNo();++i){
  590.     if(sim.GetCell(i).Velocity.y/sim.GetCell(i).SoundSpeed>mn_thres){
  591.     ys = max(ys, sim.GetMeshPoint(i).y);
  592.     }
  593.     }
  594.     cout << ys << endl;
  595.     return ys < position_;
  596.     }
  597.  
  598.     private:
  599.     const double position_;
  600.     };
  601.   */
  602. }
  603.  
  604. class CustomDiagnostic
  605. {
  606. public:
  607.  
  608.   CustomDiagnostic(double dt,
  609.            double t_start):
  610.     dt_(dt),
  611.     t_next_(t_start),
  612.     counter_(0) {}
  613.  
  614.   void operator()(const hdsim& sim)
  615.   {
  616.     if(sim.GetTime()>t_next_){
  617.       write_snapshot_to_hdf5(sim, "snapshot_"+int2str(counter_)+".h5");
  618.  
  619.       t_next_ += dt_;
  620.       ++counter_;
  621.     }
  622.   }
  623.  
  624. private:
  625.   const double dt_;
  626.   double t_next_;
  627.   int counter_;
  628. };
  629.  
  630. namespace {
  631.   /*
  632.     vector<double> calc_decade_vector(double val, size_t decades)
  633.     {
  634.     vector<double> res;
  635.     for(size_t i=0;i<decades;++i)
  636.     res.push_back(val/pow(2.0,static_cast<double>(i)+1));
  637.     return res;
  638.     }
  639.   */
  640. }
  641.  
  642. namespace {
  643.  
  644.   /*
  645.     double highest_shocked_point(const hdsim& sim)
  646.     {
  647.     const double mn_thres = 0.01;
  648.     double ys = -10;
  649.     for(int i=0;i<sim.GetCellNo();++i){
  650.     if(sim.GetCell(i).Velocity.y/sim.GetCell(i).SoundSpeed>mn_thres){
  651.     ys = max(ys,sim.GetMeshPoint(i).y);
  652.     }
  653.     }
  654.     return ys;
  655.     }
  656.   */
  657. }
  658.  
  659. namespace {
  660.   /*
  661.     class MyDiagnostic: public DiagnosticFunction
  662.     {
  663.     public:
  664.  
  665.     MyDiagnostic(double y_start,
  666.     size_t decades):
  667.     snapshot_times_(calc_decade_vector(y_start, decades)),
  668.     y_next_(y_start), counter_(0), active_(true) {}
  669.  
  670.     void operator()(const hdsim& sim)
  671.     {
  672.     if(!active_)
  673.     return;
  674.  
  675.     if(highest_shocked_point(sim)>y_next_){
  676.     write_snapshot_to_hdf5(sim,"snapshot_"+int2str(counter_)+".h5");
  677.     if(static_cast<size_t>(counter_)>snapshot_times_.size()-1)
  678.     active_ = false;
  679.     else{
  680.     y_next_ = snapshot_times_.at(static_cast<size_t>(counter_));
  681.     ++counter_;
  682.     }
  683.     }
  684.     }
  685.  
  686.     private:
  687.     const vector<double> snapshot_times_;
  688.     mutable double y_next_;
  689.     mutable int counter_;
  690.     mutable bool active_;
  691.     };
  692.   */
  693. }
  694.  
  695. namespace {
  696.   class MultipleDiagnostics: public DiagnosticFunction
  697.   {
  698.   public:
  699.  
  700.     MultipleDiagnostics(const vector<DiagnosticFunction*>& diag_list):
  701.       diag_list_(diag_list) {}
  702.  
  703.     void operator()(const hdsim& sim)
  704.     {
  705.       BOOST_FOREACH(DiagnosticFunction* df, diag_list_)
  706.     {
  707.       (*df)(sim);
  708.     }
  709.     }
  710.  
  711.   private:
  712.     const vector<DiagnosticFunction*> diag_list_;
  713.   };
  714. }
  715.  
  716. namespace {
  717.  
  718.   class ShockFrontDetector: public LocalContourCriterion
  719.   {
  720.   public:
  721.     ShockFrontDetector(const Vector2D& center,
  722.                double mn_thres):
  723.       center_(center), mn_thres_(mn_thres) {}
  724.  
  725.     std::pair<bool,Vector2D> operator()(const Edge& edge,
  726.                     const hdsim& sim) const
  727.     {
  728.       if(sim.getHydroBoundaryCondition().IsBoundary(edge,
  729.                             sim.GetTessellation()))
  730.     return std::pair<bool,Vector2D>(false,Vector2D(0,0));
  731.       if(edge.neighbors.first>sim.GetCellNo() ||
  732.      edge.neighbors.second>sim.GetCellNo())
  733.     return std::pair<bool,Vector2D>(false,Vector2D(0,0));
  734.       const std::pair<double,double> mach_numbers
  735.     (calc_mach_number(sim.GetCell(edge.neighbors.first)),
  736.      calc_mach_number(sim.GetCell(edge.neighbors.second)));
  737.       if((mach_numbers.first-mn_thres_)*(mach_numbers.second-mn_thres_)>0 ||
  738.      (mach_numbers.second-mach_numbers.first)*
  739.      (abs(sim.GetMeshPoint(edge.neighbors.second)-center_)-
  740.       abs(sim.GetMeshPoint(edge.neighbors.first)-center_))>0)
  741.     return std::pair<bool,Vector2D>(false,Vector2D(0,0));
  742.       else
  743.     return std::pair<bool,Vector2D>
  744.       (true,find_pos_on_line(sim.GetMeshPoint(edge.neighbors.first),
  745.                  mach_numbers.first,
  746.                  sim.GetMeshPoint(edge.neighbors.second),
  747.                  mach_numbers.second,
  748.                  mn_thres_));
  749.                        
  750.     }
  751.  
  752.   private:
  753.     const Vector2D center_;
  754.     const double mn_thres_;
  755.   };
  756.  
  757.   /*
  758.     class CellEdgesGetter: public Index2Member<Edge>
  759.     {
  760.     public:
  761.     CellEdgesGetter(const Tessellation& tess,
  762.     const int n):
  763.     tess_(tess), edge_indices_(tess_.GetCellEdges(n)) {}
  764.  
  765.     size_t getLength(void) const
  766.     {
  767.     return edge_indices_.size();
  768.     }
  769.  
  770.     Edge operator()(size_t i) const
  771.     {
  772.     return tess_.GetEdge(edge_indices_[i]);
  773.     }
  774.  
  775.     private:
  776.     const Tessellation& tess_;
  777.     const vector<int> edge_indices_;
  778.     };
  779.   */
  780.  
  781.   /*
  782.     class EnforceRigidWall: public Manipulate
  783.     {
  784.     public:
  785.  
  786.     EnforceRigidWall(void) {}
  787.  
  788.     void operator()(hdsim& sim)
  789.     {
  790.     const PhysicalGeometry& pg = sim.getPhysicalGeometry();
  791.     const Tessellation& tess = sim.GetTessellation();
  792.     for(size_t i=0;i<sim.GetAllCells().size();++i){
  793.     if(tess.GetMeshPoint(i).x-lower_left.x<tess.GetWidth(i)){
  794.     sim.getAllConserved()[i].Momentum.x = 0;
  795.     const double volume =
  796.     pg.calcVolume(serial_generate(CellEdgesGetter(tess,i)));
  797.     sim.GetAllCells()[i] = Conserved2Primitive
  798.     (sim.getAllConserved()[i]/volume,sim.getEos());
  799.     }
  800.     }
  801.     }
  802.     };
  803.   */
  804.  
  805.   void my_main_loop(hdsim& sim, const Constants& c)
  806.   {
  807.     //  sim.SetCfl(0.1);
  808.     const double tf = 600*c.year;
  809.     SafeTimeTermination term_cond_raw(tf,1e6);
  810.     ConsecutiveSnapshots diag1(new ConstantTimeInterval(10*units::year),
  811.                    new Rubric("snapshot_",".h5"));
  812.     WriteTime diag2("time.txt");
  813.     SequentialContour diag3(new ConstantTimeInterval(10*units::year),
  814.                 new Rubric("contour_",".h5"),
  815.                 new ShockFrontDetector(Vector2D(0,-c.offset),0.1));
  816.     WriteConserved diag4("conserved.txt");
  817.     vector<DiagnosticFunction*> diag_list;
  818.     diag_list.push_back(&diag1);
  819.     diag_list.push_back(&diag2);
  820.     diag_list.push_back(&diag3);
  821.     diag_list.push_back(&diag4);
  822.     MultipleDiagnostics diag(diag_list);
  823.  
  824.     main_loop(sim,
  825.           term_cond_raw,
  826.           &hdsim::TimeAdvance,
  827.           &diag);
  828.   }
  829. }
  830.  
  831. namespace {
  832.   void report_error(UniversalError const& eo)
  833.   {
  834.     cout << eo.GetErrorMessage() << endl;
  835.     cout.precision(14);
  836.     for(size_t i=0;i<eo.GetFields().size();++i)
  837.       cout << eo.GetFields()[i] << " = "
  838.        << eo.GetValues()[i] << endl;
  839.   }
  840. }
  841.  
  842. int main(void)
  843. {
  844.   feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
  845.   //  MPI_Init(NULL, NULL);
  846.  
  847.   try{
  848.  
  849.     Constants c;
  850.  
  851.     SimData sim_data(c);
  852.     hdsim& sim = sim_data.getSim();
  853.  
  854.     HDF5Shortcut("units.h5")
  855.       ("year",vector<double>(1,c.year))
  856.       ("solar mass",vector<double>(1,c.solar_mass))
  857.       ("parsec",vector<double>(1,c.parsec));
  858.  
  859.     write_snapshot_to_hdf5(sim,
  860.                "initial.h5");
  861.  
  862.     my_main_loop(sim,c);
  863.  
  864.     write_snapshot_to_hdf5(sim,
  865.                "final.h5");
  866.  
  867.   }
  868.   catch(const UniversalError& eo){
  869.     report_error(eo);
  870.     throw;
  871.   }
  872.  
  873.   return 0;
  874. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!