Half working 2D mpi v3, write incomplete

BasicSolver/2D-mpi-v3/src/main.c

@@ -52,20 +52,19 @@ int main (int argc, char** argv)
     initParticleTracer(&particletracer, &params, &solver);
     if (rank == 0) {
         printParameter(&params);
+        printParticleTracerParameters(&particletracer);
     }
-    printParticleTracerParameters(&particletracer);
+    
 
     initProgress(solver.te);
 
-    abort();
-
     double tau = solver.tau;
     double te  = solver.te;
     double t   = 0.0;
 
     int (*solver_generic[])() = {solve, solveRB, solveRBA};
 
-
+    
     S = getTimeStamp();
     
 
@@ -82,7 +81,7 @@ int main (int argc, char** argv)
         solver_generic[variant - 1](&solver);
         adaptUV(&solver);
 
-        //trace(&particletracer, solver.u, solver.v, t);
+        trace(&particletracer, solver.u, solver.v, t);
 
         t += solver.dt;
         
@@ -102,6 +101,8 @@ int main (int argc, char** argv)
         printf("Solution took %.2fs\n", E - S);
     }
     collectResult(&solver);
+    freeParticles(&particletracer);
+
 
     MPI_Finalize();
     return EXIT_SUCCESS;

BasicSolver/2D-mpi-v3/src/particletracing.c

@@ -20,7 +20,12 @@ static int ts = 0;
 #define IDIM 0
 #define JDIM 1
 
-static double sum(double* sizes, int size, int offset, int coord, int init)
+#define XOFFSET 0
+#define YOFFSET 1
+#define XOFFSETEND 2
+#define YOFFSETEND 3
+
+static double sumX(double* sizes, int size, int offset, int coord, int init)
 {
     double sum = 0;
 
@@ -32,13 +37,31 @@ static double sum(double* sizes, int size, int offset, int coord, int init)
     return sum;
 }
 
