Internal Energy Balance

The equation of the internal energy is

\[\pder{T}{t} = A + D_x + D_y + D_z,\]

where \(A\) is the advective terms:

\[A \equiv - \dtempadv{1} - \dtempadv{2} - \dtempadv{3},\]

while \(D_i\) is the diffusive term involving spatial differentiation in the \(i\)-th direction:

\[D_i \equiv \dtempdif{i} \,\, (\text{no summation over}\,i).\]

See the spatial discretization.

The temporal discretization for each Runge-Kutta iteration leads to

\[ \begin{align}\begin{aligned}\Delta T & = \alpha^k \Delta t \left( A^{k } + D^{k } \right) + \beta^k \Delta t \left( A^{k-1} + D^{k-1} \right),\\T^{k+1} & = T^k + \Delta T,\end{aligned}\end{align} \]

when all diffusive terms are treated explicitly, while

\[ \begin{align}\begin{aligned}\newcommand{\lap}[2]{ {#2} \gamma^k \Delta t \frac{1}{\sqrt{Pr} \sqrt{Ra}} \frac{1}{J} \dif{}{\gcs{#1}} \left( \frac{J}{\sfact{#1}} \frac{1}{\sfact{#1}} \dif{}{\gcs{#1}} \right) } \Delta T & = \alpha^k \Delta t A^{k } + \beta^k \Delta t A^{k-1} + \gamma^k \Delta t D^{k },\\T^{k+1} & = T^k + \left\{ 1 - \lap{3}{c} \right\}^{-1} \left\{ 1 - \lap{2}{c} \right\}^{-1} \left\{ 1 - \lap{1}{c} \right\}^{-1} \Delta T,\end{aligned}\end{align} \]

when all diffusive terms are treated implicitly.

Diffusive terms are sometimes partially treated implicitly (c.f., van der Poel et al., Comput. Fluids (116), 2015). When only the diffusive term in \(x\) direction is implicitly treated (the default configuration in this project), we have

\[ \begin{align}\begin{aligned}\Delta T & = \alpha^k \Delta t \left( A^{k } + D_2^{k } + D_3^{k } \right) + \beta^k \Delta t \left( A^{k-1} + D_2^{k-1} + D_3^{k-1} \right) + \gamma^k \Delta t D_1^{k },\\T^{k+1} & = T^k + \left\{ 1 - \lap{1}{c} \right\}^{-1} \Delta T.\end{aligned}\end{align} \]

First, the explicit and implicit terms are calculated and stored to the corresponding buffers:

src/fluid/predict/t.c

int compute_rhs_t (
    const domain_t * domain,
    fluid_t * fluid
) {
  if (!laplacians.is_initialised) {
    if (0 != init_laplacians(domain)) {
      return 1;
    }
  }
  const double * restrict ux = fluid->ux.data;
  const double * restrict uy = fluid->uy.data;
#if NDIMS == 3
  const double * restrict uz = fluid->uz.data;
#endif
  const double * restrict  t = fluid-> t.data;
  // buffer for explicit terms
  double * restrict srca = fluid->srct.alpha.data;
  // buffer for implicit terms
  double * restrict srcg = fluid->srct.gamma.data;
  const double diffusivity = fluid_compute_temperature_diffusivity(fluid);
  // advective contributions, always explicit
  advection_x(domain, t, ux, srca);
  advection_y(domain, t, uy, srca);
#if NDIMS == 3
  advection_z(domain, t, uz, srca);
#endif
  // diffusive contributions, can be explicit or implicit
  diffusion_x(domain, diffusivity, t, param_t_implicit_x ? srcg : srca);
  diffusion_y(domain, diffusivity, t, param_t_implicit_y ? srcg : srca);
#if NDIMS == 3
  diffusion_z(domain, diffusivity, t, param_t_implicit_z ? srcg : srca);
#endif
  return 0;
}

The buffers are used to compute \(\Delta T\):

src/fluid/predict/t.c

    // Runge-Kutta coefficients, alpha, beta, gamma
    const double coef_a = rkcoefs[rkstep].alpha;
    const double coef_b = rkcoefs[rkstep].beta ;
    const double coef_g = rkcoefs[rkstep].gamma;
    // Runge-Kutta buffers, alpha, beta, gamma
    const double * restrict srcta = fluid->srct.alpha.data;
    const double * restrict srctb = fluid->srct.beta .data;
    const double * restrict srctg = fluid->srct.gamma.data;
    const int isize = domain->mysizes[0];
    const int jsize = domain->mysizes[1];
#if NDIMS == 3
    const int ksize = domain->mysizes[2];
#endif
    double * restrict dtemp = linear_system.x1pncl;
#if NDIMS == 2
    const size_t nitems = isize * jsize;
#else
    const size_t nitems = isize * jsize * ksize;
#endif
    // compute T(new) - T(old)
    for (size_t n = 0; n < nitems; n++) {
      dtemp[n] =
        + coef_a * dt * srcta[n]
        + coef_b * dt * srctb[n]
        + coef_g * dt * srctg[n];
    }

When necessary, linear systems are solved to take care of the implicit treatments:

src/fluid/predict/t.c

if (param_t_implicit_x) {
  // prepare tri-diagonal coefficients
  tdm_info_t * tdm_info = linear_system.tdm_x;
  int size = 0;
  double * restrict tdm_l = NULL;
  double * restrict tdm_c = NULL;
  double * restrict tdm_u = NULL;
  tdm.get_size(tdm_info, &size);
  tdm.get_l(tdm_info, &tdm_l);
  tdm.get_c(tdm_info, &tdm_c);
  tdm.get_u(tdm_info, &tdm_u);
  const laplacian_t * lapx = laplacians.lapx;
  for (int i = 0; i < size; i++) {
    tdm_l[i] =    - prefactor * lapx[i].l;
    tdm_c[i] = 1. - prefactor * lapx[i].c;
    tdm_u[i] =    - prefactor * lapx[i].u;
  }
  // solve all
  tdm.solve(tdm_info, linear_system.x1pncl);
}

src/fluid/predict/t.c

if (param_t_implicit_y) {
  // make all y data accessible
  sdecomp.transpose.execute(
      linear_system.transposer_x1_to_y1,
      linear_system.x1pncl,
      linear_system.y1pncl
  );
  // prepare tri-diagonal coefficients
  tdm_info_t * tdm_info = linear_system.tdm_y;
  int size = 0;
  double * restrict tdm_l = NULL;
  double * restrict tdm_c = NULL;
  double * restrict tdm_u = NULL;
  tdm.get_size(tdm_info, &size);
  tdm.get_l(tdm_info, &tdm_l);
  tdm.get_c(tdm_info, &tdm_c);
  tdm.get_u(tdm_info, &tdm_u);
  const laplacian_t * lapy = &laplacians.lapy;
  for (int j = 0; j < size; j++) {
    tdm_l[j] =    - prefactor * (*lapy).l;
    tdm_c[j] = 1. - prefactor * (*lapy).c;
    tdm_u[j] =    - prefactor * (*lapy).u;
  }
  // solve all
  tdm.solve(tdm_info, linear_system.y1pncl);
  // recover original memory alignment
  sdecomp.transpose.execute(
      linear_system.transposer_y1_to_x1,
      linear_system.y1pncl,
      linear_system.x1pncl
  );
}

src/fluid/predict/t.c

if (param_t_implicit_z) {
  // make all z data accessible
  sdecomp.transpose.execute(
      linear_system.transposer_x1_to_z2,
      linear_system.x1pncl,
      linear_system.z2pncl
  );
  // prepare tri-diagonal coefficients
  tdm_info_t * tdm_info = linear_system.tdm_z;
  int size = 0;
  double * restrict tdm_l = NULL;
  double * restrict tdm_c = NULL;
  double * restrict tdm_u = NULL;
  tdm.get_size(tdm_info, &size);
  tdm.get_l(tdm_info, &tdm_l);
  tdm.get_c(tdm_info, &tdm_c);
  tdm.get_u(tdm_info, &tdm_u);
  const laplacian_t * lapz = &laplacians.lapz;
  for (int k = 0; k < size; k++) {
    tdm_l[k] =    - prefactor * (*lapz).l;
    tdm_c[k] = 1. - prefactor * (*lapz).c;
    tdm_u[k] =    - prefactor * (*lapz).u;
  }
  // solve all
  tdm.solve(tdm_info, linear_system.z2pncl);
  // recover original memory alignment
  sdecomp.transpose.execute(
      linear_system.transposer_z2_to_x1,
      linear_system.z2pncl,
      linear_system.x1pncl
  );
}

Finally the temperature field is updated:

src/fluid/predict/t.c

    const int isize = domain->mysizes[0];
    const int jsize = domain->mysizes[1];
#if NDIMS == 3
    const int ksize = domain->mysizes[2];
#endif
    const double * restrict dtemp = linear_system.x1pncl;
    double * restrict t = fluid->t.data;
    BEGIN
#if NDIMS == 2
      T(i, j) += dtemp[cnt];
#else
      T(i, j, k) += dtemp[cnt];
#endif
    END
    if (0 != fluid_update_boundaries_t(domain, &fluid->t)) {
      return 1;
    }