OpenFPM  5.2.0
Project that contain the implementation of distributed structures
main.cpp
Go to the documentation of this file.
1 //
2 // Created by Abhinav Singh on 15.03.20.
3 //
4 
54 // Include Vector Expression,Vector Expressions for Subset,DCPSE,Odeint header files
55 #include "config.h"
56 #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
57 #define BOOST_MPL_LIMIT_VECTOR_SIZE 40
58 #include <iostream>
59 #include "DCPSE/DCPSE_op/DCPSE_op.hpp"
60 #include "DCPSE/DCPSE_op/DCPSE_Solver.hpp"
61 #include "Operators/Vector/vector_dist_operators.hpp"
62 #include "Vector/vector_dist_subset.hpp"
63 #include "DCPSE/DCPSE_op/EqnsStruct.hpp"
64 #include "OdeIntegrators/OdeIntegrators.hpp"
66 
88 constexpr int x = 0;
89 constexpr int y = 1;
90 constexpr int POLARIZATION= 0,VELOCITY = 1, VORTICITY = 2, EXTFORCE = 3,PRESSURE = 4, STRAIN_RATE = 5, STRESS = 6, MOLFIELD = 7, DPOL = 8, DV = 9, VRHS = 10, F1 = 11, F2 = 12, F3 = 13, F4 = 14, F5 = 15, F6 = 16, V_T = 17, DIV = 18, DELMU = 19, HPB = 20, FE = 21, R = 22;
91 
92 double eta = 1.0;
93 double nu = -0.5;
94 double gama = 0.1;
95 double zeta = 0.07;
96 double Ks = 1.0;
97 double Kb = 1.0;
98 double lambda = 0.1;
99 
100 int wr_f;
101 int wr_at;
102 double V_err_eps;
103 
104 void *vectorGlobal=nullptr,*vectorGlobal_bulk=nullptr,*vectorGlobal_boundary=nullptr;
105 const openfpm::vector<std::string>
106 PropNAMES={"00-Polarization","01-Velocity","02-Vorticity","03-ExternalForce","04-Pressure","05-StrainRate","06-Stress","07-MolecularField","08-DPOL","09-DV","10-VRHS","11-f1","12-f2","13-f3","14-f4","15-f5","16-f6","17-V_T","18-DIV","19-DELMU","20-HPB","21-FrankEnergy","22-R"};
107 typedef aggregate<VectorS<2, double>,VectorS<2, double>,double[2][2],VectorS<2, double>,double,double[2][2],double[2][2],VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,double,double,double,double,double,double,VectorS<2, double>,double,double,double,double,double> Activegels;
108 typedef vector_dist_ws<2, double,Activegels> vector_type;
109 typedef vector_dist_subset<2, double, Activegels> vector_type2;
111 
138 //Functor to Compute RHS of the time derivative of the polarity
139 template<typename DX,typename DY,typename DXX,typename DXY,typename DYY>
140 struct PolarEv
141 {
142  DX &Dx;
143  DY &Dy;
144  DXX &Dxx;
145  DXY &Dxy;
146  DYY &Dyy;
147  //Constructor
148  PolarEv(DX &Dx,DY &Dy,DXX &Dxx,DXY &Dxy,DYY &Dyy):Dx(Dx),Dy(Dy),Dxx(Dxx),Dxy(Dxy),Dyy(Dyy)
149  {}
150 
151  void operator()( const state_type_2d_ofp &X , state_type_2d_ofp &dxdt , const double t ) const
152  {
153 
154  vector_type &Particles= *(vector_type *) vectorGlobal;
155  vector_type2 &Particles_bulk= *(vector_type2 *) vectorGlobal_bulk;
156 
157  auto Pol=getV<POLARIZATION>(Particles);
158  auto Pol_bulk=getV<POLARIZATION>(Particles_bulk);
159  auto h = getV<MOLFIELD>(Particles);
160  auto u = getV<STRAIN_RATE>(Particles);
161  auto dPol = getV<DPOL>(Particles);
162  auto W = getV<VORTICITY>(Particles);
163  auto delmu = getV<DELMU>(Particles);
164  auto H_p_b = getV<HPB>(Particles);
165  auto r = getV<R>(Particles);
166  auto dPol_bulk = getV<DPOL>(Particles_bulk);
167  Pol[x]=X.data.get<0>();
168  Pol[y]=X.data.get<1>();
169  Particles.ghost_get<POLARIZATION>(SKIP_LABELLING);
170  H_p_b = Pol[x] * Pol[x] + Pol[y] * Pol[y];
171 
172  h[y] = (Pol[x] * (Ks * Dyy(Pol[y]) + Kb * Dxx(Pol[y]) + (Ks - Kb) * Dxy(Pol[x])) -
173  Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y])));
174 
175  h[x] = -gama * (lambda * delmu - nu * (u[x][x] * Pol[x] * Pol[x] + u[y][y] * Pol[y] * Pol[y] + 2 * u[x][y] * Pol[x] * Pol[y]) / (H_p_b));
176 
177  dPol_bulk[x] = ((h[x] * Pol[x] - h[y] * Pol[y]) / gama + lambda * delmu * Pol[x] -
178  nu * (u[x][x] * Pol[x] + u[x][y] * Pol[y]) + W[x][x] * Pol[x] +
179  W[x][y] * Pol[y]);
180  dPol_bulk[y] = ((h[x] * Pol[y] + h[y] * Pol[x]) / gama + lambda * delmu * Pol[y] -
181  nu * (u[y][x] * Pol[x] + u[y][y] * Pol[y]) + W[y][x] * Pol[x] +
182  W[y][y] * Pol[y]);
183  dxdt.data.get<0>()=dPol[x];
184  dxdt.data.get<1>()=dPol[y];
185  }
186 };
188 
237 // Functor to calculate velocity and move particles with explicit euler
238 template<
239  typename DX,
240  typename DY,
241  typename DX_BULK,
242  typename DY_BULK,
243  typename DXX,
244  typename DXY,
245  typename DYY,
246  typename verletList_type,
247  typename verletListBulk_type>
248 struct CalcVelocity
249 {
250 
251  DX &Dx, &Bulk_Dx;
252  DY &Dy, &Bulk_Dy;
253  DXX &Dxx;
254  DXY &Dxy;
255  DYY &Dyy;
256  verletList_type& verletList;
257  verletListBulk_type& verletList_bulk;
258 
259  double t_old;
260  double rCut;
261  int ctr;
262 
263  //Constructor
264  CalcVelocity(
265  DX &Dx,
266  DY &Dy,
267  DXX &Dxx,
268  DXY &Dxy,
269  DYY &Dyy,
270  DX_BULK &Bulk_Dx,
271  DY_BULK &Bulk_Dy,
272  verletList_type& verletList,
273  verletListBulk_type& verletList_bulk,
274  double rCut
275  ):
276  Dx(Dx),
277  Dy(Dy),
278  Dxx(Dxx),
279  Dxy(Dxy),
280  Dyy(Dyy),
281  Bulk_Dx(Bulk_Dx),
282  Bulk_Dy(Bulk_Dy),
283  verletList(verletList),
284  verletList_bulk(verletList_bulk),
285  rCut(rCut)
286  {
287  t_old = 0.0;
288  ctr = 0;
289  }
290 
291  void operator() (state_type_2d_ofp &state, double t)
292  {
293 
294  double dt = t - t_old;
295  vector_type &Particles= *(vector_type *) vectorGlobal;
296  vector_type2 &Particles_bulk= *(vector_type2 *) vectorGlobal_bulk;
297  vector_type2 &Particles_boundary= *(vector_type2 *) vectorGlobal_boundary;
298  auto &v_cl = create_vcluster();
299 
300  timer tt;
301 
302  auto Pos = getV<POS_PROP>(Particles);
303  auto Pol=getV<POLARIZATION>(Particles);
304  auto V = getV<VELOCITY>(Particles);
305  auto H_p_b = getV<HPB>(Particles);
306 
307 
308  // skip in first time step
309  if (dt != 0) {
310  tt.start();
311  //normalize polarization
312  H_p_b = sqrt(Pol[x] * Pol[x] + Pol[y] * Pol[y]);
313  Pol = Pol / H_p_b;
314  Pos = Pos + (t-t_old)*V;
315  Particles.map();
316  Particles.ghost_get<POLARIZATION, EXTFORCE, DELMU>();
317  Particles_bulk.update();
318  Particles_boundary.update();
319  tt.start();
320  Particles.updateVerlet(verletList,rCut);
321  Particles_bulk.updateVerlet(verletList_bulk,rCut);
322  Dx.update(Particles);
323  Dy.update(Particles);
324  Dxy.update(Particles);
325  Dxx.update(Particles);
326  Dyy.update(Particles);
327 
328  Bulk_Dx.update(Particles_bulk);
329  Bulk_Dy.update(Particles_bulk);
330 
331  state.data.get<0>()=Pol[x];
332  state.data.get<1>()=Pol[y];
333 
334  tt.stop();
335  if (v_cl.rank() == 0) {
336  std::cout << "Updating operators took " << tt.getwct() << " seconds." << std::endl;
337  std::cout << "Time step " << ctr << " : " << t << " over." << std::endl;
338  std::cout << "----------------------------------------------------------" << std::endl;
339  }
340  ctr++;
341 
342  }
343  auto Dyx = Dxy;
344  t_old = t;
345  tt.start();
346 
347  auto & bulk = Particles_bulk.getIds();
348  auto & boundary = Particles_boundary.getIds();
349  auto Pol_bulk=getV<POLARIZATION>(Particles_bulk);
350  auto sigma = getV<STRESS>(Particles);
351  auto r = getV<R>(Particles);
352  auto h = getV<MOLFIELD>(Particles);
353  auto FranckEnergyDensity = getV<FE>(Particles);
354  auto f1 = getV<F1>(Particles);
355  auto f2 = getV<F2>(Particles);
356  auto f3 = getV<F3>(Particles);
357  auto f4 = getV<F4>(Particles);
358  auto f5 = getV<F5>(Particles);
359  auto f6 = getV<F6>(Particles);
360  auto dV = getV<DV>(Particles);
361  auto delmu = getV<DELMU>(Particles); // why is delmu a property, not a constant?
362  auto g = getV<EXTFORCE>(Particles);
363  auto P = getV<PRESSURE>(Particles);
364  auto P_bulk = getV<PRESSURE>(Particles_bulk); //Pressure only on inside
365  auto RHS = getV<VRHS>(Particles);
366  auto RHS_bulk = getV<VRHS>(Particles_bulk);
367  auto div = getV<DIV>(Particles);
368  auto V_t = getV<V_T>(Particles);
369  auto u = getV<STRAIN_RATE>(Particles);
370  auto W = getV<VORTICITY>(Particles);
371 
372  Pol_bulk[x]=state.data.get<0>();
373  Pol_bulk[y]=state.data.get<1>();
374  Particles.ghost_get<POLARIZATION>(SKIP_LABELLING);
375 
376  eq_id x_comp, y_comp;
377  x_comp.setId(0);
378  y_comp.setId(1);
379 
380  int n = 0,nmax = 300,errctr = 0, Vreset = 0;
381  double V_err = 1,V_err_eps = 5 * 1e-3, V_err_old,sum, sum1;
382  std::cout << "Calculate velocity (step t=" << t << ")" << std::endl;
383  tt.start();
384  petsc_solver<double> solverPetsc;
385  solverPetsc.setSolver(KSPGMRES);
386  solverPetsc.setPreconditioner(PCJACOBI);
387  Particles.ghost_get<POLARIZATION>(SKIP_LABELLING);
388  // calculate stress
389  sigma[x][x] =
390  -Ks * Dx(Pol[x]) * Dx(Pol[x]) - Kb * Dx(Pol[y]) * Dx(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dx(Pol[y]);
391  sigma[x][y] =
392  -Ks * Dy(Pol[y]) * Dx(Pol[y]) - Kb * Dy(Pol[x]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dx(Pol[x]);
393  sigma[y][x] =
394  -Ks * Dx(Pol[x]) * Dy(Pol[x]) - Kb * Dx(Pol[y]) * Dy(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dy(Pol[y]);
395  sigma[y][y] =
396  -Ks * Dy(Pol[y]) * Dy(Pol[y]) - Kb * Dy(Pol[x]) * Dy(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dy(Pol[x]);
397  Particles.ghost_get<STRESS>(SKIP_LABELLING);
398 
399  // if R == 0 then set to 1 to avoid division by zero for defects
400  r = Pol[x] * Pol[x] + Pol[y] * Pol[y];
401  for (int j = 0; j < bulk.size(); j++) {
402  auto p = bulk.get<0>(j);
403  Particles.getProp<R>(p) = (Particles.getProp<R>(p) == 0) ? 1 : Particles.getProp<R>(p);
404  }
405  for (int j = 0; j < boundary.size(); j++) {
406  auto p = boundary.get<0>(j);
407  Particles.getProp<R>(p) = (Particles.getProp<R>(p) == 0) ? 1 : Particles.getProp<R>(p);
408  }
409 
410  // calculate traversal molecular field (H_perpendicular)
411  h[y] = (Pol[x] * (Ks * Dyy(Pol[y]) + Kb * Dxx(Pol[y]) + (Ks - Kb) * Dxy(Pol[x])) -
412  Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y])));
413  Particles.ghost_get<MOLFIELD>(SKIP_LABELLING);
414 
415  // calulate FranckEnergyDensity
416  FranckEnergyDensity = (Ks / 2.0) *
417  ((Dx(Pol[x]) * Dx(Pol[x])) + (Dy(Pol[x]) * Dy(Pol[x])) +
418  (Dx(Pol[y]) * Dx(Pol[y])) +
419  (Dy(Pol[y]) * Dy(Pol[y]))) +
420  ((Kb - Ks) / 2.0) * ((Dx(Pol[y]) - Dy(Pol[x])) * (Dx(Pol[y]) - Dy(Pol[x])));
421  Particles.ghost_get<FE>(SKIP_LABELLING);
422 
423  // calculate preactors for LHS of Stokes Equation.
424  f1 = gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
425  f2 = 2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
426  f3 = gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
427  f4 = 2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (r);
428  f5 = 4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (r);
429  f6 = 2.0 * gama * nu * Pol[x] * Pol[y] * Pol[y] * Pol[y] / (r);
430  Particles.ghost_get<F1, F2, F3, F4, F5, F6>(SKIP_LABELLING);
431  texp_v<double> Dxf1 = Dx(f1),Dxf2 = Dx(f2),Dxf3 = Dx(f3),Dxf4 = Dx(f4),Dxf5 = Dx(f5),Dxf6 = Dx(f6),
432  Dyf1 = Dy(f1),Dyf2 = Dy(f2),Dyf3 = Dy(f3),Dyf4 = Dy(f4),Dyf5 = Dy(f5),Dyf6 = Dy(f6);
433 
434  // calculate RHS of Stokes Equation without pressure
435  dV[x] = -0.5 * Dy(h[y]) + zeta * Dx(delmu * Pol[x] * Pol[x]) + zeta * Dy(delmu * Pol[x] * Pol[y]) -
436  zeta * Dx(0.5 * delmu * (Pol[x] * Pol[x] + Pol[y] * Pol[y])) -
437  0.5 * nu * Dx(-2.0 * h[y] * Pol[x] * Pol[y])
438  - 0.5 * nu * Dy(h[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y])) - Dx(sigma[x][x]) -
439  Dy(sigma[x][y]) -
440  g[x]
441  - 0.5 * nu * Dx(-gama * lambda * delmu * (Pol[x] * Pol[x] - Pol[y] * Pol[y]))
442  - 0.5 * Dy(-2.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
443 
444  dV[y] = -0.5 * Dx(-h[y]) + zeta * Dy(delmu * Pol[y] * Pol[y]) + zeta * Dx(delmu * Pol[x] * Pol[y]) -
445  zeta * Dy(0.5 * delmu * (Pol[x] * Pol[x] + Pol[y] * Pol[y])) -
446  0.5 * nu * Dy(2.0 * h[y] * Pol[x] * Pol[y])
447  - 0.5 * nu * Dx(h[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y])) - Dx(sigma[y][x]) -
448  Dy(sigma[y][y]) -
449  g[y]
450  - 0.5 * nu * Dy(gama * lambda * delmu * (Pol[x] * Pol[x] - Pol[y] * Pol[y]))
451  - 0.5 * Dx(-2.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
452  Particles.ghost_get<DV>(SKIP_LABELLING);
453 
454  // Encode LHS of the Stokes Equations
455  auto Stokes1 = eta * (Dxx(V[x]) + Dyy(V[x]))
456  + 0.5 * nu * (Dxf1 * Dx(V[x]) + f1 * Dxx(V[x]))
457  + 0.5 * nu * (Dxf2 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
458  + 0.5 * nu * (Dxf3 * Dy(V[y]) + f3 * Dyx(V[y]))
459  + 0.5 * nu * (Dyf4 * Dx(V[x]) + f4 * Dxy(V[x]))
460  + 0.5 * nu * (Dyf5 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
461  + 0.5 * nu * (Dyf6 * Dy(V[y]) + f6 * Dyy(V[y]));
462  auto Stokes2 = eta * (Dxx(V[y]) + Dyy(V[y]))
463  - 0.5 * nu * (Dyf1 * Dx(V[x]) + f1 * Dxy(V[x]))
464  - 0.5 * nu * (Dyf2 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
465  - 0.5 * nu * (Dyf3 * Dy(V[y]) + f3 * Dyy(V[y]))
466  + 0.5 * nu * (Dxf4 * Dx(V[x]) + f4 * Dxx(V[x]))
467  + 0.5 * nu * (Dxf5 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
468  + 0.5 * nu * (Dxf6 * Dy(V[y]) + f6 * Dyx(V[y]));
469 
470  tt.stop();
471 
472  std::cout << "Init of Velocity took " << tt.getwct() << " seconds." << std::endl;
473 
474  tt.start();
475  V_err = 1;
476  n = 0;
477  errctr = 0;
478  if (Vreset == 1) {
479  P = 0;
480  Vreset = 0;
481  }
482  P=0;
483 
484  // integrate velocity
485  Particles.ghost_get<PRESSURE>(SKIP_LABELLING);
486  RHS_bulk[x] = dV[x] + Bulk_Dx(P);
487  RHS_bulk[y] = dV[y] + Bulk_Dy(P);
488  Particles.ghost_get<VRHS>(SKIP_LABELLING);
489 
490  // prepare solver
491  DCPSE_scheme<equations2d2, vector_type> Solver(Particles);
492  Solver.impose(Stokes1, bulk, RHS[0], x_comp);
493  Solver.impose(Stokes2, bulk, RHS[1], y_comp);
494  Solver.impose(V[x], boundary, 0, x_comp);
495  Solver.impose(V[y], boundary, 0, y_comp);
496  Solver.solve_with_solver(solverPetsc, V[x], V[y]);
497  Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
498  div = -(Dx(V[x]) + Dy(V[y]));
499  P_bulk = P + div;
500 
501  // approximate velocity
502  while (V_err >= V_err_eps && n <= nmax) {
503  Particles.ghost_get<PRESSURE>(SKIP_LABELLING);
504  RHS_bulk[x] = dV[x] + Bulk_Dx(P);
505  RHS_bulk[y] = dV[y] + Bulk_Dy(P);
506  Particles.ghost_get<VRHS>(SKIP_LABELLING);
507  Solver.reset_b();
508  Solver.impose_b(bulk, RHS[0], x_comp);
509  Solver.impose_b(bulk, RHS[1], y_comp);
510  Solver.impose_b(boundary, 0, x_comp);
511  Solver.impose_b(boundary, 0, y_comp);
512  Solver.solve_with_solver(solverPetsc, V[x], V[y]);
513  Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
514  div = -(Dx(V[x]) + Dy(V[y]));
515  P_bulk = P + div;
516  // calculate error
517  sum = 0;
518  sum1 = 0;
519  for (int j = 0; j < bulk.size(); j++) {
520  auto p = bulk.get<0>(j);
521  sum += (Particles.getProp<V_T>(p)[0] - Particles.getProp<VELOCITY>(p)[0]) *
522  (Particles.getProp<V_T>(p)[0] - Particles.getProp<VELOCITY>(p)[0]) +
523  (Particles.getProp<V_T>(p)[1] - Particles.getProp<VELOCITY>(p)[1]) *
524  (Particles.getProp<V_T>(p)[1] - Particles.getProp<VELOCITY>(p)[1]);
525  sum1 += Particles.getProp<VELOCITY>(p)[0] * Particles.getProp<VELOCITY>(p)[0] +
526  Particles.getProp<VELOCITY>(p)[1] * Particles.getProp<VELOCITY>(p)[1];
527  }
528  V_t = V;
529  v_cl.sum(sum);
530  v_cl.sum(sum1);
531  v_cl.execute();
532  sum = sqrt(sum);
533  sum1 = sqrt(sum1);
534  V_err_old = V_err;
535  V_err = sum / sum1;
536  if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
537  errctr++;
538  //alpha_P -= 0.1;
539  } else {
540  errctr = 0;
541  }
542  if (n > 3) {
543  if (errctr > 3) {
544  std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN DIVERGENCE" << std::endl;
545  Vreset = 1;
546  break;
547  } else {
548  Vreset = 0;
549  }
550  }
551  n++;
552 
553  }
554  tt.stop();
555 
556  Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
557  // calculate strain rate
558  u[x][x] = Dx(V[x]);
559  u[x][y] = 0.5 * (Dx(V[y]) + Dy(V[x]));
560  u[y][x] = 0.5 * (Dy(V[x]) + Dx(V[y]));
561  u[y][y] = Dy(V[y]);
562 
563  if (v_cl.rank() == 0) {
564  std::cout << "Rel l2 cgs err in V = " << V_err << " and took " << tt.getwct() << " seconds with " << n
565  << " iterations. dt for the stepper is " << dt
566  << std::endl;
567  }
568 
569  // calculate vorticity
570  W[x][x] = 0;
571  W[x][y] = 0.5 * (Dy(V[x]) - Dx(V[y]));
572  W[y][x] = 0.5 * (Dx(V[y]) - Dy(V[x]));
573  W[y][y] = 0;
574 
575  if (ctr%wr_at==0 || ctr==wr_f){
576  Particles.deleteGhost();
577  Particles.write_frame("Polar",ctr);
578  Particles.ghost_get<POLARIZATION>();
579  }
580  }
581 };
583 
595 int main(int argc, char* argv[])
596 {
597  { openfpm_init(&argc,&argv);
598  timer tt2;
599  tt2.start();
600  size_t Gd = int(std::atof(argv[1]));
601  double tf = std::atof(argv[2]);
602  double dt = tf/std::atof(argv[3]);
603  wr_f=int(std::atof(argv[3]));
604  wr_at=1;
605  V_err_eps = 5e-4;
607 
625  double boxsize = 10;
626  const size_t sz[2] = {Gd, Gd};
627  Box<2, double> box({0, 0}, {boxsize, boxsize});
628  double Lx = box.getHigh(0),Ly = box.getHigh(1);
629  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
630  double spacing = box.getHigh(0) / (sz[0] - 1),rCut = 3.9 * spacing;
631  int ord = 2;
632  Ghost<2, double> ghost(rCut);
633  auto &v_cl = create_vcluster();
634  vector_dist_ws<2, double,Activegels> Particles(0, box, bc, ghost);
635  Particles.setPropNames(PropNAMES);
636 
637  double x0=box.getLow(0), y0=box.getLow(1), x1=box.getHigh(0), y1=box.getHigh(1);
638  auto it = Particles.getGridIterator(sz);
639  while (it.isNext()) {
640  Particles.add();
641  auto key = it.get();
642  double xp = key.get(0) * it.getSpacing(0),yp = key.get(1) * it.getSpacing(1);
643  Particles.getLastPos()[x] = xp;
644  Particles.getLastPos()[y] = yp;
645  if (xp != x0 && yp != y0 && xp != x1 && yp != y1)
646  Particles.getLastSubset(0);
647  else
648  Particles.getLastSubset(1);
649  ++it;
650  }
651  Particles.map();
652  Particles.ghost_get<POLARIZATION>();
653 
654  //auto Pos = getV<POS_PROP>(Particles);
655 
656  auto Pol = getV<POLARIZATION>(Particles);
657  auto V = getV<VELOCITY>(Particles);
658  auto g = getV<EXTFORCE>(Particles);
659  auto P = getV<PRESSURE>(Particles);
660  auto delmu = getV<DELMU>(Particles);
661  auto dPol = getV<DPOL>(Particles);
662 
663  g = 0;delmu = -1.0;P = 0;V = 0;
664  auto it2 = Particles.getDomainIterator();
665  while (it2.isNext()) {
666  auto p = it2.get();
667  Point<2, double> xp = Particles.getPos(p);
668  Particles.getProp<POLARIZATION>(p)[x] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
669  Particles.getProp<POLARIZATION>(p)[y] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
670  ++it2;
671  }
672  Particles.ghost_get<POLARIZATION,EXTFORCE,DELMU>(SKIP_LABELLING);
674 
688 
689  vector_dist_subset<2, double, Activegels> Particles_bulk(Particles,0);
690  vector_dist_subset<2, double, Activegels> Particles_boundary(Particles,1);
691  auto & bulk = Particles_bulk.getIds();
692  auto & boundary = Particles_boundary.getIds();
693 
694  auto P_bulk = getV<PRESSURE>(Particles_bulk);//Pressure only on inside
695  auto Pol_bulk = getV<POLARIZATION>(Particles_bulk);;
696  auto dPol_bulk = getV<DPOL>(Particles_bulk);
697  auto dV_bulk = getV<DV>(Particles_bulk);
698  auto RHS_bulk = getV<VRHS>(Particles_bulk);
699  auto div_bulk = getV<DIV>(Particles_bulk);
700 
701  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
702  auto verletList_bulk = Particles_bulk.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
703 
704  Derivative_x<decltype(verletList)> Dx(Particles, verletList, ord,rCut);
705  Derivative_y<decltype(verletList)> Dy(Particles, verletList, ord, rCut);
706  Derivative_xy<decltype(verletList)> Dxy(Particles, verletList, ord, rCut);
707 
708  Derivative_x<decltype(verletList_bulk)> Bulk_Dx(Particles, Particles_bulk, verletList_bulk, ord, rCut);
709  Derivative_y<decltype(verletList_bulk)> Bulk_Dy(Particles, Particles_bulk, verletList_bulk, ord, rCut);
710  auto Dyx = Dxy;
711  Derivative_xx<decltype(verletList)> Dxx(Particles, verletList, ord, rCut);
712  Derivative_yy<decltype(verletList)> Dyy(Particles, verletList, ord, rCut);
713 
714  boost::numeric::odeint::runge_kutta4< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
715 
716  vectorGlobal=(void *) &Particles;
717  vectorGlobal_bulk=(void *) &Particles_bulk;
718  vectorGlobal_boundary=(void *) &Particles_boundary;
719 
720  PolarEv<
721  Derivative_x<decltype(verletList)>,
722  Derivative_y<decltype(verletList)>,
723  Derivative_xx<decltype(verletList)>,
724  Derivative_xy<decltype(verletList)>,
725  Derivative_yy<decltype(verletList)>>
726  System(Dx,Dy,Dxx,Dxy,Dyy);
727 
728  CalcVelocity<
729  Derivative_x<decltype(verletList)>,
730  Derivative_y<decltype(verletList)>,
731  Derivative_x<decltype(verletList_bulk)>,
732  Derivative_y<decltype(verletList_bulk)>,
733  Derivative_xx<decltype(verletList)>,
734  Derivative_xy<decltype(verletList)>,
735  Derivative_yy<decltype(verletList)>,
736  decltype(verletList),
737  decltype(verletList_bulk)>
738  CalcVelocityObserver(Dx,Dy,Dxx,Dxy,Dyy,Bulk_Dx,Bulk_Dy,verletList,verletList_bulk,rCut);
739 
740  state_type_2d_ofp tPol;
741  tPol.data.get<0>()=Pol[x];
742  tPol.data.get<1>()=Pol[y];
744 
745 
746  timer tt;
747  timer tt3;
748  dPol = Pol;
749  double V_err = 1, V_err_old;
750  double tim=0;
774  // intermediate time steps
775  std::vector<double> inter_times;
776 
777  size_t steps = integrate_const(rk4 , System , tPol , tim , tf , dt, CalcVelocityObserver);
778 
779 
780  std::cout << "Time steps: " << steps << std::endl;
781 
782  Pol_bulk[x]=tPol.data.get<0>();
783  Pol_bulk[y]=tPol.data.get<1>();
784 
785  Particles.deleteGhost();
786  Particles.write("Polar_Last");
787  Dx.deallocate(Particles);
788  Dy.deallocate(Particles);
789  Dxy.deallocate(Particles);
790  Dxx.deallocate(Particles);
791  Dyy.deallocate(Particles);
792  Bulk_Dx.deallocate(Particles_bulk);
793  Bulk_Dy.deallocate(Particles_bulk);
794  std::cout.precision(17);
795  tt2.stop();
796  if (v_cl.rank() == 0) {
797  std::cout << "The simulation took " << tt2.getcputime() << "(CPU) ------ " << tt2.getwct()
798  << "(Wall) Seconds.";
799  }
800  }
801  openfpm_finalize();
803 }
804 
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
Box::getHigh
__device__ __host__ T getHigh(int i) const
get the high interval of the box
Definition: Box.hpp:566
Ghost
Definition: Ghost.hpp:40
Point_test
Test structure used for several test.
Definition: Point_test.hpp:106
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:28
eq_id
Definition: eq_solve_common.hpp:15
petsc_solver
This class is able to do Matrix inversion in parallel with PETSC solvers.
Definition: petsc_solver.hpp:102
petsc_solver::setPreconditioner
void setPreconditioner(PCType type)
Set the preconditioner of the linear solver.
Definition: petsc_solver.hpp:1046
petsc_solver::setSolver
void setSolver(KSPType type)
Set the Petsc solver.
Definition: petsc_solver.hpp:918
timer
Class for cpu time benchmarking.
Definition: timer.hpp:28
timer::stop
void stop()
Stop the timer.
Definition: timer.hpp:119
timer::getcputime
double getcputime()
Return the cpu time.
Definition: timer.hpp:142
timer::start
void start()
Start the timer.
Definition: timer.hpp:90
timer::getwct
double getwct()
Return the elapsed real time.
Definition: timer.hpp:130
vector_dist_expression
Main class that encapsulate a vector properties operand to be used for expressions construction.
Definition: vector_dist_operators.hpp:898
vector_dist_subset
Definition: vector_dist_subset.hpp:84
vector_dist_subset::getIds
openfpm::vector< aggregate< int > > & getIds()
Return the ids.
Definition: vector_dist_subset.hpp:170
vector_dist_subset::updateVerlet
void updateVerlet(VerletList< dim, St, opt, Mem_type, shift< dim, St >, vPos_type > &verletList, St r_cut, bool no_se3=false)
Update an existing Verlet List.
Definition: vector_dist_subset.hpp:492
vector_dist_subset::update
void update()
Update the subset indexes.
Definition: vector_dist_subset.hpp:188
vector_dist_ws
Definition: vector_dist_subset.hpp:17
vector_dist
Distributed vector.
Definition: vector_dist.hpp:176
vector_dist::write_frame
bool write_frame(std::string out, size_t iteration, int opt=VTK_WRITER)
Output particle position and properties.
Definition: vector_dist.hpp:2612
vector_dist::getProp
auto getProp(vect_dist_key_dx vec_key) -> decltype(vPrp.template get< id >(vec_key.getKey()))
Get the property of an element.
Definition: vector_dist.hpp:655
vector_dist::updateVerlet
void updateVerlet(VerletList< dim, St, opt, Mem_type, shift< dim, St > > &verletList, St r_cut)
for each particle get the verlet list
Definition: vector_dist.hpp:1556
vector_dist::ghost_get
void ghost_get(size_t opt=WITH_POSITION)
It synchronize the properties and position of the ghost particles.
Definition: vector_dist.hpp:2203
vector_dist::map
void map(size_t opt=NONE)
It move all the particles that does not belong to the local processor to the respective processor.
Definition: vector_dist.hpp:2158
vector_dist::deleteGhost
void deleteGhost()
Delete the particles on the ghost.
Definition: vector_dist.hpp:2559
cub::int
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Definition: dispatch_radix_sort.cuh:336
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:221
eta
[Definition of the system]
Definition: eq_unit_test.cpp:72
state_type_2d_ofp
A 2d Odeint and Openfpm compatible structure.
Definition: OdeIntegrators.hpp:179
sum
It model an expression expr1 + ... exprn.
Definition: sum.hpp:93