Work in Progress 2D PT MPI v3

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

@@ -44,4 +44,18 @@ eps           0.00001   # stopping tolerance for pressure iteration
 rho           0.52  	#
 omg           1.8       # relaxation parameter for SOR iteration
 gamma         0.9       # upwind differencing factor gamma
+
+
+# Particle Tracing Data:
+# -----------------------
+
+numberOfParticles   30
+startTime           10
+injectTimePeriod    2.0
+writeTimePeriod     2.0
+
+x1                  1.0
+y1                  0.0
+x2                  1.0
+y2                  4.0
 #===============================================================================

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

@@ -44,4 +44,17 @@ eps      0.001		# stopping tolerance for pressure iteration
 rho      0.5		#
 omg      1.7		# relaxation parameter for SOR iteration
 gamma    0.9		# upwind differencing factor gamma
+
+# Particle Tracing Data:
+# -----------------------
+
+numberOfParticles   30
+startTime           0
+injectTimePeriod    2.0
+writeTimePeriod     0.5
+
+x1                  0.1
+y1                  0.9
+x2                  0.9
+y2                  0.9
 #===============================================================================

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

@@ -14,6 +14,9 @@
 #include "progress.h"
 #include "solver.h"
 #include "timing.h"
+#include "particletracing.h"
+#include "vtkWriter.h"
+
 #include <mpi.h>
 
 
@@ -27,6 +30,8 @@ int main (int argc, char** argv)
     double S, E;
     Parameter params;
     Solver solver;
+    ParticleTracer particletracer;
+
 
     MPI_Init(&argc, &argv);
     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
@@ -42,17 +47,23 @@ int main (int argc, char** argv)
     {
         variant = atoi(argv[2]);
     }
+ 
+    initSolver(&solver, &params);
+    initParticleTracer(&particletracer, &params, &solver);
     if (rank == 0) {
         printParameter(&params);
     }
-    initSolver(&solver, &params);
+    printParticleTracerParameters(&particletracer);
+
     initProgress(solver.te);
 
+    abort();
+
     double tau = solver.tau;
     double te  = solver.te;
     double t   = 0.0;
 
-    int (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
+    int (*solver_generic[])() = {solve, solveRB, solveRBA};
 
 
     S = getTimeStamp();
@@ -68,8 +79,11 @@ int main (int argc, char** argv)
         setSpecialBoundaryCondition(&solver);
         computeFG(&solver);
         computeRHS(&solver);
-        (*solver_generic[variant - 1])(&solver);
+        solver_generic[variant - 1](&solver);
         adaptUV(&solver);
+
+        //trace(&particletracer, solver.u, solver.v, t);
+
         t += solver.dt;
         
 

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

@@ -82,6 +82,16 @@ void readParameter(Parameter* param, const char* filename)
             PARSE_REAL(p_init);
             PARSE_REAL(rho);
 
+
+            /* Added new particle tracing parameters */
+            PARSE_INT(numberOfParticles);
+            PARSE_REAL(startTime);
+            PARSE_REAL(injectTimePeriod);
+            PARSE_REAL(writeTimePeriod);
+            PARSE_REAL(x1);
+            PARSE_REAL(y1);
+            PARSE_REAL(x2);
+            PARSE_REAL(y2);
         }
     }
 

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

@@ -18,6 +18,11 @@ typedef struct {
     char* name;
     int bcLeft, bcRight, bcBottom, bcTop;
     double u_init, v_init, p_init;
+
+    int numberOfParticles;
+    double startTime, injectTimePeriod, writeTimePeriod;
+
+    double x1, y1, x2, y2;
 } Parameter;
 
 void initParameter(Parameter*);

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

