main.cpp 8.75 KB
 Pietro Incardona committed Jul 30, 2016 1 2 3 4 #include "Grid/grid_dist_id.hpp" #include "data_type/aggregate.hpp" #include "timer.hpp"  Pietro Incardona committed Aug 02, 2016 5 6 /*! * \page Grid_3_gs Grid 3 Gray Scott  Pietro Incardona committed Jul 30, 2016 7  *  Pietro Incardona committed Aug 02, 2016 8  * # Solving a gray scott-system # {#e3_gs_gray_scott}  Pietro Incardona committed Jul 30, 2016 9 10 11 12  * * This example show the usage of periodic grid with ghost part given in grid units to solve * the following system of equations *  incardon committed Jun 03, 2017 13  * \f$\frac{\partial u}{\partial t} = D_u \nabla u - uv^2 + F(1-u)\f$  Pietro Incardona committed Aug 02, 2016 14 15  * *  incardon committed Jun 03, 2017 16  * \f$\frac{\partial v}{\partial t} = D_v \nabla v + uv^2 - (F + k)v\f$  Pietro Incardona committed Jul 30, 2016 17  *  Pietro Incardona committed Aug 02, 2016 18 19 20 21 22  * ## Constants and functions ## * * First we define convenient constants * * \snippet Grid/3_gray_scott/main.cpp constants  Pietro Incardona committed Jul 30, 2016 23 24 25  * */  Pietro Incardona committed Aug 02, 2016 26 27 //! \cond [constants] \endcond  Pietro Incardona committed Jul 30, 2016 28 29 30 constexpr int U = 0; constexpr int V = 1;  Pietro Incardona committed Aug 02, 2016 31 32 33 34 constexpr int x = 0; constexpr int y = 1; //! \cond [constants] \endcond  Pietro Incardona committed Jul 30, 2016 35   Pietro Incardona committed Aug 02, 2016 36 37 /*! * \page Grid_3_gs Grid 3 Gray Scott  Pietro Incardona committed Jul 30, 2016 38  *  Pietro Incardona committed Aug 02, 2016 39 40 41 42 43  * We also define an init function. This function initialize the species U and V. In the following we are going into the * detail of this function * * \snippet Grid/3_gray_scott/main.cpp init fun * \snippet Grid/3_gray_scott/main.cpp end fun  Pietro Incardona committed Jul 30, 2016 44 45  * */  Pietro Incardona committed Aug 02, 2016 46 47 48  //! \cond [init fun] \endcond  Pietro Incardona committed Jul 30, 2016 49 50 void init(grid_dist_id<2,double,aggregate > & Old, grid_dist_id<2,double,aggregate > & New, Box<2,double> & domain) {  Pietro Incardona committed Aug 02, 2016 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68  //! \cond [init fun] \endcond /*! * \page Grid_3_gs Grid 3 Gray Scott * * Here we initialize for the full domain. U and V itarating across the grid points. For the calculation * We are using 2 grids one Old and New. We initialize Old with the initial condition concentration of the * species U = 1 over all the domain and concentration of the specie V = 0 over all the domain. While the * New grid is initialized to 0 * * \snippet Grid/3_gray_scott/main.cpp init uv * */ //! \cond [init uv] \endcond auto it = Old.getDomainIterator();  Pietro Incardona committed Jul 30, 2016 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85  while (it.isNext()) { // Get the local grid key auto key = it.get(); // Old values U and V Old.template get(key) = 1.0; Old.template get(key) = 0.0; // Old values U and V New.template get(key) = 0.0; New.template get(key) = 0.0; ++it; }  Pietro Incardona committed Aug 02, 2016 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101  //! \cond [init uv] \endcond /*! * \page Grid_3_gs Grid 3 Gray Scott * * After we initialized the full grid, we create a perturbation in the domain with different values. * We do in the part of space: 1.55 < x < 1.85 and 1.55 < y < 1.85. Or more precisely on the points included * in this part of space. * * * \snippet Grid/3_gray_scott/main.cpp init per * */ //! \cond [init per] \endcond  Pietro Incardona committed Jul 30, 2016 102 103 104 105 106 107 108 109 110 111 112 113 114 115  grid_key_dx<2> start({(long int)std::floor(Old.size(0)*1.55f/domain.getHigh(0)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(1))}); grid_key_dx<2> stop ({(long int)std::ceil (Old.size(0)*1.85f/domain.getHigh(0)),(long int)std::ceil (Old.size(1)*1.85f/domain.getHigh(1))}); auto it_init = Old.getSubDomainIterator(start,stop); while (it_init.isNext()) { auto key = it_init.get(); Old.template get(key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/100.0; Old.template get(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/200.0; ++it_init; }  Pietro Incardona committed Aug 02, 2016 116 117 118  //! \cond [init per] \endcond //! \cond [end fun] \endcond  Pietro Incardona committed Jul 30, 2016 119 120 121  }  Pietro Incardona committed Aug 02, 2016 122 123 //! \cond [end fun] \endcond  Pietro Incardona committed Jul 30, 2016 124 125 126  int main(int argc, char* argv[]) {  Pietro Incardona committed Aug 02, 2016 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154  /*! * \page Grid_3_gs Grid 3 Gray Scott * * ## Initialization ## * * Initialize the library * * Create * * A 2D box that define the domain * * an array of 2 unsigned integer that will define the size of the grid on each dimension * * Periodicity of the grid * * A Ghost object that will define the extension of the ghost part for each sub-domain in grid point unit * * We also define numerical and physical parameters * * * Time stepping for the integration * * Diffusion constant for the species u * * Diffusion constant for the species v * * Number of time-steps * * Physical constant K * * Physical constant F * * \snippet Grid/3_gray_scott/main.cpp init lib * */ //! \cond [init lib] \endcond  Pietro Incardona committed Jul 30, 2016 155 156  openfpm_init(&argc,&argv);  Pietro Incardona committed Aug 02, 2016 157  // domain  Pietro Incardona committed Jul 30, 2016 158 159  Box<2,double> domain({0.0,0.0},{2.5,2.5});  Pietro Incardona committed Aug 02, 2016 160 161 162  // grid size size_t sz[2] = {128,128};  Pietro Incardona committed Jul 30, 2016 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184  // Define periodicity of the grid periodicity<2> bc = {PERIODIC,PERIODIC}; // Ghost in grid unit Ghost<2,long int> g(1); // deltaT double deltaT = 1; // Diffusion constant for specie U double du = 2*1e-5; // Diffusion constant for specie V double dv = 1*1e-5; // Number of timesteps size_t timeSteps = 15000; // K and F (Physical constant in the equation) double K = 0.055; double F = 0.03;  Pietro Incardona committed Aug 02, 2016 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199  //! \cond [init lib] \endcond /*! * \page Grid_3_gs Grid 3 Gray Scott * * Here we create 2 distributed grid in 2D Old and New. In particular because we want that * the second grid is distributed across processors in the same way we pass the decomposition * of the Old grid to the New one in the constructor with **Old.getDecomposition()**. Doing this, * we force the two grid to have the same decomposition. * * \snippet Grid/3_gray_scott/main.cpp init grid * */ //! \cond [init grid] \endcond  Pietro Incardona committed Jul 30, 2016 200 201  grid_dist_id<2, double, aggregate> Old(sz,domain,g,bc);  Pietro Incardona committed Aug 02, 2016 202 203  // New grid with the decomposition of the old grid  incardon committed Jun 03, 2017 204  grid_dist_id<2, double, aggregate> New(Old.getDecomposition(),sz,g);  Pietro Incardona committed Jul 30, 2016 205 206   Pietro Incardona committed Aug 02, 2016 207  // spacing of the grid on x and y  Pietro Incardona committed Jul 30, 2016 208 209  double spacing[2] = {Old.spacing(0),Old.spacing(1)};  Pietro Incardona committed Aug 02, 2016 210  //! \cond [init grid] \endcond  Pietro Incardona committed Jul 30, 2016 211   Pietro Incardona committed Aug 02, 2016 212 213 214 215 216 217 218 219 220 221  /*! * \page Grid_3_gs Grid 3 Gray Scott * * We use the function init to initialize U and V on the grid Old * * \snippet Grid/3_gray_scott/main.cpp init uvc * */ //! \cond [init uvc] \endcond  Pietro Incardona committed Jul 30, 2016 222   Pietro Incardona committed Aug 02, 2016 223  init(Old,New,domain);  Pietro Incardona committed Jul 30, 2016 224   Pietro Incardona committed Aug 02, 2016 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261  //! \cond [init uvc] \endcond /*! * \page Grid_3_gs Grid 3 Gray Scott * * ## Time stepping ## * * After initialization, we first synchronize the ghost part of the species U and V * for the grid, that we are going to read (Old). In the next we are going to do * 15000 time steps using Eulerian integration * * Because the update step of the Laplacian operator from \f$\frac{\partial u}{\partial t} = \nabla u + ... \f$ * discretized with eulerian time-stepping look like * * \f$\delta U_{next}(x,y) = \delta t D_u (\frac{U(x+1,y) - 2U(x,y) + U(x-1,y)}{(\delta x)^2} + \frac{U(x,y+1) - 2U(x,y) + U(x,y-1)}{(\delta y)^2}) + ... \f$ * * If \f$\delta x = \delta y \f$ we can simplify with * * \f$U_{next}(x,y) = \frac{\delta t D_u}{(\delta x)^2} (U(x+1,y) + U(x-1,y) + U(x,y-1) + U(x,y+1) -4U(x,y)) + ... \f$ (%Eq 2) * * The specie V follow the same concept while for the \f$... \f$ it simply expand into * * \f$- \delta t uv^2 - \delta t F(U - 1.0) \f$ * * and V the same concept * * * \see \ref e1_s_ghost * \see \ref e0_s_loop_gp * \see \ref e0_s_grid_coord * \see \ref e0_s_VTK_vis * * \snippet Grid/3_gray_scott/main.cpp time stepping * */ //! \cond [time stepping] \endcond  Pietro Incardona committed Jul 30, 2016 262 263  // sync the ghost  incardon committed Oct 26, 2016 264  size_t count = 0;  Pietro Incardona committed Aug 02, 2016 265  Old.template ghost_get();  Pietro Incardona committed Jul 30, 2016 266   Pietro Incardona committed Aug 02, 2016 267 268 269 270  // because we assume that spacing[x] == spacing[y] we use formula 2 // and we calculate the prefactor of Eq 2 double uFactor = deltaT * du/(spacing[x]*spacing[x]); double vFactor = deltaT * dv/(spacing[x]*spacing[x]);  Pietro Incardona committed Jul 30, 2016 271 272 273 274 275 276 277 278 279  for (size_t i = 0; i < timeSteps; ++i) { auto it = Old.getDomainIterator(); while (it.isNext()) { auto key = it.get();  Pietro Incardona committed Aug 02, 2016 280  // update based on Eq 2  Pietro Incardona committed Jul 30, 2016 281 282 283 284 285 286 287 288 289  New.get(key) = Old.get(key) + uFactor * ( Old.get(key.move(x,1)) + Old.get(key.move(x,-1)) + Old.get(key.move(y,1)) + Old.get(key.move(y,-1)) + -4.0*Old.get(key)) + - deltaT * Old.get(key) * Old.get(key) * Old.get(key) + - deltaT * F * (Old.get(key) - 1.0);  Pietro Incardona committed Aug 02, 2016 290  // update based on Eq 2  Pietro Incardona committed Jul 30, 2016 291 292 293 294 295 296 297 298 299  New.get(key) = Old.get(key) + vFactor * ( Old.get(key.move(x,1)) + Old.get(key.move(x,-1)) + Old.get(key.move(y,1)) + Old.get(key.move(y,-1)) - 4*Old.get(key)) + deltaT * Old.get(key) * Old.get(key) * Old.get(key) + - deltaT * (F+K) * Old.get(key);  Pietro Incardona committed Aug 02, 2016 300  // Next point in the grid  Pietro Incardona committed Jul 30, 2016 301 302 303  ++it; }  Pietro Incardona committed Aug 02, 2016 304 305 306 307 308  // Here we copy New into the old grid in preparation of the new step // It would be better to alternate, but using this we can show the usage // of the function copy. To note that copy work only on two grid of the same // decomposition. If you want to copy also the decomposition, or force to be // exactly the same, use Old = New  Pietro Incardona committed Jul 30, 2016 309 310  Old.copy(New);  Pietro Incardona committed Aug 02, 2016 311 312 313 314 315  // After copy we synchronize again the ghost part U and V Old.ghost_get(); // Every 100 time step we output the configuration for // visualization  Pietro Incardona committed Jul 30, 2016 316 317 318 319 320 321 322  if (i % 100 == 0) { Old.write("output",count); count++; } }  Pietro Incardona committed Aug 02, 2016 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337  //! \cond [time stepping] \endcond /*! * \page Grid_3_gs Grid 3 Gray Scott * * ## Finalize ## * * Deinitialize the library * * \snippet Grid/3_gray_scott/main.cpp finalize * */ //! \cond [finalize] \endcond  Pietro Incardona committed Jul 30, 2016 338  openfpm_finalize();  Pietro Incardona committed Aug 02, 2016 339 340  //! \cond [finalize] \endcond  Pietro Incardona committed Jul 30, 2016 341 }