Working 2D PT MPI v3 w/o vtk writer(normal writer)

BasicSolver/2D-mpi-v3/dcavity.par

@@ -32,7 +32,7 @@ jmax       100		# number of interior cells in y-direction
 # Time Data:
 # ---------
 
-te      5.0		# final time
+te      10.0		# final time
 dt     0.02	    # time stepsize
 tau     0.5		# safety factor for time stepsize control (<0 constant delt)
 
@@ -49,12 +49,12 @@ gamma    0.9		# upwind differencing factor gamma
 # -----------------------
 
 numberOfParticles   30
-startTime           0
+startTime           2
 injectTimePeriod    2.0
-writeTimePeriod     0.5
+writeTimePeriod     0.1
 
 x1                  0.1
-y1                  0.9
+y1                  0.5
 x2                  0.9
-y2                  0.9
+y2                  0.5
 #===============================================================================

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

@@ -52,10 +52,10 @@ int main (int argc, char** argv)
     initParticleTracer(&particletracer, &params, &solver);
     if (rank == 0) {
         printParameter(&params);
-        printParticleTracerParameters(&particletracer);
+        
     }
     
-
+    printParticleTracerParameters(&particletracer);
     initProgress(solver.te);
 
     double tau = solver.tau;
@@ -66,7 +66,6 @@ int main (int argc, char** argv)
 
     
     S = getTimeStamp();
-    
 
     while (t <= te) 
     {

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

@@ -12,8 +12,8 @@
 
 #include "vtkWriter.h"
 
-#define U(i, j)  u[(j) * (imax + 2) + (i)]
+#define U(i, j)   u[(j) * (imaxLocal + 2) + (i)]
-#define V(i, j)  v[(j) * (imax + 2) + (i)]
+#define V(i, j)   v[(j) * (imaxLocal + 2) + (i)]
 
 static int ts = 0;
 
@@ -49,6 +49,38 @@ static double sumY(double* sizes, int size, int offset, int coord, int init)
     return sum;
 }
 
+void printUV(ParticleTracer* particletracer, double* u, double* v)
+{
+    int imaxLocal = particletracer->imaxLocal;
+
+    for (int i = 0; i < particletracer->size; i++) {
+        if (i == particletracer->rank) {
+            printf(
+                "\n### RANK %d #######################################################\n",
+                particletracer->rank);
+                printf("\nGrid U : \n");
+            for (int j = 0; j < particletracer->jmaxLocal + 2; j++) {
+                printf("%02d: ", j);
+                for (int i = 0; i < particletracer->imaxLocal + 2; i++) {
+                    printf("%4.2f  ", u[j * (imaxLocal + 2) + i]);
+                }
+                printf("\n");
+            }
+            fflush(stdout);
+            printf("\nGrid V : \n");
+            for (int j = 0; j < particletracer->jmaxLocal + 2; j++) {
+                printf("%02d: ", j);
+                for (int i = 0; i < particletracer->imaxLocal + 2; i++) {
+                    printf("%4.2f  ", v[j * (imaxLocal + 2) + i]);
+                }
+                printf("\n");
+            }
+            fflush(stdout);
+        }
+        MPI_Barrier(MPI_COMM_WORLD);
+    }
+}
+
 void printParticles(ParticleTracer* particletracer)
 {
     for(int i = 0; i < particletracer->totalParticles; ++i)
@@ -71,8 +103,11 @@ void injectParticles(ParticleTracer* particletracer)
     {
         x = particletracer->linSpaceLine[i].x;
         y = particletracer->linSpaceLine[i].y;
-        if(x >= particletracer->xOffset && y >= particletracer->yOffset && x < particletracer->xOffsetEnd && y < particletracer->yOffsetEnd )
+        if(x >= particletracer->xOffset && y >= particletracer->yOffset && x <= particletracer->xOffsetEnd && y <= particletracer->yOffsetEnd )
         {
+            // printf("\nRank : %d\n", particletracer->rank);
+            // printf("\t%.2f >= %.2f && %.2f >= %.2f && %.2f <= %.2f && %.2f <= %.2f\n",x , particletracer->xOffset ,y , particletracer->yOffset, x , particletracer->xOffsetEnd ,y , particletracer->yOffsetEnd);
+        
             particletracer->particlePool[particletracer->pointer].x = x;
             particletracer->particlePool[particletracer->pointer].y = y;
             particletracer->particlePool[particletracer->pointer].flag = true;
@@ -155,12 +190,11 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
             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);
+            // printf("Rank : %d\n", particletracer->rank);
+            // printf("\tOld X : %.2f, translated X : %.2f, xOffset : %.2f, New X : %.2f, iCoord : %d\n\tOld Y : %.2f, translated X : %.2f, yOffset : %.2f, New Y : %.2f, jCoord : %d\n\n",xTemp, x, particletracer->xOffset, new_x, iCoord, yTemp, y, particletracer->yOffset , new_y, jCoord);
             // printf("\tU(iCoord - 1, jCoord - 1) : %.2f, U(iCoord, jCoord - 1) : %.2f, U(iCoord - 1, jCoord) : %.2f, U(iCoord, jCoord) : %.2f\n", U(iCoord - 1, jCoord - 1), U(iCoord, jCoord - 1), U(iCoord - 1, jCoord), U(iCoord, jCoord));
             // printf("\tV(iCoord - 1, jCoord - 1) : %.2f, V(iCoord, jCoord - 1) : %.2f, V(iCoord - 1, jCoord) : %.2f, V(iCoord, jCoord) : %.2f\n\n", V(iCoord - 1, jCoord - 1), V(iCoord, jCoord - 1), V(iCoord - 1, jCoord), V(iCoord, jCoord));
-            // printf("\t U N : %.2f, V N : %.2f\n\n", u_n, v_n);
+            //printf("\t U N : %.2f, V N : %.2f\n\n", u_n, v_n);
-            // 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 < particletracer->xOffset) || (new_x >= particletracer->xOffsetEnd) ||
                  (new_y < particletracer->yOffset) || (new_y >= particletracer->yOffsetEnd)))
