13
$\begingroup$

This algorithm is 3D extension of our 2D algorithm published on this page and here.
We suppose that with this code we can simulate transition from laminar to turbulent flow. In this example we compute viscous flow around cuboid at Reynolds number $Re=6250$. Note, that code has been tested up to $Re=10^6$.

dif = 1/6250; pec = .72; U0 = 1.; V0 = 0.; W0 = 0.; dn0 = 1.; n = 80;
n1 = n + 1; sm = 500; r = 20; den = ConstantArray[dn0 , {n1, n1, n1}];
u0 = ConstantArray[U0, {n1, n1, n1}];
v0 = ConstantArray[V0, {n1, n1, n1}]; w0 = 
 ConstantArray[W0, {n1, n1, n1}]; Do[u0[[i, j, k]] = 0; 
 v0[[i, j, k]] = 0; 
 w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
  Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];

periodic[n_, up_, ud_, ul_, ur_, ub_] := 
  Module[{bd = ub}, 
   Do[bd[[n + 1, i, j]] = bd[[2, i, j]]; 
    bd[[1, i, j]] = bd[[n, i, j]];, {i, 2, n}, {j, 2, n}];
   Do[bd[[i, 1, j]] = ud;
    bd[[i, n + 1, j]] = up; bd[[i, j, 1]] = ul;
    bd[[i, j, n + 1]] = ur;, {i, 1, n + 1}, {j, 1, n + 1}];
   bd];

diffuse[n_, r_, a_, c_, c0_] := 
  Module[{c1 = c}, 
   Do[Do[Do[
       Do[c1[[i, j, 
            k]] = (c0[[i, j, k]] + 
              a (c1[[i - 1, j, k]] + c1[[i + 1, j, k]] + 
                 c1[[i, j - 1, k]] + c1[[i, j + 1, k]] + 
                 c1[[i, j, k - 1]] + c1[[i, j, k + 1]]))/(1 + 
              6 a);, {k, 2, n}];, {j, 2, n}];, {i, 2, n}];
    Do[c1[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
      Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 
      1}];, {k1, 0, r}];
   c1];

advect[n_, d0_, u_, v_, w_, dt_] := 
  Module[{x, y, z, d1, dt0, i0, i1, j0, j1, k0, k1, s0, s1, t0, t1, 
    p1, p0, d000, d100, d010, d001, d110, d011, d101, d111}, 
   d1 = ConstantArray[0, {n + 1, n + 1, n + 1}]; dt0 = dt n;
   Do[Do[Do[x = i - dt0 u[[i, j, k]]; y = j - dt0 v[[i, j, k]]; 
        z = k - dt0 w[[i, j, k]];
        i0 = Which[x <= 1, 1, 1 < x < n, Floor[x], True, n];
        i1 = i0 + 1;
        j0 = Which[y <= 1, 1, 1 < y < n, Floor[y], True, n];
        j1 = j0 + 1; 
        k0 = Which[z <= 1, 1, 1 < z < n, Floor[z], True, n];
        k1 = k0 + 1; 
        d000 = (d0[[i0, j0, 
            k0]] + (x - i0) (d0[[i1, j0, k0]] - 
              d0[[i0, j0, k0]]) + (y - j0) (d0[[i0, j1, k0]] - 
              d0[[i0, j0, k0]]) + (z - k0) (d0[[i0, j0, k1]] - 
              d0[[i0, j0, k0]])); 
        d100 = (d0[[i1, j0, 
            k0]] + (x - i1) (d0[[i1, j0, k0]] - 
              d0[[i0, j0, k0]]) + (y - j0) (d0[[i1, j1, k0]] - 
              d0[[i1, j0, k0]]) + (z - k0) (d0[[i1, j0, k1]] - 
              d0[[i1, j0, k0]])); 
        d010 = (d0[[i0, j1, 
            k0]] + (x - i0) (d0[[i1, j1, k0]] - 
              d0[[i0, j1, k0]]) + (y - j1) (d0[[i0, j1, k0]] - 
              d0[[i0, j0, k0]]) + (z - k0) (d0[[i0, j1, k1]] - 
              d0[[i0, j1, k0]])); 
        d001 = (d0[[i0, j0, 
            k1]] + (x - i0) (d0[[i1, j0, k1]] - 
              d0[[i0, j0, k1]]) + (y - j0) (d0[[i0, j1, k1]] - 
              d0[[i0, j0, k1]]) + (z - k1) (d0[[i0, j0, k1]] - 
              d0[[i0, j0, k0]])); 
        d110 = (d0[[i1, j1, 
            k0]] + (x - i1) (d0[[i1, j1, k0]] - 
              d0[[i0, j1, k0]]) + (y - j1) (d0[[i1, j1, k0]] - 
              d0[[i1, j0, k0]]) + (z - k0) (d0[[i1, j1, k1]] - 
              d0[[i1, j1, k0]])); 
        d011 = (d0[[i0, j1, 
            k1]] + (x - i0) (d0[[i1, j1, k1]] - 
              d0[[i0, j1, k1]]) + (y - j1) (d0[[i0, j1, k1]] - 
              d0[[i0, j0, k1]]) + (z - k1) (d0[[i0, j1, k1]] - 
              d0[[i0, j1, k0]])); 
        d101 = (d0[[i1, j0, 
            k1]] + (x - i1) (d0[[i1, j0, k1]] - 
              d0[[i0, j0, k1]]) + (y - j0) (d0[[i1, j1, k1]] - 
              d0[[i1, j0, k1]]) + (z - k1) (d0[[i1, j0, k1]] - 
              d0[[i1, j0, k0]])); 
        d111 = (d0[[i1, j1, 
            k1]] + (x - i1) (d0[[i1, j1, k1]] - 
              d0[[i0, j1, k1]]) + (y - j1) (d0[[i1, j1, k1]] - 
              d0[[i1, j0, k1]]) + (z - k1) (d0[[i1, j1, k1]] - 
              d0[[i1, j1, k0]]));
        d1[[i, j, 
          k]] = (d000 + d100 + d010 + d001 + d110 + d011 + d101 + 
            d111)/8;, {k, 2, n}];, {j, 2, n}];, {i, 1, n + 1}]; d1];