+static double sumY(double* sizes, int size, int offset, int coord, int init)
+{
+    double sum = 0;
+
+    for (int i = init * offset; i < size && coord > 0; ++i, --coord) 
+    {
+        sum += sizes[i];
+    }
+
+    return sum;
+}
+
 void printParticles(ParticleTracer* particletracer)
 {
     for(int i = 0; i < particletracer->totalParticles; ++i)
-    {
+    {   
-        printf("Particle position X : %.2f, Y : %.2f, flag : %d\n", particletracer->particlePool[i].x,
+        printf("Rank : %d Particle position X : %.2f, Y : %.2f, flag : %d, total pt : %d, pointer : %d, xOffset : %.2f, yOffset : %.2f, xOffsetEnd : %.2f, yOffsetEnd : %.2f\n", 
-        particletracer->particlePool[i].y,
+            particletracer->rank, 
-        particletracer->particlePool[i].flag);
+            particletracer->particlePool[i].x,
+            particletracer->particlePool[i].y,
+            particletracer->particlePool[i].flag, 
+            particletracer->totalParticles, 
+            particletracer->pointer, 
+            particletracer->xOffset, particletracer->yOffset, 
+            particletracer->xOffsetEnd, particletracer->yOffsetEnd);
     }
 }
 void injectParticles(ParticleTracer* particletracer)
@@ -48,7 +71,7 @@ void injectParticles(ParticleTracer* particletracer)
     {
         x = particletracer->linSpaceLine[i].x;
         y = particletracer->linSpaceLine[i].y;
-        if(true)
+        if(x >= particletracer->xOffset && y >= particletracer->yOffset && x < particletracer->xOffsetEnd && y < particletracer->yOffsetEnd )
         {
             particletracer->particlePool[particletracer->pointer].x = x;
             particletracer->particlePool[particletracer->pointer].y = y;
@@ -61,22 +84,42 @@ void injectParticles(ParticleTracer* particletracer)
 
 void advanceParticles(ParticleTracer* particletracer, double* restrict u, double* restrict v, double time)
 {
-    int imax = particletracer->imax;
+    int imax       = particletracer->imax;
-    int jmax = particletracer->jmax;
+    int jmax       = particletracer->jmax;
+    int imaxLocal  = particletracer->imaxLocal;
+    int jmaxLocal  = particletracer->jmaxLocal;
 
-
+    double dx      = particletracer->dx;
-    double dx = particletracer->dx;
+    double dy      = particletracer->dy;
-    double dy = particletracer->dy;
 
     double xlength = particletracer->xlength;
     double ylength = particletracer->ylength;
 
+    Particle buff[particletracer->size][30];
+
+    for(int i = 0; i < particletracer->size; ++i)
+    {   
+        for(int j = 0; j < 30; ++j)
+        { 
+            buff[i][j].x = 0.0;
+            buff[i][j].y = 0.0;
+            buff[i][j].flag = false;
+        }
+    }
+    int particleBufIndex[particletracer->size];
+
+    memset(particleBufIndex, 0, sizeof(particleBufIndex));
+
+
     for(int i = 0; i < particletracer->totalParticles; ++i)
     {
         if(particletracer->particlePool[i].flag == true)
         {
-            double x = particletracer->particlePool[i].x;
+            double xTemp = particletracer->particlePool[i].x;
-            double y = particletracer->particlePool[i].y;
+            double yTemp = particletracer->particlePool[i].y;
+
+            double x   = xTemp - particletracer->xOffset;
+            double y   = yTemp - particletracer->yOffset;
 
             int iCoord = (int)(x / dx) + 1;
             int jCoord = (int)((y + 0.5 * dy) / dy) + 1;
@@ -92,7 +135,7 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
                         (x2 - x) * (y - y1) * U(iCoord - 1, jCoord) +
                         (x - x1) * (y - y1) * U(iCoord, jCoord));
 
-            double new_x = x + particletracer->dt * u_n;
+            double new_x = (x + particletracer->xOffset) + particletracer->dt * u_n;
             particletracer->particlePool[i].x = new_x;
 
             iCoord = (int)((x + 0.5 * dx) / dx) + 1;
@@ -109,7 +152,7 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
                             (x2 - x) * (y - y1) * V(iCoord - 1, jCoord) +
                             (x - x1) * (y - y1) * V(iCoord, jCoord));
 
-            double new_y = y + particletracer->dt * v_n;
+            double new_y = (y + particletracer->yOffset) + particletracer->dt * v_n;
             particletracer->particlePool[i].y = new_y;
 
             // printf("\tOld X : %.2f, New X : %.2f, iCoord : %d\n\tOld Y : %.2f, New Y : %.2f, jCoord : %d\n\n", x, new_x, iCoord, y, new_y, jCoord);
@@ -119,48 +162,110 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
             // printf("\t j-1 * (imax + 2) + i-1 = %d with element from U : %.2f", (jCoord - 1) * (200 + 2) + (iCoord - 1), u[(jCoord - 1) * (imax + 2) + (iCoord - 1)]);
             // printf("\nimax : %d, jmax : %d\n", imax, jmax);
 
-            if (((new_x < 0.0) || (new_x > xlength) || (new_y < 0.0) || (new_y > ylength)))
+            if (((new_x < particletracer->xOffset) || (new_x >= particletracer->xOffsetEnd) ||
+                 (new_y < particletracer->yOffset) || (new_y >= particletracer->yOffsetEnd)))
             {
+                //New logic to transfer particles to neighbouring ranks or discard the particle.
+
+                for(int i = 0; i < particletracer->size; ++i)
+                {
+                    if((new_x >= particletracer->offset[i + particletracer->size * XOFFSET]) && 
+                       (new_x < particletracer->offset[i + particletracer->size * XOFFSETEND]) && 
+                       (new_y >= particletracer->offset[i + particletracer->size * YOFFSET]) && 
+                       (new_y < particletracer->offset[i + particletracer->size * YOFFSETEND]) &&
+                       i != particletracer->rank)
+                    {
+                        buff[i][particleBufIndex[i]].x = new_x;
+                        buff[i][particleBufIndex[i]].y = new_y;
+                        buff[i][particleBufIndex[i]].flag = true;
+                        ++particleBufIndex[i];
+                    }
+                }
                 particletracer->particlePool[i].flag = false;
             }
         }
     }
+
+    for(int i = 0; i < particletracer->size; ++i)
+    {   
+        if(i != particletracer->rank)
+        {
+            MPI_Send(buff[i], 30, particletracer->mpi_particle, i, 0, particletracer->comm);
+        }
+    }
+    for(int i = 0; i < particletracer->size; ++i)
+    {   
+        if(i != particletracer->rank)
+        {
+            MPI_Recv(buff[i], 30, particletracer->mpi_particle, i, 0, particletracer->comm, MPI_STATUS_IGNORE);
+        }
+    }
+    for(int i = 0; i < particletracer->size; ++i)
+    {   
+        if(i != particletracer->rank)
+        {
+            for(int j = 0; j < 30; ++j)
+            {
+                if(buff[i][j].flag == true)
+                {
+                    particletracer->particlePool[particletracer->pointer].x = buff[i][j].x;
+                    particletracer->particlePool[particletracer->pointer].y = buff[i][j].y;
+                    particletracer->particlePool[particletracer->pointer].flag = true;
+                    ++(particletracer->pointer);
+                    ++(particletracer->totalParticles);
+                }
+            }
+        }
+    }
 }
 
 void freeParticles(ParticleTracer* particletracer)
 {
     free(particletracer->particlePool);
     free(particletracer->linSpaceLine);
+    free(particletracer->offset);
 }
 
 void writeParticles(ParticleTracer* particletracer)
 {
-    VtkOptions opts = { .particletracer = particletracer };
+    
-
+    if(particletracer->rank == 0)
-    char filename[50];
-    snprintf(filename, 50, "vtk_files/particles%d.vtk", ts);
-    vtkOpen(&opts, filename, ts);
-    vtkParticle(&opts, "particle");
-    vtkClose(&opts);
-
-    FILE* fp;
-    Particle* particlePool = particletracer->particlePool;
-
-    snprintf(filename, 50, "vis_files/particles_%d.dat", ts);
-    fp = fopen(filename, "w");
-
-    if (fp == NULL) {
-        printf("Error!\n");
-        exit(EXIT_FAILURE);
-    }
-
-    for (int i = 0; i < particletracer->totalParticles; ++i)
     {
-        double x = particlePool[i].x;
+        int overallTotalParticles;
-        double y = particlePool[i].y;
+        MPI_Reduce(&particletracer->totalParticles, &overallTotalParticles, 1, MPI_INT, MPI_SUM, 0, particletracer->comm);
-        fprintf(fp, "%f %f\n", x, y);
+
+        printf("Total collected particles at rank 0 : %d", overallTotalParticles);
+        
+
+        VtkOptions opts = { .particletracer = particletracer, .overallTotalParticles = overallTotalParticles };
+
+        char filename[50];
+        snprintf(filename, 50, "vtk_files/particles%d.vtk", ts);
+        vtkOpen(&opts, filename, ts);
+        vtkParticle(&opts, "particle");
+        vtkClose(&opts);
+        
+
+        
+        FILE* fp;
+        Particle* particlePool = particletracer->particlePool;
+
+        snprintf(filename, 50, "vis_files/particles_%d.dat", ts);
+        fp = fopen(filename, "w");
+
+        if (fp == NULL) {
+            printf("Error!\n");
+            exit(EXIT_FAILURE);
+        }
+
+        for (int i = 0; i < particletracer->totalParticles; ++i)
+        {
+            double x = particlePool[i].x;
+            double y = particlePool[i].y;
+            fprintf(fp, "%f %f\n", x, y);
+        }
+        fclose(fp);
     }
-    fclose(fp);
 
     ++ts;
 }
@@ -202,16 +307,27 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params, Solve
     particletracer->imaxLocal         = solver->imaxLocal;
     particletracer->jmaxLocal         = solver->jmaxLocal;
 
-    particletracer->estimatedNumParticles = (particletracer->imaxLocal * particletracer->jmaxLocal) + (4 * params->numberOfParticles);
+    particletracer->estimatedNumParticles = (particletracer->imaxLocal * particletracer->jmaxLocal);
 
     particletracer->particlePool      = malloc(sizeof(Particle) * particletracer->estimatedNumParticles);
+
+    for(int i = 0; i < particletracer->estimatedNumParticles; ++i)
+    {   
+        particletracer->particlePool[i].x = 0.0;
+        particletracer->particlePool[i].y = 0.0;
+        particletracer->particlePool[i].flag = false;
+    }
+
     particletracer->linSpaceLine      = malloc(sizeof(Particle) * particletracer->numberOfParticles);
 
     /* duplicating communication from solver */
     MPI_Comm_dup(solver->comm, &particletracer->comm);
 
-    particletracer->rank               = solver->rank;
+    particletracer->rank              = solver->rank;
-    particletracer->size               = solver->size;
+    particletracer->size              = solver->size;
+
+    particletracer->offset            = (double *)malloc(sizeof(double) * 4 * particletracer->size);
+
 
     memcpy(particletracer->dims, solver->dims, sizeof(solver->dims));
 
@@ -226,14 +342,40 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params, Solve
     particletracer->yLocal           = particletracer->jmaxLocal * particletracer->dy;
 
 
-    double xLocal[solver->size];
+    double xLocal[particletracer->size];
-    double yLocal[solver->size];
+    double yLocal[particletracer->size];
+
+    double offset[4][particletracer->size];
     
     MPI_Allgather(&particletracer->xLocal, 1, MPI_DOUBLE, xLocal, 1, MPI_DOUBLE, particletracer->comm);
     MPI_Allgather(&particletracer->yLocal, 1, MPI_DOUBLE, yLocal, 1, MPI_DOUBLE, particletracer->comm);
 
-    particletracer->xOffset = sum(xLocal, particletracer->size, particletracer->dims[JDIM], particletracer->coords[IDIM], particletracer->coords[JDIM]);
+    particletracer->xOffset     = sumX(xLocal, particletracer->size, particletracer->dims[JDIM], particletracer->coords[IDIM], particletracer->coords[JDIM]);
-    particletracer->yOffset = sum(yLocal, particletracer->size, particletracer->dims[IDIM], particletracer->coords[JDIM], particletracer->coords[IDIM]);
+    particletracer->yOffset     = sumY(yLocal, particletracer->size, particletracer->dims[JDIM], particletracer->coords[JDIM], particletracer->coords[IDIM]);
+    particletracer->xOffsetEnd  = particletracer->xOffset + particletracer->xLocal;
+    particletracer->yOffsetEnd  = particletracer->yOffset + particletracer->yLocal;
+
+    MPI_Allgather(&particletracer->xOffset, 1, MPI_DOUBLE, offset[0], 1, MPI_DOUBLE, particletracer->comm);
+    MPI_Allgather(&particletracer->yOffset, 1, MPI_DOUBLE, offset[1], 1, MPI_DOUBLE, particletracer->comm);
+    MPI_Allgather(&particletracer->xOffsetEnd, 1, MPI_DOUBLE, offset[2], 1, MPI_DOUBLE, particletracer->comm);
+    MPI_Allgather(&particletracer->yOffsetEnd, 1, MPI_DOUBLE, offset[3], 1, MPI_DOUBLE, particletracer->comm);
+
+    memcpy(particletracer->offset, offset, sizeof(offset));
+
+    
+    // if(particletracer->rank == 0)
+    // {   
+    //     for(int i = 0;i < particletracer->size; ++i)
+    //     {
+    //         printf("Rank : %d and its xOffset : %.2f, yOffset : %.2f, xLocal : %.2f, yLocal : %.2f, xOffsetEnd : %.2f, yOffsetEnd : %.2f\n", 
+    //         i, particletracer->offset[i + particletracer->size * XOFFSET],
+    //         particletracer->offset[i + particletracer->size * YOFFSET],  xLocal[i], yLocal[i], 
+    //         particletracer->offset[i + particletracer->size * XOFFSETEND], 
+    //         particletracer->offset[i + particletracer->size * YOFFSETEND]);
+    //     }
+    // }
+    
+    
 
     for (int i = 0; i < particletracer->numberOfParticles; ++i) 
     {
@@ -243,6 +385,25 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params, Solve
         particletracer->linSpaceLine[i].flag = true;
     }
 
+    // Create the mpi_particle datatype
+    MPI_Datatype mpi_particle;
+    int lengths[3] = { 1, 1, 1 };
+ 
+    MPI_Aint displacements[3];
+    Particle dummy_particle;
+    MPI_Aint base_address;
+    MPI_Get_address(&dummy_particle, &base_address);
+    MPI_Get_address(&dummy_particle.x, &displacements[0]);
+    MPI_Get_address(&dummy_particle.y, &displacements[1]);
+    MPI_Get_address(&dummy_particle.flag, &displacements[2]);
+    displacements[0] = MPI_Aint_diff(displacements[0], base_address);
+    displacements[1] = MPI_Aint_diff(displacements[1], base_address);
+    displacements[2] = MPI_Aint_diff(displacements[2], base_address);
+ 
+    MPI_Datatype types[3] = { MPI_DOUBLE, MPI_DOUBLE, MPI_C_BOOL };
+    MPI_Type_create_struct(3, lengths, displacements, types, &particletracer->mpi_particle);
+    MPI_Type_commit(&particletracer->mpi_particle);
+
 }
 
 void printParticleTracerParameters(ParticleTracer* particletracer)
