Initial working 2D-seq multigrid

BasicSolver/2D-seq-multigrid/backstep.par

@@ -32,7 +32,7 @@ jmax       45	   # number of interior cells in y-direction
 # Time Data:
 # ---------
 
-te      60.0   # final time
+te      100.0   # final time
 dt      0.02    # time stepsize
 tau     0.5     # safety factor for time stepsize control (<0 constant delt)
 
@@ -44,12 +44,13 @@ eps      0.0001   # stopping tolerance for pressure iteration
 rho      0.52
 omg      1.8       # relaxation parameter for SOR iteration
 gamma    0.9       # upwind differencing factor gamma
+levels   5         # Multigrid levels
 
 # Particle Tracing Data:
 # -----------------------
 
 numberOfParticles   200
-startTime           0
+startTime           100
 injectTimePeriod    1.0
 writeTimePeriod     0.5
 

BasicSolver/2D-seq-multigrid/canal.par

@@ -50,7 +50,7 @@ levels   5         # Multigrid levels
 # -----------------------
 
 numberOfParticles   60
-startTime           1000
+startTime           100
 injectTimePeriod    2.0
 writeTimePeriod     0.5
 

BasicSolver/2D-seq-multigrid/dcavity.par

@@ -44,6 +44,7 @@ eps      0.001		# stopping tolerance for pressure iteration
 rho      0.5
 omg      1.8		# relaxation parameter for SOR iteration
 gamma    0.9		# upwind differencing factor gamma
+levels   5         # Multigrid levels
 
 # Particle Tracing Data:
 # -----------------------

BasicSolver/2D-seq-multigrid/karman.par

@@ -44,12 +44,13 @@ eps      0.001   # stopping tolerance for pressure iteration
 rho      0.52
 omg      1.75      # relaxation parameter for SOR iteration
 gamma    0.9       # upwind differencing factor gamma
+levels   5         # Multigrid levels
 
 # Particle Tracing Data:
 # -----------------------
 
 numberOfParticles   200
-startTime           50
+startTime           201
 injectTimePeriod    1.0
 writeTimePeriod     0.5
 

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

@@ -39,6 +39,8 @@ int main(int argc, char** argv)
 
     printParameter(&params);
     initSolver(&solver, &params);
+    printf("initsolver done\n");
+
     initParticleTracer(&particletracer, &params);
     printParticleTracerParameters(&particletracer);
 
@@ -62,11 +64,12 @@ int main(int argc, char** argv)
         setBoundaryConditions(&solver);
         setSpecialBoundaryCondition(&solver);
         setObjectBoundaryCondition(&solver);
-
         computeFG(&solver);
         computeRHS(&solver);
+
         if (nt % 100 == 0) normalizePressure(&solver);
         multiGrid(&solver);
