Swarm-NG  1.1
utils.cpp
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 "common.hpp"
26 #include "utils.hpp"
27 
28 using std::max;
29 using namespace swarm;
30 using std::string;
31 
32 int number_of_disabled_systems(defaultEnsemble ens) {
33  int count_running = 0;
34  for(int i = 0; i < ens.nsys() ; i++)
35  if( ens[i].is_disabled() ) count_running++;
36  return count_running;
37 }
38 
39 
63 swarm::hostEnsemble generate_ensemble(swarm::config& cfg) {
64  int nsys = cfg.require("nsys",0);
65  int nbod = cfg.require("nbod",0);
66  double spacing_factor = cfg.optional( "spacing_factor", 1.4 );
67  double planet_mass = cfg.optional( "planet_mass" , .001 );
68  double ejection_factor = cfg.optional("ejection_factor", 1.0/sqrt(2) );
69 
70 
71  hostEnsemble ens = hostEnsemble::create( nbod, nsys );
72 
73 
74  for(unsigned int sys=0;sys<ens.nsys();++sys)
75  {
76  // set sun to unit mass and at origin
77  double mass_sun = 1.;
78  double x=0, y=0, z=0, vx=0, vy=0, vz=0;
79  ens.set_body(sys, 0, mass_sun, x, y, z, vx, vy, vz);
80 
81  // add near-Jupiter-mass planets on nearly circular orbits
82  for(unsigned int bod=1;bod<ens.nbod();++bod)
83  {
84  double rmag = pow( spacing_factor ,int(bod-1)); // semi-major axes exceeding this spacing results in systems are stable for nbody=3 and mass_planet=0.001
85  double vmag = sqrt(2*mass_sun/rmag) * ejection_factor ; // spped for uniform circular motion
86  double theta = (2.*M_PI*rand())/static_cast<double>(RAND_MAX); // randomize initial positions along ecah orbit
87  x = rmag*cos(theta); y = rmag*sin(theta); z = 0;
88  vx = -vmag*sin(theta); vy = vmag*cos(theta); vz = 0.;
89 
90  // assign body a mass, position and velocity
91  ens.set_body(sys, bod, planet_mass , x, y, z, vx, vy, vz);
92  }
93  ens[sys].set_active();
94  ens[sys].time() = 0;
95  ens[sys].id() = sys;
96  }
97  return ens;
98 }
99 
100 double find_max_energy_conservation_error(ensemble& ens, ensemble& reference_ensemble ) {
101  return energy_conservation_error_range(ens,reference_ensemble).max;
102 }
103 ensemble::range_t energy_conservation_error_range(ensemble& ens, ensemble& reference_ensemble ) {
104  std::vector<double>
105  energy_init(reference_ensemble.nsys())
106  ,energy_final(ens.nsys())
107  ,deltaE(ens.nsys());
108 
109  reference_ensemble.calc_total_energy(&energy_init[0]);
110  ens.calc_total_energy(&energy_final[0]);
111 
112  for(int sysid=0;sysid<ens.nsys();++sysid)
113  deltaE[sysid] = fabs ((energy_final[sysid]-energy_init[sysid])/energy_init[sysid] ) ;
114 
115  return ensemble::range_t::calculate( deltaE.begin(), deltaE.end() );
116 }
117 
118 
119 bool validate_configuration(config& cfg){
120  bool valid = true; // Indicates whether cfg parameters are valid
121  int nsystems = cfg.optional("nsys",1000);
122  int nbodypersystem = cfg.optional("nbod",3);
123  // WARNING: blocksize isn't being used as the block size. Trying to remove this.
124 // int bs = cfg.optional("blocksize",16);
125 
126  // Check that parameters from command line are ok
127 // if((bs<ENSEMBLE_CHUNK_SIZE)||(bs>64)) valid =false;
128 // if( bs % ENSEMBLE_CHUNK_SIZE != 0 ) valid = false;
129  if(!(nsystems>=1)||!(nsystems<=256000)) valid = false;
130  if(!(nbodypersystem>=3)||!(nbodypersystem<=MAX_NBODIES)) valid = false;
131 
132  return valid;
133 }
134 
135 void outputConfigSummary(std::ostream& o,swarm::config& cfg) {
136  o << "# Integrator:\t" << cfg["integrator"] << "\n"
137  << "# Time step\t" << cfg["time_step"] << "\n"
138  << "# Destination time\t" << cfg["destination_time"] << "\n"
139  << "# Min time step\t" << cfg["min_time_step"] << "\n"
140  << "# Max time step\t" << cfg["max_time_step"] << "\n"
141  << "# No. Systems\t" << cfg["nsys"] << "\n"
142  << "# No. Bodies\t" << cfg["nbod"] << "\n"
143 // << "# Blocksize\t" << cfg["blocksize"] << "\n"
144  << std::endl;
145 }
146 
147 // WARNING: Is it really wise to provide all of these by default in utils?
148 // This makes it easy to use a default value by accident
149 config default_config() {
150  config cfg;
151  cfg["integrator"] = "hermite_cpu"; // Set to use a GPU integrator
152  cfg["time_step"] = "0.001"; // time step
153  cfg["log_writer"] = "null";
154  cfg["nsys"] = "16";
155  cfg["nbod"] = "3";
156  return cfg;
157 }
158 std::ostream& operator << (std::ostream& o, const swarm::ensemble::range_t& r){
159  if(r.min-r.max < 1e-11)
160  return o << r.median;
161  else
162  return o << r.median << "[" << (r.min-r.median) << "," << (r.max-r.median) << "] ";
163 }
164 
165 bool compare_ensembles( swarm::ensemble& e1, swarm::ensemble &e2 , double & pos_diff, double & vel_diff, double & time_diff ) {
166  if (e1.nsys() != e2.nsys() || e1.nbod() != e2.nbod() ) return false;
167 
168  pos_diff = vel_diff = time_diff = 0;
169 
170  for(int i = 0; i < e1.nsys(); i++) {
171  for(int j = 0; j < e1.nbod() ; j++){
172 
173  // Distance between body position in ensemble e1 and e2
174  double dp = sqrt(
175  square ( e1[i][j][0].pos() - e2[i][j][0].pos() )
176  + square ( e1[i][j][1].pos() - e2[i][j][1].pos() )
177  + square ( e1[i][j][2].pos() - e2[i][j][2].pos() ) ) ;
178 
179  // Average magnitude of the position in e1 and e2
180  double ap = sqrt(
181  square ( e1[i][j][0].pos() + e2[i][j][0].pos() )
182  + square ( e1[i][j][1].pos() + e2[i][j][1].pos() )
183  + square ( e1[i][j][2].pos() + e2[i][j][2].pos() ) ) / 2.0 ;
184 
185  // Difference between body velocities in ensemble e1 and e2
186  double dv = sqrt(
187  square ( e1[i][j][0].vel() - e2[i][j][0].vel() )
188  + square ( e1[i][j][1].vel() - e2[i][j][1].vel() )
189  + square ( e1[i][j][2].vel() - e2[i][j][2].vel() ) ) ;
190 
191  // Average magnitude of the velocity in e1 and e2
192  double av = sqrt(
193  square ( e1[i][j][0].vel() + e2[i][j][0].vel() )
194  + square ( e1[i][j][1].vel() + e2[i][j][1].vel() )
195  + square ( e1[i][j][2].vel() + e2[i][j][2].vel() ) ) / 2.0 ;
196 
197  if ( dp > pos_diff ) pos_diff = dp;
198  if ( dv > vel_diff ) vel_diff = dv;
199 
200  }
201 
202  // Difference between time divided by average of the two
203  double dt = fabs(e1[i].time() - e2[i].time())
204  /(e1[i].time() + e2[i].time())*2;
205  if ( dt > time_diff ) time_diff = dt;
206 
207  }
208  return true;
209 }
210