Incompressibility constraint

The incompressibility constraint is defined at each cell center:

\[\ddiv.\]

The (left-hand-side of) incompressibility constraint is computed to find the right-hand-side term of the Poisson equation:

src/fluid/compute_potential.c
487  BEGIN
488    const double hx_xm = HXXF(i  );
489    const double hx_xp = HXXF(i+1);
490    const double hy    = HYXC(i  );
491    const double jd_xm = JDXF(i  );
492    const double jd_x0 = JDXC(i  );
493    const double jd_xp = JDXF(i+1);
494#if NDIMS == 2
495    const double ux_xm = UX(i  , j  );
496    const double ux_xp = UX(i+1, j  );
497    const double uy_ym = UY(i  , j  );
498    const double uy_yp = UY(i  , j+1);
499#else
500    const double ux_xm = UX(i  , j  , k  );
501    const double ux_xp = UX(i+1, j  , k  );
502    const double uy_ym = UY(i  , j  , k  );
503    const double uy_yp = UY(i  , j+1, k  );
504    const double uz_zm = UZ(i  , j  , k  );
505    const double uz_zp = UZ(i  , j  , k+1);
506#endif
507    const double div = 1. / jd_x0 * (
508        - jd_xm / hx_xm * ux_xm + jd_xp / hx_xp * ux_xp
509        - jd_x0 / hy    * uy_ym + jd_x0 / hy    * uy_yp
510#if NDIMS == 3
511        - jd_x0 / hz    * uz_zm + jd_x0 / hz    * uz_zp
512#endif
513    );
514    rhs[cnt] = prefactor * div;
515  END

and to monitor the maximum divergence of the flow field as a logging process:

src/logging/divergence.c
 73  BEGIN
 74    const double hx_xm = HXXF(i  );
 75    const double hx_xp = HXXF(i+1);
 76    const double hy    = HYXC(i  );
 77    const double jd_xm = JDXF(i  );
 78    const double jd_x0 = JDXC(i  );
 79    const double jd_xp = JDXF(i+1);
 80#if NDIMS == 2
 81    const double ux_xm = UX(i  , j  );
 82    const double ux_xp = UX(i+1, j  );
 83    const double uy_ym = UY(i  , j  );
 84    const double uy_yp = UY(i  , j+1);
 85#else
 86    const double ux_xm = UX(i  , j  , k  );
 87    const double ux_xp = UX(i+1, j  , k  );
 88    const double uy_ym = UY(i  , j  , k  );
 89    const double uy_yp = UY(i  , j+1, k  );
 90    const double uz_zm = UZ(i  , j  , k  );
 91    const double uz_zp = UZ(i  , j  , k+1);
 92#endif
 93    const double div = 1. / jd_x0 * (
 94        - jd_xm / hx_xm * ux_xm + jd_xp / hx_xp * ux_xp
 95        - jd_x0 / hy    * uy_ym + jd_x0 / hy    * uy_yp
 96#if NDIMS == 3
 97        - jd_x0 / hz    * uz_zm + jd_x0 / hz    * uz_zp
 98#endif
 99    );
100    // check maximum
101    divmax = fmax(divmax, fabs(div));
102  END