+
         adaptUV(&solver);
 
         /* Added function for particle tracing. Will inject and advance particles as per
@@ -84,18 +87,11 @@ int main(int argc, char** argv)
     }
     printf("Total particles : %d\n", particletracer.totalParticles);
 
-    // print(&solver, solver.p);
-
     E = getTimeStamp();
     stopProgress();
 
     freeParticles(&particletracer);
     
-    // printf("\nU : \n\n");
-    // print(&solver, solver.u);
-    // printf("\nV : \n\n");
-    // print(&solver, solver.v);
-
     printf("Solution took %.2fs\n", E - S);
     writeResult(&solver);
     return EXIT_SUCCESS;

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

@@ -120,9 +120,11 @@ void advanceParticles(ParticleTracer* particletracer,
 
 void freeParticles(ParticleTracer* particletracer)
 {
+    if (particletracer->particlePool != NULL) {
         free(particletracer->particlePool);
         free(particletracer->linSpaceLine);
     }
+}
 
 void writeParticles(ParticleTracer* particletracer)
 {
@@ -184,6 +186,7 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params)
     particletracer->imax = params->imax;
     particletracer->jmax = params->jmax;
 
+    if (params->te > params->startTime) {
         particletracer->estimatedNumParticles = ((params->te - params->startTime) + 2) *
                                                 params->numberOfParticles;
 
@@ -201,6 +204,12 @@ void initParticleTracer(ParticleTracer* particletracer, Parameter* params)
             particletracer->linSpaceLine[i].flag = true;
         }
     }
+    else
+    {
+         particletracer->particlePool = NULL;
+         particletracer->linSpaceLine = NULL;
+    }
+}
 
 void printParticleTracerParameters(ParticleTracer* particletracer)
 {

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

@@ -24,7 +24,7 @@
 #define E(i, j)    e[(j) * (imax + 2) + (i)]
 #define R(i, j)    r[(j) * (imax + 2) + (i)]
 #define oldR(i, j) oldr[(j) * (imax + 2) + (i)]
-
+#define oldE(i, j) olde[(j) * (imax + 2) + (i)]
 #define RHS(i, j)  rhs[(j) * (imax + 2) + (i)]
 
 static double distance(double i, double j, double iCenter, double jCenter)
@@ -126,10 +126,10 @@ void initSolver(Solver* solver, Parameter* params)
     solver->f   = allocate(64, size);
     solver->g   = allocate(64, size);
 
-    solver->r = malloc(levels * sizeof(int*));
+    solver->r = malloc(levels * sizeof(double*));
-    solver->e = malloc(levels * sizeof(int*));
+    solver->e = malloc(levels * sizeof(double*));
 
-    for (int j = 1; j <= levels; ++j) {
+    for (int j = 0; j < levels; ++j) {
 
         solver->r[j] = allocate(64, size);
         solver->e[j] = allocate(64, size);
@@ -144,7 +144,7 @@ void initSolver(Solver* solver, Parameter* params)
         solver->f[i]   = 0.0;
         solver->g[i]   = 0.0;
         solver->s[i]   = NONE;
-        for (int j = 1; j <= levels; ++j) {
+        for (int j = 0; j < levels; ++j) {
 
             solver->r[j][i] = 0.0;
             solver->e[j][i] = 0.0;
@@ -198,6 +198,7 @@ void initSolver(Solver* solver, Parameter* params)
         break;
     }
 
+    if (params->shape != NOSHAPE) {
         for (int j = 1; j < jmax + 1; j++) {
             for (int i = 1; i < imax + 1; i++) {
                 // if( S(i+1,j+1) == NONE && S(i+1,j) == NONE && S(i,j+1) == NONE &&
@@ -226,24 +227,25 @@ void initSolver(Solver* solver, Parameter* params)
                     S(i, j) = RIGHT; // RIGHT
                 if (S(i, j + 1) == NONE && S(i, j - 1) == LOCAL && S(i, j) == LOCAL)
                     S(i, j) = TOP; // BOTTOM
-            if (S(i - 1, j - 1) == NONE && S(i, j - 1) == NONE && S(i - 1, j) == NONE &&
+                if (S(i - 1, j - 1) == NONE && S(i, j - 1) == NONE &&
-                S(i + 1, j + 1) == LOCAL &&
+                    S(i - 1, j) == NONE && S(i + 1, j + 1) == LOCAL &&
                     (S(i, j) == LOCAL || S(i, j) == LEFT || S(i, j) == BOTTOM))
                     S(i, j) = BOTTOMLEFT; // TOPLEFT
-            if (S(i + 1, j - 1) == NONE && S(i, j - 1) == NONE && S(i + 1, j) == NONE &&
+                if (S(i + 1, j - 1) == NONE && S(i, j - 1) == NONE &&
-                S(i - 1, j + 1) == LOCAL &&
+                    S(i + 1, j) == NONE && S(i - 1, j + 1) == LOCAL &&
                     (S(i, j) == LOCAL || S(i, j) == RIGHT || S(i, j) == BOTTOM))
                     S(i, j) = BOTTOMRIGHT; // TOPRIGHT
-            if (S(i - 1, j + 1) == NONE && S(i - 1, j) == NONE && S(i, j + 1) == NONE &&
+                if (S(i - 1, j + 1) == NONE && S(i - 1, j) == NONE &&
-                S(i + 1, j - 1) == LOCAL &&
+                    S(i, j + 1) == NONE && S(i + 1, j - 1) == LOCAL &&
                     (S(i, j) == LOCAL || S(i, j) == LEFT || S(i, j) == TOP))
                     S(i, j) = TOPLEFT; // BOTTOMLEFT
-            if (S(i + 1, j + 1) == NONE && S(i + 1, j) == NONE && S(i, j + 1) == NONE &&
+                if (S(i + 1, j + 1) == NONE && S(i + 1, j) == NONE &&
-                S(i - 1, j - 1) == LOCAL &&
+                    S(i, j + 1) == NONE && S(i - 1, j - 1) == LOCAL &&
                     (S(i, j) == LOCAL || S(i, j) == RIGHT || S(i, j) == TOP))
                     S(i, j) = TOPRIGHT; // BOTTOMRIGHT
             }
         }
+    }
 
 #ifdef VERBOSE
     printConfig(solver);
@@ -649,7 +651,7 @@ void adaptUV(Solver* solver)
     }
 }
 
-void smoothRB(Solver* solver)
+double smoothRB(Solver* solver)
 {
     int imax      = solver->imax;
     int jmax      = solver->jmax;
@@ -689,12 +691,15 @@ void smoothRB(Solver* solver)
     }
 
     res = res / (double)(imax * jmax);
+    return res;
 }
 
 void multiGrid(Solver* solver)
 {
-
+    double res = 0.0;
-    if (solver->currentlevel == solver->levels) {
+    int imax = solver->imax;
+    int jmax = solver->jmax;
+    if (solver->currentlevel == (solver->levels - 1)) {
         for (int i = 0; i < 5; i++) {
             smoothRB(solver);
         }
@@ -703,11 +708,9 @@ void multiGrid(Solver* solver)
 
     for (int i = 0; i < 5; i++) {
         smoothRB(solver);
-        if (solver->currentlevel == 1) {
+        if (solver->currentlevel == 0) {
 
             double* p = solver->p;
-            int imax  = solver->imax;
-            int jmax  = solver->jmax;
 
             for (int i = 1; i < imax + 1; i++) {
                 P(i, 0)        = P(i, 1);
@@ -720,14 +723,18 @@ void multiGrid(Solver* solver)
             }
         }
     }
-
+    Solver coarseSolver = copySolver(solver);
-    Solver* coarseSolver = copySolver(solver);
 
     // restrict
     restrictMG(solver);
 
+    coarseSolver.p   = solver->e[coarseSolver.currentlevel];
+    coarseSolver.rhs = solver->r[coarseSolver.currentlevel];
+    coarseSolver.imax /= 2;
+    coarseSolver.jmax /= 2;
+
     // MGSolver on residual and error.
-    multiGrid(coarseSolver);
+    multiGrid(&coarseSolver);
 
     // prolongate
     prolongate(solver);
@@ -735,11 +742,9 @@ void multiGrid(Solver* solver)
     // correct p on finest level using residual
     correct(solver);
 
-    if (solver->currentlevel == 1) {
+    if (solver->currentlevel == 0) {
 
         double* p = solver->p;
-        int imax  = solver->imax;
-        int jmax  = solver->jmax;
 
         for (int i = 1; i < imax + 1; i++) {
             P(i, 0)        = P(i, 1);
@@ -753,16 +758,33 @@ void multiGrid(Solver* solver)
     }
 
     for (int i = 0; i < 5; i++) {
-        smoothRB(solver);
+        res = smoothRB(solver);
+        if (solver->currentlevel == 0) {
+
+            double* p = solver->p;
+
+            for (int i = 1; i < imax + 1; i++) {
+                P(i, 0)        = P(i, 1);
+                P(i, jmax + 1) = P(i, jmax);
             }
+
+            for (int j = 1; j < jmax + 1; j++) {
+                P(0, j)        = P(1, j);
+                P(imax + 1, j) = P(imax, j);
+            }
+        }
+    }
+    #ifdef VERBOSE
+        if (solver->currentlevel == 0) {printf("Residuum: %.6f\n", res); }
+#endif
 }
 
 void restrictMG(Solver* solver)
 {
     int imax     = solver->imax;
     int jmax     = solver->jmax;
-    double* r    = solver->r[solver->currentlevel];
+    double* r    = solver->r[solver->currentlevel + 1];
-    double* oldr = solver->r[solver->currentlevel - 1];
+    double* oldr = solver->r[solver->currentlevel];
 
     for (int j = 1; j < jmax + 1; j++) {
         for (int i = 1; i < imax + 1; ++i) {
@@ -780,12 +802,12 @@ void prolongate(Solver* solver)
 {
     int imax     = solver->imax;
     int jmax     = solver->jmax;
-    double* r    = solver->r[solver->currentlevel];
+    double* olde = solver->r[solver->currentlevel + 1];
-    double* oldr = solver->r[solver->currentlevel - 1];
+    double* e    = solver->r[solver->currentlevel];
 
     for (int j = 2; j < jmax + 1; j += 2) {
         for (int i = 2; i < imax + 1; i += 2) {
-            oldR(i, j) = R(i / 2, j / 2);
+            E(i, j) = oldE(i / 2, j / 2);
         }
     }
 }
@@ -803,36 +825,36 @@ void correct(Solver* solver)
     }
 }
 
-Solver* copySolver(Solver* solver)
+Solver copySolver(Solver* solver)
 {
-    Solver* newSolver;
+    Solver newSolver;
-    newSolver->problem      = solver->problem;
+    newSolver.problem      = solver->problem;
-    newSolver->bcLeft       = solver->bcLeft;
+    newSolver.bcLeft       = solver->bcLeft;
-    newSolver->bcRight      = solver->bcRight;
+    newSolver.bcRight      = solver->bcRight;
-    newSolver->bcBottom     = solver->bcBottom;
+    newSolver.bcBottom     = solver->bcBottom;
-    newSolver->bcTop        = solver->bcTop;
+    newSolver.bcTop        = solver->bcTop;
-    newSolver->imax         = solver->imax;
+    newSolver.imax         = solver->imax;
-    newSolver->jmax         = solver->jmax;
+    newSolver.jmax         = solver->jmax;
-    newSolver->xlength      = solver->xlength;
+    newSolver.xlength      = solver->xlength;
-    newSolver->ylength      = solver->ylength;
+    newSolver.ylength      = solver->ylength;
-    newSolver->dx           = solver->xlength / solver->imax;
+    newSolver.dx           = solver->xlength / solver->imax;
-    newSolver->dy           = solver->ylength / solver->jmax;
+    newSolver.dy           = solver->ylength / solver->jmax;
-    newSolver->eps          = solver->eps;
+    newSolver.eps          = solver->eps;
-    newSolver->omega        = solver->omega;
+    newSolver.omega        = solver->omega;
-    newSolver->itermax      = solver->itermax;
+    newSolver.itermax      = solver->itermax;
-    newSolver->re           = solver->re;
+    newSolver.re           = solver->re;
-    newSolver->gx           = solver->gx;
+    newSolver.gx           = solver->gx;
-    newSolver->gy           = solver->gy;
+    newSolver.gy           = solver->gy;
-    newSolver->dt           = solver->dt;
+    newSolver.dt           = solver->dt;
-    newSolver->te           = solver->te;
+    newSolver.te           = solver->te;
-    newSolver->tau          = solver->tau;
+    newSolver.tau          = solver->tau;
-    newSolver->gamma        = solver->gamma;
+    newSolver.gamma        = solver->gamma;
-    newSolver->rho          = solver->rho;
+    newSolver.rho          = solver->rho;
-    newSolver->levels       = solver->levels;
+    newSolver.levels       = solver->levels;
-    newSolver->currentlevel = solver->currentlevel;
+    newSolver.currentlevel = solver->currentlevel + 1;
 
-    newSolver->r = solver->r;
+    newSolver.r = solver->r;
-    newSolver->e = solver->e;
+    newSolver.e = solver->e;
 
     return newSolver;
 }

BasicSolver/2D-seq-multigrid/src/solver.h

@@ -48,12 +48,12 @@ typedef struct {
 
 extern void initSolver(Solver*, Parameter*);
 extern void computeRHS(Solver*);
-extern void smoothRB(Solver*);
+extern double smoothRB(Solver*);
 extern void residualRB(Solver*);
 extern void restrictMG(Solver*);
 extern void prolongate(Solver*);
 extern void correct(Solver*);
-extern Solver* copySolver(Solver*);
+extern Solver copySolver(Solver*);
 extern void multiGrid(Solver*);
 extern void normalizePressure(Solver*);
 extern void computeTimestep(Solver*);

BasicSolver/2D-seq-multigrid/vector.plot

@@ -1,7 +1,6 @@
 set terminal png size 3600,768 enhanced font ,28
 set output 'velocity.png'
-set xrange[0:7]
+
-set yrange[0:1.5]
 set size ratio -1
 set datafile separator whitespace
 

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

@@ -259,7 +259,6 @@ void solve(Solver* solver)
     double res    = 1.0;
 
     while ((res >= epssq) && (it < itermax)) {
-        res = 0.0;
 
         for (int j = 1; j < jmax + 1; j++) {
             for (int i = 1; i < imax + 1; i++) {