Momentum Balance¶
The momentum balance is
\[\pder{u_i}{t}
=
P_i
+
A_i
+
D_{i1}
+
D_{i2}
+
D_{i3},\]
where \(P_i\) and \(A_i\) are the pressure-gradient terms:
\[P_i
\equiv
-
\dmompre{i},\]
and the advective terms:
\[A_i
\equiv
-
\dmomadv{i}{1}
-
\dmomadv{i}{2}
-
\dmomadv{i}{3},\]
respectively.
\(D_{ij}\) is the diffusive term involving spatial differentiation in the \(j\)-th direction:
\[D_{ij}
\equiv
\dmomdif{j}{i}
\,\,
(\text{no summation over}\,j).\]
See the spatial discretization.
The temporal discretization for each Runge-Kutta iteration leads to
\[ \begin{align}\begin{aligned}\Delta u_i
&
=
\gamma^k \Delta t P_i^k
+
\alpha^k \Delta t \left( A_i^{k } + D_i^{k } \right)
+
\beta^k \Delta t \left( A_i^{k-1} + D_i^{k-1} \right),\\u_i^*
&
=
u_i^k
+
\Delta u_i,\end{aligned}\end{align} \]
when all advective and diffusive terms are treated explicitly, while
\[ \begin{align}\begin{aligned}\newcommand{\lap}[2]{
{#2} \gamma^k \Delta t
\frac{\sqrt{Pr}}{\sqrt{Ra}}
\frac{1}{J}
\dif{}{\gcs{#1}}
\left(
\frac{J}{\sfact{#1}}
\frac{1}{\sfact{#1}}
\dif{}{\gcs{#1}}
\right)
}
\Delta u_i
&
=
\gamma^k \Delta t P_i^k
+
\alpha^k \Delta t A_i^{k }
+
\beta^k \Delta t A_i^{k-1}
+
\gamma^k \Delta t \left( D_{i1}^k + D_{i2}^k + D_{i3}^k \right),\\u_i^*
&
=
u_i^k
+
\left\{
1
-
\lap{3}{c}
\right\}^{-1}
\left\{
1
-
\lap{2}{c}
\right\}^{-1}
\left\{
1
-
\lap{1}{c}
\right\}^{-1}
\Delta u_i,\end{aligned}\end{align} \]
when diffusive terms are treated implicitly.
Although we obtain a new velocity field, this does not satisfy the incompressibility in general, which necessitates the additional procedure below.
SMAC method¶
In addition to the incompressibility constraint, we need to somehow update the pressure field as well, which we do not have any equation such as:
\[\pder{p}{t}
=
\cdots.\]
To overcome these issues, we adopt the Simplified Marker And Cell (SMAC) method (Amsden and Harlow, J. Comput. Phys. (6), 1970), which is a two-step method.
Prediction Step¶
In the first step (prediction step), momentum equation is integrated in time, without taking care of the mass conservation, as discussed above. Basically the procedure is identical to how we handle the temperature field. First, the explicit and implicit terms are calculated and stored to the corresponding buffers:
src/fluid/predict/ux.c¶
363int compute_rhs_ux (
364 const domain_t * domain,
365 fluid_t * fluid
366) {
367 if (!laplacians.is_initialised) {
368 if (0 != init_laplacians(domain)) {
369 return 1;
370 }
371 }
372 const double * restrict ux = fluid->ux.data;
373 const double * restrict uy = fluid->uy.data;
374#if NDIMS == 3
375 const double * restrict uz = fluid->uz.data;
376#endif
377 const double * restrict p = fluid-> p.data;
378 const double * restrict t = fluid-> t.data;
379 // buffer for explicit terms
380 double * restrict srca = fluid->srcux.alpha.data;
381 // buffer for implicit terms
382 double * restrict srcg = fluid->srcux.gamma.data;
383 const double diffusivity = fluid_compute_momentum_diffusivity(fluid);
384 // advective contributions, always explicit
385 advection_x(domain, ux, srca);
386 advection_y(domain, ux, uy, srca);
387#if NDIMS == 3
388 advection_z(domain, ux, uz, srca);
389#endif
390 // pressure-gradient contribution, always implicit
391 pressure(domain, p, srcg);
392 // diffusive contributions, can be explicit or implicit
393 diffusion_x(domain, diffusivity, ux, param_m_implicit_x ? srcg : srca);
394 diffusion_y(domain, diffusivity, ux, param_m_implicit_y ? srcg : srca);
395#if NDIMS == 3
396 diffusion_z(domain, diffusivity, ux, param_m_implicit_z ? srcg : srca);
397#endif
398 // buoyancy contribution, always explicit
399 buoyancy(domain, t, srca);
400 return 0;
401}
src/fluid/predict/uy.c¶
325int compute_rhs_uy(
326 const domain_t * domain,
327 fluid_t * fluid
328) {
329 if (!laplacians.is_initialised) {
330 if (0 != init_laplacians(domain)) {
331 return 1;
332 }
333 }
334 const double * restrict ux = fluid->ux.data;
335 const double * restrict uy = fluid->uy.data;
336#if NDIMS == 3
337 const double * restrict uz = fluid->uz.data;
338#endif
339 const double * restrict p = fluid-> p.data;
340 // buffer for explicit terms
341 double * restrict srca = fluid->srcuy.alpha.data;
342 // buffer for implicit terms
343 double * restrict srcg = fluid->srcuy.gamma.data;
344 const double diffusivity = fluid_compute_momentum_diffusivity(fluid);
345 // advective contributions, always explicit
346 advection_x(domain, uy, ux, srca);
347 advection_y(domain, uy, srca);
348#if NDIMS == 3
349 advection_z(domain, uy, uz, srca);
350#endif
351 // pressure-gradient contribution, always implicit
352 pressure(domain, p, srcg);
353 // diffusive contributions, can be explicit or implicit
354 diffusion_x(domain, diffusivity, uy, param_m_implicit_x ? srcg : srca);
355 diffusion_y(domain, diffusivity, uy, param_m_implicit_y ? srcg : srca);
356#if NDIMS == 3
357 diffusion_z(domain, diffusivity, uy, param_m_implicit_z ? srcg : srca);
358#endif
359 return 0;
360}
src/fluid/predict/uz.c¶
255int compute_rhs_uz (
256 const domain_t * domain,
257 fluid_t * fluid
258) {
259 if (!laplacians.is_initialised) {
260 if (0 != init_laplacians(domain)) {
261 return 1;
262 }
263 }
264 const double * restrict ux = fluid->ux.data;
265 const double * restrict uy = fluid->uy.data;
266 const double * restrict uz = fluid->uz.data;
267 const double * restrict p = fluid-> p.data;
268 // buffer for explicit terms
269 double * restrict srca = fluid->srcuz.alpha.data;
270 // buffer for implicit terms
271 double * restrict srcg = fluid->srcuz.gamma.data;
272 const double diffusivity = fluid_compute_momentum_diffusivity(fluid);
273 // advective contributions, always explicit
274 advection_x(domain, uz, ux, srca);
275 advection_y(domain, uz, uy, srca);
276 advection_z(domain, uz, srca);
277 // pressure-gradient contribution, always implicit
278 pressure(domain, p, srcg);
279 // diffusive contributions, can be explicit or implicit
280 diffusion_x(domain, diffusivity, uz, param_m_implicit_x ? srcg : srca);
281 diffusion_y(domain, diffusivity, uz, param_m_implicit_y ? srcg : srca);
282 diffusion_z(domain, diffusivity, uz, param_m_implicit_z ? srcg : srca);
283 return 0;
284}
The stored values are used to compute \(\Delta u_i\):
src/fluid/predict/ux.c¶
436 const double coef_a = rkcoefs[rkstep].alpha;
437 const double coef_b = rkcoefs[rkstep].beta ;
438 const double coef_g = rkcoefs[rkstep].gamma;
439 const double * restrict srcuxa = fluid->srcux.alpha.data;
440 const double * restrict srcuxb = fluid->srcux.beta .data;
441 const double * restrict srcuxg = fluid->srcux.gamma.data;
442 const int isize = domain->mysizes[0];
443 const int jsize = domain->mysizes[1];
444#if NDIMS == 3
445 const int ksize = domain->mysizes[2];
446#endif
447 double * restrict dux = linear_system.x1pncl;
448#if NDIMS == 2
449 const size_t nitems = (isize - 1) * jsize;
450#else
451 const size_t nitems = (isize - 1) * jsize * ksize;
452#endif
453 for (size_t n = 0; n < nitems; n++) {
454 dux[n] =
455 + coef_a * dt * srcuxa[n]
456 + coef_b * dt * srcuxb[n]
457 + coef_g * dt * srcuxg[n];
458 }
src/fluid/predict/uy.c¶
395 const double coef_a = rkcoefs[rkstep].alpha;
396 const double coef_b = rkcoefs[rkstep].beta ;
397 const double coef_g = rkcoefs[rkstep].gamma;
398 const double * restrict srcuya = fluid->srcuy.alpha.data;
399 const double * restrict srcuyb = fluid->srcuy.beta .data;
400 const double * restrict srcuyg = fluid->srcuy.gamma.data;
401 const int isize = domain->mysizes[0];
402 const int jsize = domain->mysizes[1];
403#if NDIMS == 3
404 const int ksize = domain->mysizes[2];
405#endif
406 double * restrict duy = linear_system.x1pncl;
407#if NDIMS == 2
408 const size_t nitems = isize * jsize;
409#else
410 const size_t nitems = isize * jsize * ksize;
411#endif
412 for (size_t n = 0; n < nitems; n++) {
413 duy[n] =
414 + coef_a * dt * srcuya[n]
415 + coef_b * dt * srcuyb[n]
416 + coef_g * dt * srcuyg[n];
417 }
src/fluid/predict/uz.c¶
315const double coef_a = rkcoefs[rkstep].alpha;
316const double coef_b = rkcoefs[rkstep].beta ;
317const double coef_g = rkcoefs[rkstep].gamma;
318const double * restrict srcuza = fluid->srcuz.alpha.data;
319const double * restrict srcuzb = fluid->srcuz.beta .data;
320const double * restrict srcuzg = fluid->srcuz.gamma.data;
321const int isize = domain->mysizes[0];
322const int jsize = domain->mysizes[1];
323const int ksize = domain->mysizes[2];
324double * restrict duz = linear_system.x1pncl;
325const size_t nitems = isize * jsize * ksize;
326for (size_t n = 0; n < nitems; n++) {
327 duz[n] =
328 + coef_a * dt * srcuza[n]
329 + coef_b * dt * srcuzb[n]
330 + coef_g * dt * srcuzg[n];
331}
When necessary, linear systems are solved to take care of the implicit treatments:
src/fluid/predict/ux.c¶
465if (param_m_implicit_x) {
466 tdm_info_t * tdm_info = linear_system.tdm_x;
467 int size = 0;
468 double * restrict tdm_l = NULL;
469 double * restrict tdm_c = NULL;
470 double * restrict tdm_u = NULL;
471 tdm.get_size(tdm_info, &size);
472 tdm.get_l(tdm_info, &tdm_l);
473 tdm.get_c(tdm_info, &tdm_c);
474 tdm.get_u(tdm_info, &tdm_u);
475 const laplacian_t * lapx = laplacians.lapx;
476 for (int i = 0; i < size; i++) {
477 tdm_l[i] = - prefactor * lapx[i].l;
478 tdm_c[i] = 1. - prefactor * lapx[i].c;
479 tdm_u[i] = - prefactor * lapx[i].u;
480 }
481 tdm.solve(tdm_info, linear_system.x1pncl);
482}
src/fluid/predict/ux.c¶
484if (param_m_implicit_y) {
485 sdecomp.transpose.execute(
486 linear_system.transposer_x1_to_y1,
487 linear_system.x1pncl,
488 linear_system.y1pncl
489 );
490 tdm_info_t * tdm_info = linear_system.tdm_y;
491 int size = 0;
492 double * restrict tdm_l = NULL;
493 double * restrict tdm_c = NULL;
494 double * restrict tdm_u = NULL;
495 tdm.get_size(tdm_info, &size);
496 tdm.get_l(tdm_info, &tdm_l);
497 tdm.get_c(tdm_info, &tdm_c);
498 tdm.get_u(tdm_info, &tdm_u);
499 const laplacian_t * lapy = &laplacians.lapy;
500 for (int j = 0; j < size; j++) {
501 tdm_l[j] = - prefactor * (*lapy).l;
502 tdm_c[j] = 1. - prefactor * (*lapy).c;
503 tdm_u[j] = - prefactor * (*lapy).u;
504 }
505 tdm.solve(tdm_info, linear_system.y1pncl);
506 sdecomp.transpose.execute(
507 linear_system.transposer_y1_to_x1,
508 linear_system.y1pncl,
509 linear_system.x1pncl
510 );
511}
src/fluid/predict/ux.c¶
514if (param_m_implicit_z) {
515 sdecomp.transpose.execute(
516 linear_system.transposer_x1_to_z2,
517 linear_system.x1pncl,
518 linear_system.z2pncl
519 );
520 tdm_info_t * tdm_info = linear_system.tdm_z;
521 int size = 0;
522 double * restrict tdm_l = NULL;
523 double * restrict tdm_c = NULL;
524 double * restrict tdm_u = NULL;
525 tdm.get_size(tdm_info, &size);
526 tdm.get_l(tdm_info, &tdm_l);
527 tdm.get_c(tdm_info, &tdm_c);
528 tdm.get_u(tdm_info, &tdm_u);
529 const laplacian_t * lapz = &laplacians.lapz;
530 for (int k = 0; k < size; k++) {
531 tdm_l[k] = - prefactor * (*lapz).l;
532 tdm_c[k] = 1. - prefactor * (*lapz).c;
533 tdm_u[k] = - prefactor * (*lapz).u;
534 }
535 tdm.solve(tdm_info, linear_system.z2pncl);
536 sdecomp.transpose.execute(
537 linear_system.transposer_z2_to_x1,
538 linear_system.z2pncl,
539 linear_system.x1pncl
540 );
541}
src/fluid/predict/uy.c¶
424if (param_m_implicit_x) {
425 tdm_info_t * tdm_info = linear_system.tdm_x;
426 int size = 0;
427 double * restrict tdm_l = NULL;
428 double * restrict tdm_c = NULL;
429 double * restrict tdm_u = NULL;
430 tdm.get_size(tdm_info, &size);
431 tdm.get_l(tdm_info, &tdm_l);
432 tdm.get_c(tdm_info, &tdm_c);
433 tdm.get_u(tdm_info, &tdm_u);
434 const laplacian_t * lapx = laplacians.lapx;
435 for (int i = 0; i < size; i++) {
436 tdm_l[i] = - prefactor * lapx[i].l;
437 tdm_c[i] = 1. - prefactor * lapx[i].c;
438 tdm_u[i] = - prefactor * lapx[i].u;
439 }
440 tdm.solve(tdm_info, linear_system.x1pncl);
441}
src/fluid/predict/uy.c¶
443if (param_m_implicit_y) {
444 sdecomp.transpose.execute(
445 linear_system.transposer_x1_to_y1,
446 linear_system.x1pncl,
447 linear_system.y1pncl
448 );
449 tdm_info_t * tdm_info = linear_system.tdm_y;
450 int size = 0;
451 double * restrict tdm_l = NULL;
452 double * restrict tdm_c = NULL;
453 double * restrict tdm_u = NULL;
454 tdm.get_size(tdm_info, &size);
455 tdm.get_l(tdm_info, &tdm_l);
456 tdm.get_c(tdm_info, &tdm_c);
457 tdm.get_u(tdm_info, &tdm_u);
458 const laplacian_t * lapy = &laplacians.lapy;
459 for (int j = 0; j < size; j++) {
460 tdm_l[j] = - prefactor * (*lapy).l;
461 tdm_c[j] = 1. - prefactor * (*lapy).c;
462 tdm_u[j] = - prefactor * (*lapy).u;
463 }
464 tdm.solve(tdm_info, linear_system.y1pncl);
465 sdecomp.transpose.execute(
466 linear_system.transposer_y1_to_x1,
467 linear_system.y1pncl,
468 linear_system.x1pncl
469 );
470}
src/fluid/predict/uy.c¶
473if (param_m_implicit_z) {
474 sdecomp.transpose.execute(
475 linear_system.transposer_x1_to_z2,
476 linear_system.x1pncl,
477 linear_system.z2pncl
478 );
479 tdm_info_t * tdm_info = linear_system.tdm_z;
480 int size = 0;
481 double * restrict tdm_l = NULL;
482 double * restrict tdm_c = NULL;
483 double * restrict tdm_u = NULL;
484 tdm.get_size(tdm_info, &size);
485 tdm.get_l(tdm_info, &tdm_l);
486 tdm.get_c(tdm_info, &tdm_c);
487 tdm.get_u(tdm_info, &tdm_u);
488 const laplacian_t * lapz = &laplacians.lapz;
489 for (int k = 0; k < size; k++) {
490 tdm_l[k] = - prefactor * (*lapz).l;
491 tdm_c[k] = 1. - prefactor * (*lapz).c;
492 tdm_u[k] = - prefactor * (*lapz).u;
493 }
494 tdm.solve(tdm_info, linear_system.z2pncl);
495 sdecomp.transpose.execute(
496 linear_system.transposer_z2_to_x1,
497 linear_system.z2pncl,
498 linear_system.x1pncl
499 );
500}
src/fluid/predict/uz.c¶
338if (param_m_implicit_x) {
339 tdm_info_t * tdm_info = linear_system.tdm_x;
340 int size = 0;
341 double * restrict tdm_l = NULL;
342 double * restrict tdm_c = NULL;
343 double * restrict tdm_u = NULL;
344 tdm.get_size(tdm_info, &size);
345 tdm.get_l(tdm_info, &tdm_l);
346 tdm.get_c(tdm_info, &tdm_c);
347 tdm.get_u(tdm_info, &tdm_u);
348 const laplacian_t * lapx = laplacians.lapx;
349 for (int i = 0; i < size; i++) {
350 tdm_l[i] = - prefactor * lapx[i].l;
351 tdm_c[i] = 1. - prefactor * lapx[i].c;
352 tdm_u[i] = - prefactor * lapx[i].u;
353 }
354 tdm.solve(tdm_info, linear_system.x1pncl);
355}
src/fluid/predict/uz.c¶
357if (param_m_implicit_y) {
358 sdecomp.transpose.execute(
359 linear_system.transposer_x1_to_y1,
360 linear_system.x1pncl,
361 linear_system.y1pncl
362 );
363 tdm_info_t * tdm_info = linear_system.tdm_y;
364 int size = 0;
365 double * restrict tdm_l = NULL;
366 double * restrict tdm_c = NULL;
367 double * restrict tdm_u = NULL;
368 tdm.get_size(tdm_info, &size);
369 tdm.get_l(tdm_info, &tdm_l);
370 tdm.get_c(tdm_info, &tdm_c);
371 tdm.get_u(tdm_info, &tdm_u);
372 const laplacian_t * lapy = &laplacians.lapy;
373 for (int j = 0; j < size; j++) {
374 tdm_l[j] = - prefactor * (*lapy).l;
375 tdm_c[j] = 1. - prefactor * (*lapy).c;
376 tdm_u[j] = - prefactor * (*lapy).u;
377 }
378 tdm.solve(tdm_info, linear_system.y1pncl);
379 sdecomp.transpose.execute(
380 linear_system.transposer_y1_to_x1,
381 linear_system.y1pncl,
382 linear_system.x1pncl
383 );
384}
src/fluid/predict/uz.c¶
386if (param_m_implicit_z) {
387 sdecomp.transpose.execute(
388 linear_system.transposer_x1_to_z2,
389 linear_system.x1pncl,
390 linear_system.z2pncl
391 );
392 tdm_info_t * tdm_info = linear_system.tdm_z;
393 int size = 0;
394 double * restrict tdm_l = NULL;
395 double * restrict tdm_c = NULL;
396 double * restrict tdm_u = NULL;
397 tdm.get_size(tdm_info, &size);
398 tdm.get_l(tdm_info, &tdm_l);
399 tdm.get_c(tdm_info, &tdm_c);
400 tdm.get_u(tdm_info, &tdm_u);
401 const laplacian_t * lapz = &laplacians.lapz;
402 for (int k = 0; k < size; k++) {
403 tdm_l[k] = - prefactor * (*lapz).l;
404 tdm_c[k] = 1. - prefactor * (*lapz).c;
405 tdm_u[k] = - prefactor * (*lapz).u;
406 }
407 tdm.solve(tdm_info, linear_system.z2pncl);
408 sdecomp.transpose.execute(
409 linear_system.transposer_z2_to_x1,
410 linear_system.z2pncl,
411 linear_system.x1pncl
412 );
413}
Finally the velocity field is updated:
src/fluid/predict/ux.c¶
546 const int isize = domain->mysizes[0];
547 const int jsize = domain->mysizes[1];
548#if NDIMS == 3
549 const int ksize = domain->mysizes[2];
550#endif
551 const double * restrict dux = linear_system.x1pncl;
552 double * restrict ux = fluid->ux.data;
553 BEGIN
554#if NDIMS == 2
555 UX(i, j) += dux[cnt];
556#else
557 UX(i, j, k) += dux[cnt];
558#endif
559 END
560 if (0 != fluid_update_boundaries_ux(domain, &fluid->ux)) {
561 return 1;
562 }
src/fluid/predict/uy.c¶
505 const int isize = domain->mysizes[0];
506 const int jsize = domain->mysizes[1];
507#if NDIMS == 3
508 const int ksize = domain->mysizes[2];
509#endif
510 const double * restrict duy = linear_system.x1pncl;
511 double * restrict uy = fluid->uy.data;
512 BEGIN
513#if NDIMS == 2
514 UY(i, j) += duy[cnt];
515#else
516 UY(i, j, k) += duy[cnt];
517#endif
518 END
519 if (0 != fluid_update_boundaries_uy(domain, &fluid->uy)) {
520 return 1;
521 }
src/fluid/predict/uz.c¶
417const int isize = domain->mysizes[0];
418const int jsize = domain->mysizes[1];
419const int ksize = domain->mysizes[2];
420const double * restrict duz = linear_system.x1pncl;
421double * restrict uz = fluid->uz.data;
422BEGIN
423 UZ(i, j, k) += duz[cnt];
424END
425if (0 != fluid_update_boundaries_uz(domain, &fluid->uz)) {
426 return 1;
427}
Correction Step¶
The updated velocity field \(u_i^*\), which in general violates the incompressibility, is corrected in the second step (correction step). The idea is mathematically written as
\[u_i^{k+1}
=
u_i^*
-
\gamma^k \Delta t
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}},\]
where \(\psi\) is a scalar potential to be given. By taking the (discrete) divergence, we obtain
\[\frac{1}{J}
\dif{}{\gcs{i}}
\left(
\frac{J}{\sfact{i}}
\vel{i}^{k+1}
\right)
=
\frac{1}{J}
\dif{}{\gcs{i}}
\left(
\frac{J}{\sfact{i}}
\vel{i}^*
\right)
-
\gamma^k \Delta t
\frac{1}{J}
\dif{}{\gcs{i}}
\left(
\frac{J}{\sfact{i}}
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}}
\right).\]
By requesting the (discrete) incompressibility constraint on the new velocity field (namely the left-hand-side term to be zero), we obtain
\[\frac{1}{J}
\dif{}{\gcs{i}}
\left(
\frac{J}{\sfact{i}}
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}}
\right)
=
\frac{1}{\gamma^k \Delta t}
\frac{1}{J}
\dif{}{\gcs{i}}
\left(
\frac{J}{\sfact{i}}
\vel{i}^*
\right),\]
which is a Poisson equation with respect to \(\psi\).
The right-hand-side term is computed as follows in the code:
src/fluid/compute_potential.c¶
455static int assign_input (
456 const domain_t * domain,
457 const size_t rkstep,
458 const double dt,
459 const fluid_t * fluid,
460 double * restrict rhs
461) {
462 const int isize = domain->mysizes[0];
463 const int jsize = domain->mysizes[1];
464#if NDIMS == 3
465 const int ksize = domain->mysizes[2];
466#endif
467 const double * restrict hxxf = domain->hxxf;
468 const double hy = domain->hy;
469#if NDIMS == 3
470 const double hz = domain->hz;
471#endif
472 const double * restrict jdxf = domain->jdxf;
473 const double * restrict jdxc = domain->jdxc;
474 const double * restrict ux = fluid->ux.data;
475 const double * restrict uy = fluid->uy.data;
476#if NDIMS == 3
477 const double * restrict uz = fluid->uz.data;
478#endif
479 // normalise FFT beforehand
480#if NDIMS == 2
481 const double norm = 1. * domain->glsizes[1];
482#else
483 const double norm = 1. * domain->glsizes[1] * domain->glsizes[2];
484#endif
485 const double prefactor = 1. / (rkcoefs[rkstep].gamma * dt) / norm;
486 BEGIN
487 const double hx_xm = HXXF(i );
488 const double hx_xp = HXXF(i+1);
489 const double jd_xm = JDXF(i );
490 const double jd_x0 = JDXC(i );
491 const double jd_xp = JDXF(i+1);
492#if NDIMS == 2
493 const double ux_xm = UX(i , j );
494 const double ux_xp = UX(i+1, j );
495 const double uy_ym = UY(i , j );
496 const double uy_yp = UY(i , j+1);
497#else
498 const double ux_xm = UX(i , j , k );
499 const double ux_xp = UX(i+1, j , k );
500 const double uy_ym = UY(i , j , k );
501 const double uy_yp = UY(i , j+1, k );
502 const double uz_zm = UZ(i , j , k );
503 const double uz_zp = UZ(i , j , k+1);
504#endif
505 const double div = 1. / jd_x0 * (
506 - jd_xm / hx_xm * ux_xm + jd_xp / hx_xp * ux_xp
507 - jd_x0 / hy * uy_ym + jd_x0 / hy * uy_yp
508#if NDIMS == 3
509 - jd_x0 / hz * uz_zm + jd_x0 / hz * uz_zp
510#endif
511 );
512 rhs[cnt] = prefactor * div;
513 END
514 return 0;
515}
After solving the Poisson equation, the velocity field is corrected as follows in the code:
src/fluid/correct/ux.c¶
28int fluid_correct_velocity_ux (
29 const domain_t * domain,
30 const double prefactor,
31 fluid_t * fluid
32) {
33 const int isize = domain->mysizes[0];
34 const int jsize = domain->mysizes[1];
35#if NDIMS == 3
36 const int ksize = domain->mysizes[2];
37#endif
38 const double * restrict hxxf = domain->hxxf;
39 const double * restrict psi = fluid->psi.data;
40 double * restrict ux = fluid->ux.data;
41 BEGIN
42 const double hx = HXXF(i );
43#if NDIMS == 2
44 const double psi_xm = PSI(i-1, j );
45 const double psi_xp = PSI(i , j );
46 double * vel = &UX(i, j);
47#else
48 const double psi_xm = PSI(i-1, j , k );
49 const double psi_xp = PSI(i , j , k );
50 double * vel = &UX(i, j, k);
51#endif
52 *vel -= prefactor / hx * (
53 - psi_xm
54 + psi_xp
55 );
56 END
57 // update boundary and halo cells
58 if (0 != fluid_update_boundaries_ux(domain, &fluid->ux)) {
59 return 1;
60 }
61 return 0;
62}
src/fluid/correct/uy.c¶
27int fluid_correct_velocity_uy (
28 const domain_t * domain,
29 const double prefactor,
30 fluid_t * fluid
31) {
32 const int isize = domain->mysizes[0];
33 const int jsize = domain->mysizes[1];
34#if NDIMS == 3
35 const int ksize = domain->mysizes[2];
36#endif
37 const double hy = domain->hy;
38 const double * restrict psi = fluid->psi.data;
39 double * restrict uy = fluid->uy.data;
40 BEGIN
41#if NDIMS == 2
42 const double psi_ym = PSI(i , j-1);
43 const double psi_yp = PSI(i , j );
44 double * vel = &UY(i, j);
45#else
46 const double psi_ym = PSI(i , j-1, k );
47 const double psi_yp = PSI(i , j , k );
48 double * vel = &UY(i, j, k);
49#endif
50 *vel -= prefactor / hy * (
51 - psi_ym
52 + psi_yp
53 );
54 END
55 // update boundary and halo cells
56 if (0 != fluid_update_boundaries_uy(domain, &fluid->uy)) {
57 return 1;
58 }
59 return 0;
60}
src/fluid/correct/uz.c¶
10int fluid_correct_velocity_uz (
11 const domain_t * domain,
12 const double prefactor,
13 fluid_t * fluid
14) {
15 const int isize = domain->mysizes[0];
16 const int jsize = domain->mysizes[1];
17 const int ksize = domain->mysizes[2];
18 const double hz = domain->hz;
19 const double * restrict psi = fluid->psi.data;
20 double * restrict uz = fluid->uz.data;
21 for (int k = 1; k <= ksize; k++) {
22 for (int j = 1; j <= jsize; j++) {
23 for (int i = 1; i <= isize; i++) {
24 const double psi_zm = PSI(i , j , k-1);
25 const double psi_zp = PSI(i , j , k );
26 UZ(i, j, k) -= prefactor / hz * (
27 - psi_zm
28 + psi_zp
29 );
30 }
31 }
32 }
33 // update boundary and halo cells
34 if (0 != fluid_update_boundaries_uz(domain, &fluid->uz)) {
35 return 1;
36 }
37 return 0;
38}
Pressure and Scalar Potential¶
Finally we relate the pressure field with the scalar potential to close the system. Each Runge-Kutta step when all diffusive terms are treated explicitly is given by
\[ \begin{align}\begin{aligned}\Delta u_i
&
=
\gamma^k \Delta t P_i^k
+
\alpha^k \Delta t \left( A_i^{k } + D_{i1}^{k } + D_{i2}^{k } + D_{i3}^{k } \right)
+
\beta^k \Delta t \left( A_i^{k-1} + D_{i1}^{k-1} + D_{i2}^{k-1} + D_{i3}^{k-1} \right),\\u_i^*
&
=
u_i^k
+
\Delta u_i,\\u_i^{k+1}
&
=
u_i^*
-
\gamma^k \Delta t
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}}.\end{aligned}\end{align} \]
Summing all three steps yield
\[u_i^{k+1}
=
u_i^k
-
\gamma^k \Delta t
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}}
+
\gamma^k \Delta t P_i^k
+
\alpha^k \Delta t \left( A_i^{k } + D_{i1}^{k } + D_{i2}^{k } + D_{i3}^{k } \right)
+
\beta^k \Delta t \left( A_i^{k-1} + D_{i1}^{k-1} + D_{i2}^{k-1} + D_{i3}^{k-1} \right).\]
Since the pressure field should be treated implicitly in time, we find a requirement:
\[-
\gamma^k \Delta t
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}}
+
\gamma^k \Delta t P_i^k
=
\gamma^k \Delta t P_i^{k+1},\]
or equivalently
\[p^{k+1}
=
p^k
+
\psi.\]
Next, we consider cases where the diffusive terms are partially treated implicitly in time; for instance, the following relation holds when the diffusive terms are treated implicitly only in \(x\) direction:
\[ \begin{align}\begin{aligned}\Delta u_i
&
=
\gamma^k \Delta t P_i^k
+
\alpha^k \Delta t \left( A_i^{k } + D_{i2}^{k } + D_{i3}^{k } \right)
+
\beta^k \Delta t \left( A_i^{k-1} + D_{i2}^{k-1} + D_{i3}^{k-1} \right)
+
\gamma^k \Delta t D_{i1}^{k },\\u_i^*
&
=
u_i^k
+
\left\{
1
-
\lap{1}{c}
\right\}^{-1}
\Delta u_i,\\u_i^{k+1}
&
=
u_i^*
-
\gamma^k \Delta t
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}}.\end{aligned}\end{align} \]
Note that, in the second step, the left-hand side is \(u_i^*\) while it should be \(u_i^{k+1}\) but is unknown. By requesting the pressure-gradient term to be implicit in time and with some algebra, we obtain
\[\gamma^k \Delta t
\left\{
1
-
\lap{1}{c}
\right\}
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}}
-
\gamma^k \Delta t
P_i^k
=
-
\gamma^k \Delta t
P_i^{k+1},\]
or equivalently
\[\frac{1}{\sfact{i}}
\dif{p^{k+1}}{\gcs{i}}
=
\frac{1}{\sfact{i}}
\dif{p^k}{\gcs{i}}
+
\left\{
1
-
\lap{1}{c}
\right\}
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}}.\]
Since the spatial-differential operators are interchangeable, we simplify the relation to
\[p^{k+1}
=
p^k
+
\psi
-
\lap{1}{c}
\psi.\]
Similarly, when all diffusive terms are treated implicitly, we have
\[ \begin{align}\begin{aligned}\Delta u_i
&
=
\gamma^k \Delta t P_i^k
+
\alpha^k \Delta t A_i^{k }
+
\beta^k \Delta t A_i^{k-1}
+
\gamma^k \Delta t \left( D_{i1}^{k } + D_{i2}^{k } + D_{i3}^{k } \right),\\u_i^*
&
=
u_i^k
+
\left\{
1
-
\lap{3}{c}
\right\}^{-1}
\left\{
1
-
\lap{2}{c}
\right\}^{-1}
\left\{
1
-
\lap{1}{c}
\right\}^{-1}
\Delta u_i,\\u_i^{k+1}
&
=
u_i^*
-
\gamma^k \Delta t
\frac{1}{\sfact{i}}
\dif{\psi}{\gcs{i}},\\p^{k+1}
&
=
p^k
+
\psi
-
\lap{1}{c}
\psi
-
\lap{2}{c}
\psi
-
\lap{3}{c}
\psi.\end{aligned}\end{align} \]
Updating pressure field using the scalar potential is implemented as follows:
src/fluid/update_pressure.c¶
166int fluid_update_pressure (
167 const domain_t * domain,
168 const size_t rkstep,
169 const double dt,
170 fluid_t * fluid
171) {
172 // explicit contribution, always present
173 explicit_contribution(domain, fluid);
174 const double prefactor =
175 0.5 * rkcoefs[rkstep].gamma * dt * fluid_compute_momentum_diffusivity(fluid);
176 // additional terms if diffusive terms in the direction is treated implicitly
177 if (param_m_implicit_x) {
178 implicit_x_contribution(domain, prefactor, fluid);
179 }
180 if (param_m_implicit_y) {
181 implicit_y_contribution(domain, prefactor, fluid);
182 }
183#if NDIMS == 3
184 if (param_m_implicit_z) {
185 implicit_z_contribution(domain, prefactor, fluid);
186 }
187#endif
188 // impose boundary conditions and communicate halo cells
189 if (0 != fluid_update_boundaries_p(domain, &fluid->p)) {
190 return 1;
191 }
192 return 0;
193}
The explicit contribution, which is always present, is given here:
src/fluid/update_pressure.c¶
31static int explicit_contribution (
32 const domain_t * domain,
33 const fluid_t * fluid
34) {
35 const int isize = domain->mysizes[0];
36 const int jsize = domain->mysizes[1];
37#if NDIMS == 3
38 const int ksize = domain->mysizes[2];
39#endif
40 const double * restrict psi = fluid->psi.data;
41 double * restrict p = fluid->p.data;
42 BEGIN
43#if NDIMS == 2
44 P(i, j) += PSI(i, j);
45#else
46 P(i, j, k) += PSI(i, j, k);
47#endif
48 END
49 return 0;
50}
The implicit contributions, which is needed when the Laplace operator in the direction is implicitly treated, are given here:
src/fluid/update_pressure.c¶
53static int implicit_x_contribution (
54 const domain_t * domain,
55 const double prefactor,
56 fluid_t * fluid
57) {
58 const int isize = domain->mysizes[0];
59 const int jsize = domain->mysizes[1];
60#if NDIMS == 3
61 const int ksize = domain->mysizes[2];
62#endif
63 const double * restrict hxxf = domain->hxxf;
64 const double * restrict jdxf = domain->jdxf;
65 const double * restrict jdxc = domain->jdxc;
66 const double * restrict psi = fluid->psi.data;
67 double * restrict p = fluid->p.data;
68 BEGIN
69 const double hx_xm = HXXF(i );
70 const double hx_xp = HXXF(i+1);
71 const double jd_xm = JDXF(i );
72 const double jd_x0 = JDXC(i );
73 const double jd_xp = JDXF(i+1);
74 const double l = 1. / jd_x0 * jd_xm / hx_xm / hx_xm;
75 const double u = 1. / jd_x0 * jd_xp / hx_xp / hx_xp;
76 const double c = - l - u;
77#if NDIMS == 2
78 const double psi_xm = PSI(i-1, j );
79 const double psi_x0 = PSI(i , j );
80 const double psi_xp = PSI(i+1, j );
81 double * pre = &P(i, j);
82#else
83 const double psi_xm = PSI(i-1, j , k );
84 const double psi_x0 = PSI(i , j , k );
85 const double psi_xp = PSI(i+1, j , k );
86 double * pre = &P(i, j, k);
87#endif
88 *pre -= prefactor * (
89 + l * psi_xm
90 + c * psi_x0
91 + u * psi_xp
92 );
93 END
94 return 0;
95}
src/fluid/update_pressure.c¶
98static int implicit_y_contribution (
99 const domain_t * domain,
100 const double prefactor,
101 fluid_t * fluid
102) {
103 const int isize = domain->mysizes[0];
104 const int jsize = domain->mysizes[1];
105#if NDIMS == 3
106 const int ksize = domain->mysizes[2];
107#endif
108 const double hy = domain->hy;
109 const double * restrict psi = fluid->psi.data;
110 double * restrict p = fluid->p.data;
111 BEGIN
112 const double l = 1. / hy / hy;
113 const double u = 1. / hy / hy;
114 const double c = - l - u;
115#if NDIMS == 2
116 const double psi_ym = PSI(i , j-1);
117 const double psi_y0 = PSI(i , j );
118 const double psi_yp = PSI(i , j+1);
119 double * pre = &P(i, j);
120#else
121 const double psi_ym = PSI(i , j-1, k );
122 const double psi_y0 = PSI(i , j , k );
123 const double psi_yp = PSI(i , j+1, k );
124 double * pre = &P(i, j, k);
125#endif
126 *pre -= prefactor * (
127 + l * psi_ym
128 + c * psi_y0
129 + u * psi_yp
130 );
131 END
132 return 0;
133}
src/fluid/update_pressure.c¶
137static int implicit_z_contribution (
138 const domain_t * domain,
139 const double prefactor,
140 fluid_t * fluid
141) {
142 const int isize = domain->mysizes[0];
143 const int jsize = domain->mysizes[1];
144 const int ksize = domain->mysizes[2];
145 const double hz = domain->hz;
146 const double * restrict psi = fluid->psi.data;
147 double * restrict p = fluid->p.data;
148 BEGIN
149 const double l = 1. / hz / hz;
150 const double u = 1. / hz / hz;
151 const double c = - l - u;
152 const double psi_zm = PSI(i , j , k-1);
153 const double psi_z0 = PSI(i , j , k );
154 const double psi_zp = PSI(i , j , k+1);
155 P(i, j, k) -= prefactor * (
156 + l * psi_zm
157 + c * psi_z0
158 + u * psi_zp
159 );
160 END
161 return 0;
162}