Poisson Equation

To integrate a momentum field in time (specifically to enforce the incompressibility), we need to solve a Poisson equation:

\[\frac{1}{J} \dif{}{\gcs{1}} \left( \frac{J}{\sfact{1}} \frac{1}{\sfact{1}} \dif{p}{\gcs{1}} \right) + \frac{1}{J} \dif{}{\gcs{2}} \left( \frac{J}{\sfact{2}} \frac{1}{\sfact{2}} \dif{p}{\gcs{2}} \right) + \frac{1}{J} \dif{}{\gcs{3}} \left( \frac{J}{\sfact{3}} \frac{1}{\sfact{3}} \dif{p}{\gcs{3}} \right) = q,\]

where \(p\) and \(q\) are both defined at cell centers \(\left( \ccidx{i}, \ccidx{j}, \ccidx{k} \right)\).

Since directions except \(x\) are homogeneous (i.e., periodic boundary conditions are imposed and the grid sizes are equal), we consider the spectral representation (discrete backward Fourier transform) of \(p\):

\[\newcommand{ \wavy}{\exp \left( I \frac{2 \pi}{\ngp{2}} j m \right)} \newcommand{ \wavz}{\exp \left( I \frac{2 \pi}{\ngp{3}} k n \right)} \newcommand{\iwavy}{\exp \left( - I \frac{2 \pi}{\ngp{2}} j m^\prime \right)} \newcommand{\iwavz}{\exp \left( - I \frac{2 \pi}{\ngp{3}} k n^\prime \right)} \vat{p}{\ccidx{i}, \ccidx{j}, \ccidx{k}} = \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \wavy \wavz,\]

where \(I\) is the imaginary unit \(\sqrt{-1}\). Note that the discrete forward Fourier transform yields

\[ \begin{align}\begin{aligned}\sum_{k = 0}^{\ngp{3} - 1} \sum_{j = 0}^{\ngp{2} - 1} \vat{p}{\ccidx{i}, \ccidx{j}, \ccidx{k}} \iwavy \iwavz & = \sum_{k = 0}^{\ngp{3} - 1} \sum_{j = 0}^{\ngp{2} - 1} \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \wavy \wavz \iwavy \iwavz\\& = \sum_{k = 0}^{\ngp{3} - 1} \sum_{j = 0}^{\ngp{2} - 1} \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \delta_{m m^\prime} \delta_{n n^\prime}\\& = \ngp{2} \ngp{3} \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \delta_{m m^\prime} \delta_{n n^\prime}\\& = \ngp{2} \ngp{3} \vat{P}{\ccidx{i}, \ccidx{m^\prime}, \ccidx{n^\prime}}.\end{aligned}\end{align} \]

Since we have

\[ \begin{align}\begin{aligned}\vat{p}{\ccidx{i} \pm 1, \ccidx{j}, \ccidx{k}} & = \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i} \pm 1, \ccidx{m}, \ccidx{n}} \wavy \wavz,\\\vat{p}{\ccidx{i}, \ccidx{j} \pm 1, \ccidx{k}} & = \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \exp \left( \pm I \frac{2 \pi}{\ngp{2}} m \right) \wavy \wavz,\\\vat{p}{\ccidx{i}, \ccidx{j}, \ccidx{k} \pm 1} & = \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \exp \left( \pm I \frac{2 \pi}{\ngp{3}} n \right) \wavy \wavz,\end{aligned}\end{align} \]

the second-order derivatives in the Poisson equation can be reformulated as

