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" 77 double dt=1e-2,tf=1.0;
79 void *PointerDistGlobal, *PointerDistSubset,*PointerDistSubset2;
114 template<
typename DXX,
typename DYY>
122 RHSFunctor(DXX &Dxx,DYY &Dyy):Dxx(Dxx),Dyy(Dyy)
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);
142 C_bulk[x]=X.data.get<0>();
143 C_bulk[y]=X.data.get<1>();
147 dC_bulk[x] = 0.1*(Dxx(C[x])+Dyy(C[x]));
148 dC_bulk[y] = 0.1*(Dxx(C[y])+Dyy(C[y]));
151 dxdt.data.get<0>()=dC[x];
152 dxdt.data.get<1>()=dC[y];
185 template<
typename DXX,
typename DYY>
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);
218 Concentration_bulk[x] = X.data.get<0>();
219 Concentration_bulk[y] = X.data.get<1>();
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;
229 Particles_boundary.
update();
230 Dxx.update(Particles);
231 Dyy.update(Particles);
234 X.data.get<0>() = Concentration[x];
235 X.data.get<1>() = Concentration[y];
240 std::cout<<
"Taking a step at t="<<t<<
" with dt="<<t-t_old<<std::endl;
244 Particles.write_frame(
"PDE_sol", ctr);
261 int main(
int argc,
char *argv[])
265 openfpm_init(&argc, &argv);
280 const size_t sz[2] = {41, 41};
283 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
285 spacing[0] = 2.0 / (sz[0] - 1);
286 spacing[1] = 2.0 / (sz[1] - 1);
287 double rCut = 3.1 * spacing[0];
291 Particles.setPropNames({
"Concentration",
"Velocity",
"TempConcentration"});
293 auto it = Particles.getGridIterator(sz);
294 while (it.isNext()) {
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;
302 if (x != -1.0 && x != 1.0 && y != -1.0 && y != 1) {
303 Particles.getLastSubset(0);
305 Particles.getLastSubset(1);
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;
315 Particles.template getLastProp<0>()[0] = 0.0;
316 Particles.template getLastProp<0>()[1] = 0.0;
321 Particles.ghost_get<0>();
324 Particles.write(
"Init");
345 PointerDistGlobal = (
void *) &Particles;
346 PointerDistSubset = (
void *) &Particles_bulk;
347 PointerDistSubset2 = (
void *) &Particles_boundary;
362 Derivative_xx Dxx(Particles, 2, rCut);
363 Derivative_yy Dyy(Particles, 2, rCut);
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);
371 auto C_bulk = getV<0>(Particles_bulk);
372 auto dC_bulk = getV<2>(Particles_bulk);
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;
404 X.data.get<x>() = C[0];
405 X.data.get<y>() = C[1];
427 std::vector<double> inter_times;
431 size_t steps = boost::numeric::odeint::integrate_const(Odeint_rk4,
System, X, 0.0, tf, dt, ObserveAndUpdate);
445 std::cout <<
"No. of Time steps taken: " << steps << std::endl;
447 C_bulk[x] = X.data.get<0>();
448 C_bulk[y] = X.data.get<1>();
450 Dxx.deallocate(Particles);
451 Dyy.deallocate(Particles);
This class implement the point shape in an N-dimensional space.
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.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
void update()
Update the subset indexes.