Swarm-NG  1.1
stop_on_crossing_orbit.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 
28 #pragma once
29 
30 #include <limits>
31 
32 namespace swarm { namespace monitors {
33 
42 struct stop_on_crossing_orbit_params {
43  bool deactivate_on, log_on, verbose_on;
46  stop_on_crossing_orbit_params(const config &cfg)
47  {
48  deactivate_on = cfg.optional("deactivate_on_crossing",false);
49  log_on = cfg.optional("log_on_crossing",false);
50  verbose_on = cfg.optional("verbose_on_crossing",false);
51  }
52 };
53 
54 GPUAPI double sqr(const double& d){
55  return d*d;
56 }
57 
65 template<class log_t>
66 class stop_on_crossing_orbit {
67  public:
68  typedef stop_on_crossing_orbit_params params;
69 
70  private:
71  params _params;
72  bool condition_met;
73  ensemble::SystemRef& _sys;
74  log_t& _log;
75 
76 
77  public:
78  template<class T>
79  static GENERIC int thread_per_system(T compile_time_param){
80  return 1;
81  }
82 
83  template<class T>
84  static GENERIC int shmem_per_system(T compile_time_param) {
85  return 0;
86  }
87  GPUAPI bool is_deactivate_on() { return _params.deactivate_on; };
88  GPUAPI bool is_log_on() { return _params.log_on; };
89  GPUAPI bool is_verbose_on() { return _params.verbose_on; };
90  GPUAPI bool is_any_on() { return is_deactivate_on() || is_log_on() || is_verbose_on() ; }
91  GPUAPI bool is_condition_met () { return ( condition_met ); }
92  GPUAPI bool need_to_log_system ()
93  { return (is_log_on() && is_condition_met() ); }
94  GPUAPI bool need_to_deactivate ()
95  { return ( is_deactivate_on() && is_condition_met() ); }
96 
97  GPUAPI void log_system() { log::system(_log, _sys); }
98 
104  GPUAPI void calc_a_e(const int& b,double& a,double& e) {
105  double x,y,z,vx,vy,vz; _sys[b].get(x,y,z,vx,vy,vz);
106  // h2 = ||pos X vel||^2
107  double h2 = sqr(y*vz-z*vy) + sqr(z*vx-x*vz) + sqr(x*vy-y*vx);
108  double _GM = _sys[b].mass();
109  double r = _sys[b].distance_to_origin(), sp = _sys[b].speed();
110  double energy = sp*0.5-_GM/r;
111 
112  a = -0.5*_GM/energy;
113  double fac = 1.-h2/(_GM*a);
114  e = (fac>1.e-8) ? sqrtf(fac) : 0.;
115  }
116 
124  GPUAPI bool check_for_crossing_orbits(const int& i, const int& j) {
125 
126  double a_i, e_i, a_j, e_j;
127  calc_a_e(i, a_i, e_i);
128  calc_a_e(j, a_j, e_j);
129 
130  bool is_orbits_crossing;
131  if(a_i<=a_j)
132  is_orbits_crossing = a_i * (1. + e_i) > a_j * ( 1. - e_j ) ;
133  else
134  is_orbits_crossing = a_i * (1. + e_i) < a_j * ( 1. - e_j ) ;
135 
136  if( is_orbits_crossing && is_verbose_on())
137  lprintf(_log, "Crossing orbits detected: "
138  "sys=%d, T=%lg i=%d j=%d a_i=%lg e_i=%lg a_j=%lg e_j=%lg.\n"
139  , _sys.number(), _sys.time(),i,j, a_i,e_i, a_j, e_j);
140 
141  return is_orbits_crossing;
142  }
143 
144 #if 0
145  // GPUAPI void operator () () {
146  GPUAPI void operator () (int thread_in_system) {
147  if(!is_any_on()) return;
148 
149  if(thread_in_system==0)
150  {
151  condition_met = false;
152  // Check for crossing orbits between every pair of planets
153  for(int j = 2; j < _sys.nbod(); j++)
154  for(int i = 1; i < j; i++)
155  condition_met = condition_met || check_for_crossing_orbits(i, j);
156 
157  if(condition_met) {
158  if(is_log_on())
159  log::system(_log, _sys);
160  if(is_deactivate_on())
161  _sys.set_disabled();
162  }
163  }
164  }
165 #endif
166 
167  GPUAPI void operator () (int thread_in_system)
168  {
169  pass_one(thread_in_system);
170  pass_two(thread_in_system);
171  if(need_to_log_system() && (thread_in_system==0) )
172  log::system(_log, _sys);
173  }
174 
175 
176  GPUAPI bool pass_one (int thread_in_system)
177  {
178  condition_met = false;
179  return true;
180  }
181 
182 
183  GPUAPI int pass_two (int thread_in_system)
184  {
185  if(is_any_on() && (thread_in_system==0) )
186  {
187  // Check for crossing orbits between every pair of planets
188  for(int j = 2; j < _sys.nbod(); j++)
189  for(int i = 1; i < j; i++)
190  condition_met = condition_met || check_for_crossing_orbits(i, j);
191 
192  if(condition_met && is_deactivate_on() )
193  { _sys.set_disabled(); }
194 
195  }
196  return _sys.state();
197  }
198 
199 
200  GPUAPI stop_on_crossing_orbit(const params& p,ensemble::SystemRef& s,log_t& l)
201  :_params(p),_sys(s),_log(l){}
202 
203 };
204 
205 } } // end namespace monitors :: swarm