Correction step

Updated but still non-solenoidal velocity field \(u_i^*\) is corrected to be divergence-free by solving a Poisson equation with respect to the scalar potential \(\phi\):

\[\frac{1}{J} \dif{}{\gcs{1}} \left( \frac{J}{\sfact{1}} \frac{1}{\sfact{1}} \dif{\phi}{\gcs{1}} \right) + \frac{1}{J} \dif{}{\gcs{2}} \left( \frac{J}{\sfact{2}} \frac{1}{\sfact{2}} \dif{\phi}{\gcs{2}} \right) + \frac{1}{J} \dif{}{\gcs{3}} \left( \frac{J}{\sfact{3}} \frac{1}{\sfact{3}} \dif{\phi}{\gcs{3}} \right) = \frac{1}{\Delta t} \dif{\vel{i}^*}{\gcs{i}}.\]

The eigenvalues are computed as follows:

src/fluid/compute_potential.c

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

src/fluid/compute_potential.c

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

In the radial direction, tri-diagonal matrices are to be solved for each \(y\) and \(z\):

src/fluid/compute_potential.c

// set center diagonal components
for(size_t i = 1; i <= isize; i++){
  tdm_c[i-1] =
    - tdm_l[i-1]
    - tdm_u[i-1]
    + eval_y / HYXC(i  ) / HYXC(i  );
}
// 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);

src/fluid/compute_potential.c

// set center diagonal components
for(size_t i = 1; i <= isize; i++){
  tdm_c[i-1] =
    - tdm_l[i-1]
    - tdm_u[i-1]
    + eval_y / HYXC(i  ) / HYXC(i  )
    + 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);

After computing the scalar potential \(\phi\), the velocity field is corrected using its gradient:

\[\pder{\phi}{x_i} \approx \vec{e_{\gcs{1}}} \frac{1}{\sfact{1}} \dif{\phi}{\gcs{1}} + \vec{e_{\gcs{2}}} \frac{1}{\sfact{2}} \dif{\phi}{\gcs{2}} + \vec{e_{\gcs{3}}} \frac{1}{\sfact{3}} \dif{\phi}{\gcs{3}}.\]

src/fluid/correct/ux.c

  BEGIN
    const double hx = HXXF(i  );
#if NDIMS == 2
    const double psi_xm = PSI(i-1, j  );
    const double psi_xp = PSI(i  , j  );
    double * vel = &UX(i, j);
#else
    const double psi_xm = PSI(i-1, j  , k  );
    const double psi_xp = PSI(i  , j  , k  );
    double * vel = &UX(i, j, k);
#endif
    *vel -= prefactor / hx * (
        - psi_xm
        + psi_xp
    );
  END

src/fluid/correct/uy.c

  BEGIN
    const double hy = HYXC(i  );
#if NDIMS == 2
    const double psi_ym = PSI(i  , j-1);
    const double psi_yp = PSI(i  , j  );
    double * vel = &UY(i, j);
#else
    const double psi_ym = PSI(i  , j-1, k  );
    const double psi_yp = PSI(i  , j  , k  );
    double * vel = &UY(i, j, k);
#endif
    *vel -= prefactor / hy * (
        - psi_ym
        + psi_yp
    );
  END

src/fluid/correct/uz.c

for(int k = 1; k <= ksize; k++){
  for(int j = 1; j <= jsize; j++){
    for(int i = 1; i <= isize; i++){
      const double psi_zm = PSI(i  , j  , k-1);
      const double psi_zp = PSI(i  , j  , k  );
      UZ(i, j, k) -= prefactor / hz * (
          - psi_zm
          + psi_zp
      );
    }
  }
}