Diffusion

Wall-normal derivative

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

src/fluid/predict/ux.c

for(size_t i = 2; i <= isize; i++){
  const double l = 1. / JDXF(i  ) * JDXC(i-1) / HXXC(i-1) / HXXC(i-1);
  const double u = 1. / JDXF(i  ) * JDXC(i  ) / HXXC(i  ) / HXXC(i  );
  const double c = - l - u;
  laplacians.LAPX(i).l = l;
  laplacians.LAPX(i).c = c;
  laplacians.LAPX(i).u = u;
}

src/fluid/predict/ux.c

  BEGIN
    // NOTE: txx = "2" lxx
    src[cnt] += 2. * diffusivity * (
#if NDIMS == 2
        + LAPX(i).l * UX(i-1, j  )
        + LAPX(i).c * UX(i  , j  )
        + LAPX(i).u * UX(i+1, j  )
#else
        + LAPX(i).l * UX(i-1, j  , k  )
        + LAPX(i).c * UX(i  , j  , k  )
        + LAPX(i).u * UX(i+1, j  , k  )
#endif
    );
  END

Stream-wise derivative

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

Diagonal part

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

src/fluid/predict/ux.c

for(size_t i = 2; i <= isize; i++){
  const double l = 1. / HYXF(i  ) / HYXF(i  );
  const double u = 1. / HYXF(i  ) / HYXF(i  );
  const double c = - l - u;
  laplacians.LAPY(i).l = l;
  laplacians.LAPY(i).c = c;
  laplacians.LAPY(i).u = u;
}

src/fluid/predict/ux.c

  BEGIN
    src[cnt] += diffusivity * (
#if NDIMS == 2
        + LAPY(i).l * UX(i  , j-1)
        + LAPY(i).c * UX(i  , j  )
        + LAPY(i).u * UX(i  , j+1)
#else
        + LAPY(i).l * UX(i  , j-1, k  )
        + LAPY(i).c * UX(i  , j  , k  )
        + LAPY(i).u * UX(i  , j+1, k  )
#endif
    );
  END

Non-diagonal part

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

src/fluid/predict/ux.c

  BEGIN
    const double hy = HYXF(i  );
    const double jd = JDXF(i  );
#if NDIMS == 2
    const double lyx_ym = LYX1(i  , j  );
    const double lyx_yp = LYX1(i  , j+1);
    const double lxy_ym = LXY(i  , j  );
    const double lxy_yp = LXY(i  , j+1);
#else
    const double lyx_ym = LYX1(i  , j  , k  );
    const double lyx_yp = LYX1(i  , j+1, k  );
    const double lxy_ym = LXY(i  , j  , k  );
    const double lxy_yp = LXY(i  , j+1, k  );
#endif
    src[cnt] += diffusivity / jd * (
        - jd / hy * (lyx_ym + lxy_ym)
        + jd / hy * (lyx_yp + lxy_yp)
    );
  END

Span-wise derivative

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

Diagonal part

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

src/fluid/predict/ux.c

const double l = 1. / hz / hz;
const double u = 1. / hz / hz;
const double c = - l - u;
laplacians.lapz.l = l;
laplacians.lapz.c = c;
laplacians.lapz.u = u;

src/fluid/predict/ux.c

BEGIN
  src[cnt] += diffusivity * (
      + (*lapz).l * UX(i  , j  , k-1)
      + (*lapz).c * UX(i  , j  , k  )
      + (*lapz).u * UX(i  , j  , k+1)
  );
END

Non-diagonal part

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

src/fluid/predict/ux.c

BEGIN
  const double jd = JDXF(i  );
  const double lxz_zm = LXZ(i  , j  , k  );
  const double lxz_zp = LXZ(i  , j  , k+1);
  src[cnt] += diffusivity / jd * (
      - jd / hz * lxz_zm
      + jd / hz * lxz_zp
  );
END

Additional term

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

src/fluid/predict/ux.c

  BEGIN
    const double hx_xm = HXXF(i-1);
    const double hx_x0 = HXXF(i  );
    const double hx_xp = HXXF(i+1);
    const double jd_xm = JDXF(i-1);
    const double jd_x0 = JDXF(i  );
    const double jd_xp = JDXF(i+1);
    const double djdhx_xm = - jd_xm / hx_xm + jd_x0 / hx_x0;
    const double djdhx_xp = - jd_x0 / hx_x0 + jd_xp / hx_xp;
#if NDIMS == 2
    const double lyy_xm = LYY0(i-1, j  ) + LYY1(i-1, j  );
    const double lyy_xp = LYY0(i  , j  ) + LYY1(i  , j  );
#else
    const double lyy_xm = LYY0(i-1, j  , k  ) + LYY1(i-1, j  , k  );
    const double lyy_xp = LYY0(i  , j  , k  ) + LYY1(i  , j  , k  );
#endif
    // NOTE: tyy = "2" lyy
    src[cnt] -= 2. * diffusivity / jd_x0 * (
        + 0.5 * djdhx_xm * lyy_xm
        + 0.5 * djdhx_xp * lyy_xp
    );
  END