Swarm-NG  1.1
hermite_cpu.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 #ifdef _OPENMP
26 #include <omp.h>
27 #endif
28 
29 
30 #include "swarm/common.hpp"
31 #include "swarm/integrator.hpp"
32 #include "swarm/plugin.hpp"
33 
34 namespace swarm { namespace cpu {
45 template< class Monitor >
46 class hermite_cpu : public integrator {
47  typedef integrator base;
48  typedef Monitor monitor_t;
49  typedef typename monitor_t::params mon_params_t;
50  private:
51  double _time_step;
52  mon_params_t _mon_params;
53 
54 public:
55  hermite_cpu(const config& cfg): base(cfg),_time_step(0.001), _mon_params(cfg) {
56  _time_step = cfg.require("time_step", 0.0);
57  }
58 
59  virtual void launch_integrator() {
60  #ifdef _OPENMP
61  #pragma omp parallel for
62  #endif
63  for(int i = 0; i < _ens.nsys(); i++){
64  integrate_system(_ens[i]);
65  }
66  }
67 
69  inline static double inner_product(const double a[3],const double b[3]){
70  return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
71  }
72 
74  void calcForces(ensemble::SystemRef& sys, double acc[][3],double jerk[][3]){
75  const int nbod = sys.nbod();
76 
78  for(int b = 0; b < nbod; b++) for(int c =0; c < 3; c++)
79  acc[b][c] = 0, jerk[b][c] = 0;
80 
82  for(int i=0; i < nbod-1; i++) for(int j = i+1; j < nbod; j++) {
83 
84  double dx[3] = { sys[j][0].pos()-sys[i][0].pos(),
85  sys[j][1].pos()-sys[i][1].pos(),
86  sys[j][2].pos()-sys[i][2].pos()
87  };
88  double dv[3] = { sys[j][0].vel()-sys[i][0].vel(),
89  sys[j][1].vel()-sys[i][1].vel(),
90  sys[j][2].vel()-sys[i][2].vel()
91  };
92 
94  double r2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2] * dx[2];
95  double rinv = 1 / ( sqrt(r2) * r2 ) ;
96  double rv = inner_product(dx,dv) * 3. / r2;
97 
99  const double scalar_i = +rinv*sys[j].mass();
100  for(int c = 0; c < 3; c++) {
101  acc[i][c] += dx[c]* scalar_i;
102  jerk[i][c] += (dv[c] - dx[c] * rv) * scalar_i;
103  }
104 
106  const double scalar_j = -rinv*sys[i].mass();
107  for(int c = 0; c < 3; c++) {
108  acc[j][c] += dx[c]* scalar_j;
109  jerk[j][c] += (dv[c] - dx[c] * rv) * scalar_j;
110  }
111  }
112  }
113 
115  void integrate_system(ensemble::SystemRef sys){
116  const int nbod = sys.nbod();
117  double pre_pos[nbod][3];
118  double pre_vel[nbod][3];
119  double acc0[nbod][3];
120  double acc1[nbod][3];
121  double jerk0[nbod][3];
122  double jerk1[nbod][3];
123 
124  calcForces(sys,acc0,jerk0);
125 
126  monitor_t montest (_mon_params,sys,*_log);
127 
128 
129  for(int iter = 0 ; (iter < _max_iterations) && sys.is_active() ; iter ++ ) {
130  double h = _time_step;
131 
132  if( sys.time() + h > _destination_time ) {
133  h = _destination_time - sys.time();
134  }
135 
137  for(int b = 0; b < nbod; b++) for(int c =0; c < 3; c++) {
138  sys[b][c].pos() += h * (sys[b][c].vel()+h*0.5*(acc0[b][c]+h/3*jerk0[b][c]));
139  sys[b][c].vel() += h * (acc0[b][c]+h*0.5*jerk0[b][c]);
140  }
141 
143  for(int b = 0; b < nbod; b++) for(int c =0; c < 3; c++)
144  pre_pos[b][c] = sys[b][c].pos(), pre_vel[b][c] = sys[b][c].vel();
145 
147  {
148  calcForces(sys,acc1,jerk1);
149 
150  // Correct
151  for(int b = 0; b < nbod; b++) for(int c =0; c < 3; c++) {
152  sys[b][c].pos() = pre_pos[b][c]
153  + (.1-.25) * (acc0[b][c] - acc1[b][c]) * h * h
154  - 1/60.0 * ( 7 * jerk0[b][c] + 2 * jerk1[b][c] ) * h * h * h;
155 
156  sys[b][c].vel() = pre_vel[b][c]
157  + ( -.5 ) * (acc0[b][c] - acc1[b][c] ) * h
158  - 1/12.0 * ( 5 * jerk0[b][c] + jerk1[b][c] ) * h * h;
159  }
160  }
161 
163  {
164  calcForces(sys,acc1,jerk1);
165 
166  // Correct
167  for(int b = 0; b < nbod; b++) for(int c =0; c < 3; c++) {
168  sys[b][c].pos() = pre_pos[b][c]
169  + (.1-.25) * (acc0[b][c] - acc1[b][c]) * h * h
170  - 1/60.0 * ( 7 * jerk0[b][c] + 2 * jerk1[b][c] ) * h * h * h;
171 
172  sys[b][c].vel() = pre_vel[b][c]
173  + ( -.5 ) * (acc0[b][c] - acc1[b][c] ) * h
174  - 1/12.0 * ( 5 * jerk0[b][c] + jerk1[b][c] ) * h * h;
175  }
176  }
177 
178  for(int b = 0; b < nbod; b++) for(int c =0; c < 3; c++)
179  acc0[b][c] = acc1[b][c], jerk0[b][c] = jerk1[b][c];
180 
181  sys.time() += h;
182 
183  if( sys.is_active() ) {
184  montest(0);
185  if( sys.time() >= _destination_time )
186  sys.set_inactive();
187  }
188 
189  }
190  }
191 };
192 
193 
194 
195 } } // Close namespaces