OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
88constexpr int x = 0;
89constexpr int y = 1;
90constexpr 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
92double eta = 1.0;
93double nu = -0.5;
94double gama = 0.1;
95double zeta = 0.07;
96double Ks = 1.0;
97double Kb = 1.0;
98double lambda = 0.1;
99
100int wr_f;
101int wr_at;
102double V_err_eps;
103
104void *vectorGlobal=nullptr,*vectorGlobal_bulk=nullptr,*vectorGlobal_boundary=nullptr;
105const openfpm::vector<std::string>
106PropNAMES={"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"};
107typedef 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;
108typedef vector_dist_ws<2, double,Activegels> vector_type;
109typedef vector_dist_subset<2, double, Activegels> vector_type2;
111
138//Functor to Compute RHS of the time derivative of the polarity
139template<typename DX,typename DY,typename DXX,typename DXY,typename DYY>
140struct 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
238template<typename DX,typename DY,typename DXX,typename DXY,typename DYY>
239struct CalcVelocity
240{
241
242 DX &Dx, &Bulk_Dx;
243 DY &Dy, &Bulk_Dy;
244 DXX &Dxx;
245 DXY &Dxy;
246 DYY &Dyy;
247
248 double t_old;
249 int ctr;
250
251 //Constructor
252 CalcVelocity(DX &Dx,DY &Dy,DXX &Dxx,DXY &Dxy,DYY &Dyy,DX &Bulk_Dx,DY &Bulk_Dy):Dx(Dx),Dy(Dy),Dxx(Dxx),Dxy(Dxy),Dyy(Dyy),Bulk_Dx(Bulk_Dx),Bulk_Dy(Bulk_Dy)
253 {
254 t_old = 0.0;
255 ctr = 0;
256 }
257
258 void operator() (state_type_2d_ofp &state, double t)
259 {
260
261 double dt = t - t_old;
262 vector_type &Particles= *(vector_type *) vectorGlobal;
263 vector_type2 &Particles_bulk= *(vector_type2 *) vectorGlobal_bulk;
264 vector_type2 &Particles_boundary= *(vector_type2 *) vectorGlobal_boundary;
265 auto &v_cl = create_vcluster();
266
267 timer tt;
268
269 auto Pos = getV<PROP_POS>(Particles);
270 auto Pol=getV<POLARIZATION>(Particles);
271 auto V = getV<VELOCITY>(Particles);
272 auto H_p_b = getV<HPB>(Particles);
273
274
275 // skip in first time step
276 if (dt != 0) {
277 tt.start();
278 //normalize polarization
279 H_p_b = sqrt(Pol[x] * Pol[x] + Pol[y] * Pol[y]);
280 Pol = Pol / H_p_b;
281 Pos = Pos + (t-t_old)*V;
282 Particles.map();
283 Particles.ghost_get<POLARIZATION, EXTFORCE, DELMU>();
284 Particles_bulk.update();
285 Particles_boundary.update();
286 tt.start();
287 Dx.update(Particles);
288 Dy.update(Particles);
289 Dxy.update(Particles);
290 Dxx.update(Particles);
291 Dyy.update(Particles);
292
293 Bulk_Dx.update(Particles_bulk);
294 Bulk_Dy.update(Particles_bulk);
295
296 state.data.get<0>()=Pol[x];
297 state.data.get<1>()=Pol[y];
298
299 tt.stop();
300 if (v_cl.rank() == 0) {
301 std::cout << "Updating operators took " << tt.getwct() << " seconds." << std::endl;
302 std::cout << "Time step " << ctr << " : " << t << " over." << std::endl;
303 std::cout << "----------------------------------------------------------" << std::endl;
304 }
305 ctr++;
306
307 }
308 auto Dyx = Dxy;
309 t_old = t;
310 tt.start();
311
312 auto & bulk = Particles_bulk.getIds();
313 auto & boundary = Particles_boundary.getIds();
314 auto Pol_bulk=getV<POLARIZATION>(Particles_bulk);
315 auto sigma = getV<STRESS>(Particles);
316 auto r = getV<R>(Particles);
317 auto h = getV<MOLFIELD>(Particles);
318 auto FranckEnergyDensity = getV<FE>(Particles);
319 auto f1 = getV<F1>(Particles);
320 auto f2 = getV<F2>(Particles);
321 auto f3 = getV<F3>(Particles);
322 auto f4 = getV<F4>(Particles);
323 auto f5 = getV<F5>(Particles);
324 auto f6 = getV<F6>(Particles);
325 auto dV = getV<DV>(Particles);
326 auto delmu = getV<DELMU>(Particles); // why is delmu a property, not a constant?
327 auto g = getV<EXTFORCE>(Particles);
328 auto P = getV<PRESSURE>(Particles);
329 auto P_bulk = getV<PRESSURE>(Particles_bulk); //Pressure only on inside
330 auto RHS = getV<VRHS>(Particles);
331 auto RHS_bulk = getV<VRHS>(Particles_bulk);
332 auto div = getV<DIV>(Particles);
333 auto V_t = getV<V_T>(Particles);
334 auto u = getV<STRAIN_RATE>(Particles);
335 auto W = getV<VORTICITY>(Particles);
336
337 Pol_bulk[x]=state.data.get<0>();
338 Pol_bulk[y]=state.data.get<1>();
339 Particles.ghost_get<POLARIZATION>(SKIP_LABELLING);
340
341 eq_id x_comp, y_comp;
342 x_comp.setId(0);
343 y_comp.setId(1);
344
345 int n = 0,nmax = 300,errctr = 0, Vreset = 0;
346 double V_err = 1,V_err_eps = 5 * 1e-3, V_err_old,sum, sum1;
347 std::cout << "Calculate velocity (step t=" << t << ")" << std::endl;
348 tt.start();
349 petsc_solver<double> solverPetsc;
350 solverPetsc.setSolver(KSPGMRES);
351 solverPetsc.setPreconditioner(PCJACOBI);
352 Particles.ghost_get<POLARIZATION>(SKIP_LABELLING);
353 // calculate stress
354 sigma[x][x] =
355 -Ks * Dx(Pol[x]) * Dx(Pol[x]) - Kb * Dx(Pol[y]) * Dx(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dx(Pol[y]);
356 sigma[x][y] =
357 -Ks * Dy(Pol[y]) * Dx(Pol[y]) - Kb * Dy(Pol[x]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dx(Pol[x]);
358 sigma[y][x] =
359 -Ks * Dx(Pol[x]) * Dy(Pol[x]) - Kb * Dx(Pol[y]) * Dy(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dy(Pol[y]);
360 sigma[y][y] =
361 -Ks * Dy(Pol[y]) * Dy(Pol[y]) - Kb * Dy(Pol[x]) * Dy(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dy(Pol[x]);
362 Particles.ghost_get<STRESS>(SKIP_LABELLING);
363
364 // if R == 0 then set to 1 to avoid division by zero for defects
365 r = Pol[x] * Pol[x] + Pol[y] * Pol[y];
366 for (int j = 0; j < bulk.size(); j++) {
367 auto p = bulk.get<0>(j);
368 Particles.getProp<R>(p) = (Particles.getProp<R>(p) == 0) ? 1 : Particles.getProp<R>(p);
369 }
370 for (int j = 0; j < boundary.size(); j++) {
371 auto p = boundary.get<0>(j);
372 Particles.getProp<R>(p) = (Particles.getProp<R>(p) == 0) ? 1 : Particles.getProp<R>(p);
373 }
374
375 // calculate traversal molecular field (H_perpendicular)
376 h[y] = (Pol[x] * (Ks * Dyy(Pol[y]) + Kb * Dxx(Pol[y]) + (Ks - Kb) * Dxy(Pol[x])) -
377 Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y])));
378 Particles.ghost_get<MOLFIELD>(SKIP_LABELLING);
379
380 // calulate FranckEnergyDensity
381 FranckEnergyDensity = (Ks / 2.0) *
382 ((Dx(Pol[x]) * Dx(Pol[x])) + (Dy(Pol[x]) * Dy(Pol[x])) +
383 (Dx(Pol[y]) * Dx(Pol[y])) +
384 (Dy(Pol[y]) * Dy(Pol[y]))) +
385 ((Kb - Ks) / 2.0) * ((Dx(Pol[y]) - Dy(Pol[x])) * (Dx(Pol[y]) - Dy(Pol[x])));
386 Particles.ghost_get<FE>(SKIP_LABELLING);
387
388 // calculate preactors for LHS of Stokes Equation.
389 f1 = gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
390 f2 = 2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
391 f3 = gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (r);
392 f4 = 2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (r);
393 f5 = 4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (r);
394 f6 = 2.0 * gama * nu * Pol[x] * Pol[y] * Pol[y] * Pol[y] / (r);
395 Particles.ghost_get<F1, F2, F3, F4, F5, F6>(SKIP_LABELLING);
396 texp_v<double> Dxf1 = Dx(f1),Dxf2 = Dx(f2),Dxf3 = Dx(f3),Dxf4 = Dx(f4),Dxf5 = Dx(f5),Dxf6 = Dx(f6),
397 Dyf1 = Dy(f1),Dyf2 = Dy(f2),Dyf3 = Dy(f3),Dyf4 = Dy(f4),Dyf5 = Dy(f5),Dyf6 = Dy(f6);
398
399 // calculate RHS of Stokes Equation without pressure
400 dV[x] = -0.5 * Dy(h[y]) + zeta * Dx(delmu * Pol[x] * Pol[x]) + zeta * Dy(delmu * Pol[x] * Pol[y]) -
401 zeta * Dx(0.5 * delmu * (Pol[x] * Pol[x] + Pol[y] * Pol[y])) -
402 0.5 * nu * Dx(-2.0 * h[y] * Pol[x] * Pol[y])
403 - 0.5 * nu * Dy(h[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y])) - Dx(sigma[x][x]) -
404 Dy(sigma[x][y]) -
405 g[x]
406 - 0.5 * nu * Dx(-gama * lambda * delmu * (Pol[x] * Pol[x] - Pol[y] * Pol[y]))
407 - 0.5 * Dy(-2.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
408
409 dV[y] = -0.5 * Dx(-h[y]) + zeta * Dy(delmu * Pol[y] * Pol[y]) + zeta * Dx(delmu * Pol[x] * Pol[y]) -
410 zeta * Dy(0.5 * delmu * (Pol[x] * Pol[x] + Pol[y] * Pol[y])) -
411 0.5 * nu * Dy(2.0 * h[y] * Pol[x] * Pol[y])
412 - 0.5 * nu * Dx(h[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y])) - Dx(sigma[y][x]) -
413 Dy(sigma[y][y]) -
414 g[y]
415 - 0.5 * nu * Dy(gama * lambda * delmu * (Pol[x] * Pol[x] - Pol[y] * Pol[y]))
416 - 0.5 * Dx(-2.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
417 Particles.ghost_get<DV>(SKIP_LABELLING);
418
419 // Encode LHS of the Stokes Equations
420 auto Stokes1 = eta * (Dxx(V[x]) + Dyy(V[x]))
421 + 0.5 * nu * (Dxf1 * Dx(V[x]) + f1 * Dxx(V[x]))
422 + 0.5 * nu * (Dxf2 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
423 + 0.5 * nu * (Dxf3 * Dy(V[y]) + f3 * Dyx(V[y]))
424 + 0.5 * nu * (Dyf4 * Dx(V[x]) + f4 * Dxy(V[x]))
425 + 0.5 * nu * (Dyf5 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
426 + 0.5 * nu * (Dyf6 * Dy(V[y]) + f6 * Dyy(V[y]));
427 auto Stokes2 = eta * (Dxx(V[y]) + Dyy(V[y]))
428 - 0.5 * nu * (Dyf1 * Dx(V[x]) + f1 * Dxy(V[x]))
429 - 0.5 * nu * (Dyf2 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
430 - 0.5 * nu * (Dyf3 * Dy(V[y]) + f3 * Dyy(V[y]))
431 + 0.5 * nu * (Dxf4 * Dx(V[x]) + f4 * Dxx(V[x]))
432 + 0.5 * nu * (Dxf5 * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
433 + 0.5 * nu * (Dxf6 * Dy(V[y]) + f6 * Dyx(V[y]));
434
435 tt.stop();
436
437 std::cout << "Init of Velocity took " << tt.getwct() << " seconds." << std::endl;
438
439 tt.start();
440 V_err = 1;
441 n = 0;
442 errctr = 0;
443 if (Vreset == 1) {
444 P = 0;
445 Vreset = 0;
446 }
447 P=0;
448
449 // integrate velocity
450 Particles.ghost_get<PRESSURE>(SKIP_LABELLING);
451 RHS_bulk[x] = dV[x] + Bulk_Dx(P);
452 RHS_bulk[y] = dV[y] + Bulk_Dy(P);
453 Particles.ghost_get<VRHS>(SKIP_LABELLING);
454
455 // prepare solver
456 DCPSE_scheme<equations2d2, vector_type> Solver(Particles);
457 Solver.impose(Stokes1, bulk, RHS[0], x_comp);
458 Solver.impose(Stokes2, bulk, RHS[1], y_comp);
459 Solver.impose(V[x], boundary, 0, x_comp);
460 Solver.impose(V[y], boundary, 0, y_comp);
461 Solver.solve_with_solver(solverPetsc, V[x], V[y]);
462 Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
463 div = -(Dx(V[x]) + Dy(V[y]));
464 P_bulk = P + div;
465
466 // approximate velocity
467 while (V_err >= V_err_eps && n <= nmax) {
468 Particles.ghost_get<PRESSURE>(SKIP_LABELLING);
469 RHS_bulk[x] = dV[x] + Bulk_Dx(P);
470 RHS_bulk[y] = dV[y] + Bulk_Dy(P);
471 Particles.ghost_get<VRHS>(SKIP_LABELLING);
472 Solver.reset_b();
473 Solver.impose_b(bulk, RHS[0], x_comp);
474 Solver.impose_b(bulk, RHS[1], y_comp);
475 Solver.impose_b(boundary, 0, x_comp);
476 Solver.impose_b(boundary, 0, y_comp);
477 Solver.solve_with_solver(solverPetsc, V[x], V[y]);
478 Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
479 div = -(Dx(V[x]) + Dy(V[y]));
480 P_bulk = P + div;
481 // calculate error
482 sum = 0;
483 sum1 = 0;
484 for (int j = 0; j < bulk.size(); j++) {
485 auto p = bulk.get<0>(j);
486 sum += (Particles.getProp<V_T>(p)[0] - Particles.getProp<VELOCITY>(p)[0]) *
487 (Particles.getProp<V_T>(p)[0] - Particles.getProp<VELOCITY>(p)[0]) +
488 (Particles.getProp<V_T>(p)[1] - Particles.getProp<VELOCITY>(p)[1]) *
489 (Particles.getProp<V_T>(p)[1] - Particles.getProp<VELOCITY>(p)[1]);
490 sum1 += Particles.getProp<VELOCITY>(p)[0] * Particles.getProp<VELOCITY>(p)[0] +
491 Particles.getProp<VELOCITY>(p)[1] * Particles.getProp<VELOCITY>(p)[1];
492 }
493 V_t = V;
494 v_cl.sum(sum);
495 v_cl.sum(sum1);
496 v_cl.execute();
497 sum = sqrt(sum);
498 sum1 = sqrt(sum1);
499 V_err_old = V_err;
500 V_err = sum / sum1;
501 if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
502 errctr++;
503 //alpha_P -= 0.1;
504 } else {
505 errctr = 0;
506 }
507 if (n > 3) {
508 if (errctr > 3) {
509 std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN DIVERGENCE" << std::endl;
510 Vreset = 1;
511 break;
512 } else {
513 Vreset = 0;
514 }
515 }
516 n++;
517
518 }
519 tt.stop();
520
521 Particles.ghost_get<VELOCITY>(SKIP_LABELLING);
522 // calculate strain rate
523 u[x][x] = Dx(V[x]);
524 u[x][y] = 0.5 * (Dx(V[y]) + Dy(V[x]));
525 u[y][x] = 0.5 * (Dy(V[x]) + Dx(V[y]));
526 u[y][y] = Dy(V[y]);
527
528 if (v_cl.rank() == 0) {
529 std::cout << "Rel l2 cgs err in V = " << V_err << " and took " << tt.getwct() << " seconds with " << n
530 << " iterations. dt for the stepper is " << dt
531 << std::endl;
532 }
533
534 // calculate vorticity
535 W[x][x] = 0;
536 W[x][y] = 0.5 * (Dy(V[x]) - Dx(V[y]));
537 W[y][x] = 0.5 * (Dx(V[y]) - Dy(V[x]));
538 W[y][y] = 0;
539
540 if (ctr%wr_at==0 || ctr==wr_f){
541 Particles.deleteGhost();
542 Particles.write_frame("Polar",ctr);
543 Particles.ghost_get<POLARIZATION>();
544 }
545 }
546};
548
560int main(int argc, char* argv[])
561{
562 { openfpm_init(&argc,&argv);
563 timer tt2;
564 tt2.start();
565 size_t Gd = int(std::atof(argv[1]));
566 double tf = std::atof(argv[2]);
567 double dt = tf/std::atof(argv[3]);
568 wr_f=int(std::atof(argv[3]));
569 wr_at=1;
570 V_err_eps = 5e-4;
572
590 double boxsize = 10;
591 const size_t sz[2] = {Gd, Gd};
592 Box<2, double> box({0, 0}, {boxsize, boxsize});
593 double Lx = box.getHigh(0),Ly = box.getHigh(1);
594 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
595 double spacing = box.getHigh(0) / (sz[0] - 1),rCut = 3.9 * spacing;
596 int ord = 2;
597 Ghost<2, double> ghost(rCut);
598 auto &v_cl = create_vcluster();
599 vector_dist_ws<2, double,Activegels> Particles(0, box, bc, ghost);
600 Particles.setPropNames(PropNAMES);
601
602 double x0=box.getLow(0), y0=box.getLow(1), x1=box.getHigh(0), y1=box.getHigh(1);
603 auto it = Particles.getGridIterator(sz);
604 while (it.isNext()) {
605 Particles.add();
606 auto key = it.get();
607 double xp = key.get(0) * it.getSpacing(0),yp = key.get(1) * it.getSpacing(1);
608 Particles.getLastPos()[x] = xp;
609 Particles.getLastPos()[y] = yp;
610 if (xp != x0 && yp != y0 && xp != x1 && yp != y1)
611 Particles.getLastSubset(0);
612 else
613 Particles.getLastSubset(1);
614 ++it;
615 }
616 Particles.map();
617 Particles.ghost_get<POLARIZATION>();
618
619 //auto Pos = getV<PROP_POS>(Particles);
620
621 auto Pol = getV<POLARIZATION>(Particles);
622 auto V = getV<VELOCITY>(Particles);
623 auto g = getV<EXTFORCE>(Particles);
624 auto P = getV<PRESSURE>(Particles);
625 auto delmu = getV<DELMU>(Particles);
626 auto dPol = getV<DPOL>(Particles);
627
628 g = 0;delmu = -1.0;P = 0;V = 0;
629 auto it2 = Particles.getDomainIterator();
630 while (it2.isNext()) {
631 auto p = it2.get();
632 Point<2, double> xp = Particles.getPos(p);
633 Particles.getProp<POLARIZATION>(p)[x] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
634 Particles.getProp<POLARIZATION>(p)[y] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
635 ++it2;
636 }
637 Particles.ghost_get<POLARIZATION,EXTFORCE,DELMU>(SKIP_LABELLING);
639
653
654 vector_dist_subset<2, double, Activegels> Particles_bulk(Particles,0);
655 vector_dist_subset<2, double, Activegels> Particles_boundary(Particles,1);
656 auto & bulk = Particles_bulk.getIds();
657 auto & boundary = Particles_boundary.getIds();
658
659 auto P_bulk = getV<PRESSURE>(Particles_bulk);//Pressure only on inside
660 auto Pol_bulk = getV<POLARIZATION>(Particles_bulk);;
661 auto dPol_bulk = getV<DPOL>(Particles_bulk);
662 auto dV_bulk = getV<DV>(Particles_bulk);
663 auto RHS_bulk = getV<VRHS>(Particles_bulk);
664 auto div_bulk = getV<DIV>(Particles_bulk);
665
666 Derivative_x Dx(Particles,ord,rCut), Bulk_Dx(Particles_bulk,ord,rCut);
667 Derivative_y Dy(Particles, ord, rCut), Bulk_Dy(Particles_bulk, ord,rCut);
668 Derivative_xy Dxy(Particles, ord, rCut);
669 auto Dyx = Dxy;
670 Derivative_xx Dxx(Particles, ord, rCut);
671 Derivative_yy Dyy(Particles, ord, rCut);
672
673 boost::numeric::odeint::runge_kutta4< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
674
675 vectorGlobal=(void *) &Particles;
676 vectorGlobal_bulk=(void *) &Particles_bulk;
677 vectorGlobal_boundary=(void *) &Particles_boundary;
678
679 PolarEv<Derivative_x,Derivative_y,Derivative_xx,Derivative_xy,Derivative_yy> System(Dx,Dy,Dxx,Dxy,Dyy);
680 CalcVelocity<Derivative_x,Derivative_y,Derivative_xx,Derivative_xy,Derivative_yy> CalcVelocityObserver(Dx,Dy,Dxx,Dxy,Dyy,Bulk_Dx,Bulk_Dy);
681
682 state_type_2d_ofp tPol;
683 tPol.data.get<0>()=Pol[x];
684 tPol.data.get<1>()=Pol[y];
686
687
688 timer tt;
689 timer tt3;
690 dPol = Pol;
691 double V_err = 1, V_err_old;
692 double tim=0;
716 // intermediate time steps
717 std::vector<double> inter_times;
718
719 size_t steps = integrate_const(rk4 , System , tPol , tim , tf , dt, CalcVelocityObserver);
720
721
722 std::cout << "Time steps: " << steps << std::endl;
723
724 Pol_bulk[x]=tPol.data.get<0>();
725 Pol_bulk[y]=tPol.data.get<1>();
726
727 Particles.deleteGhost();
728 Particles.write("Polar_Last");
729 Dx.deallocate(Particles);
730 Dy.deallocate(Particles);
731 Dxy.deallocate(Particles);
732 Dxx.deallocate(Particles);
733 Dyy.deallocate(Particles);
734 Bulk_Dx.deallocate(Particles_bulk);
735 Bulk_Dy.deallocate(Particles_bulk);
736 std::cout.precision(17);
737 tt2.stop();
738 if (v_cl.rank() == 0) {
739 std::cout << "The simulation took " << tt2.getcputime() << "(CPU) ------ " << tt2.getwct()
740 << "(Wall) Seconds.";
741 }
742 }
743 openfpm_finalize();
745}
746
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Box::getHigh
__device__ __host__ T getHigh(int i) const
get the high interval of the box
Definition Box.hpp:567
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
System
System of equations.
Definition System.hpp:24
eq_id
Definition eq_solve_common.hpp:15
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition map_vector.hpp:203
petsc_solver
In case T does not match the PETSC precision compilation create a stub structure.
Definition petsc_solver.hpp:56
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:906
vector_dist_subset
Definition vector_dist_subset.hpp:84
vector_dist_subset::update
void update()
Update the subset indexes.
Definition vector_dist_subset.hpp:188
vector_dist_subset::getIds
openfpm::vector< aggregate< int > > & getIds()
Return the ids.
Definition vector_dist_subset.hpp:170
vector_dist_ws
Definition vector_dist_subset.hpp:17
vector_dist
Distributed vector.
Definition vector_dist.hpp:305
vector_dist::getProp
auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get< id >(vec_key.getKey()))
Get the property of an element.
Definition vector_dist.hpp:826
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:2481
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:2436
vector_dist::deleteGhost
void deleteGhost()
Delete the particles on the ghost.
Definition vector_dist.hpp:2820
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition aggregate.hpp:215
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:128
sum
It model an expression expr1 + ... exprn.
Definition sum.hpp:93