Corrections in enhanced solver

BasicSolver/2D-seq/canal.par

@@ -15,11 +15,11 @@ bcRight    3			#
 gx    0.0      # Body forces (e.g. gravity)
 gy    0.0      #
 
-re    100.0	   # Reynolds number
+re    500.0	   # Reynolds number
 
 u_init    1.0      # initial value for velocity in x-direction
 v_init    0.0      # initial value for velocity in y-direction
-p_init    0.0      # initial value for pressure
+p_init    1.0      # initial value for pressure
 
 # Geometry Data:
 # -------------
@@ -32,7 +32,7 @@ jmax       50	   # number of interior cells in y-direction
 # Time Data:
 # ---------
 
-te      60.0   # final time
+te      700.0   # final time
 dt      0.02    # time stepsize
 tau     0.5     # safety factor for time stepsize control (<0 constant delt)
 
@@ -49,9 +49,9 @@ gamma    0.9       # upwind differencing factor gamma
 # -----------------------
 
 numberOfParticles   60
-startTime           10
+startTime           0
 injectTimePeriod    2.0
-writeTimePeriod     1.0
+writeTimePeriod     0.5
 
 x1                  1.0
 y1                  0.0
@@ -63,10 +63,10 @@ y2                  4.0
 # Shape 0 disable, 1 Rectangle/Square, 2 Circle
 
 shape               1 
-xCenter             10.0
-yCenter             2
-xRectLength         3.0
-yRectLength         1.0
+xCenter             4.0
+yCenter             1
+xRectLength         8.0
+yRectLength         2.0
 circleRadius        1.0
 
 #===============================================================================

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

@@ -72,13 +72,13 @@ int main (int argc, char** argv)
 
         computeFG(&solver);
         computeRHS(&solver);
-        //if (nt % 100 == 0) normalizePressure(&solver);
+        if (nt % 100 == 0) normalizePressure(&solver);
         solver_generic[variant - 1](&solver);
         adaptUV(&solver);
         
 
         /* Added function for particle tracing. Will inject and advance particles as per timePeriod */
-        // trace(&particletracer, solver.u, solver.v, t);
+        trace(&particletracer, solver.u, solver.v, t);
 
         t += solver.dt;
         nt++;

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

@@ -318,6 +318,19 @@ void solveRB(Solver* solver)
         res = 0.0;
         jsw = 1;
 
+
+        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);
+        }
+
+        // setObjectPBoundaryCondition(solver);
+
         for (pass = 0; pass < 2; pass++) 
         {
             isw = jsw;
@@ -342,18 +355,6 @@ void solveRB(Solver* solver)
             jsw = 3 - jsw;
         }
 
-        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);
-        }
-
-        //setObjectPBoundaryCondition(solver);
-
         res = res / (double)(imax * jmax);
         
 #ifdef DEBUG
@@ -630,7 +631,7 @@ void setObjectBoundaryCondition(Solver* solver)
                     V(i,j)      = 0.0;
                     break;
                 case LEFT:
-                    U(i,j)      = 0.0;
+                    U(i-1,j)      = 0.0;
                     V(i,j)      = -V(i-1,j);
                     V(i,j-1)    = -V(i-1,j-1);
                     break;
@@ -643,13 +644,13 @@ void setObjectBoundaryCondition(Solver* solver)
                     U(i,j)      = -U(i,j+1);
                     U(i-1,j)    = 0.0;
                     V(i,j)      = 0.0;
-                    V(i,j+1)    = -V(i-1,j+1);
+                    V(i,j-1)    = -V(i-1,j-1);
                     break;
                 case TOPRIGHT:
                     U(i,j)      = 0.0;
-                    U(i+1,j)    = -U(i+1,j+1);
+                    U(i-1,j)    = -U(i-1,j+1);
                     V(i,j)      = 0.0;
-                    V(i,j+1)    = -V(i+1,j+1);
+                    V(i,j-1)    = -V(i+1,j-1);
                     break;
                 case BOTTOMLEFT:
                     U(i,j)      = -U(i,j-1);
@@ -659,7 +660,7 @@ void setObjectBoundaryCondition(Solver* solver)
                     break;
                 case BOTTOMRIGHT:
                     U(i,j)      = 0.0;
-                    U(i+1,j)    = -U(i+1,j-1);
+                    U(i-1,j)    = -U(i-1,j-1);
                     V(i,j)      = -V(i,j+1);
                     V(i,j-1)    = 0.0;
                     break;

BasicSolver/2D-seq/vis_files/animate.plot

@@ -3,7 +3,7 @@ set term gif animate delay 50
 set output "trace.gif"
 set xrange [0:30]
 set yrange [0:4]
-do for [ts=0:600] {
+do for [ts=0:1340] {
     plot "particles_".ts.".dat" with points pointtype 7
 }
 unset output