Swarm-NG  1.1
hermite_integrator.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/gpu/bppt.hpp"
27 
28 namespace swarm { namespace gpu { namespace bppt {
29 
34 template< class Monitor , template<class T> class Gravitation >
35 class hermite: public integrator {
36  typedef integrator base;
37  typedef Monitor monitor_t;
38  typedef typename monitor_t::params mon_params_t;
39  private:
40  double _time_step;
41  mon_params_t _mon_params;
42 
43 public:
44  hermite(const config& cfg): base(cfg),_time_step(0.001), _mon_params(cfg) {
45  _time_step = cfg.require("time_step", 0.0);
46  }
47 
48  virtual void launch_integrator() {
49  launch_templatized_integrator(this);
50  }
51 
53  GPUAPI void convert_internal_to_std_coord() {}
55  GPUAPI void convert_std_to_internal_coord() {}
56 
57  template<class T>
58  __device__ void kernel(T compile_time_param){
59 
60  if(sysid()>=_dens.nsys()) return;
61  // References to Ensemble and Shared Memory
62  typedef Gravitation<T> Grav;
63  ensemble::SystemRef sys = _dens[sysid()];
64  typedef typename Grav::shared_data grav_t;
65  Grav calcForces(sys,*( (grav_t*) system_shared_data_pointer(this,compile_time_param) ) );
66 
67  // Local variables
68  const int nbod = T::n;
69  // Body number
70  const int b = thread_body_idx(nbod);
71  // Component number
72  const int c = thread_component_idx(nbod);
73 
74  // local variables
75  monitor_t montest(_mon_params,sys,*_log) ;
76 
77 
78  // local information per component per body
79  double pos = 0.0, vel = 0.0 , acc0 = 0.0, jerk0 = 0.0;
80  if( (b < nbod) && (c < 3) )
81  { pos = sys[b][c].pos(); vel = sys[b][c].vel(); }
82 
83 
84 // if( thread_in_system()==0 ) {
85  montest( thread_in_system() );
86 // }
87 
89 
91  calcForces(thread_in_system(),b,c,pos,vel,acc0,jerk0);
92 
93  for(int iter = 0 ; (iter < _max_iterations) && sys.is_active() ; iter ++ )
94  {
95  double h = _time_step;
96 
97  if( sys.time() + h > _destination_time ) {
98  h = _destination_time - sys.time();
99  }
100 
101 
104  //calcForces(thread_in_system(),b,c,pos,vel,acc0,jerk0);
105 
106  // Predict
107  pos = pos + h*(vel+(h*0.5)*(acc0+(h/3.0)*jerk0));
108  vel = vel + h*(acc0+(h*0.5)*jerk0);
109 
110  double pre_pos = pos, pre_vel = vel;
111 
112  double acc1,jerk1;
113  {
114  // Evaluation
115  calcForces(thread_in_system(),b,c,pos,vel,acc1,jerk1);
116 
117  // Correct
118 #if 0 // OLD
119  pos = pre_pos + (0.1-0.25) * (acc0 - acc1) * h * h - 1.0/60.0 * ( 7.0 * jerk0 + 2.0 * jerk1 ) * h * h * h;
120  vel = pre_vel + ( -0.5 ) * (acc0 - acc1 ) * h - 1.0/12.0 * ( 5.0 * jerk0 + jerk1 ) * h * h;
121 #else
122  pos = pre_pos + ( (0.1-0.25) * (acc0 - acc1) - 1.0/60.0 * ( 7.0 * jerk0 + 2.0 * jerk1 ) * h) * h * h;
123  vel = pre_vel + (( -0.5 ) * (acc0 - acc1 ) - 1.0/12.0 * ( 5.0 * jerk0 + jerk1 ) * h )* h ;
124 #endif
125  }
126  {
127  // Evaluation
128  calcForces(thread_in_system(),b,c,pos,vel,acc1,jerk1);
129 
130  // Correct
131 #if 0 // OLD
132  pos = pre_pos + (0.1-0.25) * (acc0 - acc1) * h * h - 1.0/60.0 * ( 7.0 * jerk0 + 2.0 * jerk1 ) * h * h * h;
133  vel = pre_vel + ( -0.5 ) * (acc0 - acc1 ) * h - 1.0/12.0 * ( 5.0 * jerk0 + jerk1 ) * h * h;
134 #else
135  pos = pre_pos + ((0.1-0.25) * (acc0 - acc1) - 1.0/60.0 * ( 7.0 * jerk0 + 2.0 * jerk1 ) * h )* h * h ;
136  vel = pre_vel + (( -0.5 ) * (acc0 - acc1 ) - 1.0/12.0 * ( 5.0 * jerk0 + jerk1 ) * h ) * h ;
137 #endif
138  }
139  acc0 = acc1, jerk0 = jerk1;
140 
142  if( (b < nbod) && (c < 3) )
143  { sys[b][c].pos() = pos; sys[b][c].vel() = vel; }
144  if( thread_in_system()==0 )
145  sys.time() += h;
146  __syncthreads();
147  montest( thread_in_system() );
148  __syncthreads();
149  if( sys.is_active() && thread_in_system()==0 ) {
150  if( sys.time() >= _destination_time )
151  { sys.set_inactive(); }
152  }
153 
154  __syncthreads();
155 
156 
157  }
158 
159  }
160 
161 
162 };
163 
164 } } } // end namespace bppt :: integrators :: swarm