\[ \begin{align}\begin{aligned}\newcommand{\coefl}{C_l} \newcommand{\coefu}{C_u} \frac{1}{J} \dif{}{\gcs{1}} \left( \frac{J}{\sfact{1}} \frac{1}{\sfact{1}} \dif{p}{\gcs{1}} \right) & = \coefl \vat{p}{\ccidx{i} - 1, \ccidx{j}, \ccidx{k}} - \left( \coefl + \coefu \right) \vat{p}{\ccidx{i}, \ccidx{j}, \ccidx{k}} + \coefu \vat{p}{\ccidx{i} + 1, \ccidx{j}, \ccidx{k}}\\& = \coefl \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i} - 1, \ccidx{m}, \ccidx{n}} \wavy \wavz - \left( \coefl + \coefu \right) \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \wavy \wavz + \coefu \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i} + 1, \ccidx{m}, \ccidx{n}} \wavy \wavz,\\\frac{1}{J} \dif{}{\gcs{2}} \left( \frac{J}{\sfact{2}} \frac{1}{\sfact{2}} \dif{p}{\gcs{2}} \right) & = \frac{1}{\sfact{2}} \frac{1}{\sfact{2}} \left( \vat{p}{\ccidx{i}, \ccidx{j} - 1, \ccidx{k}} - 2 \vat{p}{\ccidx{i}, \ccidx{j} , \ccidx{k}} + \vat{p}{\ccidx{i}, \ccidx{j} + 1, \ccidx{k}} \right)\\& = \frac{1}{\sfact{2}} \frac{1}{\sfact{2}} \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \wavy \wavz \left\{ \exp \left( I \frac{2 \pi}{\ngp{2}} m \right) - 2 + \exp \left( - I \frac{2 \pi}{\ngp{2}} m \right) \right\}\\& = - \frac{1}{\sfact{2}} \frac{1}{\sfact{2}} \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} 4 \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \wavy \wavz \sin^2 \left( \frac{\pi}{\ngp{2}} m \right),\\\frac{1}{J} \dif{}{\gcs{3}} \left( \frac{J}{\sfact{3}} \frac{1}{\sfact{3}} \dif{p}{\gcs{3}} \right) & = - \frac{1}{\sfact{3}} \frac{1}{\sfact{3}} \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} 4 \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \wavy \wavz \sin^2 \left( \frac{\pi}{\ngp{3}} n \right),\end{aligned}\end{align} \]

where

\[ \begin{align}\begin{aligned}\coefl & \equiv \frac{ 1 }{ \vat{J}{\ccidx{i}} } \frac{ \vat{J}{\cmidx{i}} }{ \vat{\sfact{1}}{\cmidx{i}} } \frac{ 1 }{ \vat{\sfact{1}}{\cmidx{i}} },\\\coefu & \equiv \frac{ 1 }{ \vat{J}{\ccidx{i}} } \frac{ \vat{J}{\cpidx{i}} }{ \vat{\sfact{1}}{\cpidx{i}} } \frac{ 1 }{ \vat{\sfact{1}}{\cpidx{i}} },\end{aligned}\end{align} \]

which are computed here:

src/fluid/compute_potential.c

double * tdm_l = NULL;
double * tdm_u = NULL;
tdm.get_l(*tdm_info, &tdm_l);
tdm.get_u(*tdm_info, &tdm_u);
const double * hxxf = domain->hxxf;
const double * jdxf = domain->jdxf;
const double * jdxc = domain->jdxc;
for (int i = 1; i <= (int)tdm_sizes[0]; i++) {
  // N.B. iterated from i = 1 to use macros
  tdm_l[i-1] = 1. / JDXC(i  ) * JDXF(i  ) / HXXF(i  ) / HXXF(i  );
  tdm_u[i-1] = 1. / JDXC(i  ) * JDXF(i+1) / HXXF(i+1) / HXXF(i+1);
}

Consequently, we obtain

\[ \begin{align}\begin{aligned}& \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \coefl \vat{P}{\ccidx{i} - 1, \ccidx{m}, \ccidx{n}} \wavy \wavz\\+ & \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \left\{ - \left( \coefl + \coefu \right) - 4 \frac{1}{\sfact{2}} \frac{1}{\sfact{2}} \sin^2 \left( \frac{\pi}{\ngp{2}} m \right) - 4 \frac{1}{\sfact{3}} \frac{1}{\sfact{3}} \sin^2 \left( \frac{\pi}{\ngp{3}} n \right) \right\} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}} \wavy \wavz\\+ & \sum_{n = 0}^{\ngp{3} - 1} \sum_{m = 0}^{\ngp{2} - 1} \coefu \vat{P}{\ccidx{i} + 1, \ccidx{m}, \ccidx{n}} \wavy \wavz\\= & \vat{q}{\ccidx{i}, \ccidx{j}, \ccidx{k}}.\end{aligned}\end{align} \]

