Swarm-NG  1.1
midpoint.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/swarmplugin.h"
26 
27 namespace swarm { namespace gpu { namespace bppt {
28 
33 struct MidpointPropagatorParams {
34  double time_step;
36  MidpointPropagatorParams(const config& cfg){
37  time_step = cfg.require("time_step", 0.0);
38  }
39 };
40 
52 template<class T,class Gravitation>
53 struct MidpointPropagator {
54  typedef MidpointPropagatorParams params;
55  const static int nbod = T::n;
56 
57  params _params;
58 
59 
61  ensemble::SystemRef& sys;
62  Gravitation& calcForces;
63  int b;
64  int c;
65  int ij;
66 // bool body_component_grid;
67 // bool first_thread_in_system;
68  double max_timestep;
69 
71  GPUAPI MidpointPropagator(const params& p,ensemble::SystemRef& s,
72  Gravitation& calc)
73  :_params(p),sys(s),calcForces(calc){}
74 
76  static GENERIC int thread_per_system(){
77  return nbod * 3;
78  }
80  static GENERIC int shmem_per_system() {
81  return 0;
82  }
83 
84  GPUAPI void init() { }
85 
86  GPUAPI void shutdown() { }
87 
89  GPUAPI void convert_internal_to_std_coord() {}
91  GPUAPI void convert_std_to_internal_coord() {}
92 
93  __device__ bool is_in_body_component_grid()
94 // { return body_component_grid; }
95  { return ((b < T::n) && (c < 3)); }
96 
97  __device__ bool is_in_body_component_grid_no_star()
98 // { return ( body_component_grid && (b!=0) ); }
99  { return ( (b!=0) && (b < T::n) && (c < 3) ); }
100 
101  __device__ bool is_first_thread_in_system()
102 // { return first_thread_in_system; }
103  { return (thread_in_system()==0); }
104 
106  GPUAPI void advance(){
107  double H = min( max_timestep , _params.time_step );
108  double pos = 0, vel = 0;
109 
110  if( is_in_body_component_grid() )
111  pos = sys[b][c].pos() , vel = sys[b][c].vel();
112 
113 
115 
117  const int n = 4;
118  double h = H / n;
119 
120  double p_i , p_im1, p_im2;
121  double v_i, v_im1, v_im2;
122  double a_im1;
123 
124  // Step 0
125  p_i = pos;
126  v_i = vel;
127 
128  // Step 1
129  p_im1 = p_i;
130  v_im1 = v_i;
131 
132  a_im1 = calcForces.acc(ij,b,c,p_im1,v_im1);
133 
134  p_i = p_im1 + h * v_im1;
135  v_i = v_im1 + h * a_im1;
136 
137  // Step 2 .. n
138  for(int i = 2; i <= n; i++){
139  p_im2 = p_im1;
140  p_im1 = p_i;
141  v_im2 = v_im1;
142  v_im1 = v_i;
143 
144  a_im1 = calcForces.acc(ij,b,c,p_im1,v_im1);
145 
146  p_i = p_im2 + 2.0 * h * v_im1;
147  v_i = v_im2 + 2.0 * h * a_im1;
148  }
149  double a_i = calcForces.acc(ij,b,c,p_i,v_i);
150 
151  pos = 0.5 * ( p_i + p_im1 + h * v_i );
152  vel = 0.5 * ( v_i + v_im1 + h * a_i );
153 
154 
155  // Finalize the step
156  if( is_in_body_component_grid() )
157  sys[b][c].pos() = pos , sys[b][c].vel() = vel;
158  if( is_first_thread_in_system() )
159  sys.time() += H;
160  }
161 };
162 
163 
164 } } } // End namespace bppt :: gpu :: swarm