OpenFPM  5.2.0
Project that contain the implementation of distributed structures
main.cpp
1 //
2 // Created by Abhinav Singh on 15.03.20.
3 //
4 
5 
61 // Include Vector Expression,Vector Expressions for Subset,DCPSE,Odeint header files
62 #include "config.h"
63 #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
64 #define BOOST_MPL_LIMIT_VECTOR_SIZE 40
65 #include <iostream>
66 #include "DCPSE/DCPSE_op/DCPSE_op.hpp"
67 #include "DCPSE/DCPSE_op/DCPSE_Solver.hpp"
68 #include "Operators/Vector/vector_dist_operators.hpp"
69 #include "Vector/vector_dist_subset.hpp"
70 #include "DCPSE/DCPSE_op/EqnsStruct.hpp"
71 #include "OdeIntegrators/OdeIntegrators.hpp"
73 
95 constexpr int x = 0;
96 constexpr int y = 1;
97 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;
98 
99 double eta = 1.0;
100 double nu = -0.5;
101 double gama = 0.1;
102 double zeta = 0.07;
103 double Ks = 1.0;
104 double Kb = 1.0;
105 double lambda = 0.1;
106 
107 int wr_f;
108 int wr_at;
109 double V_err_eps;
110 
111 void *vectorGlobal=nullptr,*vectorGlobal_bulk=nullptr,*vectorGlobal_boundary=nullptr;
112 const openfpm::vector<std::string>
113 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"};
114 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;
115 typedef vector_dist_ws<2, double,Activegels> vector_type;
116 typedef vector_dist_subset<2, double, Activegels> vector_type2;
118 
145 //Functor to Compute RHS of the time derivative of the polarity
146 template<typename DX,typename DY,typename DXX,typename DXY,typename DYY>
147 struct PolarEv
148 {
149  DX &Dx;
150  DY &Dy;
151  DXX &Dxx;
152  DXY &Dxy;
153  DYY &Dyy;
154  //Constructor
155  PolarEv(DX &Dx,DY &Dy,DXX &Dxx,DXY &Dxy,DYY &Dyy):Dx(Dx),Dy(Dy),Dxx(Dxx),Dxy(Dxy),Dyy(Dyy)
156  {}
157 
158  void operator()( const state_type_2d_ofp &X , state_type_2d_ofp &dxdt , const double t ) const
159  {
160 
161  vector_type &Particles= *(vector_type *) vectorGlobal;
162  vector_type2 &Particles_bulk= *(vector_type2 *) vectorGlobal_bulk;
163 
164  auto Pol=getV<POLARIZATION>(Particles);
165  auto Pol_bulk=getV<POLARIZATION>(Particles_bulk);
166  auto h = getV<MOLFIELD>(Particles);
167  auto u = getV<STRAIN_RATE>(Particles);
168  auto dPol = getV<DPOL>(Particles);
169  auto W = getV<VORTICITY>(Particles);
170  auto delmu = getV<DELMU>(Particles);
171  auto H_p_b = getV<HPB>(Particles);
172  auto r = getV<R>(Particles);
173  auto dPol_bulk = getV<DPOL>(Particles_bulk);
174  Pol[x]=X.data.get<0>();
175  Pol[y]=X.data.get<1>();
176  Particles.ghost_get<POLARIZATION>(SKIP_LABELLING);
177  H_p_b = Pol[x] * Pol[x] + Pol[y] * Pol[y];
178 
179  h[y] = (Pol[x] * (Ks * Dyy(Pol[y]) + Kb * Dxx(Pol[y]) + (Ks - Kb) * Dxy(Pol[x])) -
180  Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y])));
181 
182  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));
183 
184  dPol_bulk[x] = ((h[x] * Pol[x] - h[y] * Pol[y]) / gama + lambda * delmu * Pol[x] -
185  nu * (u[x][x] * Pol[x] + u[x][y] * Pol[y]) + W[x][x] * Pol[x] +
186  W[x][y] * Pol[y]);
187  dPol_bulk[y] = ((h[x] * Pol[y] + h[y] * Pol[x]) / gama + lambda * delmu * Pol[y] -
188  nu * (u[y][x] * Pol[x] + u[y][y] * Pol[y]) + W[y][x] * Pol[x] +
189  W[y][y] * Pol[y]);
190  dxdt.data.get<0>()=dPol[x];
191  dxdt.data.get<1>()=dPol[y];
192  }
193 };
195 
244 // Functor to calculate velocity and move particles with explicit euler
245 template<
246  typename DX,
247  typename DY,
248  typename DX_BULK,
249  typename DY_BULK,
250  typename DXX,
251  typename DXY,
252  typename DYY,
253  typename verletList_type,
254  typename verletListBulk_type>
255 struct CalcVelocity
256 {
257 
258  DX &Dx, &Bulk_Dx;
259  DY &Dy, &Bulk_Dy;
260  DXX &Dxx;
261  DXY &Dxy;
262  DYY &Dyy;
263  verletList_type& verletList;
264  verletListBulk_type& verletList_bulk;
265 
266  double t_old;
267  double rCut;
268  int ctr;
269 
270  //Constructor
271  CalcVelocity(
272  DX &Dx,
273  DY &Dy,
274  DXX &Dxx,
275  DXY &Dxy,
276  DYY &Dyy,
277  DX_BULK &Bulk_Dx,
278  DY_BULK &Bulk_Dy,
279  verletList_type& verletList,
280  verletListBulk_type& verletList_bulk,
281  double rCut
282  ):
283  Dx(Dx),
284  Dy(Dy),
285  Dxx(Dxx),
286  Dxy(Dxy),
287  Dyy(Dyy),
288  Bulk_Dx(Bulk_Dx),
289  Bulk_Dy(Bulk_Dy),
290  verletList(verletList),
291  verletList_bulk(verletList_bulk),
292  rCut(rCut)
293  {
294  t_old = 0.0;
295  ctr = 0;
296  }
297 
298  void operator() (state_type_2d_ofp &state, double t)
299  {
300 
301  double dt = t - t_old;
302  vector_type &Particles= *(vector_type *) vectorGlobal;
303  vector_type2 &Particles_bulk= *(vector_type2 *) vectorGlobal_bulk;
304  vector_type2 &Particles_boundary= *(vector_type2 *) vectorGlobal_boundary;
305  auto &v_cl = create_vcluster();
306 
307  timer tt;
308 
309  auto Pos = getV<POS_PROP>(Particles);
310  auto Pol=getV<POLARIZATION>(Particles);
311  auto V = getV<VELOCITY>(Particles);
312  auto H_p_b = getV<HPB>(Particles);
313 
314 
315  // skip in first time step
316  if (dt != 0) {
317  tt.start();
318  //normalize polarization
319  H_p_b = sqrt(Pol[x] * Pol[x] + Pol[y] * Pol[y]);
320  Pol = Pol / H_p_b;
321  Pos = Pos + (t-t_old)*V;
322  Particles.map();
323  Particles.ghost_get<POLARIZATION, EXTFORCE, DELMU>();
324  Particles_bulk.update();
325  Particles_boundary.update();
326  tt.start();
327  Particles.updateVerlet(verletList,rCut);
328  Particles_bulk.updateVerlet(verletList_bulk,rCut);
329  Dx.update(Particles);
330  Dy.update(Particles);
331  Dxy.update(Particles);
332  Dxx.update(Particles);
333  Dyy.update(Particles);
334 
335  Bulk_Dx.update(Particles_bulk);
336  Bulk_Dy.update(Particles_bulk);
337 
338  state.data.get<0>()=Pol[x];
339  state.data.get<1>()=Pol[y];
340 
341  tt.stop();
342  if (v_cl.rank() == 0) {
343  std::cout << "Updating operators took " << tt.getwct() << " seconds." << std::endl;
344  std::cout << "Time step " << ctr << " : " << t << " over." << std::endl;
345  std::cout << "----------------------------------------------------------" << std::endl;
346  }
347  ctr++;
348 
349  }
350  auto Dyx = Dxy;
351  t_old = t;
352  tt.start();
353 
354  auto & bulk = Particles_bulk.getIds();
355  auto & boundary = Particles_boundary.getIds();
356  auto Pol_bulk=getV<POLARIZATION>(Particles_bulk);
357  auto sigma = getV<STRESS>(Particles);
358  auto r = getV<R>(Particles);
359  auto h = getV<MOLFIELD>(Particles);
360  auto FranckEnergyDensity = getV<FE>(Particles);
361  auto f1 = getV<F1>(Particles);
362  auto f2 = getV<F2>(Particles);
363  auto f3 = getV<F3>(Particles);
364  auto f4 = getV<F4>(Particles);
365  auto f5 = getV<F5>(Particles);
366  auto f6 = getV<F6>(Particles);
367  auto dV = getV<DV>(Particles);
368  auto delmu = getV<DELMU>(Particles); // why is delmu a property, not a constant?
369  auto g = getV<EXTFORCE>(Particles);
370  auto P = getV<PRESSURE>(Particles);
371  auto P_bulk = getV<PRESSURE>(Particles_bulk); //Pressure only on inside
372  auto RHS = getV<VRHS>(Particles);
373  auto RHS_bulk = getV<VRHS>(Particles_bulk);
374  auto div = getV<DIV>(Particles);
375  auto V_t = getV<V_T>(Particles);
376  auto u = getV<STRAIN_RATE>(Particles);
377  auto W = getV<VORTICITY>(Particles);
378 
379  Pol_bulk[x]=state.data.get<0>();
380  Pol_bulk[y]=state.data.get<1>();
381  Particles.ghost_get<POLARIZATION>(SKIP_LABELLING);
382 
383  eq_id x_comp, y_comp;
384  x_comp.setId(0);
385  y_comp.setId(1);
386 
387  int n = 0,nmax = 300,errctr = 0, Vreset = 0;
388  double V_err = 1,V_err_eps = 5 * 1e-3, V_err_old,sum, sum1;
389  std::cout << "Calculate velocity (step t=" << t << ")" << std::endl;
390  tt.start();
391  petsc_solver<double> solverPetsc;
392  solverPetsc.setSolver(KSPGMRES);
393  solverPetsc.setPreconditioner(PCJACOBI);
394  Particles.ghost_get<POLARIZATION>(SKIP_LABELLING);
395  // calculate stress
396  sigma[x][x] =
397  -Ks * Dx(Pol[x]) * Dx(Pol[x]) - Kb * Dx(Pol[y]) * Dx(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dx(Pol[y]);
398  sigma[x][y] =
399  -Ks * Dy(Pol[y]) * Dx(Pol[y]) - Kb * Dy(Pol[x]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dx(Pol[x]);
400  sigma[y][x] =
401  -Ks * Dx(Pol[x]) * Dy(Pol[x]) - Kb * Dx(Pol[y]) * Dy(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dy(Pol[y]);
402  sigma[y][y] =
403  -Ks * Dy(Pol[y]) * Dy(Pol[y]) - Kb * Dy(Pol[x]) * Dy(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dy(Pol[x]);
404  Particles.ghost_get<STRESS>(SKIP_LABELLING);
405 
406  // if R == 0 then set to 1 to avoid division by zero for defects
407  r = Pol[x] * Pol[x] + Pol[y] * Pol[y];
408  for (int j = 0; j < bulk.size(); j++) {
409  auto p = bulk.get<0>(j);
410  Particles.getProp<R>(p) = (Particles.getProp<R>(p) == 0) ? 1 : Particles.getProp<R>(p);
411  }
412  for (int j = 0; j < boundary.size(); j++) {
413  auto p = boundary.get<0>(j);
414  Particles.getProp<R>(p) = (Particles.getProp<R>(p) == 0) ? 1 : Particles.getProp<R>(p);
415  }
416 
417  // calculate traversal molecular field (H_perpendicular)
418  h[y] = (Pol[x] * (Ks * Dyy(Pol[y]) + Kb * Dxx(Pol[y]) + (Ks - Kb) * Dxy(Pol[x])) -
419  Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y])));
420  Particles.ghost_get<MOLFIELD>(SKIP_LABELLING);
421 
422  // calulate FranckEnergyDensity
423  FranckEnergyDensity = (Ks / 2.0) *
424  ((Dx(Pol[x]) * Dx(Pol[x])) + (Dy(Pol[x]) * Dy(Pol[x])) +
425  (Dx(Pol[y]) * Dx(Pol[y])) +
426  (Dy(Pol[y]) * Dy(Pol[y]))) +
427  ((Kb - Ks) / 2.0) * ((Dx(Pol[y]) - Dy(Pol[x])) * (Dx(Pol[y]) - Dy(Pol[x])));
428  Particles.ghost_get<FE>(SKIP_LABELLING);
429 
430  // calculate preactors for LHS of Stokes Equation.
431  f1 = gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
432  f2 = 2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
433  f3 = gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
434  f4 = 2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (r);
435  f5 = 4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (r);
436  f6 = 2.0 * gama * nu * Pol[x] * Pol[y] * Pol[y] * Pol[y] / (r);
437  Particles.ghost_get<F1, F2, F3, F4, F5, F6>(SKIP_LABELLING);
438  texp_v<double> Dxf1 = Dx(f1),Dxf2 = Dx(f2),Dxf3 = Dx(f3),Dxf4 = Dx(f4),Dxf5 = Dx(f5),Dxf6 = Dx(f6),
439  Dyf1 = Dy(f1),Dyf2 = Dy(f2),Dyf3 = Dy(f3),Dyf4 = Dy(f4),Dyf5 = Dy(f5),Dyf6 = Dy(f6);
440 
441  // calculate RHS of Stokes Equation without pressure
442  dV[x] = -0.5 * Dy(h[y]) + zeta * Dx(delmu * Pol[x] * Pol[x]) + zeta * Dy(delmu * Pol[x] * Pol[y]) -
443  zeta * Dx(0.5 * delmu * (Pol[x] * Pol[x] + Pol[y] * Pol[y])) -
444  0.5 * nu * Dx(-2.0 * h[y] * Pol[x] * Pol[y])
445  - 0.5 * nu * Dy(h[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y])) - Dx(sigma[x][x]) -
446  Dy(sigma[x][y]) -
447  g[x]
448  - 0.5 * nu * Dx(-gama * lambda * delmu * (Pol[x] * Pol[x] - Pol[y] * Pol[y]))
449  - 0.5 * Dy(-2.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
450 
451  dV[y] = -0.5 * Dx(-h[y]) + zeta * Dy(delmu * Pol[y] * Pol[y]) + zeta * Dx(delmu * Pol[x] * Pol[y]) -
452  zeta * Dy(0.5 * delmu * (Pol[x] * Pol[x] + Pol[y] * Pol[y])) -
453  0.5 * nu * Dy(2.0 * h[y] * Pol[x] * Pol[y])
454  - 0.5 * nu * Dx(h[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y])) - Dx(sigma[y][x]) -
455  Dy(sigma[y][y]) -
456  g[y]
457  - 0.5 * nu * Dy(gama * lambda * delmu * (Pol[x] * Pol[x] - Pol[y] * Pol[y]))
458  - 0.5 * Dx(-2.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
459  Particles.ghost_get<DV>(SKIP_LABELLING);
460 
461  // Encode LHS of the Stokes Equations
462  auto Stokes1 = eta * (Dxx(V[x]) + Dyy(V[x]))
463  + 0.5 * nu * (Dxf1 * Dx(V[x]) + f1 * Dxx(V[x]))
464  + 0.5 * nu * (Dxf2 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
465  + 0.5 * nu * (Dxf3 * Dy(V[y]) + f3 * Dyx(V[y]))
466  + 0.5 * nu * (Dyf4 * Dx(V[x]) + f4 * Dxy(V[x]))
467  + 0.5 * nu * (Dyf5 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
468  + 0.5 * nu * (Dyf6 * Dy(V[y]) + f6 * Dyy(V[y]));
469  auto Stokes2 = eta * (Dxx(V[y]) + Dyy(V[y]))
470  - 0.5 * nu * (Dyf1 * Dx(V[x]) + f1 * Dxy(V[x]))
471  - 0.5 * nu * (Dyf2 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
472  - 0.5 * nu * (Dyf3 * Dy(V[y]) + f3 * Dyy(V[y]))
473  + 0.5 * nu * (Dxf4 * Dx(V[x]) + f4 * Dxx(V[x]))
474  + 0.5 * nu * (Dxf5 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
475  + 0.5 * nu * (Dxf6 * Dy(V[y]) + f6 * Dyx(V[y]));
476 
477  tt.stop();
478 
479  std::cout << "Init of Velocity took " << tt.getwct() << " seconds." << std::endl;
480 
481  tt.start();
482  V_err = 1;
483  n = 0;
484  errctr = 0;
485  if (Vreset == 1) {
486  P = 0;
487  Vreset = 0;
488  }
489  P=0;
490 
491  // integrate velocity
492  Particles.ghost_get<PRESSURE>(SKIP_LABELLING);
493  RHS_bulk[x] = dV[x] + Bulk_Dx(P);
494  RHS_bulk[y] = dV[y] + Bulk_Dy(P);
495  Particles.ghost_get<VRHS>(SKIP_LABELLING);
496 
497  // prepare solver
498  DCPSE_scheme<equations2d2, vector_type> Solver(Particles);
499  Solver.impose(Stokes1, bulk, RHS[0], x_comp);
500  Solver.impose(Stokes2, bulk, RHS[1], y_comp);
501  Solver.impose(V[x], boundary, 0, x_comp);
502  Solver.impose(V[y], boundary, 0, y_comp);
503  Solver.solve_with_solver(solverPetsc, V[x], V[y]);
504  Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
505  div = -(Dx(V[x]) + Dy(V[y]));
506  P_bulk = P + div;
507 
508  // approximate velocity
509  while (V_err >= V_err_eps && n <= nmax) {
510  Particles.ghost_get<PRESSURE>(SKIP_LABELLING);
511  RHS_bulk[x] = dV[x] + Bulk_Dx(P);
512  RHS_bulk[y] = dV[y] + Bulk_Dy(P);
513  Particles.ghost_get<VRHS>(SKIP_LABELLING);
514  Solver.reset_b();
515  Solver.impose_b(bulk, RHS[0], x_comp);
516  Solver.impose_b(bulk, RHS[1], y_comp);
517  Solver.impose_b(boundary, 0, x_comp);
518  Solver.impose_b(boundary, 0, y_comp);
519  Solver.solve_with_solver(solverPetsc, V[x], V[y]);
520  Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
521  div = -(Dx(V[x]) + Dy(V[y]));
522  P_bulk = P + div;
523  // calculate error
524  sum = 0;
525  sum1 = 0;
526  for (int j = 0; j < bulk.size(); j++) {
527  auto p = bulk.get<0>(j);
528  sum += (Particles.getProp<V_T>(p)[0] - Particles.getProp<VELOCITY>(p)[0]) *
529  (Particles.getProp<V_T>(p)[0] - Particles.getProp<VELOCITY>(p)[0]) +
530  (Particles.getProp<V_T>(p)[1] - Particles.getProp<VELOCITY>(p)[1]) *
531  (Particles.getProp<V_T>(p)[1] - Particles.getProp<VELOCITY>(p)[1]);
532  sum1 += Particles.getProp<VELOCITY>(p)[0] * Particles.getProp<VELOCITY>(p)[0] +
533  Particles.getProp<VELOCITY>(p)[1] * Particles.getProp<VELOCITY>(p)[1];
534  }
535  V_t = V;
536  v_cl.sum(sum);
537  v_cl.sum(sum1);
538  v_cl.execute();
539  sum = sqrt(sum);
540  sum1 = sqrt(sum1);
541  V_err_old = V_err;
542  V_err = sum / sum1;
543  if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
544  errctr++;
545  //alpha_P -= 0.1;
546  } else {
547  errctr = 0;
548  }
549  if (n > 3) {
550  if (errctr > 3) {
551  std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN DIVERGENCE" << std::endl;
552  Vreset = 1;
553  break;
554  } else {
555  Vreset = 0;
556  }
557  }
558  n++;
559 
560  }
561  tt.stop();
562 
563  Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
564  // calculate strain rate
565  u[x][x] = Dx(V[x]);
566  u[x][y] = 0.5 * (Dx(V[y]) + Dy(V[x]));
567  u[y][x] = 0.5 * (Dy(V[x]) + Dx(V[y]));
568  u[y][y] = Dy(V[y]);
569 
570  if (v_cl.rank() == 0) {
571  std::cout << "Rel l2 cgs err in V = " << V_err << " and took " << tt.getwct() << " seconds with " << n
572  << " iterations. dt for the stepper is " << dt
573  << std::endl;
574  }
575 
576  // calculate vorticity
577  W[x][x] = 0;
578  W[x][y] = 0.5 * (Dy(V[x]) - Dx(V[y]));
579  W[y][x] = 0.5 * (Dx(V[y]) - Dy(V[x]));
580  W[y][y] = 0;
581 
582  if (ctr%wr_at==0 || ctr==wr_f){
583  Particles.deleteGhost();
584  Particles.write_frame("Polar",ctr);
585  Particles.ghost_get<POLARIZATION>();
586  }
587  }
588 };
590 
602 int main(int argc, char* argv[])
603 {
604  { openfpm_init(&argc,&argv);
605  timer tt2;
606  tt2.start();
607  size_t Gd = int(std::atof(argv[1]));
608  double tf = std::atof(argv[2]);
609  double dt = tf/std::atof(argv[3]);
610  wr_f=int(std::atof(argv[3]));
611  wr_at=1;
612  V_err_eps = 5e-4;
614 
632  double boxsize = 10;
633  const size_t sz[2] = {Gd, Gd};
634  Box<2, double> box({0, 0}, {boxsize, boxsize});
635  double Lx = box.getHigh(0),Ly = box.getHigh(1);
636  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
637  double spacing = box.getHigh(0) / (sz[0] - 1),rCut = 3.9 * spacing;
638  int ord = 2;
639  Ghost<2, double> ghost(rCut);
640  auto &v_cl = create_vcluster();
641  vector_dist_ws<2, double,Activegels> Particles(0, box, bc, ghost);
642  Particles.setPropNames(PropNAMES);
643 
644  double x0=box.getLow(0), y0=box.getLow(1), x1=box.getHigh(0), y1=box.getHigh(1);
645  auto it = Particles.getGridIterator(sz);
646  while (it.isNext()) {
647  Particles.add();
648  auto key = it.get();
649  double xp = key.get(0) * it.getSpacing(0),yp = key.get(1) * it.getSpacing(1);
650  Particles.getLastPos()[x] = xp;
651  Particles.getLastPos()[y] = yp;
652  if (xp != x0 && yp != y0 && xp != x1 && yp != y1)
653  Particles.getLastSubset(0);
654  else
655  Particles.getLastSubset(1);
656  ++it;
657  }
658  Particles.map();
659  Particles.ghost_get<POLARIZATION>();
660 
661  //auto Pos = getV<POS_PROP>(Particles);
662 
663  auto Pol = getV<POLARIZATION>(Particles);
664  auto V = getV<VELOCITY>(Particles);
665  auto g = getV<EXTFORCE>(Particles);
666  auto P = getV<PRESSURE>(Particles);
667  auto delmu = getV<DELMU>(Particles);
668  auto dPol = getV<DPOL>(Particles);
669 
670  g = 0;delmu = -1.0;P = 0;V = 0;
671  auto it2 = Particles.getDomainIterator();
672  while (it2.isNext()) {
673  auto p = it2.get();
674  Point<2, double> xp = Particles.getPos(p);
675  Particles.getProp<POLARIZATION>(p)[x] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
676  Particles.getProp<POLARIZATION>(p)[y] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
677  ++it2;
678  }
679  Particles.ghost_get<POLARIZATION,EXTFORCE,DELMU>(SKIP_LABELLING);
681 
695 
696  vector_dist_subset<2, double, Activegels> Particles_bulk(Particles,0);
697  vector_dist_subset<2, double, Activegels> Particles_boundary(Particles,1);
698  auto & bulk = Particles_bulk.getIds();
699  auto & boundary = Particles_boundary.getIds();
700 
701  auto P_bulk = getV<PRESSURE>(Particles_bulk);//Pressure only on inside
702  auto Pol_bulk = getV<POLARIZATION>(Particles_bulk);;
703  auto dPol_bulk = getV<DPOL>(Particles_bulk);
704  auto dV_bulk = getV<DV>(Particles_bulk);
705  auto RHS_bulk = getV<VRHS>(Particles_bulk);
706  auto div_bulk = getV<DIV>(Particles_bulk);
707 
708  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
709  auto verletList_bulk = Particles_bulk.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
710 
711  Derivative_x<decltype(verletList)> Dx(Particles, verletList, ord,rCut);
712  Derivative_y<decltype(verletList)> Dy(Particles, verletList, ord, rCut);
713  Derivative_xy<decltype(verletList)> Dxy(Particles, verletList, ord, rCut);
714 
715  Derivative_x<decltype(verletList_bulk)> Bulk_Dx(Particles, Particles_bulk, verletList_bulk, ord, rCut);
716  Derivative_y<decltype(verletList_bulk)> Bulk_Dy(Particles, Particles_bulk, verletList_bulk, ord, rCut);
717  auto Dyx = Dxy;
718  Derivative_xx<decltype(verletList)> Dxx(Particles, verletList, ord, rCut);
719  Derivative_yy<decltype(verletList)> Dyy(Particles, verletList, ord, rCut);
720 
721  boost::numeric::odeint::runge_kutta4< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
722 
723  vectorGlobal=(void *) &Particles;
724  vectorGlobal_bulk=(void *) &Particles_bulk;
725  vectorGlobal_boundary=(void *) &Particles_boundary;
726 
727  PolarEv<
728  Derivative_x<decltype(verletList)>,
729  Derivative_y<decltype(verletList)>,
730  Derivative_xx<decltype(verletList)>,
731  Derivative_xy<decltype(verletList)>,
732  Derivative_yy<decltype(verletList)>>
733  System(Dx,Dy,Dxx,Dxy,Dyy);
734 
735  CalcVelocity<
736  Derivative_x<decltype(verletList)>,
737  Derivative_y<decltype(verletList)>,
738  Derivative_x<decltype(verletList_bulk)>,
739  Derivative_y<decltype(verletList_bulk)>,
740  Derivative_xx<decltype(verletList)>,
741  Derivative_xy<decltype(verletList)>,
742  Derivative_yy<decltype(verletList)>,
743  decltype(verletList),
744  decltype(verletList_bulk)>
745  CalcVelocityObserver(Dx,Dy,Dxx,Dxy,Dyy,Bulk_Dx,Bulk_Dy,verletList,verletList_bulk,rCut);
746 
747  state_type_2d_ofp tPol;
748  tPol.data.get<0>()=Pol[x];
749  tPol.data.get<1>()=Pol[y];
751 
752 
753  timer tt;
754  timer tt3;
755  dPol = Pol;
756  double V_err = 1, V_err_old;
757  double tim=0;
781  // intermediate time steps
782  std::vector<double> inter_times;
783 
784  size_t steps = integrate_const(rk4 , System , tPol , tim , tf , dt, CalcVelocityObserver);
785 
786 
787  std::cout << "Time steps: " << steps << std::endl;
788 
789  Pol_bulk[x]=tPol.data.get<0>();
790  Pol_bulk[y]=tPol.data.get<1>();
791 
792  Particles.deleteGhost();
793  Particles.write("Polar_Last");
794  Dx.deallocate(Particles);
795  Dy.deallocate(Particles);
796  Dxy.deallocate(Particles);
797  Dxx.deallocate(Particles);
798  Dyy.deallocate(Particles);
799  Bulk_Dx.deallocate(Particles_bulk);
800  Bulk_Dy.deallocate(Particles_bulk);
801  std::cout.precision(17);
802  tt2.stop();
803  if (v_cl.rank() == 0) {
804  std::cout << "The simulation took " << tt2.getcputime() << "(CPU) ------ " << tt2.getwct()
805  << "(Wall) Seconds.";
806  }
807  }
808  openfpm_finalize();
810 }
811 
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