By applying the discrete forward Fourier transform, we obtain

\[ \begin{align}\begin{aligned}& \coefl \vat{P}{\ccidx{i} - 1, \ccidx{m}, \ccidx{n}}\\+ & \left\{ - \left( \coefl + \coefu \right) - 4 \frac{1}{\sfact{2}} \frac{1}{\sfact{2}} \sin^2 \left( \frac{\pi}{\ngp{2}} m \right) - 4 \frac{1}{\sfact{3}} \frac{1}{\sfact{3}} \sin^2 \left( \frac{\pi}{\ngp{3}} n \right) \right\} \vat{P}{\ccidx{i}, \ccidx{m}, \ccidx{n}}\\+ & \coefu \vat{P}{\ccidx{i} + 1, \ccidx{m}, \ccidx{n}}\\= & \frac{1}{\ngp{2}} \frac{1}{\ngp{3}} \sum_{k = 0}^{\ngp{3} - 1} \sum_{j = 0}^{\ngp{2} - 1} \vat{q}{\ccidx{i}, \ccidx{j}, \ccidx{k}} \exp \left( - I \frac{2 \pi}{\ngp{2}} j m \right) \exp \left( - I \frac{2 \pi}{\ngp{3}} k n \right).\end{aligned}\end{align} \]

The resulting tri-diagonal matrix for each \(m\) and \(n\) is solved here:

src/fluid/compute_potential.c

#if NDIMS == 2
  for (size_t j = 0; j < jsize; j++) {
    const double eval_y = eval_ys[j];
    // set center diagonal components
    for (size_t i = 0; i < isize; i++) {
      tdm_c[i] =
        - tdm_l[i]
        - tdm_u[i]
        + eval_y / hy / hy;
    }
    // boundary treatment (Neumann boundary condition)
    tdm_c[        0] += tdm_l[        0];
    tdm_c[isize - 1] += tdm_u[isize - 1];
    // solve
    tdm.solve(tdm_info, rhs + j * isize);
  }
#else
  for (size_t k = 0; k < ksize; k++) {
    const double eval_z = eval_zs[k];
    for (size_t j = 0; j < jsize; j++) {
      const double eval_y = eval_ys[j];
      // set center diagonal components
      for (size_t i = 0; i < isize; i++) {
        tdm_c[i] =
          - tdm_l[i]
          - tdm_u[i]
          + eval_y / hy / hy
          + eval_z / hz / hz;
      }
      // boundary treatment (Neumann boundary condition)
      tdm_c[        0] += tdm_l[        0];
      tdm_c[isize - 1] += tdm_u[isize - 1];
      // solve
      tdm.solve(tdm_info, rhs + (k * jsize + j) * isize);
    }
  }
#endif

The sinusoidal functions (wave numbers) are computed here:

src/fluid/compute_potential.c

double * restrict * evals = &poisson_solver->eval_ys;
size_t mysize = 0;
size_t myoffs = 0;
sdecomp.get_pencil_mysize(info, pencil, SDECOMP_YDIR, c_gl_sizes[SDECOMP_YDIR], &mysize);
sdecomp.get_pencil_offset(info, pencil, SDECOMP_YDIR, c_gl_sizes[SDECOMP_YDIR], &myoffs);
*evals = memory_calloc(mysize, sizeof(double));
for (size_t local_n = 0; local_n < mysize; local_n++) {
  const size_t global_n = local_n + myoffs;
  (*evals)[local_n] =
    - 4. * pow(
      sin( g_pi * global_n / r_gl_sizes[SDECOMP_YDIR] ),
      2.
    );
}

src/fluid/compute_potential.c

double * restrict * evals = &poisson_solver->eval_zs;
size_t mysize = 0;
size_t myoffs = 0;
sdecomp.get_pencil_mysize(info, pencil, SDECOMP_ZDIR, c_gl_sizes[SDECOMP_ZDIR], &mysize);
sdecomp.get_pencil_offset(info, pencil, SDECOMP_ZDIR, c_gl_sizes[SDECOMP_ZDIR], &myoffs);
*evals = memory_calloc(mysize, sizeof(double));
for (size_t local_n = 0; local_n < mysize; local_n++) {
  const size_t global_n = local_n + myoffs;
  (*evals)[local_n] =
    - 4. * pow(
      sin( g_pi * global_n / r_gl_sizes[SDECOMP_ZDIR] ),
      2.
    );
}

