Diffusion

Wall-normal derivative

\[\dmomdif{1}{2} = \mu \frac{1}{J} \dif{}{\gcs{1}} \left( \frac{J}{\sfact{1}} \vgt{1}{2} \right) + \mu \frac{1}{J} \dif{}{\gcs{1}} \left( \frac{J}{\sfact{1}} \vgt{2}{1} \right)\]

Diagonal part

\[\mu \frac{1}{J} \dif{}{\gcs{1}} \left( \frac{J}{\sfact{1}} \vgt{1}{2} \right) = \mu \frac{1}{J} \dif{}{\gcs{1}} \left( \frac{J}{\sfact{1}} \dlxy \right)\]
src/fluid/predict/uy.c
49for(size_t i = 1; i <= isize; i++){
50  const double l = 1. / JDXC(i  ) * JDXF(i  ) / HXXF(i  ) / HXXF(i  );
51  const double u = 1. / JDXC(i  ) * JDXF(i+1) / HXXF(i+1) / HXXF(i+1);
52  const double c = - l - u;
53  laplacians.LAPX(i).l = l;
54  laplacians.LAPX(i).c = c;
55  laplacians.LAPX(i).u = u;
56}
src/fluid/predict/uy.c
298  BEGIN
299    src[cnt] += diffusivity * (
300#if NDIMS == 2
301        + LAPX(i).l * UY(i-1, j  )
302        + LAPX(i).c * UY(i  , j  )
303        + LAPX(i).u * UY(i+1, j  )
304#else
305        + LAPX(i).l * UY(i-1, j  , k  )
306        + LAPX(i).c * UY(i  , j  , k  )
307        + LAPX(i).u * UY(i+1, j  , k  )
308#endif
309    );
310  END

Non-diagonal part

\[\mu \frac{1}{J} \dif{}{\gcs{1}} \left( \frac{J}{\sfact{1}} \vgt{2}{1} \right)\]
src/fluid/predict/uy.c
330  BEGIN
331    const double hx_xm = HXXF(i  );
332    const double hx_xp = HXXF(i+1);
333    const double jd_xm = JDXF(i  );
334    const double jd_x0 = JDXC(i  );
335    const double jd_xp = JDXF(i+1);
336#if NDIMS == 2
337    const double lyx_xm = LYX0(i  , j  ) + LYX1(i  , j  );
338    const double lyx_xp = LYX0(i+1, j  ) + LYX1(i+1, j  );
339#else
340    const double lyx_xm = LYX0(i  , j  , k  ) + LYX1(i  , j  , k  );
341    const double lyx_xp = LYX0(i+1, j  , k  ) + LYX1(i+1, j  , k  );
342#endif
343    src[cnt] += diffusivity / jd_x0 * (
344        - jd_xm / hx_xm * lyx_xm
345        + jd_xp / hx_xp * lyx_xp
346    );
347  END

Stream-wise derivative

\[\dmomdif{2}{2} = 2 \mu \frac{1}{J} \dif{}{\gcs{2}} \left( \frac{J}{\sfact{2}} \vgt{2}{2} \right)\]

Diagonal part

\[2 \mu \frac{1}{J} \dif{}{\gcs{2}} \left( \frac{J}{\sfact{2}} \vgt{2}{2}^{\prime} \right) = 2 \mu \frac{1}{J} \dif{}{\gcs{2}} \left( \frac{J}{\sfact{2}} \frac{1}{\sfact{2}} \dif{\vel{2}}{\gcs{2}} \right)\]
src/fluid/predict/uy.c
64for(size_t i = 1; i <= isize; i++){
65  const double l = 1. / HYXC(i  ) / HYXC(i  );
66  const double u = 1. / HYXC(i  ) / HYXC(i  );
67  const double c = - l - u;
68  laplacians.LAPY(i).l = l;
69  laplacians.LAPY(i).c = c;
70  laplacians.LAPY(i).u = u;
71}
src/fluid/predict/uy.c
364  BEGIN
365    // NOTE: tyy = "2" lyy
366    src[cnt] += 2. * diffusivity * (
367#if NDIMS == 2
368        + LAPY(i).l * UY(i  , j-1)
369        + LAPY(i).c * UY(i  , j  )
370        + LAPY(i).u * UY(i  , j+1)
371#else
372        + LAPY(i).l * UY(i  , j-1, k  )
373        + LAPY(i).c * UY(i  , j  , k  )
374        + LAPY(i).u * UY(i  , j+1, k  )
375#endif
376    );
377  END

