arepo_interp.cpp
1 #include "arepo_interp.hpp"
2 #include "../common/hydrodynamics.hpp"
3 #include <cmath>
4 
5 using namespace std;
6 
8  eos_(eos) {}
9 
10 namespace {
11  Primitive Derivative
12  (vector<double> const& grid,
13  vector<Primitive> const& cells,
14  size_t idx)
15  {
16  return (cells[size_t(idx)]-cells[size_t(idx-1)])/
17  (0.5*(grid[size_t(idx+1)]-grid[size_t(idx-1)]));
18  }
19 
20  Primitive SymmetricDerivative
21  (vector<double> const& grid,
22  vector<Primitive> const& cells,
23  size_t idx)
24  {
25  const double xr = 0.5*(grid[size_t(idx+2)]+grid[size_t(idx+1)]);
26  const double xl = 0.5*(grid[size_t(idx)]+grid[size_t(idx-1)]);
27  return (cells[size_t(idx+1)]-cells[size_t(idx-1)])/(xr-xl);
28  }
29 
30  bool effectively_zero(double x)
31  {
32  return abs(x)<1e-14;
33  }
34 
35  double sgn(double x)
36  {
37  if(x>0)
38  return 1;
39  else if(x<0)
40  return -1;
41  else if(effectively_zero(x))
42  return 0;
43  else
44  throw UniversalError("NaN detected in sgn");
45  }
46 
47  double my_min(double x, double y, double z)
48  {
49  return min(min(x,y),z);
50  }
51 
52  double minmod(double x, double y, double z)
53  {
54  return 0.125*(sgn(x)+sgn(y))*
55  (sgn(y)+sgn(z))*
56  (sgn(z)+sgn(x))*
57  my_min(abs(x),abs(y),abs(z));
58  }
59 
60  Primitive minmod(Primitive const& p1,
61  Primitive const& p2,
62  Primitive const& p3)
63  {
64  Primitive res;
65  for(int i=0;i<res.GetVarNo();++i)
66  res[i] = minmod(p1[i],p2[i],p3[i]);
67  return res;
68  }
69 }
70 
71 Primitive ArepoInterp::operator()
72  (vector<double> const& vp,
73  vector<Primitive> const& hv,
74  double v_i,
75  size_t idx, int dir, double dt) const
76 {
77  if(dir==0){
78  if(idx==1)
79  return hv[0];
80  else if(idx==vp.size()-1)
81  return hv[vp.size()-2];
82  else if(idx==0)
83  throw UniversalError("Error in PLM1D::InterpState\n Target of interpolation lies outside grid, on the left side");
84  else{
85  const Primitive slope_left = Derivative(vp,hv,idx-1);
86  const Primitive slope_center = Derivative(vp,hv,idx);
87  const Primitive slope_symmetric = SymmetricDerivative(vp,hv,idx-1);
88  const Primitive slope = minmod(slope_left,slope_center,2*slope_symmetric);
89  const double dx = vp[size_t(idx)]-vp[size_t(idx-1)];
90  const double edge_density =
91  hv[size_t(idx-1)].Density
92  +slope.Density*dx/2
93  -(dt/2)*slope.Density*(hv[size_t(idx-1)].Velocity.x-v_i)
94  -(dt/2)*slope.Velocity.x*hv[size_t(idx-1)].Density;
95  const double g = eos_.getAdiabaticIndex();
96  const double edge_pressure =
97  hv[size_t(idx-1)].Pressure
98  +slope.Pressure*dx/2
99  -(dt/2)*slope.Velocity.x*g*hv[size_t(idx-1)].Pressure
100  -(dt/2)*slope.Pressure*(hv[size_t(idx-1)].Velocity.x-v_i);
101  const double edge_xvelocity =
102  hv[size_t(idx-1)].Velocity.x
103  +(dx/2)*slope.Velocity.x
104  -(dt/2)*slope.Velocity.x*(hv[size_t(idx-1)].Velocity.x-v_i)
105  -(dt/2)*slope.Pressure/hv[size_t(idx-1)].Density;
106  const double edge_yvelocity =
107  hv[size_t(idx-1)].Velocity.y
108  +(dx/2)*slope.Velocity.y
109  -(dt/2)*slope.Velocity.y*(hv[size_t(idx-1)].Velocity.x-v_i);
110  return CalcPrimitive
111  (edge_density,
112  edge_pressure,
113  Vector2D(edge_xvelocity,
114  edge_yvelocity),
115  eos_);
116  }
117  }
118  else if(dir==1){
119  if(idx==vp.size()-2)
120  return hv[vp.size()-2];
121  else if(idx==0)
122  return hv[0];
123  else if(idx==vp.size()-1)
124  throw UniversalError("Error in PLM1D::InterpState\n Target of interpolation lies outside grid, on the right side");
125  else{
126  const Primitive slope_right = Derivative(vp,hv,idx+1);
127  const Primitive slope_center = Derivative(vp,hv,idx);
128  const Primitive slope_symmetric = SymmetricDerivative(vp,hv,idx);
129  const Primitive slope = minmod(slope_center,slope_right,2*slope_symmetric);
130  const double dx = vp[size_t(idx+1)]-vp[size_t(idx)];
131  const double edge_density =
132  hv[size_t(idx)].Density
133  -slope.Density*dx/2
134  -(dt/2)*slope.Density*(hv[size_t(idx)].Velocity.x-v_i)
135  -(dt/2)*slope.Velocity.x*hv[size_t(idx)].Density;
136  const double g = eos_.getAdiabaticIndex();
137  const double edge_pressure =
138  hv[size_t(idx)].Pressure
139  -slope.Pressure*dx/2
140  -(dt/2)*slope.Velocity.x*g*hv[size_t(idx)].Pressure
141  -(dt/2)*slope.Pressure*(hv[size_t(idx)].Velocity.x-v_i);
142  const double edge_xvelocity =
143  hv[size_t(idx)].Velocity.x
144  -(dx/2)*slope.Velocity.x
145  -(dt/2)*slope.Velocity.x*(hv[size_t(idx)].Velocity.x-v_i)
146  -(dt/2)*slope.Pressure/hv[size_t(idx)].Density;
147  const double edge_yvelocity =
148  hv[size_t(idx)].Velocity.y
149  -(dx/2)*slope.Velocity.y
150  -(dt/2)*slope.Velocity.y*(hv[size_t(idx)].Velocity.x-v_i);
151  return CalcPrimitive
152  (edge_density,
153  edge_pressure,
154  Vector2D(edge_xvelocity,
155  edge_yvelocity),
156  eos_);
157  }
158  }
159  else
160  throw UniversalError("Unknown direction in plm1d");
161 }
double Density
Density.
Ideal gas equation of state.
Definition: ideal_gas.hpp:13
int GetVarNo(void) const
Returns the numbers of members.
Container for error reports.
Primitive CalcPrimitive(double density, double pressure, Vector2D const &velocity, EquationOfState const &eos)
Calculates the primitive variables.
Interpolation described in arepo paper.
Vector2D Velocity
Velocity.
double getAdiabaticIndex(void) const
Returns the adiabatic index.
Definition: ideal_gas.cpp:8
double y
Component in the y direction.
Definition: geometry.hpp:48
ArepoInterp(IdealGas const &eos)
Class constructor.
Definition: arepo_interp.cpp:7
double min(vector< double > const &v)
Returns the minimal term in a vector.
Definition: utils.cpp:44
double Pressure
Pressure.
double abs(Vector3D const &v)
Norm of a vector.
Definition: Vector3D.cpp:44
Primitive hydrodynamic variables.
2D Mathematical vector
Definition: geometry.hpp:15
double x
Component in the x direction.
Definition: geometry.hpp:45