@@ -170,9 +204,9 @@ void advanceParticles(ParticleTracer* particletracer, double* restrict u, double
                 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_x <= particletracer->offset[i + particletracer->size * XOFFSETEND]) && 
                        (new_y >= particletracer->offset[i + particletracer->size * YOFFSET]) && 
-                       (new_y < particletracer->offset[i + particletracer->size * YOFFSETEND]) &&
+                       (new_y <= particletracer->offset[i + particletracer->size * YOFFSETEND]) &&
                        i != particletracer->rank)
                     {
                         buff[i][particleBufIndex[i]].x = new_x;
@@ -228,27 +262,26 @@ void freeParticles(ParticleTracer* particletracer)
 
 void writeParticles(ParticleTracer* particletracer)
 {
-    
+    int collectedBuffIndex[particletracer->size];
+
+    MPI_Gather(&particletracer->totalParticles, 1, MPI_INT, collectedBuffIndex, 1, MPI_INT, 0, particletracer->comm);
+
+    if(particletracer->rank != 0)
+    {
+        Particle buff[particletracer->totalParticles];
+        for(int i = 0; i < particletracer->totalParticles; ++i)
+        {   
+            buff[i].x = particletracer->particlePool[i].x;
+            buff[i].y = particletracer->particlePool[i].y;
+            buff[i].flag =  particletracer->particlePool[i].flag;
+            //printf("Rank : %d sending to rank 0 X : %.2f, Y : %.2f with totalpt : %d\n", particletracer->rank, buff[i].x, buff[i].y, particletracer->totalParticles);
+        }
+        MPI_Send(buff, particletracer->totalParticles, particletracer->mpi_particle, 0, 1, particletracer->comm);
+    }
     if(particletracer->rank == 0)
     {
-        int overallTotalParticles;
-        MPI_Reduce(&particletracer->totalParticles, &overallTotalParticles, 1, MPI_INT, MPI_SUM, 0, particletracer->comm);
-
-        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");
@@ -258,10 +291,25 @@ void writeParticles(ParticleTracer* particletracer)
             exit(EXIT_FAILURE);
         }
 
+        for (int i = 1; i < particletracer->size; ++i)
+        {
+            Particle recvBuff[collectedBuffIndex[i]];
+            MPI_Recv(&recvBuff, collectedBuffIndex[i], particletracer->mpi_particle, i, 1, particletracer->comm, MPI_STATUS_IGNORE);
+
+            for (int j = 0; j < collectedBuffIndex[i]; ++j)
+            {
+                double x = recvBuff[j].x;
+                double y = recvBuff[j].y;
+                fprintf(fp, "%f %f\n", x, y);
+                //printf("Rank : 0 receiving from rank %d X : %.2f, Y : %.2f with totalpt : %d\n", i, x, y, particletracer->totalParticles);
+
+            }
+        
+        }
         for (int i = 0; i < particletracer->totalParticles; ++i)
         {
-            double x = particlePool[i].x;
+            double x = particletracer->particlePool[i].x;
-            double y = particlePool[i].y;
+            double y = particletracer->particlePool[i].y;
             fprintf(fp, "%f %f\n", x, y);
         }
         fclose(fp);
@@ -355,6 +403,8 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params, Solve
     particletracer->xOffsetEnd  = particletracer->xOffset + particletracer->xLocal;
     particletracer->yOffsetEnd  = particletracer->yOffset + particletracer->yLocal;
 
+    printf("Rank : %d, xOffset : %.2f, yOffset : %.2f, xOffsetEnd : %.2f, yOffsetEnd : %.2f\n", particletracer->rank, particletracer->xOffset, particletracer->yOffset, particletracer->xOffsetEnd, particletracer->yOffsetEnd);
+
     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);
@@ -383,8 +433,14 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params, Solve
         particletracer->linSpaceLine[i].x    = spacing * particletracer->x1 + (1.0 - spacing) * particletracer->x2;
         particletracer->linSpaceLine[i].y    = spacing * particletracer->y1 + (1.0 - spacing) * particletracer->y2;
         particletracer->linSpaceLine[i].flag = true;
+
+        //if(particletracer->rank == 1) printf("\nRank : %d, with linspace X : %.2f and linspace Y : %.2f\n", particletracer->rank, particletracer->linSpaceLine[i].x , particletracer->linSpaceLine[i].y);
     }
 
