Swarm-NG  1.1
mvs.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 #include "swarm/common.hpp"
26 #include "swarm/swarmplugin.h"
27 #include "keplerian.hpp"
28 
29 namespace swarm {
30 
31 namespace gpu {
32 namespace bppt {
33 
38 struct MVSPropagatorParams {
39  double time_step;
41  MVSPropagatorParams(const config& cfg){
42  time_step = cfg.require("time_step", 0.0);
43  }
44 };
45 
51 template<class T,class Gravitation>
52 struct MVSPropagator {
53  typedef MVSPropagatorParams params;
54  static const int nbod = T::n;
55 
56  params _params;
57 
58 
60  ensemble::SystemRef& sys;
61  Gravitation& calcForces;
62  int b;
63  int c;
64  int ij;
65 
66  double sqrtGM;
67  double max_timestep;
68 
69  double acc_bc;
70 
72  GPUAPI MVSPropagator(const params& p,ensemble::SystemRef& s,
73  Gravitation& calc)
74  :_params(p),sys(s),calcForces(calc){}
75 
76  __device__ bool is_in_body_component_grid()
77  { return ((b < nbod) && (c < 3)); }
78 
79  __device__ bool is_in_body_component_grid_no_star()
80  { return ( (b!=0) && (b < nbod) && (c < 3) ); }
81 
82  __device__ bool is_first_thread_in_system()
83  { return (thread_in_system()==0); }
84 
85  static GENERIC int thread_per_system(){
86  return nbod * 3;
87  }
88 
89  static GENERIC int shmem_per_system() {
90  return 0;
91  }
92 
93 
94 
99  GPUAPI void init() {
100  sqrtGM = sqrt(sys[0].mass());
101  convert_std_to_helio_pos_bary_vel_coord();
102  __syncthreads();
103  acc_bc = calcForces.acc_planets(ij,b,c);
104  }
105 
107  GPUAPI void shutdown() {
108  convert_helio_pos_bary_vel_to_std_coord ();
109  }
110 
112  GPUAPI void convert_std_to_helio_pos_bary_vel_coord()
113  { convert_std_to_helio_pos_bary_vel_coord_without_shared(); }
114 
117  GPUAPI void convert_std_to_helio_pos_bary_vel_coord_with_shared() {
118  double sump = 0., sumv = 0., mtot = 0.;
119  if( is_in_body_component_grid() )
120  {
121 
122  if( b==0 )
123  {
124  calcForces.shared[1][c].acc() = sys[0][c].pos();
126  for(int j=0;j<nbod;++j) {
127  const double mj = sys[j].mass();
128  mtot += mj;
129  sump += mj*sys[j][c].pos();
130  sumv += mj*sys[j][c].vel();
131  }
132  sumv /= mtot;
133  // There should be a better way to access shared memory
134  // What if we need to change coordinates
135  calcForces.shared[0][c].acc() = sumv;
136  }
137 
138  }
139  __syncthreads();
140 
141  if( is_in_body_component_grid() )
142  {
143  if(b==0) // For sun
144  {
145  sys[b][c].vel() = sumv;
146  sys[b][c].pos() = sump/mtot;
147  }
148  else // For planets
149  {
150  sys[b][c].vel() -= calcForces.shared[0][c].acc(); // really sumv from shared;
151  sys[b][c].pos() -= calcForces.shared[1][c].acc(); // really pc0 = original sys[0][c].pos() from shared
152  }
153  }
154 
155  }
156 
157 
159  GPUAPI void convert_std_to_helio_pos_bary_vel_coord_without_shared() {
160  double pc0;
161  double sump = 0., sumv = 0., mtot = 0.;
162  if( is_in_body_component_grid() )
163  {
164  pc0 = sys[0][c].pos();
165  // Find Center of mass and momentum
166  for(int j=0;j<nbod;++j) {
167  const double mj = sys[j].mass();
168  mtot += mj;
169  sump += mj*sys[j][c].pos();
170  sumv += mj*sys[j][c].vel();
171  }
172  sumv /= mtot;
173  }
174  __syncthreads();
175 
176  if( is_in_body_component_grid() )
177  {
178  if(b==0) // For sun
179  {
180  sys[b][c].vel() = sumv;
181  sys[b][c].pos() = sump/mtot;
182  }
183  else // For planets
184  {
185  sys[b][c].vel() -= sumv;
186  sys[b][c].pos() -= pc0;
187  }
188  }
189 
190  }
191 
192 
193 
195  GPUAPI void convert_helio_pos_bary_vel_to_std_coord ()
196  {
197  double sump = 0., sumv = 0., mtot;
198  double m0, pc0, vc0;
199  if ( is_in_body_component_grid() )
200  {
201  m0 = sys[0].mass();
202  pc0 = sys[0][c].pos();
203  vc0 = sys[0][c].vel();
204  mtot = m0;
205  for(int j=1;j<nbod;++j)
206  {
207  const double mj = sys[j].mass();
208  mtot += mj;
209  sump += mj*sys[j][c].pos();
210  sumv += mj*sys[j][c].vel();
211  }
212  sump /= mtot;
213  }
214  __syncthreads();
215 
216  if ( is_in_body_component_grid() )
217  {
218  if(b==0) // For sun only
219  {
220  sys[b][c].pos() -= sump;
221  sys[b][c].vel() -= sumv/m0;
222  }
223  else
224  {
225  sys[b][c].pos() += pc0 - sump;
226  sys[b][c].vel() += vc0;
227  }
228  }
229 
230  }
231 
233  GPUAPI void convert_internal_to_std_coord()
234  { convert_helio_pos_bary_vel_to_std_coord (); }
235 
237  GPUAPI void convert_std_to_internal_coord()
238  { convert_std_to_helio_pos_bary_vel_coord_without_shared(); }
239 
240 
242  GPUAPI void drift_step(const double hby2)
243  {
244  if(b==0)
245  {
246  sys[b][c].pos() += hby2*sys[b][c].vel();
247  }
248  else
249  {
250  double mv = 0;
251  for(int j=1;j<nbod;++j)
252  mv += sys[j].mass() * sys[j][c].vel();
253 
254  sys[b][c].pos() += mv*hby2/sys[0].mass();
255  }
256  }
257 
258 
260  GPUAPI void advance()
261  {
262  double hby2 = 0.5 * min( max_timestep , _params.time_step );
263 
264  // Step 1
265  if ( is_in_body_component_grid() )
266  drift_step(hby2);
267 
268  // Step 2: Kick Step
269  if( is_in_body_component_grid_no_star() )
270  sys[b][c].vel() += hby2 * acc_bc;
271 
272  __syncthreads();
273 
274  // 3: Kepler Drift Step (Keplerian orbit about sun/central body)
275  if( (ij>0) && (ij<nbod) )
276  drift_kepler( sys[ij][0].pos(),sys[ij][1].pos(),sys[ij][2].pos(),sys[ij][0].vel(),sys[ij][1].vel(),sys[ij][2].vel(),sqrtGM, 2.0*hby2 );
277  __syncthreads();
278 
279  // TODO: check for close encounters here
280  acc_bc = calcForces.acc_planets(ij,b,c);
281 
282  // Step 4: Kick Step
283  if( is_in_body_component_grid_no_star() )
284  sys[b][c].vel() += hby2 * acc_bc;
285  __syncthreads();
286 
287  // Step 5
288  if ( is_in_body_component_grid() )
289  drift_step(hby2);
290 
291  if( is_first_thread_in_system() )
292  sys.time() += 2.0*hby2;
293  }
294 };
295 
296 
297 
298 
299 }
300 }
301 }
302