@@ -254,9 +415,9 @@ void printParticleTracerParameters(ParticleTracer* particletracer)
     printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : %.2f, x2 : %.2f, y2 : %.2f\n", particletracer->x1, particletracer->y1, particletracer->x2, particletracer->y2);
     printf("\tPointer : %d, TotalParticles : %d\n", particletracer->pointer, particletracer->totalParticles);
     printf("\tdt : %.2f, dx : %.2f, dy : %.2f\n", particletracer->dt, particletracer->dx, particletracer->dy);
-    printf("\tcoord[0] : %d, coord[1] : %d\n", particletracer->coords[IDIM], particletracer->coords[JDIM]);
+    // printf("\tcoord[0] : %d, coord[1] : %d\n", particletracer->coords[IDIM], particletracer->coords[JDIM]);
-    printf("\txOffset : %.2f, yOffset : %.2f\n", particletracer->xOffset, particletracer->yOffset);
+    // printf("\txOffset : %.2f, yOffset : %.2f\n", particletracer->xOffset, particletracer->yOffset);
-    printf("\txLocal : %.2f, yLocal : %.2f\n", particletracer->xLocal, particletracer->yLocal);
+    // printf("\txLocal : %.2f, yLocal : %.2f\n", particletracer->xLocal, particletracer->yLocal);
 }   
 
 void trace(ParticleTracer* particletracer, double* u, double* v, double time)
