Partially working vtk implementation

BasicSolver/2D-seq/animate.plot

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

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

@@ -15,6 +15,7 @@
 #include "solver.h"
 #include "particletracing.h"
 #include "timing.h"
+#include "vtkWriter.h"
 
 enum VARIANT { SOR = 1, RB, RBA };
 
@@ -54,7 +55,7 @@ int main (int argc, char** argv)
     double t   = 0.0;
     int nt     = 0;
 
-    void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
+    void (*solver_generic[])() = {solve, solveRB, solveRBA};
 
     S = getTimeStamp();
 
@@ -66,7 +67,7 @@ int main (int argc, char** argv)
         computeFG(&solver);
         computeRHS(&solver);
         if (nt % 100 == 0) normalizePressure(&solver);
-        (*solver_generic[variant - 1])(&solver);
+        solver_generic[variant - 1](&solver);
         adaptUV(&solver);
 
         /* Added function for particle tracing. Will inject and advance particles as per timePeriod */
@@ -75,13 +76,12 @@ int main (int argc, char** argv)
         t += solver.dt;
         nt++;
 
-    /*
+    
-    #ifdef VERBOSE
+#ifdef VERBOSE
-            printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
+        printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
-    #else
+#else
-            printProgress(t);
+        printProgress(t);
-    #endif
+#endif
-    */
     }
     printf("Total particles : %d\n", particletracer.totalParticles);
     E = getTimeStamp();

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

@@ -10,29 +10,11 @@
 #include <stdlib.h>
 #include <string.h>
 
-#include "allocate.h"
+#include "vtkWriter.h"
-#include "parameter.h"
-#include "solver.h"
-#include "particletracing.h"
-#include "util.h"
 
 #define U(i, j)  u[(j) * (imax + 2) + (i)]
 #define V(i, j)  v[(j) * (imax + 2) + (i)]
 
-static void print(Solver* solver, double* grid)
-{
-    int imax = solver->imax;
-
-    for (int j = 0; j < solver->jmax + 2; j++) {
-        printf("%02d: ", j);
-        for (int i = 0; i < solver->imax + 2; i++) {
-            printf("%3.2f  ", grid[j * (imax + 2) + i]);
-        }
-        printf("\n");
-    }
-    fflush(stdout);
-}
-
 static int ts = 0;
 void printParticles(ParticleTracer* particletracer)
 {
@@ -55,7 +37,7 @@ void injectParticles(ParticleTracer* particletracer)
     }
 }
 
-void advanceParticles(ParticleTracer* particletracer, Solver* solver, double* restrict u, double* restrict v, double time)
+void advanceParticles(ParticleTracer* particletracer, double* restrict u, double* restrict v, double time)
 {
     int imax = particletracer->imax;
     int jmax = particletracer->jmax;
@@ -129,30 +111,15 @@ void freeParticles(ParticleTracer* particletracer)
     free(particletracer->linSpaceLine);
 }
 
-void writeParticles(ParticleTracer* particletracer)
+void writeParticles(ParticleTracer* particletracer, Solver* solver)
 {
-    FILE* fp;
+    VtkOptions opts = { .solver = solver, .particletracer = particletracer };
-    Particle* particlePool = particletracer->particlePool;
 
     char filename[50];
-    snprintf(filename, 50, "particles_%d.dat", ts++);
+    snprintf(filename, 50, "vtk_files/particles_%d.vtk", ts++);
-    fp = fopen(filename, "w");
+    vtkOpen(&opts, filename);
-
+    vtkParticle(&opts, "particle");
-    if (fp == NULL) {
+    vtkClose(&opts);
-        printf("Error!\n");
-        exit(EXIT_FAILURE);
-    }
-
-    for (int i = 0; i < particletracer->totalParticles; ++i)
-    {
-         if(particlePool[i].flag == true)
-        {
-            double x = particlePool[i].x;
-            double y = particlePool[i].y;
-            fprintf(fp, "%f %f\n", x, y);
-        }
-    }
-    fclose(fp);
 }
 
 void initParticleTracer(ParticleTracer* particletracer, Parameter* params)
@@ -221,11 +188,11 @@ void trace(ParticleTracer* particletracer, Solver* solver, double* u, double* v,
         }
         if( (time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod)
         {
-            writeParticles(particletracer);
+            writeParticles(particletracer, solver);
             particletracer->lastWriteTime = time;
 
         }
-        advanceParticles(particletracer, solver, u, v, time);
+        advanceParticles(particletracer, u, v, time);
         particletracer->lastUpdateTime = time;   
     }
 }

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

@@ -6,7 +6,10 @@
  */
 #ifndef __PARTICLETRACING_H_
 #define __PARTICLETRACING_H_
+#include "allocate.h"
 #include "parameter.h"
+#include "solver.h"
+#include "particletracing.h"
 #include <stdbool.h>
 
 typedef enum COORD { X = 0, Y, NCOORD } COORD;