+    
+
+
+
     // Create the mpi_particle datatype
     MPI_Datatype mpi_particle;
     int lengths[3] = { 1, 1, 1 };
@@ -415,16 +471,16 @@ 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("\txOffsetEnd : %.2f, yOffsetEnd : %.2f\n", particletracer->xOffsetEnd, particletracer->yOffsetEnd);
+    printf("\txLocal : %.2f, yLocal : %.2f\n", particletracer->xLocal, particletracer->yLocal);
 }   
 
-void trace(ParticleTracer* particletracer, double* u, double* v, double time)
+void trace(ParticleTracer* particletracer, double* restrict u, double* restrict v, double time)
 {
     if (time >= particletracer->startTime)
     {
-        //printParticles(particletracer);
         if( (time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod)
         {
             injectParticles(particletracer);
@@ -432,10 +488,8 @@ void trace(ParticleTracer* particletracer, double* u, double* v, double time)
         }
         if( (time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod)
         {
-            //writeParticles(particletracer);
+            writeParticles(particletracer);
-            printParticles(particletracer);
             particletracer->lastWriteTime = time;
-
         }
         advanceParticles(particletracer, u, v, time);
         compress(particletracer);
@@ -468,7 +522,7 @@ void compress(ParticleTracer* particletracer)
     }
 
     particletracer->totalParticles = totalParticles;
-    particletracer->pointer        = totalParticles + 1;
+    particletracer->pointer        = totalParticles;
 
-    memcpy(particletracer->particlePool, tempPool, sizeof(tempPool));
+    memcpy(particletracer->particlePool, tempPool, totalParticles * sizeof(Particle));
 }

BasicSolver/2D-mpi-v3/vis_files/animate.plot

@@ -0,0 +1,9 @@
+unset border; unset tics; unset key;
+set term gif animate delay 50
+set output "trace.gif"
+set xrange [0:30]
+set yrange [0:4]
+do for [ts=0:200] {
+    plot "particles_".ts.".dat" with points pointtype 7
+}
+unset output

BasicSolver/2D-seq/src/particletracing.c

@@ -205,6 +205,7 @@ void trace(ParticleTracer* particletracer, double* u, double* v, double time)
         if( (time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod)
         {
             injectParticles(particletracer);
+            printf("Rank : %d particles after particles injected !", particletracer->rank);
             particletracer->lastInjectTime = time;
         }
         if( (time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod)