main.cpp

#include "Grid/grid_dist_id.hpp"
#include "data_type/aggregate.hpp"
#include "timer.hpp"


constexpr int U = 0;
constexpr int V = 1;

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



void init(grid_dist_id<2,double,aggregate<double,double> > & Old, grid_dist_id<2,double,aggregate<double,double> > & New, Box<2,double> & domain)
{



    auto it = Old.getDomainIterator();

    while (it.isNext())
    {
        // Get the local grid key
        auto key = it.get();

        // Old values U and V
        Old.template get<U>(key) = 1.0;
        Old.template get<V>(key) = 0.0;

        // Old values U and V
        New.template get<U>(key) = 0.0;
        New.template get<V>(key) = 0.0;

        ++it;
    }



    grid_key_dx<2> start({(long int)std::floor(Old.size(0)*1.55f/domain.getHigh(0)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(1))});
    grid_key_dx<2> stop ({(long int)std::ceil (Old.size(0)*1.85f/domain.getHigh(0)),(long int)std::ceil (Old.size(1)*1.85f/domain.getHigh(1))});
    auto it_init = Old.getSubDomainIterator(start,stop);

    while (it_init.isNext())
    {
        auto key = it_init.get();

        Old.template get<U>(key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/100.0;
        Old.template get<V>(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/200.0;

        ++it_init;
    }



}



int main(int argc, char* argv[])
{

    openfpm_init(&argc,&argv);

    // domain
    Box<2,double> domain({0.0,0.0},{2.5,2.5});
    
    // grid size
    size_t sz[2] = {128,128};

    // Define periodicity of the grid
    periodicity<2> bc = {PERIODIC,PERIODIC};
    
    // Ghost in grid unit
    Ghost<2,long int> g(1);
    
    // deltaT
    double deltaT = 1;

    // Diffusion constant for specie U
    double du = 2*1e-5;

    // Diffusion constant for specie V
    double dv = 1*1e-5;

    // Number of timesteps
    size_t timeSteps = 15000;

    // K and F (Physical constant in the equation)
    double K = 0.055;
    double F = 0.03;



    grid_dist_id<2, double, aggregate<double,double>> Old(sz,domain,g,bc);

    // New grid with the decomposition of the old grid
    grid_dist_id<2, double, aggregate<double,double>> New(Old.getDecomposition(),sz,g);

    
    // spacing of the grid on x and y
    double spacing[2] = {Old.spacing(0),Old.spacing(1)};



    init(Old,New,domain);



    // sync the ghost
    size_t count = 0;
    Old.template ghost_get<U,V>();

    // because we assume that spacing[x] == spacing[y] we use formula 2
    // and we calculate the prefactor of Eq 2
    double uFactor = deltaT * du/(spacing[x]*spacing[x]);
    double vFactor = deltaT * dv/(spacing[x]*spacing[x]);

    for (size_t i = 0; i < timeSteps; ++i)
    {
        auto it = Old.getDomainIterator();

        while (it.isNext())
        {
            auto key = it.get();

            // update based on Eq 2
            New.get<U>(key) = Old.get<U>(key) + uFactor * (
                                        Old.get<U>(key.move(x,1)) +
                                        Old.get<U>(key.move(x,-1)) +
                                        Old.get<U>(key.move(y,1)) +
                                        Old.get<U>(key.move(y,-1)) +
                                        -4.0*Old.get<U>(key)) +
                                        - deltaT * Old.get<U>(key) * Old.get<V>(key) * Old.get<V>(key) +
                                        - deltaT * F * (Old.get<U>(key) - 1.0);

            // update based on Eq 2
            New.get<V>(key) = Old.get<V>(key) + vFactor * (
                                        Old.get<V>(key.move(x,1)) +
                                        Old.get<V>(key.move(x,-1)) +
                                        Old.get<V>(key.move(y,1)) +
                                        Old.get<V>(key.move(y,-1)) -
                                        4*Old.get<V>(key)) +
                                        deltaT * Old.get<U>(key) * Old.get<V>(key) * Old.get<V>(key) +
                                        - deltaT * (F+K) * Old.get<V>(key);

            // Next point in the grid
            ++it;
        }

        // Here we copy New into the old grid in preparation of the new step
        // It would be better to alternate, but using this we can show the usage
        // of the function copy. To note that copy work only on two grid of the same
        // decomposition. If you want to copy also the decomposition, or force to be
        // exactly the same, use Old = New
        Old.copy(New);

        // After copy we synchronize again the ghost part U and V
        Old.ghost_get<U,V>();

        // Every 100 time step we output the configuration for
        // visualization
        if (i % 100 == 0)
        {
            Old.write_frame("output",count);
            count++;
        }
    }
    


    openfpm_finalize();

}