Merge branch 'FD_solver' into 'develop'

Fixing odeint Examples

See merge request openfpm/openfpm_pdata!9

example/Numerics/OdeInt/main.cpp

@@ -272,6 +272,8 @@ int main(int argc, char* argv[]) {
     // Now we initialize the grid with a filled circle. Outside the circle, the value of Phi_0 will be -1, inside +1.
     //Now we construct the subsets based on the subset number.
     dist_vector_subset_type Particles_bulk(Particles, 0);
+    dist_vector_subset_type Particles_boundary(Particles, 1);
+
     //We cast the global pointers to Particles and Particles_bulk as expected by the RHS functor.
     PointerDistGlobal = (void *) &Particles;
     PointerDistSubset = (void *) &Particles_bulk;
@@ -360,13 +362,15 @@ int main(int argc, char* argv[]) {
     //! @cond [OdeintT] @endcond
     // Now we will create a time loop for controlling the step size ourselves but call odeint to do the 4 stages of RK4.
     int ctr = 0;
-    double t = 0, tf = 1e-1, dt = 1e-2;
+    double t = 0, tf = 1, dt = 1e-2;
     while (t < tf) {
         //computing the velocity at the current position and at time step t (only in the bulk, so boundary remains 0).
         V_bulk[x] = -Pos[y] * exp(-10.0 * (Pos[x] * Pos[x] + Pos[y] * Pos[y]));
         V_bulk[y] = Pos[x] * exp(-10.0 * (Pos[x] * Pos[x] + Pos[y] * Pos[y]));
         //Observing the state
+        Particles.deleteGhost();
         Particles.write_frame("PDE_Sol", ctr);
+        Particles.ghost_get<0>();
         //calling rk4 with the function System, the state_type, current time t and the stepsize dt. It computes one step of the RK4.
         Odeint_rk4.do_step(System, X, t, dt);
         //Copying back the step, updating only the bulk values.
@@ -379,8 +383,12 @@ int main(int argc, char* argv[]) {
         Particles.ghost_get<0>();
         //We update the subset and operators as the particles moved.
         Particles_bulk.update();
+        Particles_boundary.update();
         Dxx.update(Particles);
         Dyy.update(Particles);
+        //Reinitialzing as distributed size can change.
+        X.data.get<0>()=C[x];
+        X.data.get<1>()=C[y];
         ctr++;
         t += dt;
     } //time loop end
@@ -398,9 +406,4 @@ int main(int argc, char* argv[]) {
  * ## Full code ## {#odeint_c_full}
  *
  * @include example/Numerics/OdeInt/main.cpp
- */
-
-
-
-
-
+ */
\ No newline at end of file

example/Numerics/OdeInt/main2.cpp

@@ -74,9 +74,9 @@
 constexpr int x = 0;
 constexpr int y = 1;
 
-double dt;
+double dt=1e-2,tf=1.0;
 
-void *PointerDistGlobal, *PointerDistSubset;
+void *PointerDistGlobal, *PointerDistSubset,*PointerDistSubset2;
 
 typedef aggregate<VectorS<2, double>, VectorS<2, double>, VectorS<2, double>> Property_type;
 typedef vector_dist_ws<2, double, Property_type> dist_vector_type;
@@ -203,6 +203,7 @@ struct ObserverFunctor {
         //Casting the pointers to OpenFPM vector distributions
         dist_vector_type &Particles = *(dist_vector_type *) PointerDistGlobal;
         dist_vector_subset_type &Particles_bulk = *(dist_vector_subset_type *) PointerDistSubset;
+        dist_vector_subset_type &Particles_boundary = *(dist_vector_subset_type *) PointerDistSubset2;
 
         //Aliasing the position and properties.
         auto Pos = getV<PROP_POS>(Particles);
@@ -225,6 +226,7 @@ struct ObserverFunctor {
             Particles.ghost_get<0>();
             //Updating the subset and operators based on new positions
             Particles_bulk.update();
+            Particles_boundary.update();
             Dxx.update(Particles);
             Dyy.update(Particles);
 
@@ -232,9 +234,11 @@ struct ObserverFunctor {
             X.data.get<0>() = Concentration[x];
             X.data.get<1>() = Concentration[y];
             //counting the step number
-            ctr++;
-        }
 
+        }
+        ctr++;
+        std::cout<<"Taking a step at t="<<t<<" with dt="<<t-t_old<<std::endl;
+        t_old=t;
         //Writing the data after every step
         Particles.deleteGhost();
         Particles.write_frame("PDE_sol", ctr);
@@ -335,9 +339,12 @@ int main(int argc, char *argv[])
     // Now we initialize the grid with a filled circle. Outside the circle, the value of Phi_0 will be -1, inside +1.
     //Now we construct the subsets based on the subset number.
     dist_vector_subset_type Particles_bulk(Particles, 0);
+    dist_vector_subset_type Particles_boundary(Particles, 1);
+
     //We cast the global pointers to Particles and Particles_bulk as expected by the RHS functor.
     PointerDistGlobal = (void *) &Particles;
     PointerDistSubset = (void *) &Particles_bulk;
+    PointerDistSubset2 = (void *) &Particles_boundary;
     //! @cond [Pointer2Init] @endcond
 
     /**
@@ -419,7 +426,6 @@ int main(int argc, char *argv[])
     */
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     //! @cond [OdeintTCall] @endcond
-    double t = 0, tf = 1e-1, dt = 1e-2;
     std::vector<double> inter_times; // vector to store intermediate time steps taken by odeint.
 
     size_t steps = boost::numeric::odeint::integrate_const(Odeint_rk4, System, X, 0.0, tf, dt, ObserveAndUpdate);