Linear System¶
To treat the diffusive terms implicitly, we need to solve a linear system:
\[\left\{
1
-
C
\frac{1}{J}
\dif{}{\gcs{i}}
\left(
\frac{J}{\sfact{i}}
\frac{1}{\sfact{i}}
\dif{}{\gcs{i}}
\right)
\right\}
p
=
q,\]
where \(C\) is a coefficient (e.g., coming from the diffusivity and time-step sizes).
Laplace Operators¶
Note that linear systems involve discrete Laplace operators:
\[\frac{1}{J}
\dif{}{\gcs{i}}
\left(
\frac{J}{\sfact{i}}
\frac{1}{\sfact{i}}
\dif{}{\gcs{i}}
\right),\]
which are pre-computed as follows.
Wall-normal velocity:
src/fluid/predict/ux.c¶
34const int isize = domain->glsizes[0];
35const double * hxxc = domain->hxxc;
36const double * jdxf = domain->jdxf;
37const double * jdxc = domain->jdxc;
38laplacians.lapx = memory_calloc(isize - 1, sizeof(laplacian_t));
39for (int i = 2; i <= isize; i++) {
40 const double l = 1. / JDXF(i ) * JDXC(i-1) / HXXC(i-1) / HXXC(i-1);
41 const double u = 1. / JDXF(i ) * JDXC(i ) / HXXC(i ) / HXXC(i );
42 const double c = - l - u;
43 laplacians.LAPX(i).l = l;
44 laplacians.LAPX(i).c = c;
45 laplacians.LAPX(i).u = u;
46}
src/fluid/predict/ux.c¶
50const double hy = domain->hy;
51const double l = 1. / hy / hy;
52const double u = 1. / hy / hy;
53const double c = - l - u;
54laplacians.lapy.l = l;
55laplacians.lapy.c = c;
56laplacians.lapy.u = u;
src/fluid/predict/ux.c¶
61const double hz = domain->hz;
62const double l = 1. / hz / hz;
63const double u = 1. / hz / hz;
64const double c = - l - u;
65laplacians.lapz.l = l;
66laplacians.lapz.c = c;
67laplacians.lapz.u = u;
Stream-wise velocity:
src/fluid/predict/uy.c¶
33const int isize = domain->glsizes[0];
34const double * hxxf = domain->hxxf;
35const double * jdxf = domain->jdxf;
36const double * jdxc = domain->jdxc;
37laplacians.lapx = memory_calloc(isize, sizeof(laplacian_t));
38for (int i = 1; i <= isize; i++) {
39 const double l = 1. / JDXC(i ) * JDXF(i ) / HXXF(i ) / HXXF(i );
40 const double u = 1. / JDXC(i ) * JDXF(i+1) / HXXF(i+1) / HXXF(i+1);
41 const double c = - l - u;
42 laplacians.LAPX(i).l = l;
43 laplacians.LAPX(i).c = c;
44 laplacians.LAPX(i).u = u;
45}
src/fluid/predict/uy.c¶
49const double hy = domain->hy;
50const double l = 1. / hy / hy;
51const double u = 1. / hy / hy;
52const double c = - l - u;
53laplacians.lapy.l = l;
54laplacians.lapy.c = c;
55laplacians.lapy.u = u;
src/fluid/predict/uy.c¶
60const double hz = domain->hz;
61const double l = 1. / hz / hz;
62const double u = 1. / hz / hz;
63const double c = - l - u;
64laplacians.lapz.l = l;
65laplacians.lapz.c = c;
66laplacians.lapz.u = u;
Span-wise velocity:
src/fluid/predict/uz.c¶
31const int isize = domain->glsizes[0];
32const double * hxxf = domain->hxxf;
33const double * jdxf = domain->jdxf;
34const double * jdxc = domain->jdxc;
35laplacians.lapx = memory_calloc(isize, sizeof(laplacian_t));
36for (int i = 1; i <= isize; i++) {
37 const double l = 1. / JDXC(i ) * JDXF(i ) / HXXF(i ) / HXXF(i );
38 const double u = 1. / JDXC(i ) * JDXF(i+1) / HXXF(i+1) / HXXF(i+1);
39 const double c = - l - u;
40 laplacians.LAPX(i).l = l;
41 laplacians.LAPX(i).c = c;
42 laplacians.LAPX(i).u = u;
43}
src/fluid/predict/uz.c¶
47const double hy = domain->hy;
48const double l = 1. / hy / hy;
49const double u = 1. / hy / hy;
50const double c = - l - u;
51laplacians.lapy.l = l;
52laplacians.lapy.c = c;
53laplacians.lapy.u = u;
src/fluid/predict/uz.c¶
57const double hz = domain->hz;
58const double l = 1. / hz / hz;
59const double u = 1. / hz / hz;
60const double c = - l - u;
61laplacians.lapz.l = l;
62laplacians.lapz.c = c;
63laplacians.lapz.u = u;
Temperature:
src/fluid/predict/t.c¶
33const int isize = domain->glsizes[0];
34const double * hxxf = domain->hxxf;
35const double * jdxf = domain->jdxf;
36const double * jdxc = domain->jdxc;
37laplacians.lapx = memory_calloc(isize, sizeof(laplacian_t));
38for (int i = 1; i <= isize; i++) {
39 const double l = 1. / JDXC(i ) * JDXF(i ) / HXXF(i ) / HXXF(i );
40 const double u = 1. / JDXC(i ) * JDXF(i+1) / HXXF(i+1) / HXXF(i+1);
41 const double c = - l - u;
42 laplacians.LAPX(i).l = l;
43 laplacians.LAPX(i).c = c;
44 laplacians.LAPX(i).u = u;
45}
src/fluid/predict/t.c¶
49const double hy = domain->hy;
50const double l = 1. / hy / hy;
51const double u = 1. / hy / hy;
52const double c = - l - u;
53laplacians.lapy.l = l;
54laplacians.lapy.c = c;
55laplacians.lapy.u = u;
src/fluid/predict/t.c¶
60const double hz = domain->hz;
61const double l = 1. / hz / hz;
62const double u = 1. / hz / hz;
63const double c = - l - u;
64laplacians.lapz.l = l;
65laplacians.lapz.c = c;
66laplacians.lapz.u = u;
Solving Linear Systems¶
Wall-normal velocity:
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}
Stream-wise velocity:
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}
Span-wise velocity:
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}
Temperature:
src/fluid/predict/t.c¶
390if (param_t_implicit_x) {
391 // prepare tri-diagonal coefficients
392 tdm_info_t * tdm_info = linear_system.tdm_x;
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 * lapx = laplacians.lapx;
402 for (int i = 0; i < size; i++) {
403 tdm_l[i] = - prefactor * lapx[i].l;
404 tdm_c[i] = 1. - prefactor * lapx[i].c;
405 tdm_u[i] = - prefactor * lapx[i].u;
406 }
407 // solve all
408 tdm.solve(tdm_info, linear_system.x1pncl);
409}
src/fluid/predict/t.c¶
411if (param_t_implicit_y) {
412 // make all y data accessible
413 sdecomp.transpose.execute(
414 linear_system.transposer_x1_to_y1,
415 linear_system.x1pncl,
416 linear_system.y1pncl
417 );
418 // prepare tri-diagonal coefficients
419 tdm_info_t * tdm_info = linear_system.tdm_y;
420 int size = 0;
421 double * restrict tdm_l = NULL;
422 double * restrict tdm_c = NULL;
423 double * restrict tdm_u = NULL;
424 tdm.get_size(tdm_info, &size);
425 tdm.get_l(tdm_info, &tdm_l);
426 tdm.get_c(tdm_info, &tdm_c);
427 tdm.get_u(tdm_info, &tdm_u);
428 const laplacian_t * lapy = &laplacians.lapy;
429 for (int j = 0; j < size; j++) {
430 tdm_l[j] = - prefactor * (*lapy).l;
431 tdm_c[j] = 1. - prefactor * (*lapy).c;
432 tdm_u[j] = - prefactor * (*lapy).u;
433 }
434 // solve all
435 tdm.solve(tdm_info, linear_system.y1pncl);
436 // recover original memory alignment
437 sdecomp.transpose.execute(
438 linear_system.transposer_y1_to_x1,
439 linear_system.y1pncl,
440 linear_system.x1pncl
441 );
442}
src/fluid/predict/t.c¶
445if (param_t_implicit_z) {
446 // make all z data accessible
447 sdecomp.transpose.execute(
448 linear_system.transposer_x1_to_z2,
449 linear_system.x1pncl,
450 linear_system.z2pncl
451 );
452 // prepare tri-diagonal coefficients
453 tdm_info_t * tdm_info = linear_system.tdm_z;
454 int size = 0;
455 double * restrict tdm_l = NULL;
456 double * restrict tdm_c = NULL;
457 double * restrict tdm_u = NULL;
458 tdm.get_size(tdm_info, &size);
459 tdm.get_l(tdm_info, &tdm_l);
460 tdm.get_c(tdm_info, &tdm_c);
461 tdm.get_u(tdm_info, &tdm_u);
462 const laplacian_t * lapz = &laplacians.lapz;
463 for (int k = 0; k < size; k++) {
464 tdm_l[k] = - prefactor * (*lapz).l;
465 tdm_c[k] = 1. - prefactor * (*lapz).c;
466 tdm_u[k] = - prefactor * (*lapz).u;
467 }
468 // solve all
469 tdm.solve(tdm_info, linear_system.z2pncl);
470 // recover original memory alignment
471 sdecomp.transpose.execute(
472 linear_system.transposer_z2_to_x1,
473 linear_system.z2pncl,
474 linear_system.x1pncl
475 );
476}
We need to solve tri-diagonal matrix.