Non-diagonal part

\[2 \mu \frac{1}{J} \dif{}{\gcs{2}} \left( \frac{J}{\sfact{2}} \vgt{2}{2}^{\prime\prime} \right)\]
src/fluid/predict/uy.c
395  BEGIN
396    // NOTE: tyy = "2" lyy
397    const double hy = HYXC(i  );
398    const double jd = JDXC(i  );
399    src[cnt] += 2. * diffusivity / jd * (
400#if NDIMS == 2
401        - jd / hy * LYY1(i  , j-1)
402        + jd / hy * LYY1(i  , j  )
403#else
404        - jd / hy * LYY1(i  , j-1, k  )
405        + jd / hy * LYY1(i  , j  , k  )
406#endif
407    );
408  END

Span-wise derivative

\[\dmomdif{3}{2} = \mu \frac{1}{J} \dif{}{\gcs{3}} \left( \frac{J}{\sfact{3}} \vgt{3}{2} \right) + \mu \frac{1}{J} \dif{}{\gcs{3}} \left( \frac{J}{\sfact{3}} \vgt{2}{3} \right)\]

Diagonal part

\[\mu \frac{1}{J} \dif{}{\gcs{3}} \left( \frac{J}{\sfact{3}} \dlzy \right)\]
src/fluid/predict/uy.c
78const double l = 1. / hz / hz;
79const double u = 1. / hz / hz;
80const double c = - l - u;
81laplacians.lapz.l = l;
82laplacians.lapz.c = c;
83laplacians.lapz.u = u;
src/fluid/predict/uy.c
424BEGIN
425  src[cnt] += diffusivity * (
426      + (*lapz).l * UY(i  , j  , k-1)
427      + (*lapz).c * UY(i  , j  , k  )
428      + (*lapz).u * UY(i  , j  , k+1)
429  );
430END

Non-diagonal part

\[\mu \frac{1}{J} \dif{}{\gcs{3}} \left( \frac{J}{\sfact{3}} \vgt{2}{3} \right)\]
src/fluid/predict/uy.c
448BEGIN
449  const double jd = JDXC(i  );
450  const double lyz_zm = LYZ(i  , j  , k  );
451  const double lyz_zp = LYZ(i  , j  , k+1);
452  src[cnt] += diffusivity / jd * (
453      - jd / hz * lyz_zm
454      + jd / hz * lyz_zp
455  );
456END

Additional term

\[\dmomdify = \mu \frac{1}{J} \ave{ \dif{ \left( \frac{J}{\sfact{1}} \right) }{\gcs{1}} \left( \vgt{2}{1} + \vgt{1}{2} \right) }{\gcs{1}}\]
src/fluid/predict/uy.c
477  BEGIN
478    const double hx_xm = HXXC(i-1);
479    const double hx_x0 = HXXC(i  );
480    const double hx_xp = HXXC(i+1);
481    const double jd_xm = JDXC(i-1);
482    const double jd_x0 = JDXC(i  );
483    const double jd_xp = JDXC(i+1);
484    const double djdhx_xm = - jd_xm / hx_xm + jd_x0 / hx_x0;
485    const double djdhx_xp = - jd_x0 / hx_x0 + jd_xp / hx_xp;
486#if NDIMS == 2
487    const double lyx_xm = LYX0(i  , j  ) + LYX1(i  , j  );
488    const double lyx_xp = LYX0(i+1, j  ) + LYX1(i+1, j  );
489    const double lxy_xm = LXY(i  , j  );
490    const double lxy_xp = LXY(i+1, j  );
491#else
492    const double lyx_xm = LYX0(i  , j  , k  ) + LYX1(i  , j  , k  );
493    const double lyx_xp = LYX0(i+1, j  , k  ) + LYX1(i+1, j  , k  );
494    const double lxy_xm = LXY(i  , j  , k  );
495    const double lxy_xp = LXY(i+1, j  , k  );
496#endif
497    src[cnt] += diffusivity / jd_x0 * (
498        + 0.5 * djdhx_xm * (lyx_xm + lxy_xm)
499        + 0.5 * djdhx_xp * (lyx_xp + lxy_xp)
500    );
501  END