main.cpp

//
// Created by Abhinav Singh @absingh on 26/04/2021
//
// Include Vector Expression,Vector Expressions for Subset,DCPSE,Odeint header files
#include "Operators/Vector/vector_dist_operators.hpp"
#include "Vector/vector_dist_subset.hpp"
#include "DCPSE/DCPSE_op/DCPSE_op.hpp"
#include "OdeIntegrators/OdeIntegrators.hpp"



constexpr int x = 0;
constexpr int y = 1;

double dt;

void *PointerDistGlobal, *PointerDistSubset;

typedef aggregate<VectorS<2, double>, VectorS<2, double>, VectorS<2, double>> Property_type;
typedef vector_dist_ws<2, double, Property_type> dist_vector_type;
typedef vector_dist_subset<2, double, Property_type> dist_vector_subset_type;

template<typename DXX,typename DYY>
struct RHSFunctor
{
    //Intializing the operators
    DXX &Dxx;
    DYY &Dyy;
    //Constructor
    RHSFunctor(DXX &Dxx,DYY &Dyy):Dxx(Dxx),Dyy(Dyy)
    {}

    void operator()( const state_type_2d_ofp &X , state_type_2d_ofp &dxdt , const double t ) const
    {
        //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;

        //Aliasing the properties.
        auto C = getV<0>(Particles);
        auto C_bulk = getV<0>(Particles_bulk);
        auto dC = getV<2>(Particles);
        auto dC_bulk = getV<2>(Particles_bulk);

        //Since the state vector does not have enough information to compute derivates, we copy them back to the Particle Distribution.
        //In this way we can compute all the things that are required to be computed in each stage of the ode integrator.
        // 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).

        //These expressions only update the bulk values of C.
        C_bulk[x]=X.data.get<0>();
        C_bulk[y]=X.data.get<1>();
        Particles.ghost_get<0>();

        // We do the RHS computations for the Laplacian (Updating bulk only).
        dC_bulk[x] = 0.1*(Dxx(C[x])+Dyy(C[x]));
        dC_bulk[y] = 0.1*(Dxx(C[y])+Dyy(C[y]));

        //We copy back to the dxdt state_type for Odeint
        dxdt.data.get<0>()=dC[x];
        dxdt.data.get<1>()=dC[y];
    }
};



int main(int argc, char* argv[]) {
    //  initialize library
    openfpm_init(&argc, &argv);


    const size_t sz[2] = {41, 41};
    Box<2, double> box({-1, -1}, {1.0, 1.0});
    size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
    double spacing[2];
    spacing[0] = 2.0 / (sz[0] - 1);
    spacing[1] = 2.0 / (sz[1] - 1);
    double rCut = 3.1 * spacing[0];
    Ghost<2, double> ghost(rCut);

    dist_vector_type Particles(0, box, bc, ghost);
    Particles.setPropNames({"Concentration", "Velocity", "TempConcentration"});

    auto it = Particles.getGridIterator(sz);
    while (it.isNext()) {
        Particles.add();
        auto key = it.get();
        double x = -1.0 + key.get(0) * spacing[0];
        Particles.getLastPos()[0] = x;
        double y = -1.0 + key.get(1) * spacing[1];
        Particles.getLastPos()[1] = y;
        //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.
        if (x != -1.0 && x != 1.0 && y != -1.0 && y != 1) {
            Particles.getLastSubset(0);
        } else {
            Particles.getLastSubset(1);
        }
        // Here fill the Initial value of the concentration.
        if (x == 0.0 && y > -0.5 && y < 0) {
            Particles.template getLastProp<0>()[0] = 1.0;
            Particles.template getLastProp<0>()[1] = 0.0;
        } else if (x == 0.0 && y > 0 && y < 0.5) {
            Particles.template getLastProp<0>()[0] = 0.0;
            Particles.template getLastProp<0>()[1] = 1.0;
        } else {
            Particles.template getLastProp<0>()[0] = 0.0;
            Particles.template getLastProp<0>()[1] = 0.0;
        }
        ++it;
    }
    Particles.map();
    Particles.ghost_get<0>();

    //We write the particles to check if the initialization is correct.
    Particles.write("Init");

    // 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;

    //We create the DCPSE Based operators for discretization of the operators.
    Derivative_xx Dxx(Particles, 2, rCut);
    Derivative_yy Dyy(Particles, 2, rCut);
    //We create aliases for referring to property and and positions.
    auto Pos = getV<PROP_POS>(Particles);
    auto C = getV<0>(Particles);
    auto V = getV<1>(Particles);
    auto dC = getV<2>(Particles);

    //We create aliases for referring to the subset properties.
    auto C_bulk = getV<0>(Particles_bulk);
    auto V_bulk = getV<1>(Particles_bulk);
    auto dC_bulk = getV<2>(Particles_bulk);


    //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.
    // 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)
    boost::numeric::odeint::runge_kutta4<state_type_2d_ofp, double, state_type_2d_ofp, double, boost::numeric::odeint::vector_space_algebra_ofp> Odeint_rk4;

    //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.
    RHSFunctor<Derivative_xx, Derivative_yy> System(Dxx, Dyy);

    //Furhter, odeint needs data in a state type "state_type_2d_ofp", we create one and fill in the initial condition.
    state_type_2d_ofp X;
    //Since we created a 2d state_type we initialize the two fields in the object data using the method get.
    X.data.get<x>() = C[x];
    X.data.get<y>() = C[y];

    // 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 = 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.
        C_bulk[x] = X.data.get<0>();
        C_bulk[y] = X.data.get<1>();
        //We advect the particles with the velocity using an Euler step.
        Pos = Pos + dt * V;
        //We need to map and get the ghost, as we moved the particles.
        Particles.map();
        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

    //Deallocating the operators
    Dxx.deallocate(Particles);
    Dyy.deallocate(Particles);
    openfpm_finalize(); // Finalize openFPM library
    return 0;
} //main end