diff --git a/BasicSolver/2D-mpi-v3/dcavity.par b/BasicSolver/2D-mpi-v3/dcavity.par index 3ada9c1..4b56dc3 100644 --- a/BasicSolver/2D-mpi-v3/dcavity.par +++ b/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 #=============================================================================== diff --git a/BasicSolver/2D-mpi-v3/src/main.c b/BasicSolver/2D-mpi-v3/src/main.c index b7f495a..1eef68d 100644 --- a/BasicSolver/2D-mpi-v3/src/main.c +++ b/BasicSolver/2D-mpi-v3/src/main.c @@ -52,10 +52,10 @@ int main (int argc, char** argv) initParticleTracer(&particletracer, ¶ms, &solver); if (rank == 0) { printParameter(¶ms); - 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) { diff --git a/BasicSolver/2D-mpi-v3/src/particletracing.c b/BasicSolver/2D-mpi-v3/src/particletracing.c index c8c67fc..4a0493a 100644 --- a/BasicSolver/2D-mpi-v3/src/particletracing.c +++ b/BasicSolver/2D-mpi-v3/src/particletracing.c @@ -12,8 +12,8 @@ #include "vtkWriter.h" -#define U(i, j) u[(j) * (imax + 2) + (i)] -#define V(i, j) v[(j) * (imax + 2) + (i)] +#define U(i, j) u[(j) * (imaxLocal + 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 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); + //printf("\t U N : %.2f, V N : %.2f\n\n", u_n, v_n); 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 y = particlePool[i].y; + double x = particletracer->particlePool[i].x; + 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("\txOffset : %.2f, yOffset : %.2f\n", particletracer->xOffset, particletracer->yOffset); - // printf("\txLocal : %.2f, yLocal : %.2f\n", particletracer->xLocal, particletracer->yLocal); + 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("\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); - printParticles(particletracer); + writeParticles(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)); } \ No newline at end of file diff --git a/BasicSolver/2D-mpi-v3/velocity.png b/BasicSolver/2D-mpi-v3/velocity.png index 2685be6..674b7a2 100644 Binary files a/BasicSolver/2D-mpi-v3/velocity.png and b/BasicSolver/2D-mpi-v3/velocity.png differ diff --git a/BasicSolver/2D-mpi-v3/vis_files/animate.plot b/BasicSolver/2D-mpi-v3/vis_files/animate.plot new file mode 100644 index 0000000..aa80137 --- /dev/null +++ b/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 diff --git a/BasicSolver/2D-mpi-v3/vis_files/trace.gif b/BasicSolver/2D-mpi-v3/vis_files/trace.gif new file mode 100644 index 0000000..2a8d407 Binary files /dev/null and b/BasicSolver/2D-mpi-v3/vis_files/trace.gif differ diff --git a/BasicSolver/2D-seq/src/particletracing.c b/BasicSolver/2D-seq/src/particletracing.c index 0ca1023..b1dd11e 100644 --- a/BasicSolver/2D-seq/src/particletracing.c +++ b/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)