OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
main2.cpp
1 //
2 // Created by Abhinav Singh @absingh on 26/04/2021
3 //
41 
43 // Include Vector Expression,Vector Expressions for Subset,DCPSE,Odeint header files
44 #include "Operators/Vector/vector_dist_operators.hpp"
45 #include "Vector/vector_dist_subset.hpp"
46 #include "DCPSE/DCPSE_op/DCPSE_op.hpp"
47 #include "OdeIntegrators/OdeIntegrators.hpp"
49 
50 
72 
74 constexpr int x = 0;
75 constexpr int y = 1;
76 
77 double dt=1e-2,tf=1.0;
78 
79 void *PointerDistGlobal, *PointerDistSubset,*PointerDistSubset2;
80 
85 
114 template<typename DXX,typename DYY>
116 struct RHSFunctor
117 {
118  //Intializing the operators
119  DXX &Dxx;
120  DYY &Dyy;
121  //Constructor
122  RHSFunctor(DXX &Dxx,DYY &Dyy):Dxx(Dxx),Dyy(Dyy)
123  {}
124 
125  void operator()( const state_type_2d_ofp &X , state_type_2d_ofp &dxdt , const double t ) const
126  {
127  //Casting the pointers to OpenFPM vector distributions
128  dist_vector_type &Particles= *(dist_vector_type *) PointerDistGlobal;
129  dist_vector_subset_type &Particles_bulk= *(dist_vector_subset_type *) PointerDistSubset;
130 
131  //Aliasing the properties.
132  auto C = getV<0>(Particles);
133  auto C_bulk = getV<0>(Particles_bulk);
134  auto dC = getV<2>(Particles);
135  auto dC_bulk = getV<2>(Particles_bulk);
136 
137  //Since the state vector does not have enough information to compute derivates, we copy them back to the Particle Distribution.
138  //In this way we can compute all the things that are required to be computed in each stage of the ode integrator.
139  // Note that we are going to use bulk expressions as needed. Whenever we don't want to update the boundaries (Due to the boundary conditions).
140 
141  //These expressions only update the bulk values of C.
142  C_bulk[x]=X.data.get<0>();
143  C_bulk[y]=X.data.get<1>();
144  Particles.ghost_get<0>();
145 
146  // We do the RHS computations for the Laplacian (Updating bulk only).
147  dC_bulk[x] = 0.1*(Dxx(C[x])+Dyy(C[x]));
148  dC_bulk[y] = 0.1*(Dxx(C[y])+Dyy(C[y]));
149 
150  //We copy back to the dxdt state_type for Odeint
151  dxdt.data.get<0>()=dC[x];
152  dxdt.data.get<1>()=dC[y];
153  }
154 };
156 
157 
185 template<typename DXX,typename DYY>
188 
189  DXX &Dxx;
190  DYY &Dyy;
191  int ctr;
192  double t_old;
193 
194  //Constructor
195  ObserverFunctor(DXX &Dxx, DYY &Dyy) : Dxx(Dxx), Dyy(Dyy) {
196  //a counter for counting the np. of steps
197  ctr = 0;
198  //Starting with t=0, we compute the step size take by t-t_old. So for the first observed step size is what we provide. Which will be 0-(-dt)=dt.
199  t_old = -dt;
200  }
201 
202  void operator()(state_type_2d_ofp &X, double t) {
203  //Casting the pointers to OpenFPM vector distributions
204  dist_vector_type &Particles = *(dist_vector_type *) PointerDistGlobal;
205  dist_vector_subset_type &Particles_bulk = *(dist_vector_subset_type *) PointerDistSubset;
206  dist_vector_subset_type &Particles_boundary = *(dist_vector_subset_type *) PointerDistSubset2;
207 
208  //Aliasing the position and properties.
209  auto Pos = getV<PROP_POS>(Particles);
210  auto Concentration = getV<0>(Particles);
211  auto Velocity = getV<1>(Particles);
212  auto Concentration_bulk = getV<0>(Particles_bulk);
213  auto Velocity_bulk = getV<1>(Particles_bulk);
214 
215  //We would like to move the particles after t=0. (Remember, odeint calls the observer before taking the step.)
216  if (t != 0) {
217  //copy back the state after the time step into data structure. This is required as we are going to move the particles and the distributed state can be resized correctly (by copying back after map). Also this expression preserves the boundary condition on concentration.
218  Concentration_bulk[x] = X.data.get<0>();
219  Concentration_bulk[y] = X.data.get<1>();
220  //computing the velocity and move the particles
221  Velocity_bulk[x] = -Pos[y] * exp(-10.0 * (Pos[x] * Pos[x] + Pos[y] * Pos[y]));
222  Velocity_bulk[y] = Pos[x] * exp(-10.0 * (Pos[x] * Pos[x] + Pos[y] * Pos[y]));
223  Pos = Pos + dt * Velocity;
224  //Map and ghost_get is required after moving particles.
225  Particles.map();
226  Particles.ghost_get<0>();
227  //Updating the subset and operators based on new positions
228  Particles_bulk.update();
229  Particles_boundary.update();
230  Dxx.update(Particles);
231  Dyy.update(Particles);
232 
233  //Since we did Map, we assign the Odeint state again.
234  X.data.get<0>() = Concentration[x];
235  X.data.get<1>() = Concentration[y];
236  //counting the step number
237 
238  }
239  ctr++;
240  std::cout<<"Taking a step at t="<<t<<" with dt="<<t-t_old<<std::endl;
241  t_old=t;
242  //Writing the data after every step
243  Particles.deleteGhost();
244  Particles.write_frame("PDE_sol", ctr);
245  Particles.ghost_get<0>();
246 
247  }
248 };
250 
261 int main(int argc, char *argv[])
263 {
264  // initialize library
265  openfpm_init(&argc, &argv);
267 
280  const size_t sz[2] = {41, 41};
282  Box<2, double> box({-1, -1}, {1.0, 1.0});
283  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
284  double spacing[2];
285  spacing[0] = 2.0 / (sz[0] - 1);
286  spacing[1] = 2.0 / (sz[1] - 1);
287  double rCut = 3.1 * spacing[0];
288  Ghost<2, double> ghost(rCut);
289 
290  dist_vector_type Particles(0, box, bc, ghost);
291  Particles.setPropNames({"Concentration", "Velocity", "TempConcentration"});
292 
293  auto it = Particles.getGridIterator(sz);
294  while (it.isNext()) {
295  Particles.add();
296  auto key = it.get();
297  double x = -1.0 + key.get(0) * spacing[0];
298  Particles.getLastPos()[0] = x;
299  double y = -1.0 + key.get(1) * spacing[1];
300  Particles.getLastPos()[1] = y;
301  //We detect Boundaries. By labelling the subset to be 0 for bulk and 1 for the boundary. Subsets will be constructed later based on the label.
302  if (x != -1.0 && x != 1.0 && y != -1.0 && y != 1) {
303  Particles.getLastSubset(0);
304  } else {
305  Particles.getLastSubset(1);
306  }
307  // Here fill the Initial value of the concentration.
308  if (x == 0.0 && y > -0.5 && y < 0) {
309  Particles.template getLastProp<0>()[0] = 1.0;
310  Particles.template getLastProp<0>()[1] = 0.0;
311  } else if (x == 0.0 && y > 0 && y < 0.5) {
312  Particles.template getLastProp<0>()[0] = 0.0;
313  Particles.template getLastProp<0>()[1] = 1.0;
314  } else {
315  Particles.template getLastProp<0>()[0] = 0.0;
316  Particles.template getLastProp<0>()[1] = 0.0;
317  }
318  ++it;
319  }
320  Particles.map();
321  Particles.ghost_get<0>();
322 
323  //We write the particles to check if the initialization is correct.
324  Particles.write("Init");
326 
338  // Now we initialize the grid with a filled circle. Outside the circle, the value of Phi_0 will be -1, inside +1.
340  //Now we construct the subsets based on the subset number.
341  dist_vector_subset_type Particles_bulk(Particles, 0);
342  dist_vector_subset_type Particles_boundary(Particles, 1);
343 
344  //We cast the global pointers to Particles and Particles_bulk as expected by the RHS functor.
345  PointerDistGlobal = (void *) &Particles;
346  PointerDistSubset = (void *) &Particles_bulk;
347  PointerDistSubset2 = (void *) &Particles_boundary;
349 
360  //We create the DCPSE Based operators for discretization of the operators.
362  Derivative_xx Dxx(Particles, 2, rCut);
363  Derivative_yy Dyy(Particles, 2, rCut);
364  //We create aliases for referring to property and and positions.
365  auto Pos = getV<PROP_POS>(Particles);
366  auto C = getV<0>(Particles);
367  auto V = getV<1>(Particles);
368  auto dC = getV<2>(Particles);
369 
370  //We create aliases for referring to the subset properties.
371  auto C_bulk = getV<0>(Particles_bulk);
372  auto dC_bulk = getV<2>(Particles_bulk);
374 
389  //Now we create a odeint stepper object (RK4). Since we are in 2d, we are going to use "state_type_2d_ofp". Which is a structure or state_type compatible with odeint. We further pass all the parameters including "boost::numeric::odeint::vector_space_algebra_ofp",which tell odeint to use openfpm algebra.
391  // The template parameters are: state_type_2d_ofp (state type of X), double (type of the value inside the state), state_type_2d_ofp (state type of DxDt), double (type of the time), boost::numeric::odeint::vector_space_algebra_ofp (our algebra)
392  boost::numeric::odeint::runge_kutta4<state_type_2d_ofp, double, state_type_2d_ofp, double, boost::numeric::odeint::vector_space_algebra_ofp> Odeint_rk4;
393 
394  //The method Odeint_rk4 from Odeint, requires system (a function which computes RHS of the PDE), an instance of the Compute RHS functor. We create the System with the correct types and parameteres for the operators as declared before.
396 
397  //Since we are using Odeint to control the time steps, we also create a observer instance. Which also updates the position via an euler step for moving thr particles.
398  ObserverFunctor<Derivative_xx, Derivative_yy> ObserveAndUpdate(Dxx, Dyy);
399 
400 
401  //Furhter, odeint needs data in a state type "state_type_2d_ofp", we create one and fill in the initial condition.
403  //Since we created a 2d state_type we initialize the two fields in the object data using the method get.
404  X.data.get<x>() = C[0];
405  X.data.get<y>() = C[1];
407 
427  std::vector<double> inter_times; // vector to store intermediate time steps taken by odeint.
430 
431  size_t steps = boost::numeric::odeint::integrate_const(Odeint_rk4, System, X, 0.0, tf, dt, ObserveAndUpdate);
432  /* From Odeint: (Please refer to Odeint on more details about these steppers.)
433  runge_kutta_dopri5 is maybe the best default stepper. It has step size control as well as dense-output functionality. Simple create a dense-output stepper by make_dense_output( 1.0e-6 , 1.0e-5 , runge_kutta_dopri5< state_type >() ).
434  runge_kutta4 is a good stepper for constant step sizes. It is widely used and very well known. If you need to create artificial time series this stepper should be the first choice.
435  'runge_kutta_fehlberg78' is similar to the 'runge_kutta4' with the advantage that it has higher precision. It can also be used with step size control.
436  adams_bashforth_moulton is very well suited for ODEs where the r.h.s. is expensive (in terms of computation time). It will calculate the system function only once during each step.*/
437  /*Below are listed examples of using different steppers and their usage.*/
438  //size_t steps = integrate_adaptive(controlled_stepper_type() , System , X , t , tf , dt, ObserveAndUpdate);
439  //size_t steps = boost::numeric::odeint::integrate_adaptive( boost::numeric::odeint::make_controlled( 1.0e-7 , 1.0e-7 ,boost::numeric::odeint::runge_kutta_cash_karp54< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp>() ) , System , X , t , tf , dt, ObserveAndUpdate );
440  //size_t steps = boost::numeric::odeint::integrate_const( boost::numeric::odeint::make_dense_output(1.0e-7, 1.0e-7, boost::numeric::odeint::runge_kutta_dopri5< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp >() ) , System , X , t , tf , dt, ObserveAndUpdate );
441  //size_t steps = boost::numeric::odeint::integrate_adaptive( boost::numeric::odeint::make_dense_output( 1.0e-8 , 1.0e-8 , boost::numeric::odeint::runge_kutta_dopri5< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp >() ) , System , X , t , tf , dt, ObserveAndUpdate );
442  //size_t steps = boost::numeric::odeint::integrate_adaptive( boost::numeric::odeint::make_controlled( 1.0e-7 , 1.0e-7 , boost::numeric::odeint::runge_kutta_dopri5< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp >() ) , System , X , t , tf , dt, ObserveAndUpdate );
443 
444 
445  std::cout << "No. of Time steps taken: " << steps << std::endl;
446  //Copying back the final solution and outputting the number of steps taken by Odeint.
447  C_bulk[x] = X.data.get<0>();
448  C_bulk[y] = X.data.get<1>();
449  //Deallocating the operators
450  Dxx.deallocate(Particles);
451  Dyy.deallocate(Particles);
452  openfpm_finalize(); // Finalize openFPM library
453  return 0;
454 } //main end
456 
System of equations.
Definition: System.hpp:23
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:27
Definition: Ghost.hpp:39
void deleteGhost()
Delete the particles on the ghost.
void map(size_t opt=NONE)
It move all the particles that does not belong to the local processor to the respective processor.
A 2d Odeint and Openfpm compatible structure.
void ghost_get(size_t opt=WITH_POSITION)
It synchronize the properties and position of the ghost particles.
This class represent an N-dimensional box.
Definition: Box.hpp:60
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:214
void update()
Update the subset indexes.