project[n_, r_, u0_, v0_, w0_, u_, v_, w_] := 
  Module[{ux = u, vy = v, wz = w, div, p}, 
   p = ConstantArray[0, {n + 1, n + 1, n + 1}];
   div = ConstantArray[0, {n + 1, n + 1, n + 1}];
   ux = ConstantArray[0, {n + 1, n + 1, n + 1}];
   vy = ConstantArray[0, {n + 1, n + 1, n + 1}]; 
   wz = ConstantArray[0, {n + 1, n + 1, n + 1}];
   Do[div[[i, j, 
       k]] = .5/
        n (u0[[i + 1, j, k]] - u0[[i - 1, j, k]] + v0[[i, 1 + j, k]] -
          v0[[i, j - 1, k]] + w0[[i, j, k + 1]] - 
         w0[[i, j, k - 1]]);, {i, 2, n}, {j, 2, n}, {k, 2, n}]; 
   Do[div[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
     Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}]; 
   div = periodic[n, 0, 0, 0, 0, div];
   Do[Do[Do[
       Do[p[[i, j, 
            k]] = (div[[i, j, 
               k]] + (p[[i - 1, j, k]] + p[[i + 1, j, k]] + 
                p[[i, j - 1, k]] + p[[i, j + 1, k]] + 
                p[[i, j, k - 1]] + p[[i, j, k + 1]]))/6;, {k, 2, 
          n}];, {j, 2, n}], {i, 2, n}];, {k1, 0, r}]; 
   Do[p[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
     Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}]; 
   p = periodic[n, 0, 0, 0, 0, p];
   Do[ux[[i, j, k]] = 
     u0[[i, j, k]] + .5 n (p[[i + 1, j, k]] - p[[i - 1, j, k]]);
    vy[[i, j, k]] = 
     v0[[i, j, k]] + .5 n (p[[i, j + 1, k]] - p[[i, j - 1, k]]); 
    wz[[i, j, k]] = 
     w0[[i, j, k]] + .5 n (p[[i, j, k + 1]] - p[[i, j, k - 1]]);, {i, 
     2, n}, {j, 2, n}, {k, 2, n}]; 
   Do[ux[[i, j, k]] = 0; vy[[i, j, k]] = 0; 
    wz[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
     Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}]; {ux,
     vy, wz}];