The overall procedure can be found here:

src/fluid/compute_potential.c

int fluid_compute_potential (
    const domain_t * domain,
    const size_t rkstep,
    const double dt,
    fluid_t * fluid
) {
  static poisson_solver_t poisson_solver = {
    .is_initialised = false,
  };
  // initialise Poisson solver
  if (!poisson_solver.is_initialised) {
    if (0 != init_poisson_solver(domain, &poisson_solver)) {
      // failed to initialise Poisson solver
      return 1;
    }
  }
  // compute right-hand side of Poisson equation
  // assigned to buf0
  assign_input(domain, rkstep, dt, fluid, poisson_solver.buf0);
  // transpose real x1pencil to y1pencil
  // from buf0 to buf1
  sdecomp.transpose.execute(
      poisson_solver.r_transposer_x1_to_y1,
      poisson_solver.buf0,
      poisson_solver.buf1
  );
  // project y to wave space
  // f(x, y)    -> f(x, k_y)
  // f(x, y, z) -> f(x, k_y, z)
  // from buf1 to buf0
  fftw_execute(poisson_solver.fftw_plan_y[0]);
#if NDIMS == 2
  // transpose complex y1pencil to x1pencil
  // from buf0 to buf1
  sdecomp.transpose.execute(
      poisson_solver.c_transposer_y1_to_x1,
      poisson_solver.buf0,
      poisson_solver.buf1
  );
#else
  // transpose complex y1pencil to z1pencil
  // from buf0 to buf1
  sdecomp.transpose.execute(
      poisson_solver.c_transposer_y1_to_z1,
      poisson_solver.buf0,
      poisson_solver.buf1
  );
  // project z to wave space
  // f(x, k_y, z) -> f(x, k_y, k_z)
  // from buf1 to buf0
  fftw_execute(poisson_solver.fftw_plan_z[0]);
  // transpose complex z1pencil to x2pencil
  // from buf0 to buf1
  sdecomp.transpose.execute(
      poisson_solver.c_transposer_z1_to_x2,
      poisson_solver.buf0,
      poisson_solver.buf1
  );
#endif
  // solve linear systems in x direction
  solve_linear_systems(domain, &poisson_solver);
#if NDIMS == 2
  // transpose complex x1pencil to y1pencil
  // from buf1 to buf0
  sdecomp.transpose.execute(
      poisson_solver.c_transposer_x1_to_y1,
      poisson_solver.buf1,
      poisson_solver.buf0
  );
#else
  // transpose complex x2pencil to z1pencil
  // from buf1 to buf0
  sdecomp.transpose.execute(
      poisson_solver.c_transposer_x2_to_z1,
      poisson_solver.buf1,
      poisson_solver.buf0
  );
  // project z to physical space
  // f(x, k_y, k_z) -> f(x, k_y, z)
  // from buf0 to buf1
  fftw_execute(poisson_solver.fftw_plan_z[1]);
  // transpose complex z1pencil to y1pencil
  // from buf1 to buf0
  sdecomp.transpose.execute(
      poisson_solver.c_transposer_z1_to_y1,
      poisson_solver.buf1,
      poisson_solver.buf0
  );
#endif
  // project y to physical space
  // f(x, k_y)    -> f(x, y)
  // f(x, k_y, z) -> f(x, y, z)
  // from buf0 to buf1
  fftw_execute(poisson_solver.fftw_plan_y[1]);
  // transpose real y1pencil to x1pencil
  // from buf1 to buf0
  sdecomp.transpose.execute(
      poisson_solver.r_transposer_y1_to_x1,
      poisson_solver.buf1,
      poisson_solver.buf0
  );
  // copy result to the original buffer
  extract_output(domain, poisson_solver.buf0, fluid);
  return 0;
}