From e8909dcb0a10eb8e063d64fa5fcadb8a93c5e8dd Mon Sep 17 00:00:00 2001 From: absingh <absingh@mpi-cbg.de> Date: Mon, 3 May 2021 19:32:09 +0200 Subject: [PATCH] Fixing new Examples --- example/Numerics/OdeInt/main.cpp | 17 ++++++++++------- example/Numerics/OdeInt/main2.cpp | 16 +++++++++++----- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/example/Numerics/OdeInt/main.cpp b/example/Numerics/OdeInt/main.cpp index 9e33edae7..39908317d 100644 --- a/example/Numerics/OdeInt/main.cpp +++ b/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 diff --git a/example/Numerics/OdeInt/main2.cpp b/example/Numerics/OdeInt/main2.cpp index 395070eb7..235507819 100644 --- a/example/Numerics/OdeInt/main2.cpp +++ b/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); -- GitLab