hdsim.cpp
1 #include <cmath>
2 #include <algorithm>
3 #include "hdsim.hpp"
4 #include "../../misc/universal_error.hpp"
5 #include "../common/hydrodynamic_variables.hpp"
6 #include "../common/hydrodynamics.hpp"
7 #include "../../misc/utils.hpp"
8 #include "../../misc/lazy_list.hpp"
9 #include "spdlog/spdlog.h"
10 
11 using namespace std;
12 
13 // Diagnostics
14 
15 double hdsim1D::GetTime(void) const
16 {
17  return time_;
18 }
19 
20 int hdsim1D::GetCycle(void) const
21 {
22  return cycle_;
23 }
24 
26 {
27  return ss_;
28 }
29 
30 // External functions
31 
32 namespace {
33 
34  double GetVolume
35  (const vector<double>& v,
36  const PhysicalGeometry1D& pg,
37  size_t i)
38  {
39  return pg.calcVolume(v.at(i+1))
40  - pg.calcVolume(v.at(i));
41  }
42 
43  vector<Extensive> calc_extensives
44  (const PhysicalGeometry1D& pg,
45  const SimulationState1D& ss,
46  const EquationOfState& eos)
47  {
48  const vector<double>& vertices = ss.getVertices();
49  const vector<ComputationalCell>& cells = ss.getCells();
50  vector<Extensive> res(ss.getCells().size());
51  for(size_t i=0;i<res.size();++i){
52  const double volume = GetVolume(vertices,
53  pg,
54  i);
55  res.at(i).mass = cells.at(i).density*volume;
56  res.at(i).momentum = res.at(i).mass*cells.at(i).velocity;
57  const double kinetic_specific_energy =
58  0.5*pow(abs(cells.at(i).velocity),2);
59  const double thermal_specific_energy =
60  eos.dp2e(cells.at(i).density, cells.at(i).pressure);
61  res.at(i).energy = res.at(i).mass*
62  (kinetic_specific_energy+thermal_specific_energy);
63  for(size_t j=0;j<cells.at(0).tracers.size();++j){
64  res.at(i).tracers.push_back
65  (cells.at(i).tracers.at(j)*res.at(i).mass);
66  }
67  }
68 
69  return res;
70  }
71 }
72 
74 (const PhysicalGeometry1D& pg,
75  const SimulationState1D& ss,
76  const EquationOfState& eos,
77  const VertexMotion& vm,
78  const SourceTerm1D& force,
79  const TimeStepFunction1D& tsf,
80  const FluxCalculator1D& fc,
81  const ExtensiveUpdater1D& eu,
82  const CellUpdater1D& cu):
83  pg_(pg),
84  ss_(ss),
85  eos_(eos),
86  extensives_(calc_extensives(pg,ss,eos)),
87  vm_(vm),
88  force_(force),
89  tsf_(tsf),
90  fc_(fc),
91  eu_(eu),
92  cu_(cu),
93  time_(0),
94  cycle_(0),
95  tracers_intensive_(vector<vector<double> >()),
96  tracers_extensive_(vector<vector<double> >())
97 {
98  spdlog::debug("hdsim1D initialisation completed");
99 }
100 
101 namespace {
102  vector<double> CalcVertexVelocities
103  (const SimulationState1D& state,
104  VertexMotion const& vm)
105  {
106  vector<double> res(state.getVertices().size());
107  for(size_t i = 0; i<state.getVertices().size();i++)
108  res[i] = vm(i, state.getVertices(), state.getCells());
109 
110  return res;
111  }
112 
113  vector<double> calc_new_vertices
114  (const vector<double>& vv_,
115  double dt,
116  const vector<double>& vertices)
117  {
118  vector<double> res = vertices;
119  for(size_t i=0;i<res.size();++i)
120  res.at(i) += dt*vv_.at(i);
121  return res;
122  }
123 }
124 
125 namespace {
126 
127  /*
128  Extensive conserved2extensive
129  (const Conserved& c)
130  {
131  Extensive res;
132  res.mass = c.Mass;
133  res.momentum = c.Momentum;
134  res.energy = c.Energy;
135  return res;
136  }
137  */
138 
139  void force_contribution
140  (const SimulationState1D& state,
141  const SourceTerm1D& force,
142  double t,
143  double dt,
144  vector<Extensive>& extensive)
145  {
146  for(size_t i=0;i<extensive.size();++i)
147  extensive[i] +=
148  dt*force(state, i, t, dt);
149  }
150 }
151 
153 {
154  spdlog::debug("begin time advance 1o iteration {0}, virtual time {1}",
155  cycle_, time_);
156 
157  const vector<double> _VertexVelocity = CalcVertexVelocities
158  (ss_, vm_);
159 
160  spdlog::debug("Vertex velocity calculated");
161 
162  const double dt = tsf_(ss_,eos_);
163 
164  spdlog::debug("Time step calculated");
165 
166  const vector<Extensive> fluxes =
167  fc_(ss_, _VertexVelocity, eos_, dt);
168 
169  spdlog::debug("Fluxes calculated");
170 
171  eu_
172  (fluxes,
173  pg_,
174  ss_,
175  dt,
176  extensives_);
177 
178  spdlog::debug("Extensives updated");
179 
180  force_contribution(ss_,
181  force_, time_, dt,
182  extensives_);
183 
184  spdlog::debug("source term calculated");
185 
186  ss_.updateVertices(calc_new_vertices(_VertexVelocity,
187  dt,
188  ss_.getVertices()));
189 
190  spdlog::debug("Vertices updated");
191 
192  ss_.updateCells(cu_(pg_,
193  extensives_,
194  ss_,
195  eos_));
196 
197  spdlog::debug("cells updated");
198 
199  time_ += dt;
200  cycle_++;
201 
202  spdlog::debug("time advance interation finished");
203 }
204 
206 {
207  const vector<double> mid_vertex_velocities =
208  CalcVertexVelocities(ss_, vm_);
209 
210  const double dt = tsf_(ss_, eos_);
211 
212  const vector<Extensive> mid_fluxes =
213  fc_(ss_, mid_vertex_velocities, eos_, dt);
214 
215  vector<Extensive> mid_extensive = extensives_;
216 
217  eu_(mid_fluxes,
218  pg_,
219  ss_,
220  dt/2,
221  mid_extensive);
222 
223  force_contribution(ss_,
224  force_, time_, dt/2,
225  mid_extensive);
226 
227  SimulationState1D mid_state = ss_;
228  mid_state.updateVertices
229  (calc_new_vertices(mid_vertex_velocities,
230  dt,
231  ss_.getVertices()));
232 
233  /*
234  const vector<Conserved> mid_intesive =
235  UpdateConservedIntensive
236  (extensive2conserved(mid_extensive),
237  mid_state.getVertices(),
238  pg_);
239  */
240  mid_state.updateCells
241  (cu_(pg_,
242  mid_extensive,
243  mid_state,
244  eos_));
245 
246  const vector<double> _VertexVelocity = CalcVertexVelocities
247  (mid_state, vm_);
248 
249  const vector<Extensive> fluxes =
250  fc_(mid_state, _VertexVelocity, eos_, dt);
251 
252  eu_(fluxes,
253  pg_,
254  ss_,
255  dt,
256  extensives_);
257 
258  force_contribution
259  (mid_state,
260  force_, time_, dt,
261  extensives_);
262 
263  ss_.updateVertices(calc_new_vertices(_VertexVelocity,
264  dt,
265  ss_.getVertices()));
266 
267  ss_.updateCells(cu_(pg_,
268  extensives_,
269  ss_,
270  eos_));
271 
272  time_ += dt;
273  ++cycle_;
274 }
Base class for a flux calculator.
const vector< double > & getVertices(void) const
Access to vertices.
Base class for physical geometry.
double GetTime(void) const
Returns the time of the simulation.
Definition: hdsim.cpp:15
Base class for vertex motion.
int GetCycle(void) const
Returns the number of times time advance was called.
Definition: hdsim.cpp:20
void updateCells(const vector< ComputationalCell > &cells)
Updates hydro cellls.
hdsim1D(const PhysicalGeometry1D &pg, const SimulationState1D &ss, const EquationOfState &eos, const VertexMotion &vm, const SourceTerm1D &force, const TimeStepFunction1D &tsf, const FluxCalculator1D &fc, const ExtensiveUpdater1D &eu, const CellUpdater1D &cu)
Class constructor.
Definition: hdsim.cpp:74
void TimeAdvance2(void)
Second order time advance.
Definition: hdsim.cpp:205
void TimeAdvance(void)
Advances the simulation in time.
Definition: hdsim.cpp:152
const vector< ComputationalCell > & getCells(void) const
Access to hydro cells.
Base class for equation of state.
Abstract class for external forces.
Abstract class for the cell update scheme.
Method for updating the extensive variables.
virtual double dp2e(double d, double p, tvector const &tracers=tvector(), vector< string > const &tracernames=vector< string >()) const =0
Calculates the thermal energy per unit mass.
virtual double calcVolume(double radius) const =0
Calculates the volume.
Base class for a time step calculator.
double abs(Vector3D const &v)
Norm of a vector.
Definition: Vector3D.cpp:44
const SimulationState1D & getState(void) const
Access to computational domain and hydro cells.
Definition: hdsim.cpp:25
Package for computational domain and hydro cells.
void updateVertices(const vector< double > &vertices)
Updates positions of vertices.
One dimensional hydrodynamic simulation.