Swarm-NG  1.1
mvs_cpu.hpp
Go to the documentation of this file.
1 /*************************************************************************
2  * Copyright (C) 2011 by Saleh Dindar and the Swarm-NG Development Team *
3  * *
4  * This program is free software; you can redistribute it and/or modify *
5  * it under the terms of the GNU General Public License as published by *
6  * the Free Software Foundation; either version 3 of the License. *
7  * *
8  * This program is distributed in the hope that it will be useful, *
9  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
11  * GNU General Public License for more details. *
12  * *
13  * You should have received a copy of the GNU General Public License *
14  * along with this program; if not, write to the *
15  * Free Software Foundation, Inc., *
16  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
17  ************************************************************************/
18 
25 #ifndef H_MVS_CPU
26 #define H_MVS_CPU
27 
28 #include "swarm/common.hpp"
29 #include "swarm/integrator.hpp"
30 #include "swarm/plugin.hpp"
31 
33 #define ASSUME_PROPAGATOR_USES_STD_COORDINATES 0
34 
35 namespace swarm { namespace cpu {
51 template< class Monitor >
52 class mvs_cpu : public integrator {
53  typedef integrator base;
54  typedef Monitor monitor_t;
55  typedef typename monitor_t::params mon_params_t;
56  private:
57  double _time_step;
58  mon_params_t _mon_params;
59 
60  // included here so as to avoid namespace conflicts between CPU and OMP integrators
61 #include "../propagators/keplerian.hpp"
62 
63 public:
64  mvs_cpu(const config& cfg): base(cfg),_time_step(0.001), _mon_params(cfg) {
65  _time_step = cfg.require("time_step", 0.0);
66  }
67 
69  virtual void launch_integrator() {
70  for(int i = 0; i < _ens.nsys(); i++){
71  integrate_system(_ens[i]);
72  }
73  }
74 
76  inline static double inner_product(const double a[3],const double b[3]){
77  return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
78  }
79 
81  void calcForces(ensemble::SystemRef& sys, double acc[][3]){
82  const int nbod = sys.nbod();
83 
84  // Clear acc
85  for(int b = 0; b < nbod; b++) for(int c =0; c < 3; c++)
86  acc[b][c] = 0.;
87 
88  // Loop through all pairs
89  for(int i=0; i < nbod-1; i++) for(int j = i+1; j < nbod; j++) {
90 
91  double dx[3] = { sys[j][0].pos()-sys[i][0].pos(),
92  sys[j][1].pos()-sys[i][1].pos(),
93  sys[j][2].pos()-sys[i][2].pos()
94  };
95  double dv[3] = { sys[j][0].vel()-sys[i][0].vel(),
96  sys[j][1].vel()-sys[i][1].vel(),
97  sys[j][2].vel()-sys[i][2].vel()
98  };
99 
100  // Calculated the magnitude
101  double r2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2] * dx[2];
102  double rinv = 1. / ( sqrt(r2) * r2 ) ;
103  double rv = inner_product(dx,dv) * 3. / r2;
104 
105  // Update acc for i
106  const double scalar_i = +rinv*sys[j].mass();
107  for(int c = 0; c < 3; c++) {
108  acc[i][c] += dx[c]* scalar_i;
109  }
110 
111  // Update acc for j
112  const double scalar_j = -rinv*sys[i].mass();
113  for(int c = 0; c < 3; c++) {
114  acc[j][c] += dx[c]* scalar_j;
115  }
116  }
117  }
118 
119 
120  // MVS specific routines
121 
123  void convert_internal_to_std_coord(ensemble::SystemRef sys)
124  { convert_helio_pos_bary_vel_to_std_coord(sys); }
125 
127  GPUAPI void convert_std_to_internal_coord(ensemble::SystemRef sys)
128  { convert_std_to_helio_pos_bary_vel_coord(sys); }
129 
130 
132  void convert_std_to_helio_pos_bary_vel_coord(ensemble::SystemRef sys) {
133  const int nbod = sys.nbod();
134  double pc0;
135  double sump = 0., sumv = 0., mtot = 0.;
136  for(int c=0;c<3;++c)
137  {
138  pc0 = sys[0][c].pos();
139  // Find Center of mass and momentum
140  for(int j=0;j<nbod;++j) {
141  const double mj = sys[j].mass();
142  mtot += mj;
143  sump += mj*sys[j][c].pos();
144  sumv += mj*sys[j][c].vel();
145  }
146  sumv /= mtot;
147  int b = 0; // For sun
148  sys[b][c].vel() = sumv;
149  sys[b][c].pos() = sump/mtot;
150 
151  for(b=1;b<nbod;++b) // For planets
152  {
153  sys[b][c].vel() -= sumv;
154  sys[b][c].pos() -= pc0;
155  }
156  }
157  }
158 
160  void convert_helio_pos_bary_vel_to_std_coord (ensemble::SystemRef sys)
161  {
162  const int nbod = sys.nbod();
163  double sump = 0., sumv = 0., mtot;
164  double m0 = sys[0].mass();
165  double pc0, vc0;
166 
167  for(int c=0;c<3;++c)
168  {
169  pc0 = sys[0][c].pos();
170  vc0 = sys[0][c].vel();
171  mtot = m0;
172  for(int j=1;j<nbod;++j)
173  {
174  const double mj = sys[j].mass();
175  mtot += mj;
176  sump += mj*sys[j][c].pos();
177  sumv += mj*sys[j][c].vel();
178  }
179  sump /= mtot;
180 
181  int b = 0; // For sun only
182  sys[b][c].pos() -= sump;
183  sys[b][c].vel() -= sumv/m0;
184 
185  for(int b=1;b<nbod;++b) // For planets
186  {
187  sys[b][c].pos() += pc0 - sump;
188  sys[b][c].vel() += vc0;
189  }
190 
191  }
192  }
193 
194 
196  void drift_step(ensemble::SystemRef sys, const double hby2)
197  {
198  const double hby2m = hby2/sys[0].mass();
199  const int nbod = sys.nbod();
200  for(int c=0;c<3;++c)
201  {
202  int b = 0;
203  sys[b][c].pos() += hby2*sys[b][c].vel();
204  double mv = 0;
205  for(int j=1;j<nbod;++j)
206  mv += sys[j].mass() * sys[j][c].vel();
207  for(b=1;b<nbod;++b)
208  {
209  sys[b][c].pos() += mv*hby2m;
210  }
211  }
212  }
213 
215  void integrate_system(ensemble::SystemRef sys){
216  const int nbod = sys.nbod();
217  double acc[nbod][3];
218 
219  // Setting up Monitor
220  monitor_t montest(_mon_params,sys,*_log) ;
221 
222  // begin init();
223  const double sqrtGM = sqrt(sys[0].mass());
224  convert_std_to_helio_pos_bary_vel_coord(sys);
225  calcForces(sys,acc);
226  // end init()
227 
228  for(int iter = 0 ; (iter < _max_iterations) && sys.is_active() ; iter ++ )
229  {
230 
231  // begin advance();
232  double hby2 = 0.5 * std::min( _destination_time - sys.time() , _time_step );
233 
234  // Step 1
235  drift_step(sys,hby2);
236 
237  // Step 2: Kick Step
238  for(int b=1;b<nbod;++b)
239  for(int c=0;c<3;++c)
240  sys[b][c].vel() += hby2 * acc[b][c];
241 
242  // __syncthreads();
243 
244  // 3: Kepler Drift Step (Keplerian orbit about sun/central body)
245  for(int b=1;b<nbod;++b)
246  drift_kepler( sys[b][0].pos(),sys[b][1].pos(),sys[b][2].pos(),sys[b][0].vel(),sys[b][1].vel(),sys[b][2].vel(),sqrtGM, 2.0*hby2 );
247  // __syncthreads();
248 
249  // TODO: check for close encounters here
250  calcForces(sys,acc);
251 
252  // Step 4: Kick Step
253  for(int b=1;b<nbod;++b)
254  for(int c=0;c<3;++c)
255  sys[b][c].vel() += hby2 * acc[b][c];
256  // __syncthreads();
257 
258  // Step 5
259  drift_step(sys,hby2);
260 
261  sys.time() += 2.0*hby2;
262 
263  // end advance
264 
265  const int thread_in_system_for_monitor = 0;
266 #if ASSUME_PROPAGATOR_USES_STD_COORDINATES
267  montest( thread_in_system_for_monitor );
268 #else
269  bool using_std_coord = false;
270  bool needs_std_coord = montest.pass_one( thread_in_system_for_monitor );
271 
272  if(needs_std_coord)
273  {
274  convert_internal_to_std_coord(sys);
275  using_std_coord = true;
276  }
277 
278  // __syncthreads();
279  int new_state = montest.pass_two ( thread_in_system_for_monitor );
280 
281  if( montest.need_to_log_system() )
282  { log::system(*_log, sys); }
283 
284  // __syncthreads();
285  if(using_std_coord)
286  {
287  convert_std_to_internal_coord(sys);
288  using_std_coord = false;
289  }
290 #endif
291  // __syncthreads();
292 
293  if( sys.is_active() )
294  {
295  if( sys.time() >= _destination_time )
296  { sys.set_inactive(); }
297  }
298  // __syncthreads();
299 
300  }
301 
302  // shutdown();
303  convert_helio_pos_bary_vel_to_std_coord (sys);
304  }
305 };
306 
307 
308 
309  } } // Close namespaces
310 
311 
312 #endif