Swarm-NG  1.1
log_transit.hpp
Go to the documentation of this file.
1 /*************************************************************************
2  * Copyright (C) 2011 by Eric Ford 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 
28 #pragma once
29 
30 #include "swarm/gpu/gravitation_accjerk.hpp"
31 
32 #define USE_WORKS 1
33 
34 namespace swarm {
35  namespace monitors {
36 
42 struct log_transit_params {
43  double tol, dt;
44  int nbod, num_max_iter;
45  bool log_on_transit, log_on_occultation;
46 
48  log_transit_params(const config &cfg)
49  {
50  nbod = cfg.require<int>("nbod");
51  num_max_iter = cfg.optional("num_max_transit_iter", 2);
52  tol = cfg.optional("log_transit_tol", 2.e-8);
53  // WARNING assumes Hermite Fixed
54  dt = cfg.require("time_step", 0.0);
55  log_on_transit = cfg.optional("log_transits", true);
56  log_on_occultation = cfg.optional("log_occultations", false);
57 
58  }
59 };
60 
73 template<class log_t>
74 class log_transit {
75  public:
76  typedef log_transit_params params;
77 
78  private:
79  params _params;
80  bool condition_met;
81 
82  ensemble::SystemRef& _sys;
83 
84  log_t& _log;
85 
86  public:
87  template<class T>
88  static GENERIC int thread_per_system(T compile_time_param){
89  return 1;
90  }
91 
92  template<class T>
93  static GENERIC int shmem_per_system(T compile_time_param) {
94  return 0;
95  }
96  GPUAPI bool is_deactivate_on() { return false; };
97  GPUAPI bool is_log_on() { return _params.tol!=0.; };
98  GPUAPI bool is_verbose_on() { return false; };
99  GPUAPI bool is_any_on() { return is_deactivate_on() || is_log_on() || is_verbose_on() ; }
100  GPUAPI bool is_condition_met () { return ( condition_met ); }
101  GPUAPI bool need_to_log_system ()
102  { return (is_log_on() && is_condition_met() ); }
103  GPUAPI bool need_to_deactivate ()
104  { return ( is_deactivate_on() && is_condition_met() ); }
105 
106  GPUAPI void log_system() { log::system(_log, _sys); }
107 
108  GPUAPI void operator () (const int thread_in_system)
109  {
110  pass_one(thread_in_system);
111  pass_two(thread_in_system);
112  if(need_to_log_system() && (thread_in_system==0) )
113  log_system();
114  }
115 
116  GPUAPI bool pass_one (const int thread_in_system)
117  {
118  condition_met = false;
119  return true;
120  }
121 
122 
123  GPUAPI double square(const double x) { return x*x; }
124 
125 
127  template<int nbod>
128  GPUAPI void calc_transit_time(const int& thread_in_system, const int& i,const int& j, const double& dt, double dx[2], double dv[2], const double& b2begin, double& db2dt, const double& pos_step_end, const double& vel_step_end, double& dt_min_b2,double& b,double & vproj)
129  {
130  extern __shared__ char shared_mem[];
131  typedef swarm::compile_time_params_t<nbod> par_t;
132  typedef swarm::gpu::bppt::GravitationAccJerk<par_t> calcForces_t;
133  typedef typename calcForces_t::shared_data grav_t;
134 
136  calcForces_t calcForces(_sys,*( (grav_t*) (&shared_mem[(swarm::gpu::bppt::sysid_in_block()%SHMEM_CHUNK_SIZE)*sizeof(double)+(swarm::gpu::bppt::sysid_in_block()/SHMEM_CHUNK_SIZE)*SHMEM_CHUNK_SIZE*(nbod*(nbod-1)*3*sizeof(double))]) ) );
137 
138  const int bid = swarm::gpu::bppt::thread_body_idx(nbod);
139  const int cid = swarm::gpu::bppt::thread_component_idx(nbod);
140  double pos = pos_step_end, vel = vel_step_end;
141  double acc, jerk;
142 #if USING_INTEGRATOR_THAT_OVER_WRITES_SHARED_DATA || 1
143  calcForces(thread_in_system,bid,cid,pos,vel,acc,jerk);
144  __syncthreads();
145  /*
146  if(thread_in_system < calcForces_t::pair_count)
147  calcForces.calc_pair(thread_in_system);
148  __syncthreads(); */
149 #endif
150  double acc_ix, acc_jx, acc_iy, acc_jy, jerk_ix, jerk_jx, jerk_iy, jerk_jy;
151  calcForces.sum(i,0,acc_ix,jerk_ix);
152  calcForces.sum(i,1,acc_iy,jerk_iy);
153  calcForces.sum(j,0,acc_jx,jerk_jx);
154  calcForces.sum(j,1,acc_jy,jerk_jy);
155  double da[2] = { acc_jx - acc_ix, acc_jy- acc_iy };
156  double dj[2] = { jerk_jx-jerk_ix, jerk_jy-jerk_iy };
157  double d2b2dt2 = 2.*(dv[0]*dv[0]+dv[1]*dv[1]+dx[0]*da[0]+dx[1]*da[1]);
158  double d3b2dt3 = 6.*(dv[0]*da[0]+dv[1]*da[1])+2.*(dx[0]*dj[0]+dx[1]*dj[1]);
159  double dt_try = -db2dt/d2b2dt2; // find vertex of parabola
160  dt_try = -db2dt/(d2b2dt2+0.5*dt_try*d3b2dt3);
161 
162  // printf("begin: b2=%g dx[0]=%g dx[1]=%g dv[0]=%g dv[1]=%g dt_try=%g\n",b2begin,dx[0],dx[1],dv[0],dv[1],dt_try);
163  // if((dt_try<0.) || (dt_try>dt) ) // ignore if already past or vertex past next step
164  if((dt_try<=-0.5*dt) || (dt_try>0.5*dt) ) // ignore if already past or vertex past next step
165  { dt_min_b2 = dt_try; b=-1.; vproj = -1.; return; }
166  double dt_cum = dt_try;
167  double b2;
168  for(int iter=0;iter<_params.num_max_iter;++iter)
169  {
170  // take trial step towards transit mid-point
171  pos = pos + dt_try*(vel+(dt_try*0.5)*(acc+(dt_try/3.0)*jerk));
172  vel = vel + dt_try*(acc+(dt_try*0.5)*jerk);
173  calcForces(thread_in_system,bid,cid,pos,vel,acc,jerk);
174  dx[0] = _sys[j][0].pos()-_sys[i][0].pos();
175  dx[1] = _sys[j][1].pos()-_sys[i][1].pos();
176  dv[0] = _sys[j][0].vel()-_sys[i][0].vel();
177  dv[1] = _sys[j][1].vel()-_sys[i][1].vel();
178  b2 = dx[0]*dx[0]+dx[1]*dx[1];
179  db2dt = 2.*(dx[0]*dv[0]+dx[1]*dv[1]);
180  calcForces.sum(i,0,acc_ix,jerk_ix);
181  calcForces.sum(i,1,acc_iy,jerk_iy);
182  calcForces.sum(j,0,acc_jx,jerk_jx);
183  calcForces.sum(j,1,acc_jy,jerk_jy);
184  da[0] = acc_jx - acc_ix; da[1] = acc_jy- acc_iy;
185  dj[0] = jerk_jx-jerk_ix, dj[1] = jerk_jy-jerk_iy;
186  d2b2dt2 = 2.*(dv[0]*dv[0]+dv[1]*dv[1]+dx[0]*da[0]+dx[1]*da[1]);
187  d3b2dt3 = 6.*(dv[0]*da[0]+dv[1]*da[1])+2.*(dx[0]*dj[0]+dx[1]*dj[1]);
188  double dt_try2 = -db2dt/d2b2dt2; // find vertex of parabola
189  dt_try2 = -db2dt/(d2b2dt2+0.5*dt_try2*d3b2dt3);
190  // dt_try += dt_try2;
191  dt_cum += dt_try2;
192  dt_try = dt_try2;
193  if(dt_try2>_params.tol) continue;
194  b2 += dt_try2*(db2dt+0.5*dt_try2*(d2b2dt2+dt_try2*d3b2dt3/3.));
195  dv[0] += dt_try2*(da[0]+0.5*dt_try2*dj[0]);
196  dv[1] += dt_try2*(da[1]+0.5*dt_try2*dj[1]);
197  if(dt_try2<=_params.tol) break;
198  }
199  dt_try = dt_cum;
200  // if((dt_try<0.) || (dt_try>dt) ) // ignore if already past or vertex past next step
201  if((dt_try<=-0.5*dt) || (dt_try>0.5*dt) ) // ignore if already past or vertex past next step
202  { dt_min_b2 = dt_try; b=-1.; vproj = -1.; return; }
203 
204  dt_min_b2 = 0.;
205  double min_b2 = b2begin;
206  if(b2<min_b2) { min_b2 = b2; dt_min_b2 = dt_try; }
207  b = (min_b2>=0) ? sqrt(min_b2) : -sqrt(-min_b2);
208  vproj = sqrt(dv[0]*dv[0]+dv[1]*dv[1]);
209  }
210 
211 
213  template<int nbod>
214  GPUAPI void calc_transit_time_works(const int& i,const int& j, const double& dt, double dx[2], double dv[2], const double& b2begin, double& db2dt,double& dt_min_b2,double& b,double & vproj)
215  {
216  extern __shared__ char shared_mem[];
217  typedef swarm::compile_time_params_t<nbod> par_t;
218  typedef swarm::gpu::bppt::GravitationAccJerk<par_t> calcForces_t;
219  typedef typename calcForces_t::shared_data grav_t;
220  calcForces_t calcForces(_sys,*( (grav_t*) (&shared_mem[(swarm::gpu::bppt::sysid_in_block()%SHMEM_CHUNK_SIZE)*sizeof(double)+(swarm::gpu::bppt::sysid_in_block()/SHMEM_CHUNK_SIZE)*SHMEM_CHUNK_SIZE*(nbod*(nbod-1)*3*sizeof(double))]) ) );
221 
222  // if integrator could have overwritten data in shared memory, need new call to calcforces
223  double acc_ix, acc_jx, acc_iy, acc_jy, jerk_ix, jerk_jx, jerk_iy, jerk_jy;
224  calcForces.sum(i,0,acc_ix,jerk_ix);
225  calcForces.sum(i,1,acc_iy,jerk_iy);
226  calcForces.sum(j,0,acc_jx,jerk_jx);
227  calcForces.sum(j,1,acc_jy,jerk_jy);
228  double da[2] = { acc_jx - acc_ix, acc_jy- acc_iy };
229  double dj[2] = { jerk_jx-jerk_ix, jerk_jy-jerk_iy };
230  double d2b2dt2 = 2.*(dv[0]*dv[0]+dv[1]*dv[1]+dx[0]*da[0]+dx[1]*da[1]);
231  double d3b2dt3 = 6.*(dv[0]*da[0]+dv[1]*da[1])+2.*(dx[0]*dj[0]+dx[1]*dj[1]);
232  double dt_try = -db2dt/d2b2dt2; // find vertex of parabola
233  dt_try = -db2dt/(d2b2dt2+0.5*dt_try*d3b2dt3);
234 
235  // printf("begin: i=%d j=%d b2=%g dx[0]=%g dx[1]=%g dv[0]=%g dv[1]=%g dt_try=%g (works)\n",i,j,b2begin,dx[0],dx[1],dv[0],dv[1],dt_try);
236  // if((dt_try<0.) || (dt_try>dt) ) // ignore if already past or vertex past next step
237  if((dt_try<-0.5*dt) || (dt_try>0.5*dt) ) // ignore if already past or vertex past next step
238  { dt_min_b2 = dt_try; b=-1.; vproj = -1.; return; }
239 
240  dt_min_b2 = 0.;
241  double min_b2 = b2begin;
242  // double b2 = b2begin + dt_try*(db2dt + dt_try*0.5*(d2b2dt2+dt_try*d3b2dt3/3.));
243  dx[0] += dt_try*(dv[0]+dt_try*0.5*(da[0]+dt_try*dj[0]/3.));
244  dx[1] += dt_try*(dv[1]+dt_try*0.5*(da[1]+dt_try*dj[1]/3.));
245  double b2 = (dx[0]*dx[0]+dx[1]*dx[1]);
246  if(b2<min_b2) { min_b2 = b2; dt_min_b2 = dt_try; }
247  // if(min_b2>square((radi+radj)*(1.+_params.tol))) return false;
248  b = (min_b2>0) ? sqrt(min_b2) : -sqrt(-min_b2);
249  dv[0] += dt_min_b2*(da[0]+0.5*dt_min_b2*dj[0]);
250  dv[1] += dt_min_b2*(da[1]+0.5*dt_min_b2*dj[1]);
251  vproj = sqrt(dv[0]*dv[0]+dv[1]*dv[1]);
252  // printf("end: i=%d j=%d b2=%g dx[0]=%g dx[1]=%g dv[0]=%g dv[1]=%g dt_try=%g\n",i,j,b2,dx[0],dx[1],dv[0],dv[1],dt_try);
253  }
254 
255 
256  GPUAPI bool check_in_transit(const int& thread_in_system, const int& i, const int& j, const double dt)
257  {
258  double depth = _sys[j][2].pos()-_sys[i][2].pos();
259  if(!_params.log_on_occultation && (depth < 0.)) return false;
260  if(!_params.log_on_transit && (depth > 0.)) return false;
261 
262  double dx[2] = { _sys[j][0].pos()-_sys[i][0].pos(), _sys[j][1].pos()-_sys[i][1].pos() };
263  double dv[2] = { _sys[j][0].vel()-_sys[i][0].vel(), _sys[j][1].vel()-_sys[i][1].vel() };
264  double b2begin = dx[0]*dx[0]+dx[1]*dx[1];
265  double db2dt = 2.*(dx[0]*dv[0]+dx[1]*dv[1]);
266 
267  double radi = (_sys.num_body_attributes()>=1) ? _sys[i].attribute(0) : 0.;
268  double radj = (_sys.num_body_attributes()>=1) ? _sys[j].attribute(0) : 0.;
269  // if( min(b2begin,b2begin+dt*db2dt)>square((radi+radj)*(1.+_params.tol)) )
270  if( min(b2begin-0.5*dt*db2dt,b2begin+0.5*dt*db2dt)>square((radi+radj)*(1.+_params.tol)) )
271  return false;
272 
273  // store values from actual integration to return to after finding transit time
274  const int nbod = _params.nbod;
275  const int bid = swarm::gpu::bppt::thread_body_idx(nbod);
276  const int cid = swarm::gpu::bppt::thread_component_idx(nbod);
277  double pos_step_end, vel_step_end;
278  if( (bid < nbod) && (cid < 3) ) {
279  pos_step_end = _sys[bid][cid].pos();
280  vel_step_end = _sys[bid][cid].vel();
281  }
282 
283  // WARNING: hard coded to assume Hermite Fixed, nbod=3..8 and Gravitation
284  double dt_min_b2, b, vproj;
285  if((nbod==3) && (nbod<=MAX_NBODIES))
286  {
287  const int nbod_hardwired = 3;
288 #if USE_WORKS
289  if(thread_in_system==0)
290  calc_transit_time_works<nbod_hardwired> (i,j,dt,dx,dv,b2begin,db2dt,dt_min_b2,b,vproj);
291 #else
292  calc_transit_time<nbod_hardwired> (thread_in_system,i,j,dt,dx,dv,b2begin,db2dt,pos_step_end,vel_step_end,dt_min_b2,b,vproj);
293 #endif
294  }
295 #if 1 // Can be set to zero to reduce compile times when testing
296  else if((nbod==4) && (nbod<=MAX_NBODIES))
297  {
298  const int nbod_hardwired = 4;
299 #if USE_WORKS
300  if(thread_in_system==0)
301  calc_transit_time_works<nbod_hardwired> (i,j,dt,dx,dv,b2begin,db2dt,dt_min_b2,b,vproj);
302 #else
303  calc_transit_time<nbod_hardwired> (thread_in_system,i,j,dt,dx,dv,b2begin,db2dt,pos_step_end,vel_step_end,dt_min_b2,b,vproj);
304 #endif
305  }
306  else if(nbod==5)
307  {
308  const int nbod_hardwired = 5;
309 #if USE_WORKS
310  if(thread_in_system==0)
311  calc_transit_time_works<nbod_hardwired> (i,j,dt,dx,dv,b2begin,db2dt,dt_min_b2,b,vproj);
312 #else
313  calc_transit_time<nbod_hardwired> (thread_in_system,i,j,dt,dx,dv,b2begin,db2dt,pos_step_end,vel_step_end,dt_min_b2,b,vproj);
314 #endif
315  }
316  else if((nbod==6) && (nbod<=MAX_NBODIES))
317  {
318  const int nbod_hardwired = 6;
319 #if USE_WORKS
320  if(thread_in_system==0)
321  calc_transit_time_works<nbod_hardwired> (i,j,dt,dx,dv,b2begin,db2dt,dt_min_b2,b,vproj);
322 #else
323  calc_transit_time<nbod_hardwired> (thread_in_system,i,j,dt,dx,dv,b2begin,db2dt,pos_step_end,vel_step_end,dt_min_b2,b,vproj);
324 #endif
325  }
326  else if((nbod==7) && (nbod<=MAX_NBODIES))
327  {
328  const int nbod_hardwired = 7;
329 #if USE_WORKS
330  if(thread_in_system==0)
331  calc_transit_time_works<nbod_hardwired> (i,j,dt,dx,dv,b2begin,db2dt,dt_min_b2,b,vproj);
332 #else
333  calc_transit_time<nbod_hardwired> (thread_in_system,i,j,dt,dx,dv,b2begin,db2dt,pos_step_end,vel_step_end,dt_min_b2,b,vproj);
334 #endif
335  }
336  else if((nbod==8) && (nbod<=MAX_NBODIES))
337  {
338  const int nbod_hardwired = 8;
339 #if USE_WORKS
340  if(thread_in_system==0)
341  calc_transit_time_works<nbod_hardwired> (i,j,dt,dx,dv,b2begin,db2dt,dt_min_b2,b,vproj);
342 #else
343  calc_transit_time<nbod_hardwired> (thread_in_system,i,j,dt,dx,dv,b2begin,db2dt,pos_step_end,vel_step_end,dt_min_b2,b,vproj);
344 #endif
345  }
346 #endif
347 #if 0
348  if(vproj<0.)
349  { // ignore if already past or vertex past next step
350 #if 0
351  if( (bid < nbod) && (cid < 3) ) {
352  _sys[bid][cid].pos() = pos_step_end;
353  _sys[bid][cid].vel() = vel_step_end;
354  }
355 #endif
356  return false;
357  }
358 #endif
359  if((thread_in_system==0) && (vproj>=0.) )
360  {
361  if(_sys.num_body_attributes()>=1)
362  {
363  if (depth>0.) { b /= radi; vproj /= radi; }
364  else { b /= radj; vproj /= radj; }
365  }
366 
367  if((dt_min_b2>=-0.5*dt)&&(dt_min_b2<0.5*dt))
368  {
369  double tout = _sys.time()+dt_min_b2;
370  int event_id = (depth>0.) ? log::EVT_TRANSIT : log::EVT_OCCULTATION;
371  log::event(_log,event_id,tout,_sys.id(),j,b,vproj);
372  // printf("loging: event_id=%d time=%g dt_min_b2=%g time_tr=%g sys=%d bod=%d b=%g v=%g\n",event_id,_sys.time(),dt_min_b2,tout,_sys.id(),j,b,vproj);
373  }
374  }
375 
376  // restore values from main integration
377  if( (bid < nbod) && (cid < 3) ) {
378  _sys[bid][cid].pos() = pos_step_end;
379  _sys[bid][cid].vel() = vel_step_end;
380  }
381 
382  if((thread_in_system==0) && (vproj>=0.) )
383  return true;
384  else
385  return false;
386  }
387 
388 
389  GPUAPI int pass_two (const int thread_in_system)
390  {
391  if(is_any_on())
392  {
393  // Chcek for close encounters
394  for(int b = 1; b < _sys.nbod(); b++)
395  {
396  condition_met = condition_met ||
397  check_in_transit(thread_in_system,0,b,_params.dt);
398 
399  }
400  }
401  return condition_met;
402  }
403 
404 
405 
406  GPUAPI log_transit(const params& p,ensemble::SystemRef& s,log_t& l)
407  :_params(p),_sys(s),_log(l) {}
408 
409 };
410 
411 }
412 
413 
414 }