@@ -35,9 +38,9 @@ typedef struct {
 
 void initParticleTracer(ParticleTracer*, Parameter*);
 void injectParticles(ParticleTracer*);
-void advanceParticles(ParticleTracer*, Solver*, double* , double*, double);
+void advanceParticles(ParticleTracer*, double* , double*, double);
 void freeParticles(ParticleTracer*);
-void writeParticles(ParticleTracer*);
+void writeParticles(ParticleTracer*, Solver*);
 void printParticleTracerParameters(ParticleTracer*);
 void printParticles(ParticleTracer*);
 void trace(ParticleTracer*, Solver*, double* , double* , double );

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

@@ -150,9 +150,6 @@ void solve(Solver* solver)
     int it        = 0;
     double res    = 1.0;
 
-    // printf("\nBefore P : \n");
-    // print(solver, p);
-
     while ((res >= epssq) && (it < itermax)) {
         res = 0.0;
 
@@ -179,20 +176,15 @@ void solve(Solver* solver)
         }
 
         res = res / (double)(imax * jmax);
-        /*
+        
 #ifdef DEBUG
         printf("%d Residuum: %e\n", it, res);
 #endif
-*/
         it++;
     }
-    // printf("\nAfter P : \n");
-    // print(solver, p);
-/*
 #ifdef VERBOSE
     printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
 #endif
-*/
 }
 
 void solveRB(Solver* solver)
@@ -213,9 +205,6 @@ void solveRB(Solver* solver)
     double res    = 1.0;
     int pass, jsw, isw;
 
-    // printf("\nBefore P : \n");
-    // print(solver, p);
-
     while ((res >= epssq) && (it < itermax)) {
         res = 0.0;
         jsw = 1;
@@ -249,20 +238,16 @@ void solveRB(Solver* solver)
         }
 
         res = res / (double)(imax * jmax);
-        /*
+        
 #ifdef DEBUG
         printf("%d Residuum: %e\n", it, res);
 #endif
-*/
+
         it++;
     }
-    // printf("\nAfter P : \n");
-    // print(solver, p);
-/*
 #ifdef VERBOSE
     printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
 #endif
-*/
 }
 
 void solveRBA(Solver* solver)
@@ -479,12 +464,6 @@ void setBoundaryConditions(Solver* solver)
     case PERIODIC:
         break;
     }
-
-    // printf("\nBoundary condition U : \n");
-    // print(solver, u);
-
-    // printf("\nBoundary condition V : \n");
-    // print(solver, v);
 }
 
 void setSpecialBoundaryCondition(Solver* solver)
@@ -507,10 +486,6 @@ void setSpecialBoundaryCondition(Solver* solver)
             U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
         }
     }
-
-    // printf("\nSpecial Boundary condition U : \n");
-    // print(solver, u);
-
 }
 
 void computeFG(Solver* solver)
@@ -584,13 +559,6 @@ void computeFG(Solver* solver)
         G(i, 0)    = V(i, 0);
         G(i, jmax) = V(i, jmax);
     } 
-    
-    // printf("\nComputed F : \n");
-    // print(solver, f);
-
-    // printf("\nComputed G : \n");
-    // print(solver, g);
-    
 }
 
 void adaptUV(Solver* solver)
@@ -611,17 +579,6 @@ void adaptUV(Solver* solver)
             V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
         }
     }
-
-    // printf("\nApadted U : \n");
-    // print(solver, u);
-
-    // printf("\nAdapted V : \n");
-    // print(solver, v);
-
-    // char *b;
-    // b='\0';
-    // printf("\nDEBUGGER : ");
-    // scanf("%s", &b);
 }
 
 void writeResult(Solver* solver)

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

@@ -0,0 +1,112 @@
+/*
+ * 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)
+{
+    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 STRUCTURED_POINTS\n");
+    fprintf(o->fh, "DIMENSIONS %d %d 0\n", o->solver->imax, o->solver->jmax);
+    fprintf(o->fh,
+        "ORIGIN %f %f 0.0\n",
+        o->solver->dx * 0.5,
+        o->solver->dy * 0.5);
+    fprintf(o->fh, "SPACING %f %f 0\n", o->solver->dx, o->solver->dy);
+    fprintf(o->fh, "POINT_DATA %d\n", o->particletracer->totalParticles);
+}
+
+void vtkOpen(VtkOptions* o, char* problem)
+{
+    o->fh = fopen(problem, "w");
+    writeHeader(o);
+
+    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, "VECTORS %s float\n", name);
+    fprintf(o->fh, "LOOKUP_TABLE default\n");
+
+
+    for (int i = 0; i < o->particletracer->totalParticles; ++i)
+    {
+        if(particlePool[i].flag == true)
+        {
+            double x = particlePool[i].x;
+            double y = particlePool[i].y;
+            fprintf(o->fh, "%f %f 0.0\n", x, y);
+        }
+    }
+
+    /*
+    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-seq/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);
+extern void vtkParticle(VtkOptions* opts, char* name);
+extern void vtkClose(VtkOptions* opts);
+#endif // __VTKWRITER_H_

BasicSolver/3D-seq/src/main.c

@@ -100,7 +100,7 @@ int main(int argc, char** argv)
 
     timeStart = getTimeStamp();
 
-    void (*solver_generic[])(solver) = {solve, solveRB, solveRBA};
+    void (*solver_generic[])() = {solve, solveRB, solveRBA};
 
     while (t <= te) {
         if (tau > 0.0) computeTimestep(&s);
@@ -108,7 +108,7 @@ int main(int argc, char** argv)
         setSpecialBoundaryCondition(&s);
         computeFG(&s);
         computeRHS(&s);
-        (*solver_generic[variant - 1])(&s);
+        solver_generic[variant - 1](&s);
         adaptUV(&s);
         t += s.dt;