@@ -272,6 +433,7 @@ void trace(ParticleTracer* particletracer, double* u, double* v, double time)
         if( (time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod)
         {
             //writeParticles(particletracer);
+            printParticles(particletracer);
             particletracer->lastWriteTime = time;
 
         }
@@ -285,11 +447,18 @@ void compress(ParticleTracer* particletracer)
 {
     Particle* memPool = particletracer->particlePool;
     Particle tempPool[particletracer->totalParticles];
+
+    for(int i = 0; i < particletracer->totalParticles; ++i)
+    {   
+        tempPool[i].x = 0.0;
+        tempPool[i].y = 0.0;
+        tempPool[i].flag = false;
+    }
     int totalParticles = 0;
 
     for(int i=0; i < particletracer->totalParticles; ++i)
     {
-        if(memPool[i].flag == 1)
+        if(memPool[i].flag == true)
         {
             tempPool[totalParticles].x    = memPool[i].x;
             tempPool[totalParticles].y    = memPool[i].y;
@@ -301,5 +470,5 @@ void compress(ParticleTracer* particletracer)
     particletracer->totalParticles = totalParticles;
     particletracer->pointer        = totalParticles + 1;
 
-    memcpy(particletracer->particlePool, tempPool, totalParticles*sizeof(Particle));
+    memcpy(particletracer->particlePool, tempPool, sizeof(tempPool));
 }

BasicSolver/2D-mpi-v3/src/particletracing.h

@@ -39,11 +39,16 @@ typedef struct {
     double x1, y1, x2, y2;
 
     MPI_Comm comm;
+    MPI_Datatype mpi_particle;
+    
     int rank, size;
     int iNeighbours[NDIMS], jNeighbours[NDIMS];
     int coords[NDIMS], dims[NDIMS];
 
-    double xLocal, yLocal, xOffset, yOffset;
+    double xLocal, yLocal, xOffset, yOffset, xOffsetEnd, yOffsetEnd;
+
+    double* offset;
+
 } ParticleTracer;
 
 void initParticleTracer(ParticleTracer*, Parameter*, Solver* );

BasicSolver/2D-mpi-v3/src/vtkWriter.c

@@ -61,19 +61,17 @@ void vtkParticle(VtkOptions* o, char* name)
 {
     Particle* particlePool = o->particletracer->particlePool;
 
-    int imax = o->solver->imax;
-    int jmax = o->solver->jmax;
-
     if (o->fh == NULL) {
         printf("vtkWriter not initialize! Call vtkOpen first!\n");
         exit(EXIT_FAILURE);
     }
 
-    fprintf(o->fh, "POINTS %d float\n", o->particletracer->totalParticles);
+    fprintf(o->fh, "POINTS %d float\n", o->overallTotalParticles);
 
 
-    for (int i = 0; i < o->particletracer->totalParticles; ++i)
+    for (int i = 0; i < o->particletracer->size; ++i)
     {
+        
         double x = particlePool[i].x;
         double y = particlePool[i].y;
         fprintf(o->fh, "%.2f %.2f 0.0\n", x, y);

BasicSolver/2D-mpi-v3/src/vtkWriter.h

@@ -16,7 +16,7 @@ typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
 typedef struct VtkOptions {
     VtkFormat fmt;
     FILE* fh;
-    Solver* solver;
+    int overallTotalParticles;
     ParticleTracer* particletracer;
 } VtkOptions;