diff --git a/BasicSolver/2D-mpi-v3/src/main.c b/BasicSolver/2D-mpi-v3/src/main.c index 25c235f..b7f495a 100644 --- a/BasicSolver/2D-mpi-v3/src/main.c +++ b/BasicSolver/2D-mpi-v3/src/main.c @@ -52,20 +52,19 @@ int main (int argc, char** argv) initParticleTracer(&particletracer, ¶ms, &solver); if (rank == 0) { printParameter(¶ms); + 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; diff --git a/BasicSolver/2D-mpi-v3/src/particletracing.c b/BasicSolver/2D-mpi-v3/src/particletracing.c index ba8b6e1..c8c67fc 100644 --- a/BasicSolver/2D-mpi-v3/src/particletracing.c +++ b/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, - particletracer->particlePool[i].y, - particletracer->particlePool[i].flag); + { + 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->rank, + 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 jmax = particletracer->jmax; + int imax = particletracer->imax; + int jmax = particletracer->jmax; + int imaxLocal = particletracer->imaxLocal; + int jmaxLocal = particletracer->jmaxLocal; - - double dx = particletracer->dx; - double dy = particletracer->dy; + double dx = particletracer->dx; + 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 y = particletracer->particlePool[i].y; + double xTemp = particletracer->particlePool[i].x; + 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 }; - - 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) + + if(particletracer->rank == 0) { - double x = particlePool[i].x; - double y = particlePool[i].y; - fprintf(fp, "%f %f\n", x, y); + 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"); + + 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->size = solver->size; + particletracer->rank = solver->rank; + 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 yLocal[solver->size]; + double xLocal[particletracer->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->yOffset = sum(yLocal, particletracer->size, particletracer->dims[IDIM], particletracer->coords[JDIM], particletracer->coords[IDIM]); + particletracer->xOffset = sumX(xLocal, particletracer->size, particletracer->dims[JDIM], particletracer->coords[IDIM], particletracer->coords[JDIM]); + 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("\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("\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)); } \ No newline at end of file diff --git a/BasicSolver/2D-mpi-v3/src/particletracing.h b/BasicSolver/2D-mpi-v3/src/particletracing.h index 6a4dbd2..f44dbf3 100644 --- a/BasicSolver/2D-mpi-v3/src/particletracing.h +++ b/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* ); diff --git a/BasicSolver/2D-mpi-v3/src/vtkWriter.c b/BasicSolver/2D-mpi-v3/src/vtkWriter.c index 54942b7..47cf604 100644 --- a/BasicSolver/2D-mpi-v3/src/vtkWriter.c +++ b/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); diff --git a/BasicSolver/2D-mpi-v3/src/vtkWriter.h b/BasicSolver/2D-mpi-v3/src/vtkWriter.h index 22f724b..6132ed3 100644 --- a/BasicSolver/2D-mpi-v3/src/vtkWriter.h +++ b/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;