onestep[n_, step_, r_, a_, uin_, vin_, win_, dt_, c_] := 
 Module[{u1, v1, w1, u, v, w, u0, v0, w0},
  u0 = ConstantArray[0., {n + 1, n + 1, n + 1}];
  v0 = ConstantArray[0., {n + 1, n + 1, n + 1}]; 
  w0 = ConstantArray[0., {n + 1, n + 1, n + 1}];
  u = ConstantArray[0., {n + 1, n + 1, n + 1}];
  v = ConstantArray[0., {n + 1, n + 1, n + 1}]; 
  w = ConstantArray[0., {n + 1, n + 1, n + 1}];
  u1 = ConstantArray[0., {n + 1, n + 1, n + 1}];
  v1 = ConstantArray[0., {n + 1, n + 1, n + 1}]; 
  w1 = ConstantArray[0., {n + 1, n + 1, n + 1}]; u0 = uin; v0 = vin; 
  w0 = win; 
  Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0; 
   w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
    Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
  u0 = advect[n, u0, u0, v0, w0, dt]; 
  v0 = advect[n, v0, u0, v0, w0, dt]; 
  w0 = advect[n, w0, u0, v0, w0, dt]; 
  Do[u0[[i, j, k]] = 0; 
   v0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
    Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
  u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0]; 
  w0 = periodic[n, 0, 0, 0, 0, w0];
  u0 = diffuse[n, r, a, c, u0]; v0 = diffuse[n, r, a, c, v0]; 
  w0 = diffuse[n, r, a, c, w0]; 
  Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0; 
   w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
    Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
  u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0]; 
  w0 = periodic[n, 0, 0, 0, 0, w0];
  {u0, v0, w0} = project[n, r, u0, v0, w0, u, v, w]; 
  Do[u0[[i, j, k]] = 0; v0[[i, j, k]] = 0; 
   w0[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, Round[n/2 - 5], 
    Round[n/2 + 5], 1}, {k, Round[n/2 - 5], Round[n/2 + 5], 1}];
  u0 = periodic[n, 0, 0, 0, 0, u0]; v0 = periodic[n, 0, 0, 0, 0, v0]; 
  w0 = periodic[n, 0, 0, 0, 0, w0]; {u0, v0, w0}]