@@ -0,0 +1,305 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "vtkWriter.h"
+
+#define U(i, j)  u[(j) * (imax + 2) + (i)]
+#define V(i, j)  v[(j) * (imax + 2) + (i)]
+
+static int ts = 0;
+
+#define IDIM 0
+#define JDIM 1
+
+static double sum(double* sizes, int size, int offset, int coord, int init)
+{
+    double sum = 0;
+
+    for (int i = init; i < size && coord > 0; i+=offset, --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);
+    }
+}
+void injectParticles(ParticleTracer* particletracer)
+{
+    double x, y;
+    for(int i = 0; i < particletracer->numberOfParticles; ++i)
+    {
+        x = particletracer->linSpaceLine[i].x;
+        y = particletracer->linSpaceLine[i].y;
+        if(true)
+        {
+            particletracer->particlePool[particletracer->pointer].x = x;
+            particletracer->particlePool[particletracer->pointer].y = y;
+            particletracer->particlePool[particletracer->pointer].flag = true;
+            ++(particletracer->pointer);
+            ++(particletracer->totalParticles);
+        }
+    }
+}
+
+void advanceParticles(ParticleTracer* particletracer, double* restrict u, double* restrict v, double time)
+{
+    int imax = particletracer->imax;
+    int jmax = particletracer->jmax;
+
+
+    double dx = particletracer->dx;
+    double dy = particletracer->dy;
+
+    double xlength = particletracer->xlength;
+    double ylength = particletracer->ylength;
+
+    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;
+
+            int iCoord = (int)(x / dx) + 1;
+            int jCoord = (int)((y + 0.5 * dy) / dy) + 1;
+
+            double x1 = (double)(iCoord - 1) * dx;
+            double y1 = ((double)(jCoord - 1) - 0.5) * dy;
+            double x2 = (double)iCoord * dx;
+            double y2 = ((double)jCoord - 0.5) * dy;
+
+            double u_n = (1.0 / (dx * dy)) *
+                    ((x2 - x) * (y2 - y) * U(iCoord - 1, jCoord - 1) +
+                        (x - x1) * (y2 - y) * U(iCoord, jCoord - 1) +
+                        (x2 - x) * (y - y1) * U(iCoord - 1, jCoord) +
+                        (x - x1) * (y - y1) * U(iCoord, jCoord));
+
+            double new_x = x + particletracer->dt * u_n;
+            particletracer->particlePool[i].x = new_x;
+
+            iCoord = (int)((x + 0.5 * dx) / dx) + 1;
+            jCoord = (int)(y / dy) + 1;
+
+            x1 = ((double)(iCoord - 1) - 0.5) * dx;
+            y1 = (double)(jCoord - 1) * dy;
+            x2 = ((double)iCoord - 0.5) * dx;
+            y2 = (double)jCoord * dy;
+
+            double v_n = (1.0 / (dx * dy)) *
+                        ((x2 - x) * (y2 - y) * V(iCoord - 1, jCoord - 1) +
+                            (x - x1) * (y2 - y) * V(iCoord, jCoord - 1) +
+                            (x2 - x) * (y - y1) * V(iCoord - 1, jCoord) +
+                            (x - x1) * (y - y1) * V(iCoord, jCoord));
+
+            double new_y = y + 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("\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);
+
+            if (((new_x < 0.0) || (new_x > xlength) || (new_y < 0.0) || (new_y > ylength)))
+            {
+                particletracer->particlePool[i].flag = false;
+            }
+        }
+    }
+}
+
+void freeParticles(ParticleTracer* particletracer)
+{
+    free(particletracer->particlePool);
+    free(particletracer->linSpaceLine);
+}
+
+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)
+    {
+        double x = particlePool[i].x;
+        double y = particlePool[i].y;
+        fprintf(fp, "%f %f\n", x, y);
+    }
+    fclose(fp);
+
+    ++ts;
+}
+
+void initParticleTracer(ParticleTracer* particletracer, Parameter* params, Solver* solver)
+{
+    int dims[NDIMS]    = { 0, 0 };
+    int periods[NDIMS] = { 0, 0 };
+
+    /* initializing local properties from params */
+
+    particletracer->numberOfParticles = params->numberOfParticles;
+    particletracer->startTime         = params->startTime;
+    particletracer->injectTimePeriod  = params->injectTimePeriod;
+    particletracer->writeTimePeriod   = params->writeTimePeriod;
+
+    particletracer->dt                = params->dt;
+    particletracer->dx                = params->xlength / params->imax;
+    particletracer->dy                = params->ylength / params->jmax;
+
+    particletracer->xlength           = params->xlength;
+    particletracer->ylength           = params->ylength;
+    
+    particletracer->x1                = params->x1;
+    particletracer->y1                = params->y1;
+    particletracer->x2                = params->x2;
+    particletracer->y2                = params->y2;
+
+    particletracer->lastInjectTime    = params->startTime;
+    particletracer->lastUpdateTime    = params->startTime;
+    particletracer->lastWriteTime     = params->startTime;
+
+    particletracer->pointer           = 0;
+    particletracer->totalParticles    = 0;
+
+    particletracer->imax              = params->imax;
+    particletracer->jmax              = params->jmax;
+
+    particletracer->imaxLocal         = solver->imaxLocal;
+    particletracer->jmaxLocal         = solver->jmaxLocal;
+
+    particletracer->estimatedNumParticles = (particletracer->imaxLocal * particletracer->jmaxLocal) + (4 * params->numberOfParticles);
+
+    particletracer->particlePool      = malloc(sizeof(Particle) * particletracer->estimatedNumParticles);
+    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;
+
+    memcpy(particletracer->dims, solver->dims, sizeof(solver->dims));
+
+    memcpy(particletracer->coords, solver->coords, sizeof(solver->coords));
+
+    memcpy(particletracer->iNeighbours, solver->iNeighbours, sizeof(solver->iNeighbours));
+
+    memcpy(particletracer->jNeighbours, solver->jNeighbours, sizeof(solver->jNeighbours));
+
+
+    particletracer->xLocal           = particletracer->imaxLocal * particletracer->dx;
+    particletracer->yLocal           = particletracer->jmaxLocal * particletracer->dy;
+
+
+    double xLocal[solver->size];
+    double yLocal[solver->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]);
+
+    for (int i = 0; i < particletracer->numberOfParticles; ++i) 
+    {
+        double spacing                       = (double)i / (double)(particletracer->numberOfParticles - 1);
+        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;
+    }
+
+}
+
+void printParticleTracerParameters(ParticleTracer* particletracer)
+{
+    printf("Particle Tracing data:\n");
+    printf("Rank : %d\n", particletracer->rank);
+    printf("\tNumber of particles : %d being injected for every period of %.2f\n", particletracer->numberOfParticles, particletracer->injectTimePeriod);
+    printf("\tstartTime : %.2f\n", particletracer->startTime);
+    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);
+}   
+
+void trace(ParticleTracer* particletracer, double* u, double* v, double time)
+{
+    if (time >= particletracer->startTime)
+    {
+        //printParticles(particletracer);
+        if( (time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod)
+        {
+            injectParticles(particletracer);
+            particletracer->lastInjectTime = time;
+        }
+        if( (time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod)
+        {
+            //writeParticles(particletracer);
+            particletracer->lastWriteTime = time;
+
+        }
+        advanceParticles(particletracer, u, v, time);
+        compress(particletracer);
+        particletracer->lastUpdateTime = time;   
+    }
+}
+
+void compress(ParticleTracer* particletracer)
+{
+    Particle* memPool = particletracer->particlePool;
+    Particle tempPool[particletracer->totalParticles];
+    int totalParticles = 0;
+
+    for(int i=0; i < particletracer->totalParticles; ++i)
+    {
+        if(memPool[i].flag == 1)
+        {
+            tempPool[totalParticles].x    = memPool[i].x;
+            tempPool[totalParticles].y    = memPool[i].y;
+            tempPool[totalParticles].flag = memPool[i].flag;
+            ++totalParticles;
+        }
+    }
+
+    particletracer->totalParticles = totalParticles;
+    particletracer->pointer        = totalParticles + 1;
+
+    memcpy(particletracer->particlePool, tempPool, totalParticles*sizeof(Particle));
+}

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

@@ -0,0 +1,58 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __PARTICLETRACING_H_
+#define __PARTICLETRACING_H_
+#include "allocate.h"
+#include "parameter.h"
+#include "solver.h"
+#include "particletracing.h"
+#include <stdbool.h>
+#include <mpi.h>
+
+#define NDIMS 2
+
+typedef enum COORD { X = 0, Y, NCOORD } COORD;
+typedef struct
+{
+    double x, y;
+    bool flag;
+} Particle;
+
+typedef struct {
+    int numberOfParticles, totalParticles;
+    double startTime, injectTimePeriod, writeTimePeriod, lastInjectTime, lastUpdateTime, lastWriteTime;
+
+    int estimatedNumParticles;
+
+    double dx, dy, dt;
+    Particle* linSpaceLine;
+    Particle* particlePool;
+
+    int pointer;
+    
+    double imax, jmax, xlength, ylength, imaxLocal, jmaxLocal;
+
+    double x1, y1, x2, y2;
+
+    MPI_Comm comm;
+    int rank, size;
+    int iNeighbours[NDIMS], jNeighbours[NDIMS];
+    int coords[NDIMS], dims[NDIMS];
+
+    double xLocal, yLocal, xOffset, yOffset;
+} ParticleTracer;
+
+void initParticleTracer(ParticleTracer*, Parameter*, Solver* );
+void injectParticles(ParticleTracer*);
+void advanceParticles(ParticleTracer*, double* , double*, double);
+void freeParticles(ParticleTracer*);
+void writeParticles(ParticleTracer*);
+void printParticleTracerParameters(ParticleTracer*);
+void printParticles(ParticleTracer*);
+void trace(ParticleTracer*, double* , double* , double );
+void compress(ParticleTracer* );
+#endif

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

@@ -0,0 +1,128 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "vtkWriter.h"
+
+static float floatSwap(float f)
+{
+    union {
+        float f;
+        char b[4];
+    } dat1, dat2;
+
+    dat1.f    = f;
+    dat2.b[0] = dat1.b[3];
+    dat2.b[1] = dat1.b[2];
+    dat2.b[2] = dat1.b[1];
+    dat2.b[3] = dat1.b[0];
+    return dat2.f;
+}
+
+static void writeHeader(VtkOptions* o, int ts)
+{
+    fprintf(o->fh, "# vtk DataFile Version 3.0\n");
+    fprintf(o->fh, "PAMPI cfd solver particle tracing file\n");
+    if (o->fmt == ASCII) {
+        fprintf(o->fh, "ASCII\n");
+    } else if (o->fmt == BINARY) {
+        fprintf(o->fh, "BINARY\n");
+    }
+
+    fprintf(o->fh, "DATASET UNSTRUCTURED_GRID\n");
+    fprintf(o->fh, "FIELD FieldData 2\n");
+    fprintf(o->fh, "TIME 1 1 double\n");
+    fprintf(o->fh, "%d\n", ts);
+    fprintf(o->fh, "CYCLE 1 1 int\n");
+    fprintf(o->fh, "1\n");
+}
+
+void vtkOpen(VtkOptions* o, char* problem, int ts)
+{
+    o->fh = fopen(problem, "w");
+
+    if (o->fh == NULL) {
+        printf("vtkWriter not initialize! Call vtkOpen first!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    writeHeader(o, ts);
+
+    printf("Writing VTK output for %s\n", problem);
+}
+
+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);
+
+
+    for (int i = 0; i < o->particletracer->totalParticles; ++i)
+    {
+        double x = particlePool[i].x;
+        double y = particlePool[i].y;
+        fprintf(o->fh, "%.2f %.2f 0.0\n", x, y);
+    }
+
+    fprintf(o->fh, "CELLS %d %d\n", o->particletracer->totalParticles, 2 * o->particletracer->totalParticles);
+
+
+    for (int i = 0; i < o->particletracer->totalParticles; ++i)
+    {
+        fprintf(o->fh, "1 %d\n", i);
+    }
+
+    fprintf(o->fh, "CELL_TYPES %d\n", o->particletracer->totalParticles);
+
+
+    for (int i = 0; i < o->particletracer->totalParticles; ++i)
+    {
+        fprintf(o->fh, "1\n");
+    }
+
+    /*
+    for (int k = 0; k < kmax; k++) {
+        for (int j = 0; j < jmax; j++) {
+            for (int i = 0; i < imax; i++) {
+                if (o->fmt == ASCII) {
+                    fprintf(o->fh,
+                        "%f %f %f\n",
+                        G(vec.u, i, j, k),
+                        G(vec.v, i, j, k),
+                        G(vec.w, i, j, k));
+                } else if (o->fmt == BINARY) {
+                    fwrite((float[3]) { floatSwap(G(vec.u, i, j, k)),
+                               floatSwap(G(vec.v, i, j, k)),
+                               floatSwap(G(vec.w, i, j, k)) },
+                        sizeof(float),
+                        3,
+                        o->fh);
+                }
+            }
+        }
+    }
+    if (o->fmt == BINARY) fprintf(o->fh, "\n");
+
+    */
+}
+
+void vtkClose(VtkOptions* o)
+{
+    fclose(o->fh);
+    o->fh = NULL;
+}

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

@@ -0,0 +1,26 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __VTKWRITER_H_
+#define __VTKWRITER_H_
+#include <stdio.h>
+
+#include "particletracing.h"
+#include "solver.h"
+
+typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
+
+typedef struct VtkOptions {
+    VtkFormat fmt;
+    FILE* fh;
+    Solver* solver;
+    ParticleTracer* particletracer;
+} VtkOptions;
+
+extern void vtkOpen(VtkOptions* opts, char* filename, int ts);
+extern void vtkParticle(VtkOptions* opts, char* name);
+extern void vtkClose(VtkOptions* opts);
+#endif // __VTKWRITER_H_