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;
    }