cf = With[{cg = Compile`GetElement, hp = HoldPattern, 
     dv = DownValues}, 
    Hold@Compile[{{u0argu, _Real, 3}, {v0argu, _Real, 
              3}, {w0argu, _Real, 3}, {denargu, _Real, 
              3}, {sm, _Integer}, {n, _Integer}, {r, _Integer}, dif, 
             pec}, Module[{u0 = u0argu, v0 = v0argu, w0 = w0argu, uu, 
              vv, ww, dd, den = denargu, 
              c = Table[0., {n + 1}, {n + 1}, {n + 1}], dt = 40./n^2, 
              a, dnup = den[[1, n + 1, 1]], dnd = den[[1, 1, 1]], 
              dnl = den[[1, 1, 1]], dnr = den[[1, 1, n + 1]]}, 
             a = dt dif n n;
             
             uu = vv = 
               ww = dd = 
                 Table[0., {sm + 1}, {n + 1}, {n + 1}, {n + 1}];
             
             Do[{u0, v0, w0} = 
               onestep[n, step, r, a, u0, v0, w0, dt, c];
              uu[[step + 1]] = u0;
              vv[[step + 1]] = v0; ww[[step + 1]] = v0;
              den = diffuse[n, r, a/pec, c, den]; 
              Do[den[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, 
            
                   Round[n/2 - 5], Round[n/2 + 5], 1}, {k, 
                Round[n/2 - 5], Round[n/2 + 5], 1}];
              den = periodic[n, dnup, dnd, dnl, dnr, den];
              den = advect[n, den, u0, v0, w0, dt]; 
              Do[den[[i, j, k]] = 0;, {i, 10, 20, 1}, {j, 
                Round[n/2 - 5], Round[n/2 + 5], 1}, {k, 
                Round[n/2 - 5], Round[n/2 + 5], 1}];
              den = periodic[n, dnup, dnd, dnl, dnr, den];
              dd[[step + 1]] = den;, {step, 0, sm}]; {uu, vv, ww, 
              dd}], CompilationTarget -> "C", 
            RuntimeOptions -> "Speed"] /. dv@onestep /. 
         Flatten[dv /@ {advect, diffuse, periodic, project}] /. 
        hp@ConstantArray[c_, {i_, j_, kc_}] :> 
         Table[0., {i}, {j}, {kc}] /. hp@Part[a__] :> cg[a] /. 
      hp[cg[a__] = rhs_] :> (Part[a] = rhs) // 
     ReleaseHold]; 



rst = cf[u0, v0, w0, den, sm, n, r, dif, pec]; 

Visualization

Do[lst11[s] = 
  Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[1, s, i, j, k]]}, {i,
     n1}, {j, n1}, {k, n1}];
 lst12[s] = 
  Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[2, s, i, j, k]]}, {i,
     n1}, {j, n1}, {k, n1}]; 
 lst13[s] = 
  Table[{{(i - 1)/n, (j - 1)/n, (k - 1)/n}, rst[[3, s, i, j, k]]}, {i,
     n1}, {j, n1}, {k, n1}];, {s, 25, sm, 25}]

Do[su1[i] = 
  Interpolation[Flatten[lst11[i], 2], InterpolationOrder -> 3]; 
 sv2[i] = Interpolation[Flatten[lst12[i], 2], 
   InterpolationOrder -> 3]; 
 sw3[i] = Interpolation[Flatten[lst13[i], 2], 
   InterpolationOrder -> 3];, {i, 25, sm, 25}]

Table[Show[
  DensityPlot3D[
   Norm[{su1[i][x, y, z], sv2[i][x, y, z], sw3[i][x, y, z]}], {x, 0, 
    1}, {y, 0.44, 1}, {z, 0, 1}, BoxRatios -> Automatic, 
   ColorFunction -> "BlueGreenYellow", PlotLabel -> i], 
  VectorPlot3D[{su1[i][x, y, z], sv2[i][x, y, z], 
    sw3[i][x, y, z]}, {x, 0, 1}, {y, 0, .5}, {z, 0, 1}, 
   VectorPoints -> Fine, VectorMarkers -> "Arrow"], 
  Graphics3D[{{Blue, 
     Cuboid[{10, Round[n/2 - 5] - 1/2, 
        Round[n/2 - 5] - 1/2}/(n + 1), {20, Round[n/2 + 5] - 1/2, 
        Round[n/2 + 5] - 1/2}/(n + 1)]}}]], {i, 50, sm, 50}]

Figure 1

Animation

Do[lstu1[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, 
      rst[[1, k, i, Round[(n + 1)/2], j]]}, {i, n1}, {j, n1}], 1]; 
  lstw1[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, 
      rst[[2, k, i, Round[(n + 1)/2], j]]}, {i, n1}, {j, n1}], 
    1];, {k, sm}];
Do[Uvel1[i] = Interpolation[lstu1[i], InterpolationOrder -> 3];, {i, 
   1, sm}];
Do[Wvel1[i] = Interpolation[lstw1[i], InterpolationOrder -> 3];, {i, 
   1, sm}];
frame = Table[
   Show[DensityPlot[
     Norm[{Uvel[m][x, y], Vvel[m][x, y]}], {x, 0, 1}, {y, 0, 1}, 
     PlotRange -> All, ColorFunction -> "BlueGreenYellow", 
     Frame -> False, ImageSize -> Tiny, PlotLabel -> m, 
     PlotPoints -> 50], 
    Graphics[{Gray, 
      Rectangle[{10, Round[n/2 - 5] - 1/2}/(n + 1), {20, 
         Round[n/2 + 5] - 1/2}/(n + 1)]}]], {m, 5, sm, 3}];

Figure 2

The question is about code improvement. How can we define parameter r to solve Laplace and Poison equations in this code? Note, r is number of iterations in diffuse and project module (we use Gauss-Seidel relaxation to solve Laplace and Poison's equations).

Update 1. We can reduce r from r=20 as above to r=5 as below. Globally there is no difference in two animations.

Figure 3

But if we compare velocity in some point like $x=0.65, y=0.5$ then we have good agreement for $r=5$ (red points) and $r=20$ (gray points) in laminar flow at $t<0.5$, but a big difference in the turbulent flow at $t>0.5$. Figure 4

$\endgroup$
16
  • 1
    $\begingroup$ What is the meaning of r? Reynolds number? I think it is not a good idea to name a parameter r because it is nearly impossible to find it in the code. How do you want to improve your code: in readability (sure you can), speed, memory use? Your code looks like Fortran in Mathematica. Why? You may gain a lot by using Mma expressive power. $\endgroup$ Dec 25 '21 at 15:06
  • $\begingroup$ @PierreALBARÈDE r is number of iterations in diffuse and project module (we use Gauss-Seidel relaxation to solve Laplace and Poison's equations). Actually this code optimized for C and basically it coming from C code published in this paper dgp.toronto.edu/public_user/stam/reality/Research/pdf/ns.pdf $\endgroup$ Dec 25 '21 at 15:15
  • $\begingroup$ Don't you have some convergence criterion on this r? How did you choose r=20? $\endgroup$ Dec 25 '21 at 15:28
  • $\begingroup$ Er… so you mean you're looking for a systematic way to determine the proper value of r? $\endgroup$
    – xzczd
    Dec 25 '21 at 15:29
  • $\begingroup$ @PierreALBARÈDE For real time simulations in NVIDIA applications they usually use $r \ge 10$ for 2D and $r \ge 4$ for 3D flow. My choose r=20 is based on several tests and comparisons with other solvers. But, I think, that r=20 is too high for this kind of problems. $\endgroup$ Dec 25 '21 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.