PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Functions
rhs.h File Reference
#include "variables.h"
#include "logging.h"
#include "Metric.h"
#include "BodyForces.h"
Include dependency graph for rhs.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

PetscErrorCode Viscous (UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
 Computes the viscous contribution to the contravariant momentum RHS.
 
PetscErrorCode Convection (UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
 Computes the convective contribution to the contravariant momentum RHS.
 
PetscErrorCode ComputeBodyForces (UserCtx *user, Vec Rct)
 General dispatcher for applying all active body forces (momentum sources).
 
PetscErrorCode ComputeRHS (UserCtx *user, Vec Rhs)
 Computes the Right-Hand Side (RHS) of the momentum equations.
 
PetscErrorCode ComputeEulerianDiffusivity (UserCtx *user)
 Computes the effective diffusivity scalar field (Gamma_eff) on the Eulerian grid.
 
PetscErrorCode ComputeEulerianDiffusivityGradient (UserCtx *user)
 Computes the Eulerian gradient of the effective diffusivity field.
 

Function Documentation

◆ Viscous()

PetscErrorCode Viscous ( UserCtx *  user,
Vec  Ucont,
Vec  Ucat,
Vec  Visc 
)

Computes the viscous contribution to the contravariant momentum RHS.

This routine evaluates diffusive fluxes on the curvilinear grid and writes the resulting term into Visc. The caller is responsible for providing compatible vectors and for assembling any additional source terms afterwards.

Parameters
[in]userBlock-level solver context containing metrics and model parameters.
[in]UcontContravariant velocity field used by the discretization.
[in]UcatCartesian velocity field used for derivative evaluation.
[out]ViscOutput vector receiving the viscous RHS contribution.
Returns
PetscErrorCode 0 on success.

Computes the viscous contribution to the contravariant momentum RHS.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/rhs.h.

See also
Viscous()

Definition at line 519 of file rhs.c.

520{
521
522 Vec Csi = user->lCsi, Eta = user->lEta, Zet = user->lZet;
523
524 Cmpnts ***ucont, ***ucat;
525
526 Cmpnts ***csi, ***eta, ***zet;
527 Cmpnts ***icsi, ***ieta, ***izet;
528 Cmpnts ***jcsi, ***jeta, ***jzet;
529 Cmpnts ***kcsi, ***keta, ***kzet;
530
531 PetscReal ***nvert;
532
533 DM da = user->da, fda = user->fda;
534 DMDALocalInfo info;
535 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
536 PetscInt mx, my, mz; // Dimensions in three directions
537 PetscInt i, j, k;
538 Vec Fp1, Fp2, Fp3;
539 Cmpnts ***fp1, ***fp2, ***fp3;
540 Cmpnts ***visc;
541 PetscReal ***aj, ***iaj, ***jaj, ***kaj;
542
543 PetscInt lxs, lxe, lys, lye, lzs, lze;
544
545 PetscReal ajc;
546
547 PetscReal dudc, dude, dudz, dvdc, dvde, dvdz, dwdc, dwde, dwdz;
548 PetscReal csi0, csi1, csi2, eta0, eta1, eta2, zet0, zet1, zet2;
549 PetscReal g11, g21, g31;
550 PetscReal r11, r21, r31, r12, r22, r32, r13, r23, r33;
551
552 PetscScalar solid,innerblank;
553
554 // --- CONTEXT ACQUISITION BLOCK ---
555 // Get the master simulation context from the UserCtx.
556 SimCtx *simCtx = user->simCtx;
557
558 // Create local variables to mirror the legacy globals for minimal code changes.
559 const LESModelType les = simCtx->les;
560 const PetscInt rans = simCtx->rans;
561 const PetscInt ti = simCtx->step; // Assuming simCtx->step is the new integer time counter
562 const PetscReal ren = simCtx->ren;
563 const PetscInt clark = simCtx->clark;
564 solid = 0.5;
565 innerblank = 7.;
566
567 PROFILE_FUNCTION_BEGIN;
568
569 DMDAVecGetArray(fda, Ucont, &ucont);
570 DMDAVecGetArray(fda, Ucat, &ucat);
571 DMDAVecGetArray(fda, Visc, &visc);
572
573 DMDAVecGetArray(fda, Csi, &csi);
574 DMDAVecGetArray(fda, Eta, &eta);
575 DMDAVecGetArray(fda, Zet, &zet);
576
577 DMDAVecGetArray(fda, user->lICsi, &icsi);
578 DMDAVecGetArray(fda, user->lIEta, &ieta);
579 DMDAVecGetArray(fda, user->lIZet, &izet);
580
581 DMDAVecGetArray(fda, user->lJCsi, &jcsi);
582 DMDAVecGetArray(fda, user->lJEta, &jeta);
583 DMDAVecGetArray(fda, user->lJZet, &jzet);
584
585 DMDAVecGetArray(fda, user->lKCsi, &kcsi);
586 DMDAVecGetArray(fda, user->lKEta, &keta);
587 DMDAVecGetArray(fda, user->lKZet, &kzet);
588
589 DMDAVecGetArray(da, user->lNvert, &nvert);
590
591 VecDuplicate(Ucont, &Fp1);
592 VecDuplicate(Ucont, &Fp2);
593 VecDuplicate(Ucont, &Fp3);
594
595 DMDAVecGetArray(fda, Fp1, &fp1);
596 DMDAVecGetArray(fda, Fp2, &fp2);
597 DMDAVecGetArray(fda, Fp3, &fp3);
598
599 DMDAVecGetArray(da, user->lAj, &aj);
600
601 DMDAGetLocalInfo(da, &info);
602
603 mx = info.mx; my = info.my; mz = info.mz;
604 xs = info.xs; xe = xs + info.xm;
605 ys = info.ys; ye = ys + info.ym;
606 zs = info.zs; ze = zs + info.zm;
607
608 /* First we calculate the flux on cell surfaces. Stored on the upper integer
609 node. For example, along i direction, the flux are stored at node 0:mx-2*/
610 lxs = xs; lxe = xe;
611 lys = ys; lye = ye;
612 lzs = zs; lze = ze;
613
614 if (xs==0) lxs = xs+1;
615 if (ys==0) lys = ys+1;
616 if (zs==0) lzs = zs+1;
617
618
619 if (xe==mx) lxe=xe-1;
620 if (ye==my) lye=ye-1;
621 if (ze==mz) lze=ze-1;
622
623 VecSet(Visc,0.0);
624
625 PetscReal ***lnu_t;
626
627 if(les) {
628 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
629 } else if (rans) {
630
631 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
632 }
633
634 /* The visc flux on each surface center is stored at previous integer node */
635
636 DMDAVecGetArray(da, user->lIAj, &iaj);
637 /* for (k=zs; k<ze; k++) { */
638/* for (j=ys; j<ye; j++) { */
639/* for (i=xs; i<xe; i++) { */
640/* if (i==1 && (j==0 ||j==1 || j==2) && (k==21 || k==22|| k==20)) */
641/* PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d u is %.15le v is %.15le w is %.15le \n",i,j,k,ucat[k][j][i].x,ucat[k][j][i].y,ucat[k][j][i].z ); */
642/* } */
643/* } */
644/* } */
645 // i direction
646 for (k=lzs; k<lze; k++) {
647 for (j=lys; j<lye; j++) {
648 for (i=lxs-1; i<lxe; i++) {
649
650 dudc = ucat[k][j][i+1].x - ucat[k][j][i].x;
651 dvdc = ucat[k][j][i+1].y - ucat[k][j][i].y;
652 dwdc = ucat[k][j][i+1].z - ucat[k][j][i].z;
653
654 if ((nvert[k][j+1][i ]> solid && nvert[k][j+1][i ]<innerblank) ||
655 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
656 dude = (ucat[k][j ][i+1].x + ucat[k][j ][i].x -
657 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.5;
658 dvde = (ucat[k][j ][i+1].y + ucat[k][j ][i].y -
659 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.5;
660 dwde = (ucat[k][j ][i+1].z + ucat[k][j ][i].z -
661 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.5;
662 }
663 else if ((nvert[k][j-1][i ]> solid && nvert[k][j-1][i ]<innerblank) ||
664 (nvert[k][j-1][i+1]> solid && nvert[k][j-1][i+1]<innerblank)) {
665 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
666 ucat[k][j ][i+1].x - ucat[k][j ][i].x) * 0.5;
667 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
668 ucat[k][j ][i+1].y - ucat[k][j ][i].y) * 0.5;
669 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
670 ucat[k][j ][i+1].z - ucat[k][j ][i].z) * 0.5;
671 }
672 else {
673 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
674 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.25;
675 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
676 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.25;
677 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
678 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.25;
679 }
680
681 if ((nvert[k+1][j][i ]> solid && nvert[k+1][j][i ]<innerblank)||
682 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
683 dudz = (ucat[k ][j][i+1].x + ucat[k ][j][i].x -
684 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.5;
685 dvdz = (ucat[k ][j][i+1].y + ucat[k ][j][i].y -
686 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.5;
687 dwdz = (ucat[k ][j][i+1].z + ucat[k ][j][i].z -
688 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.5;
689 }
690 else if ((nvert[k-1][j][i ]> solid && nvert[k-1][j][i ]<innerblank) ||
691 (nvert[k-1][j][i+1]> solid && nvert[k-1][j][i+1]<innerblank)) {
692
693 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
694 ucat[k ][j][i+1].x - ucat[k ][j][i].x) * 0.5;
695 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
696 ucat[k ][j][i+1].y - ucat[k ][j][i].y) * 0.5;
697 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
698 ucat[k ][j][i+1].z - ucat[k ][j][i].z) * 0.5;
699 }
700 else {
701 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
702 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.25;
703 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
704 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.25;
705 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
706 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.25;
707 }
708
709 csi0 = icsi[k][j][i].x;
710 csi1 = icsi[k][j][i].y;
711 csi2 = icsi[k][j][i].z;
712
713 eta0 = ieta[k][j][i].x;
714 eta1 = ieta[k][j][i].y;
715 eta2 = ieta[k][j][i].z;
716
717 zet0 = izet[k][j][i].x;
718 zet1 = izet[k][j][i].y;
719 zet2 = izet[k][j][i].z;
720
721 g11 = csi0 * csi0 + csi1 * csi1 + csi2 * csi2;
722 g21 = eta0 * csi0 + eta1 * csi1 + eta2 * csi2;
723 g31 = zet0 * csi0 + zet1 * csi1 + zet2 * csi2;
724
725 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
726 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
727 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
728
729 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
730 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
731 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
732
733 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
734 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
735 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
736
737 ajc = iaj[k][j][i];
738
739 double nu = 1./ren, nu_t=0;
740
741 if( les || (rans && ti>0) ) {
742 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j][i+1]) ), 2.0) * Sabs;
743 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j][i+1]);
744 if ( (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == WALL && i==0) || (user->boundary_faces[BC_FACE_POS_X].mathematical_type == WALL && i==mx-2) ) nu_t=0;
745 fp1[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * csi0 + r21 * csi1 + r31 * csi2) * ajc * (nu_t);
746 fp1[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * csi0 + r22 * csi1 + r32 * csi2) * ajc * (nu_t);
747 fp1[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * csi0 + r23 * csi1 + r33 * csi2) * ajc * (nu_t);
748 }
749 else {
750 fp1[k][j][i].x = 0;
751 fp1[k][j][i].y = 0;
752 fp1[k][j][i].z = 0;
753 }
754
755 fp1[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * csi0 + r21 * csi1 + r31 * csi2 ) * ajc * (nu);
756 fp1[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * csi0 + r22 * csi1 + r32 * csi2 ) * ajc * (nu);
757 fp1[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * csi0 + r23 * csi1 + r33 * csi2 ) * ajc * (nu);
758
759
760 if(clark) {
761 double dc, de, dz;
762 ComputeCellCharacteristicLengthScale (ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
763 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
764
765 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
766 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
767 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
768 double t21 = t12;
769 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
770 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
771 double t31 = t13;
772 double t32 = t23;
773 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
774
775 fp1[k][j][i].x -= ( t11 * csi0 + t12 * csi1 + t13 * csi2 ) / 12.;
776 fp1[k][j][i].y -= ( t21 * csi0 + t22 * csi1 + t23 * csi2 ) / 12.;
777 fp1[k][j][i].z -= ( t31 * csi0 + t32 * csi1 + t33 * csi2 ) / 12.;
778 }
779
780 }
781 }
782 }
783 DMDAVecRestoreArray(da, user->lIAj, &iaj);
784
785
786 // j direction
787 DMDAVecGetArray(da, user->lJAj, &jaj);
788 for (k=lzs; k<lze; k++) {
789 for (j=lys-1; j<lye; j++) {
790 for (i=lxs; i<lxe; i++) {
791
792 if ((nvert[k][j ][i+1]> solid && nvert[k][j ][i+1]<innerblank)||
793 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
794 dudc = (ucat[k][j+1][i ].x + ucat[k][j][i ].x -
795 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.5;
796 dvdc = (ucat[k][j+1][i ].y + ucat[k][j][i ].y -
797 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.5;
798 dwdc = (ucat[k][j+1][i ].z + ucat[k][j][i ].z -
799 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.5;
800 }
801 else if ((nvert[k][j ][i-1]> solid && nvert[k][j ][i-1]<innerblank) ||
802 (nvert[k][j+1][i-1]> solid && nvert[k][j+1][i-1]<innerblank)) {
803 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
804 ucat[k][j+1][i ].x - ucat[k][j][i ].x) * 0.5;
805 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
806 ucat[k][j+1][i ].y - ucat[k][j][i ].y) * 0.5;
807 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
808 ucat[k][j+1][i ].z - ucat[k][j][i ].z) * 0.5;
809 }
810 else {
811 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
812 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.25;
813 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
814 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.25;
815 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
816 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.25;
817 }
818
819 dude = ucat[k][j+1][i].x - ucat[k][j][i].x;
820 dvde = ucat[k][j+1][i].y - ucat[k][j][i].y;
821 dwde = ucat[k][j+1][i].z - ucat[k][j][i].z;
822
823 if ((nvert[k+1][j ][i]> solid && nvert[k+1][j ][i]<innerblank)||
824 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
825 dudz = (ucat[k ][j+1][i].x + ucat[k ][j][i].x -
826 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.5;
827 dvdz = (ucat[k ][j+1][i].y + ucat[k ][j][i].y -
828 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.5;
829 dwdz = (ucat[k ][j+1][i].z + ucat[k ][j][i].z -
830 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.5;
831 }
832 else if ((nvert[k-1][j ][i]> solid && nvert[k-1][j ][i]<innerblank)||
833 (nvert[k-1][j+1][i]> solid && nvert[k-1][j+1][i]<innerblank)) {
834 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
835 ucat[k ][j+1][i].x - ucat[k ][j][i].x) * 0.5;
836 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
837 ucat[k ][j+1][i].y - ucat[k ][j][i].y) * 0.5;
838 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
839 ucat[k ][j+1][i].z - ucat[k ][j][i].z) * 0.5;
840 }
841 else {
842 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
843 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.25;
844 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
845 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.25;
846 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
847 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.25;
848 }
849
850 csi0 = jcsi[k][j][i].x;
851 csi1 = jcsi[k][j][i].y;
852 csi2 = jcsi[k][j][i].z;
853
854 eta0 = jeta[k][j][i].x;
855 eta1 = jeta[k][j][i].y;
856 eta2 = jeta[k][j][i].z;
857
858 zet0 = jzet[k][j][i].x;
859 zet1 = jzet[k][j][i].y;
860 zet2 = jzet[k][j][i].z;
861
862
863 g11 = csi0 * eta0 + csi1 * eta1 + csi2 * eta2;
864 g21 = eta0 * eta0 + eta1 * eta1 + eta2 * eta2;
865 g31 = zet0 * eta0 + zet1 * eta1 + zet2 * eta2;
866
867 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
868 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
869 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
870
871 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
872 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
873 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
874
875 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
876 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
877 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
878
879 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d dvdc is %.15le dvde is %.15le dvdz is %.15le \n",i,j,k,dvdc,dvde,dvdz);
880 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d dwdc is %.15le dwde is %.15le dwdz is %.15le \n",i,j,k,dwdc,dwde,dwdz);
881 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d jcsi is %.15le jeta is %.15le jzet is %.15le \n",i,j,k,jcsi[k][j][i].z,jeta[k][j][i].z,jzet[k][j][i].z);
882 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d r13 is %.15le r23 is %.15le r33 is %.15le \n",i,j,k,r13,r23,r33);
883
884
885
886 ajc = jaj[k][j][i];
887
888 double nu = 1./ren, nu_t = 0;
889
890 if( les || (rans && ti>0) ) {
891 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j+1][i]) ), 2.0) * Sabs;
892 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j+1][i]);
893 if ( (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == WALL && j==0) || (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == WALL && j==my-2) ) nu_t=0;
894
895 fp2[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * eta0 + r21 * eta1 + r31 * eta2) * ajc * (nu_t);
896 fp2[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * eta0 + r22 * eta1 + r32 * eta2) * ajc * (nu_t);
897 fp2[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * eta0 + r23 * eta1 + r33 * eta2) * ajc * (nu_t);
898 }
899 else {
900 fp2[k][j][i].x = 0;
901 fp2[k][j][i].y = 0;
902 fp2[k][j][i].z = 0;
903 }
904
905 fp2[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * eta0 + r21 * eta1 + r31 * eta2 ) * ajc * (nu);
906 fp2[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * eta0 + r22 * eta1 + r32 * eta2 ) * ajc * (nu);
907 fp2[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * eta0 + r23 * eta1 + r33 * eta2 ) * ajc * (nu);
908
909 if(clark) {
910 double dc, de, dz;
911 ComputeCellCharacteristicLengthScale(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
912 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
913
914 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
915 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
916 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
917 double t21 = t12;
918 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
919 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
920 double t31 = t13;
921 double t32 = t23;
922 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
923
924 fp2[k][j][i].x -= ( t11 * eta0 + t12 * eta1 + t13 * eta2 ) / 12.;
925 fp2[k][j][i].y -= ( t21 * eta0 + t22 * eta1 + t23 * eta2 ) / 12.;
926 fp2[k][j][i].z -= ( t31 * eta0 + t32 * eta1 + t33 * eta2 ) / 12.;
927 }
928 }
929 }
930 }
931
932 DMDAVecRestoreArray(da, user->lJAj, &jaj);
933 // k direction
934
935 DMDAVecGetArray(da, user->lKAj, &kaj);
936 for (k=lzs-1; k<lze; k++) {
937 for (j=lys; j<lye; j++) {
938 for (i=lxs; i<lxe; i++) {
939 if ((nvert[k ][j][i+1]> solid && nvert[k ][j][i+1]<innerblank)||
940 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
941 dudc = (ucat[k+1][j][i ].x + ucat[k][j][i ].x -
942 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.5;
943 dvdc = (ucat[k+1][j][i ].y + ucat[k][j][i ].y -
944 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.5;
945 dwdc = (ucat[k+1][j][i ].z + ucat[k][j][i ].z -
946 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.5;
947 }
948 else if ((nvert[k ][j][i-1]> solid && nvert[k ][j][i-1]<innerblank) ||
949 (nvert[k+1][j][i-1]> solid && nvert[k+1][j][i-1]<innerblank)) {
950 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
951 ucat[k+1][j][i ].x - ucat[k][j][i ].x) * 0.5;
952 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
953 ucat[k+1][j][i ].y - ucat[k][j][i ].y) * 0.5;
954 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
955 ucat[k+1][j][i ].z - ucat[k][j][i ].z) * 0.5;
956 }
957 else {
958 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
959 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.25;
960 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
961 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.25;
962 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
963 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.25;
964 }
965
966 if ((nvert[k ][j+1][i]> solid && nvert[k ][j+1][i]<innerblank)||
967 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
968 dude = (ucat[k+1][j ][i].x + ucat[k][j ][i].x -
969 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.5;
970 dvde = (ucat[k+1][j ][i].y + ucat[k][j ][i].y -
971 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.5;
972 dwde = (ucat[k+1][j ][i].z + ucat[k][j ][i].z -
973 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.5;
974 }
975 else if ((nvert[k ][j-1][i]> solid && nvert[k ][j-1][i]<innerblank) ||
976 (nvert[k+1][j-1][i]> solid && nvert[k+1][j-1][i]<innerblank)){
977 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
978 ucat[k+1][j ][i].x - ucat[k][j ][i].x) * 0.5;
979 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
980 ucat[k+1][j ][i].y - ucat[k][j ][i].y) * 0.5;
981 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
982 ucat[k+1][j ][i].z - ucat[k][j ][i].z) * 0.5;
983 }
984 else {
985 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
986 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.25;
987 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
988 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.25;
989 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
990 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.25;
991 }
992
993 dudz = ucat[k+1][j][i].x - ucat[k][j][i].x;
994 dvdz = ucat[k+1][j][i].y - ucat[k][j][i].y;
995 dwdz = ucat[k+1][j][i].z - ucat[k][j][i].z;
996
997
998 csi0 = kcsi[k][j][i].x;
999 csi1 = kcsi[k][j][i].y;
1000 csi2 = kcsi[k][j][i].z;
1001
1002 eta0 = keta[k][j][i].x;
1003 eta1 = keta[k][j][i].y;
1004 eta2 = keta[k][j][i].z;
1005
1006 zet0 = kzet[k][j][i].x;
1007 zet1 = kzet[k][j][i].y;
1008 zet2 = kzet[k][j][i].z;
1009
1010
1011 g11 = csi0 * zet0 + csi1 * zet1 + csi2 * zet2;
1012 g21 = eta0 * zet0 + eta1 * zet1 + eta2 * zet2;
1013 g31 = zet0 * zet0 + zet1 * zet1 + zet2 * zet2;
1014
1015 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
1016 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
1017 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
1018
1019 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
1020 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
1021 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
1022
1023 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
1024 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
1025 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
1026
1027 ajc = kaj[k][j][i];
1028
1029 double nu = 1./ren, nu_t =0;
1030
1031 if( les || (rans && ti>0) ) {
1032 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k+1][j][i]) ), 2.0) * Sabs;
1033 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k+1][j][i]);
1034 if ( (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == WALL && k==0) || (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == WALL && k==mz-2) ) nu_t=0;
1035
1036 fp3[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu_t);
1037 fp3[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu_t);
1038 fp3[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu_t);
1039 }
1040 else {
1041 fp3[k][j][i].x = 0;
1042 fp3[k][j][i].y = 0;
1043 fp3[k][j][i].z = 0;
1044 }
1045 fp3[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu);//
1046 fp3[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu);//
1047 fp3[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu);//
1048
1049 if(clark) {
1050 double dc, de, dz;
1051 ComputeCellCharacteristicLengthScale(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
1052 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
1053
1054 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
1055 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
1056 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
1057 double t21 = t12;
1058 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
1059 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
1060 double t31 = t13;
1061 double t32 = t23;
1062 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
1063
1064 fp3[k][j][i].x -= ( t11 * zet0 + t12 * zet1 + t13 * zet2 ) / 12.;
1065 fp3[k][j][i].y -= ( t21 * zet0 + t22 * zet1 + t23 * zet2 ) / 12.;
1066 fp3[k][j][i].z -= ( t31 * zet0 + t32 * zet1 + t33 * zet2 ) / 12.;
1067 }
1068 }
1069 }
1070 }
1071
1072 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1073
1074 for (k=lzs; k<lze; k++) {
1075 for (j=lys; j<lye; j++) {
1076 for (i=lxs; i<lxe; i++) {
1077 visc[k][j][i].x =
1078 (fp1[k][j][i].x - fp1[k][j][i-1].x +
1079 fp2[k][j][i].x - fp2[k][j-1][i].x +
1080 fp3[k][j][i].x - fp3[k-1][j][i].x);
1081
1082 visc[k][j][i].y =
1083 (fp1[k][j][i].y - fp1[k][j][i-1].y +
1084 fp2[k][j][i].y - fp2[k][j-1][i].y +
1085 fp3[k][j][i].y - fp3[k-1][j][i].y);
1086
1087 visc[k][j][i].z =
1088 (fp1[k][j][i].z - fp1[k][j][i-1].z +
1089 fp2[k][j][i].z - fp2[k][j-1][i].z +
1090 fp3[k][j][i].z - fp3[k-1][j][i].z);
1091
1092 }
1093 }
1094 }
1095/* for (k=zs; k<ze; k++) { */
1096/* for (j=ys; j<ye; j++) { */
1097/* for (i=xs; i<xe; i++) { */
1098/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp1.z is %.15le \n",i,j,k,fp1[k][j][i].z); */
1099/* if (i==0 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp1.z is %.15le \n",i,j,k,fp1[k][j][i].z); */
1100/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp2.z is %.15le \n",i,j,k,fp2[k][j][i].z); */
1101/* if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp2.z is %.15le \n",i,j,k,fp2[k][j][i].z); */
1102/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp3.z is %.15le \n",i,j,k,fp3[k][j][i].z); */
1103/* if (i==1 && j==1 && k==20) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp3.z is %.15le \n",i,j,k,fp3[k][j][i].z); */
1104
1105/* } */
1106/* } */
1107/* } */
1108 DMDAVecRestoreArray(fda, Ucont, &ucont);
1109 DMDAVecRestoreArray(fda, Ucat, &ucat);
1110 DMDAVecRestoreArray(fda, Visc, &visc);
1111
1112 DMDAVecRestoreArray(fda, Csi, &csi);
1113 DMDAVecRestoreArray(fda, Eta, &eta);
1114 DMDAVecRestoreArray(fda, Zet, &zet);
1115
1116 DMDAVecRestoreArray(fda, Fp1, &fp1);
1117 DMDAVecRestoreArray(fda, Fp2, &fp2);
1118 DMDAVecRestoreArray(fda, Fp3, &fp3);
1119
1120 DMDAVecRestoreArray(da, user->lAj, &aj);
1121
1122 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1123 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1124 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1125
1126 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1127 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1128 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1129
1130 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1131 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1132 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1133
1134 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1135
1136 if(les) {
1137 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1138 } else if (rans) {
1139
1140 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1141 }
1142
1143
1144 VecDestroy(&Fp1);
1145 VecDestroy(&Fp2);
1146 VecDestroy(&Fp3);
1147
1148
1149 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Viscous terms calculated .\n");
1150
1151 PROFILE_FUNCTION_END;
1152
1153 return(0);
1154}
ComputeCellCharacteristicLengthScale
PetscErrorCode ComputeCellCharacteristicLengthScale(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
Computes characteristic length scales (dx, dy, dz) for a curvilinear cell.
Definition Metric.c:283
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:199
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:779
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:31
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:770
LESModelType
LESModelType
Identifies the six logical faces of a structured computational block.
Definition variables.h:488
SimCtx::clark
PetscInt clark
Definition variables.h:733
WALL
@ WALL
Definition variables.h:254
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:829
UserCtx::lIEta
Vec lIEta
Definition variables.h:860
UserCtx::lIZet
Vec lIZet
Definition variables.h:860
UserCtx::lNvert
Vec lNvert
Definition variables.h:837
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:814
SimCtx::rans
PetscInt rans
Definition variables.h:732
UserCtx::lZet
Vec lZet
Definition variables.h:858
SimCtx::ren
PetscReal ren
Definition variables.h:691
UserCtx::lIAj
Vec lIAj
Definition variables.h:860
UserCtx::lKEta
Vec lKEta
Definition variables.h:862
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:861
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::lKZet
Vec lKZet
Definition variables.h:862
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:865
UserCtx::lJEta
Vec lJEta
Definition variables.h:861
UserCtx::lCsi
Vec lCsi
Definition variables.h:858
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:862
UserCtx::lJZet
Vec lJZet
Definition variables.h:861
SimCtx::step
PetscInt step
Definition variables.h:651
UserCtx::lAj
Vec lAj
Definition variables.h:858
UserCtx::lICsi
Vec lICsi
Definition variables.h:860
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:858
SimCtx::les
PetscInt les
Definition variables.h:732
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:336
UserCtx::lJAj
Vec lJAj
Definition variables.h:861
UserCtx::lKAj
Vec lKAj
Definition variables.h:862
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:245
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:247
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:246
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:247
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:245
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:246
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:643
nu_t
double nu_t(double yplus)
Computes turbulent eddy viscosity ratio (ν_t / ν)
Definition wallfunction.c:303
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Convection()

PetscErrorCode Convection ( UserCtx *  user,
Vec  Ucont,
Vec  Ucat,
Vec  Conv 
)

Computes the convective contribution to the contravariant momentum RHS.

This routine evaluates the advection operator on the current velocity state and stores the contribution in Conv for subsequent combination with viscous and body-force terms.

Parameters
[in]userBlock-level solver context containing metrics and numerics settings.
[in]UcontContravariant velocity field used in face-normal flux construction.
[in]UcatCartesian velocity field used by the convective stencil.
[out]ConvOutput vector receiving the convection RHS contribution.
Returns
PetscErrorCode 0 on success.

Computes the convective contribution to the contravariant momentum RHS.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/rhs.h.

See also
Convection()

Definition at line 13 of file rhs.c.

14{
15
16 // --- CONTEXT ACQUISITION BLOCK ---
17 // Get the master simulation context from the UserCtx.
18 SimCtx *simCtx = user->simCtx;
19
20 // Create local variables to mirror the legacy globals for minimal code changes.
21 const LESModelType les = simCtx->les;
22 const PetscInt central = simCtx->central; // Get this from SimCtx now
23 // --- END CONTEXT ACQUISITION BLOCK ---
24
25 Cmpnts ***ucont, ***ucat;
26 DM da = user->da, fda = user->fda;
27 DMDALocalInfo info;
28 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
29 PetscInt mx, my, mz; // Dimensions in three directions
30 PetscInt i, j, k;
31 Vec Fp1, Fp2, Fp3;
32 Cmpnts ***fp1, ***fp2, ***fp3;
33 Cmpnts ***conv;
34
35 PetscReal ucon, up, um;
36 PetscReal coef = 0.125, innerblank=7.;
37
38 PetscInt lxs, lxe, lys, lye, lzs, lze, gxs, gxe, gys, gye, gzs,gze;
39
40 PetscReal ***nvert,***aj;
41
42 PROFILE_FUNCTION_BEGIN;
43
44 DMDAGetLocalInfo(da, &info);
45 mx = info.mx; my = info.my; mz = info.mz;
46 xs = info.xs; xe = xs + info.xm;
47 ys = info.ys; ye = ys + info.ym;
48 zs = info.zs; ze = zs + info.zm;
49 gxs = info.gxs; gxe = gxs + info.gxm;
50 gys = info.gys; gye = gys + info.gym;
51 gzs = info.gzs; gze = gzs + info.gzm;
52
53 DMDAVecGetArray(fda, Ucont, &ucont);
54 DMDAVecGetArray(fda, Ucat, &ucat);
55 DMDAVecGetArray(fda, Conv, &conv);
56 DMDAVecGetArray(da, user->lAj, &aj);
57
58 VecDuplicate(Ucont, &Fp1);
59 VecDuplicate(Ucont, &Fp2);
60 VecDuplicate(Ucont, &Fp3);
61
62 DMDAVecGetArray(fda, Fp1, &fp1);
63 DMDAVecGetArray(fda, Fp2, &fp2);
64 DMDAVecGetArray(fda, Fp3, &fp3);
65
66 DMDAVecGetArray(da, user->lNvert, &nvert);
67
68
69 /* We have two different sets of node: 1. grid node, the physical points
70 where grid lines intercross; 2. storage node, where we store variables.
71 All node without explicitly specified as "grid node" refers to
72 storage node.
73
74 The integer node is defined at cell center while half node refers to
75 the actual grid node. (The reason to choose this arrangement is we need
76 ghost node, which is half node away from boundaries, to specify boundary
77 conditions. By using this storage arrangement, the actual storage need
78 is (IM+1) * (JM + 1) * (KM+1) where IM, JM, & KM refer to the number of
79 grid nodes along i, j, k directions.)
80
81 DA, the data structure used to define the storage of 3D arrays, is defined
82 as mx * my * mz. mx = IM+1, my = JM+1, mz = KM+1.
83
84 Staggered grid arrangement is used in this solver.
85 Pressure is stored at interger node (hence the cell center) and volume
86 fluxes defined on the center of each surface of a given control volume
87 is stored on the cloest upper integer node. */
88
89 /* First we calculate the flux on cell surfaces. Stored on the upper integer
90 node. For example, along i direction, the flux are stored at node 0:mx-2*/
91
92 lxs = xs; lxe = xe;
93 lys = ys; lye = ye;
94 lzs = zs; lze = ze;
95
96 if (xs==0) lxs = xs+1;
97 if (ys==0) lys = ys+1;
98 if (zs==0) lzs = zs+1;
99
100 if (xe==mx) lxe=xe-1;
101 if (ye==my) lye=ye-1;
102 if (ze==mz) lze=ze-1;
103
104 //Mohsen Sep 2012//
105/* First update the computational ghost points velocity for periodic boundary conditions
106 just for this subroutine because of Quick scheme for velocity deravatives */
107 /* for (k=zs; k<ze; k++) { */
108/* for (j=ys; j<ye; j++) { */
109/* for (i=xs; i<xe; i++) { */
110/* if (i==1 && (j==1) && (k==0 || k==1)) */
111/* PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d u is %.15le v is %.15le w is %.15le \n",i,j,k,ucat[k][j][i].x,ucat[k][j][i].y,ucat[k][j][i].z ); */
112/* } */
113/* } */
114/* } */
115 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
116 if (xs==0){
117 i=xs-1;
118 for (k=gzs; k<gze; k++) {
119 for (j=gys; j<gye; j++) {
120 ucat[k][j][i].x=ucat[k][j][i-2].x;
121 ucat[k][j][i].y=ucat[k][j][i-2].y;
122 ucat[k][j][i].z=ucat[k][j][i-2].z;
123 nvert[k][j][i]=nvert[k][j][i-2];
124 }
125 }
126 }
127 if (xe==mx){
128 i=mx;
129 for (k=gzs; k<gze; k++) {
130 for (j=gys; j<gye; j++) {
131 ucat[k][j][i].x=ucat[k][j][i+2].x;
132 ucat[k][j][i].y=ucat[k][j][i+2].y;
133 ucat[k][j][i].z=ucat[k][j][i+2].z;
134 nvert[k][j][i]=nvert[k][j][i+2];
135 }
136 }
137 }
138 }
139 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
140 if (ys==0){
141 j=ys-1;
142 for (k=gzs; k<gze; k++) {
143 for (i=gxs; i<gxe; i++) {
144 ucat[k][j][i].x=ucat[k][j-2][i].x;
145 ucat[k][j][i].y=ucat[k][j-2][i].y;
146 ucat[k][j][i].z=ucat[k][j-2][i].z;
147 nvert[k][j][i]=nvert[k][j-2][i];
148 }
149 }
150 }
151 if (ye==my){
152 j=my;
153 for (k=gzs; k<gze; k++) {
154 for (i=gxs; i<gxe; i++) {
155 ucat[k][j][i].x=ucat[k][j+2][i].x;
156 ucat[k][j][i].y=ucat[k][j+2][i].y;
157 ucat[k][j][i].z=ucat[k][j+2][i].z;
158 nvert[k][j][i]=nvert[k][j+2][i];
159 }
160 }
161 }
162 }
163 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
164 if (zs==0){
165 k=zs-1;
166 for (j=gys; j<gye; j++) {
167 for (i=gxs; i<gxe; i++) {
168 ucat[k][j][i].x=ucat[k-2][j][i].x;
169 ucat[k][j][i].y=ucat[k-2][j][i].y;
170 ucat[k][j][i].z=ucat[k-2][j][i].z;
171 nvert[k][j][i]=nvert[k-2][j][i];
172 }
173 }
174 }
175 if (ze==mz){
176 k=mz;
177 for (j=gys; j<gye; j++) {
178 for (i=gxs; i<gxe; i++) {
179 ucat[k][j][i].x=ucat[k+2][j][i].x;
180 ucat[k][j][i].y=ucat[k+2][j][i].y;
181 ucat[k][j][i].z=ucat[k+2][j][i].z;
182 nvert[k][j][i]=nvert[k+2][j][i];
183 }
184 }
185 }
186 }
187
188 VecSet(Conv, 0.0);
189
190 /* Calculating the convective terms on cell centers.
191 First calcualte the contribution from i direction
192 The flux is evaluated by QUICK scheme */
193
194 for (k=lzs; k<lze; k++){
195 for (j=lys; j<lye; j++){
196 for (i=lxs-1; i<lxe; i++){
197
198
199 ucon = ucont[k][j][i].x * 0.5;
200
201 up = ucon + fabs(ucon);
202 um = ucon - fabs(ucon);
203
204 if (i>0 && i<mx-2 &&
205 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
206 (nvert[k][j][i-1] < 0.1 || nvert[k][j][i-1]>innerblank)) { // interial nodes
207 if ((les || central)) {
208 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
209 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
210 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
211
212 } else {
213 fp1[k][j][i].x =
214 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
215 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
216 fp1[k][j][i].y =
217 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
218 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
219 fp1[k][j][i].z =
220 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
221 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
222 }
223 }
224 else if ((les || central) && (i==0 || i==mx-2) &&
225 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
226 (nvert[k][j][i ] < 0.1 || nvert[k][j][i ]>innerblank))
227 {
228 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
229 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
230 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
231 }
232 else if (i==0 ||(nvert[k][j][i-1] > 0.1) ) {
233 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==0 && (nvert[k][j][i-1]<0.1 && nvert[k][j][i+1]<0.1)){//Mohsen Feb 12
234 fp1[k][j][i].x =
235 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
236 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
237 fp1[k][j][i].y =
238 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
239 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
240 fp1[k][j][i].z =
241 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
242 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
243 }else{
244 fp1[k][j][i].x =
245 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
246 up * (coef * (-ucat[k][j][i ].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
247 fp1[k][j][i].y =
248 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
249 up * (coef * (-ucat[k][j][i ].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
250 fp1[k][j][i].z =
251 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
252 up * (coef * (-ucat[k][j][i ].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
253 }
254 }
255 else if (i==mx-2 ||(nvert[k][j][i+1]) > 0.1) {
256 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2 &&(nvert[k][j][i-1]<0.1 && nvert[k][j][i+1]<0.1)){//Mohsen Feb 12
257 fp1[k][j][i].x =
258 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
259 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
260 fp1[k][j][i].y =
261 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
262 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
263 fp1[k][j][i].z =
264 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
265 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
266 }else{
267 fp1[k][j][i].x =
268 um * (coef * (-ucat[k][j][i+1].x -2. * ucat[k][j][i+1].x +3. * ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
269 up * (coef * (-ucat[k][j][i-1].x -2. * ucat[k][j][i ].x +3. * ucat[k][j][i+1].x) +ucat[k][j][i ].x);
270 fp1[k][j][i].y =
271 um * (coef * (-ucat[k][j][i+1].y -2. * ucat[k][j][i+1].y +3. * ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
272 up * (coef * (-ucat[k][j][i-1].y -2. * ucat[k][j][i ].y +3. * ucat[k][j][i+1].y) +ucat[k][j][i ].y);
273 fp1[k][j][i].z =
274 um * (coef * (-ucat[k][j][i+1].z -2. * ucat[k][j][i+1].z +3. * ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
275 up * (coef * (-ucat[k][j][i-1].z -2. * ucat[k][j][i ].z +3. * ucat[k][j][i+1].z) +ucat[k][j][i ].z);
276 }
277 }
278 }
279 }
280 }
281
282 /* j direction */
283 for (k=lzs; k<lze; k++) {
284 for(j=lys-1; j<lye; j++) {
285 for(i=lxs; i<lxe; i++) {
286 ucon = ucont[k][j][i].y * 0.5;
287
288 up = ucon + fabs(ucon);
289 um = ucon - fabs(ucon);
290
291 if (j>0 && j<my-2 &&
292 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i] > innerblank) &&
293 (nvert[k][j-1][i] < 0.1 || nvert[k][j-1][i] > innerblank)) {
294 if ((les || central)) {
295 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
296 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
297 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
298
299 } else {
300 fp2[k][j][i].x =
301 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
302 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
303 fp2[k][j][i].y =
304 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
305 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
306 fp2[k][j][i].z =
307 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
308 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
309 }
310 }
311 else if ((les || central) && (j==0 || i==my-2) &&
312 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i]>innerblank) &&
313 (nvert[k][j ][i] < 0.1 || nvert[k][j ][i]>innerblank))
314 {
315 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
316 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
317 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
318 }
319 else if (j==0 || (nvert[k][j-1][i]) > 0.1) {
320 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==0 && (nvert[k][j-1][i]<0.1 && nvert[k][j+1][i]<0.1 )){//Mohsen Feb 12 //
321 fp2[k][j][i].x =
322 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
323 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
324 fp2[k][j][i].y =
325 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
326 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
327 fp2[k][j][i].z =
328 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
329 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
330 }else{
331 fp2[k][j][i].x =
332 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
333 up * (coef * (-ucat[k][j ][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
334 fp2[k][j][i].y =
335 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
336 up * (coef * (-ucat[k][j ][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
337 fp2[k][j][i].z =
338 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
339 up * (coef * (-ucat[k][j ][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
340 }
341 }
342 else if (j==my-2 ||(nvert[k][j+1][i]) > 0.1) {
343 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2 && (nvert[k][j-1][i]<0.1 && nvert[k][j+1][i]<0.1 )){//Mohsen Feb 12//
344 fp2[k][j][i].x =
345 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
346 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
347 fp2[k][j][i].y =
348 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
349 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
350 fp2[k][j][i].z =
351 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
352 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
353 }else{
354 fp2[k][j][i].x =
355 um * (coef * (-ucat[k][j+1][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
356 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
357 fp2[k][j][i].y =
358 um * (coef * (-ucat[k][j+1][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
359 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j][i ].y);
360 fp2[k][j][i].z =
361 um * (coef * (-ucat[k][j+1][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
362 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j][i ].z);
363 }
364 }
365 }
366 }
367 }
368
369
370 /* k direction */
371 for (k=lzs-1; k<lze; k++) {
372 for(j=lys; j<lye; j++) {
373 for(i=lxs; i<lxe; i++) {
374 ucon = ucont[k][j][i].z * 0.5;
375
376 up = ucon + fabs(ucon);
377 um = ucon - fabs(ucon);
378
379 if (k>0 && k<mz-2 &&
380 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i] > innerblank) &&
381 (nvert[k-1][j][i] < 0.1 || nvert[k-1][j][i] > innerblank)) {
382 if ((les || central)) {
383 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
384 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
385 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
386
387 } else {
388 fp3[k][j][i].x =
389 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
390 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
391 fp3[k][j][i].y =
392 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
393 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
394 fp3[k][j][i].z =
395 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
396 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
397 }
398 }
399 else if ((les || central) && (k==0 || k==mz-2) &&
400 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i]>innerblank) &&
401 (nvert[k ][j][i] < 0.1 || nvert[k ][j][i]>innerblank))
402 {
403 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
404 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
405 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
406 }
407 else if (k<mz-2 && (k==0 ||(nvert[k-1][j][i]) > 0.1)) {
408 if(user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==0 && (nvert[k-1][j][i]<0.1 && nvert[k+1][j][i]<0.1)){//Mohsen Feb 12//
409 fp3[k][j][i].x =
410 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
411 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
412 fp3[k][j][i].y =
413 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
414 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
415 fp3[k][j][i].z =
416 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
417 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
418 }else{
419 fp3[k][j][i].x =
420 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
421 up * (coef * (-ucat[k ][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k][j][i ].x);
422 fp3[k][j][i].y =
423 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
424 up * (coef * (-ucat[k ][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k][j][i ].y);
425 fp3[k][j][i].z =
426 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
427 up * (coef * (-ucat[k ][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k][j][i ].z);
428 }
429 }
430 else if (k>0 && (k==mz-2 ||(nvert[k+1][j][i]) > 0.1)) {
431 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2 && (nvert[k-1][j][i]<0.1 && nvert[k+1][j][i]<0.1)){//Mohsen Feb 12//
432 fp3[k][j][i].x =
433 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
434 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
435 fp3[k][j][i].y =
436 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
437 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
438 fp3[k][j][i].z =
439 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
440 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
441 }else{
442 fp3[k][j][i].x =
443 um * (coef * (-ucat[k+1][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
444 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k][j][i ].x);
445 fp3[k][j][i].y =
446 um * (coef * (-ucat[k+1][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
447 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k][j][i ].y);
448 fp3[k][j][i].z =
449 um * (coef * (-ucat[k+1][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
450 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k][j][i ].z);
451 }
452 }
453 }
454 }
455 }
456
457 /* Calculate the convective terms under cartesian coordinates */
458
459 for (k=lzs; k<lze; k++) {
460 for (j=lys; j<lye; j++) {
461 for (i=lxs; i<lxe; i++) {
462 conv[k][j][i].x =
463 fp1[k][j][i].x - fp1[k][j][i-1].x +
464 fp2[k][j][i].x - fp2[k][j-1][i].x +
465 fp3[k][j][i].x - fp3[k-1][j][i].x;
466
467 conv[k][j][i].y =
468 fp1[k][j][i].y - fp1[k][j][i-1].y +
469 fp2[k][j][i].y - fp2[k][j-1][i].y +
470 fp3[k][j][i].y - fp3[k-1][j][i].y;
471
472 conv[k][j][i].z =
473 fp1[k][j][i].z - fp1[k][j][i-1].z +
474 fp2[k][j][i].z - fp2[k][j-1][i].z +
475 fp3[k][j][i].z - fp3[k-1][j][i].z;
476 }
477 }
478 }
479 /* for (k=zs; k<ze; k++) { */
480/* for (j=ys; j<ye; j++) { */
481/* for (i=xs; i<xe; i++) { */
482/* if (i==1 && (j==1) && (k==1 || k==21 || k==22|| k==200)) */
483/* PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d conv.y is %.15le conv.z is %.15le \n",i,j,k,conv[k][j][i].y,conv[k][j][i].z); */
484/* } */
485/* } */
486/* } */
487
488 DMDAVecRestoreArray(fda, Ucont, &ucont);
489 DMDAVecRestoreArray(fda, Ucat, &ucat);
490 DMDAVecRestoreArray(fda, Conv, &conv);
491 DMDAVecRestoreArray(da, user->lAj, &aj);
492
493 DMDAVecRestoreArray(fda, Fp1, &fp1);
494 DMDAVecRestoreArray(fda, Fp2, &fp2);
495 DMDAVecRestoreArray(fda, Fp3, &fp3);
496 DMDAVecRestoreArray(da, user->lNvert, &nvert);
497
498 VecDestroy(&Fp1);
499 VecDestroy(&Fp2);
500 VecDestroy(&Fp3);
501
502
503 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Convective term calculated .\n");
504
505 PROFILE_FUNCTION_END;
506 return (0);
507}
PERIODIC
@ PERIODIC
Definition variables.h:260
SimCtx::central
PetscInt central
Definition variables.h:689
Here is the caller graph for this function:

◆ ComputeBodyForces()

PetscErrorCode ComputeBodyForces ( UserCtx *  user,
Vec  Rct 
)

General dispatcher for applying all active body forces (momentum sources).

This function serves as a central hub for adding momentum source terms to the contravariant right-hand-side (Rct) of the momentum equations. It is called once per RHS evaluation (e.g., once per Runge-Kutta stage).

It introspects the simulation configuration to determine which, if any, body forces are active and calls their specific implementation functions.

Parameters
userThe UserCtx containing the simulation state for a single block.
RctThe PETSc Vec for the contravariant RHS, modified in place.
Returns
PetscErrorCode 0 on success.

General dispatcher for applying all active body forces (momentum sources).

Local to this translation unit.

Definition at line 1162 of file rhs.c.

1163{
1164 PetscErrorCode ierr;
1165 PetscFunctionBeginUser;
1166
1167 // --- 1. Apply momentum source for driven channel/pipe flows ---
1168 // This function will internally check if a driven flow BC is active.
1169 ierr = ComputeDrivenChannelFlowSource(user, Rct); CHKERRQ(ierr);
1170
1171 // --- 2. (Future Extension) Apply gravitational force ---
1172 // if (user->simCtx->gravityEnabled) {
1173 // ierr = ApplyGravitationalForce(user, Rhs); CHKERRQ(ierr);
1174 // }
1175
1176 PetscFunctionReturn(0);
1177}
ComputeDrivenChannelFlowSource
PetscErrorCode ComputeDrivenChannelFlowSource(UserCtx *user, Vec Rct)
Applies a momentum source term to drive flow in a periodic channel or pipe.
Definition BodyForces.c:14
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeRHS()

PetscErrorCode ComputeRHS ( UserCtx *  user,
Vec  Rhs 
)
extern

Computes the Right-Hand Side (RHS) of the momentum equations.

This function calculates the contribution of the convective and diffusive terms. It is called by the momentum solvers (e.g., RungeKutta).

Parameters
userThe UserCtx for a single block.
RhsThe PETSc Vec where the RHS result will be stored.
Returns
PetscErrorCode 0 on success.

Computes the Right-Hand Side (RHS) of the momentum equations.

Local to this translation unit.

Definition at line 1185 of file rhs.c.

1186{
1187 PetscErrorCode ierr;
1188 SimCtx *simCtx = user->simCtx;
1189 DM da = user->da, fda = user->fda;
1190 DMDALocalInfo info = user->info;
1191 PetscInt i,j,k;
1192 // --- Local Grid Indices and Parameters ---
1193 PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
1194 PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
1195 PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
1196 PetscInt lxs = (xs==0) ? xs+1 : xs;
1197 PetscInt lys = (ys==0) ? ys+1 : ys;
1198 PetscInt lzs = (zs==0) ? zs+1 : zs;
1199 PetscInt lxe = (xe==mx) ? xe-1 : xe;
1200 PetscInt lye = (ye==my) ? ye-1 : ye;
1201 PetscInt lze = (ze==mz) ? ze-1 : ze;
1202
1203 // --- Array Pointers ---
1204 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1205 PetscReal ***p, ***iaj, ***jaj, ***kaj, ***aj, ***nvert;
1206 Cmpnts ***rhs, ***rc, ***rct;
1207
1208 // --- Temporary Vectors ---
1209 Vec Conv, Visc, Rc, Rct;
1210
1211 PetscFunctionBeginUser;
1212 PROFILE_FUNCTION_BEGIN;
1213 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Computing RHS (FormFunction1)...\n",
1214 simCtx->rank, user->_this);
1215
1216 // --- Get all necessary array pointers ---
1217 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1218 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1219 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1220 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1221 ierr = DMDAVecGetArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1222 ierr = DMDAVecGetArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1223 ierr = DMDAVecGetArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1224 ierr = DMDAVecGetArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1225 ierr = DMDAVecGetArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1226 ierr = DMDAVecGetArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1227 ierr = DMDAVecGetArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1228 ierr = DMDAVecGetArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1229 ierr = DMDAVecGetArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1230 ierr = DMDAVecGetArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1231 ierr = DMDAVecGetArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1232 ierr = DMDAVecGetArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1233 ierr = DMDAVecGetArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1234 ierr = DMDAVecGetArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1235 ierr = DMDAVecGetArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1236
1237 // --- Create temporary work vectors ---
1238 ierr = VecDuplicate(user->lUcont, &Rc); CHKERRQ(ierr);
1239 ierr = VecDuplicate(Rc, &Rct); CHKERRQ(ierr);
1240 ierr = VecDuplicate(Rct, &Conv); CHKERRQ(ierr);
1241 ierr = VecDuplicate(Rct, &Visc); CHKERRQ(ierr);
1242
1243 // ========================================================================
1244 // CORE LOGIC (UNCHANGED FROM LEGACY CODE)
1245 // ========================================================================
1246
1247 // 1. Obtain Cartesian velocity from Contravariant velocity
1248 ierr = Contra2Cart(user); CHKERRQ(ierr);
1249 ierr = UpdateLocalGhosts(user,"Ucat"); CHKERRQ(ierr);
1250
1251 // 2. Compute Convective term
1252 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating convective terms...\n");
1253 if (simCtx->moveframe || simCtx->rotateframe) {
1254 // ierr = Convection_MV(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1255 } else {
1256 ierr = Convection(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1257 }
1258
1259 // 3. Compute Viscous term
1260 if (simCtx->invicid) {
1261 ierr = VecSet(Visc, 0.0); CHKERRQ(ierr);
1262 } else {
1263 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating viscous terms...\n");
1264 ierr = Viscous(user, user->lUcont, user->lUcat, Visc); CHKERRQ(ierr);
1265 }
1266
1267 // 4. Combine terms to get Cartesian RHS: Rc = Visc - Conv
1268 ierr = VecWAXPY(Rc, -1.0, Conv, Visc); CHKERRQ(ierr);
1269
1270 // 5. Convert Cartesian RHS (Rc) to Contravariant RHS (Rct)
1271 LOG_ALLOW(LOCAL, LOG_DEBUG, " Converting Cartesian RHS to Contravariant RHS...\n");
1272 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1273 ierr = DMDAVecGetArray(fda, Rc, &rc); CHKERRQ(ierr);
1274
1275 for (k = lzs; k < lze; k++) {
1276 for (j = lys; j < lye; j++) {
1277 for (i = lxs; i < lxe; i++) {
1278 rct[k][j][i].x = aj[k][j][i] *
1279 (0.5 * (csi[k][j][i].x + csi[k][j][i-1].x) * rc[k][j][i].x +
1280 0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
1281 0.5 * (csi[k][j][i].z + csi[k][j][i-1].z) * rc[k][j][i].z);
1282 rct[k][j][i].y = aj[k][j][i] *
1283 (0.5 * (eta[k][j][i].x + eta[k][j-1][i].x) * rc[k][j][i].x +
1284 0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
1285 0.5 * (eta[k][j][i].z + eta[k][j-1][i].z) * rc[k][j][i].z);
1286 rct[k][j][i].z = aj[k][j][i] *
1287 (0.5 * (zet[k][j][i].x + zet[k-1][j][i].x) * rc[k][j][i].x +
1288 0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
1289 0.5 * (zet[k][j][i].z + zet[k-1][j][i].z) * rc[k][j][i].z);
1290 }
1291 }
1292 }
1293 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1294 ierr = DMDAVecRestoreArray(fda, Rc, &rc); CHKERRQ(ierr);
1295
1296 PetscBarrier(NULL);
1297
1298 DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct);
1299 DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct);
1300
1301 // Compute and Add Body Force term if applicable.
1302 ierr = ComputeBodyForces(user,Rct);
1303
1304 DMDAVecGetArray(fda, Rct, &rct);
1305
1306
1307 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
1308 if (xs==0){
1309 i=xs;
1310 for (k=lzs; k<lze; k++) {
1311 for (j=lys; j<lye; j++) {
1312 rct[k][j][i].x=rct[k][j][i-2].x;
1313 }
1314 }
1315 }
1316
1317 if (xe==mx){
1318 i=mx-1;
1319 for (k=lzs; k<lze; k++) {
1320 for (j=lys; j<lye; j++) {
1321 rct[k][j][i].x=rct[k][j][i+2].x;
1322 }
1323 }
1324 }
1325 }
1326
1327 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
1328 if (ys==0){
1329 j=ys;
1330 for (k=lzs; k<lze; k++) {
1331 for (i=lxs; i<lxe; i++) {
1332 rct[k][j][i].y=rct[k][j-2][i].y;
1333 }
1334 }
1335 }
1336
1337 if (ye==my){
1338 j=my-1;
1339 for (k=lzs; k<lze; k++) {
1340 for (i=lxs; i<lxe; i++) {
1341 rct[k][j][i].y=rct[k][j+2][i].y;
1342 }
1343 }
1344 }
1345 }
1346 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
1347 if (zs==0){
1348 k=zs;
1349 for (j=lys; j<lye; j++) {
1350 for (i=lxs; i<lxe; i++) {
1351 rct[k][j][i].z=rct[k-2][j][i].z;
1352 }
1353 }
1354 }
1355 if (ze==mz){
1356 k=mz-1;
1357 for (j=lys; j<lye; j++) {
1358 for (i=lxs; i<lxe; i++) {
1359 rct[k][j][i].z=rct[k+2][j][i].z;
1360 }
1361 }
1362 }
1363 }
1364
1365 // 6. Add Pressure Gradient Term and Finalize RHS
1366 // This involves calculating pressure derivatives (dpdc, dpde, dpdz) and using
1367 // them to adjust the contravariant RHS. The full stencil logic is preserved.
1368 LOG_ALLOW(LOCAL, LOG_DEBUG, " Adding pressure gradient term to RHS...\n");
1369
1370 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1371
1372 ierr = DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1373 ierr = DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1374
1375 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1376
1377 for (k = lzs; k < lze; k++) {
1378 for (j = lys; j < lye; j++) {
1379 for (i = lxs; i < lxe; i++) {
1380 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1381 dpdc = p[k][j][i+1] - p[k][j][i];
1382
1383 if ((j==my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1384 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1385 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1386 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1387 }
1388 }
1389 else if ((j==my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1390 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1391 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1392 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1393 }
1394 }
1395 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1396 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1397 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1398 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1399 }
1400 }
1401 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1402 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1403 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1404 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1405 }
1406 }
1407 else {
1408 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1409 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1410 }
1411
1412 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1413 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1414 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1415 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1416 }
1417 }
1418 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1419 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1420 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1421 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1422 }
1423 }
1424 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1425 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1426 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1427 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1428 }
1429 }
1430 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1431 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1432 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1433 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1434 }
1435 }
1436 else {
1437 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1438 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1439 }
1440
1441 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1442
1443
1444 rhs[k][j][i].x -=
1445 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1446 icsi[k][j][i].y * icsi[k][j][i].y +
1447 icsi[k][j][i].z * icsi[k][j][i].z)+
1448 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1449 ieta[k][j][i].y * icsi[k][j][i].y +
1450 ieta[k][j][i].z * icsi[k][j][i].z)+
1451 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1452 izet[k][j][i].y * icsi[k][j][i].y +
1453 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1454
1455 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1456 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1457 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1458 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1459 }
1460 }
1461 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1462 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1463 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1464 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1465 }
1466 }
1467 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1468 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1469 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1470 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1471 }
1472 }
1473 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1474 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1475 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1476 p[k][j][i] - p[k][j+1][i]) * 0.5;
1477 }
1478 }
1479 else {
1480 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1481 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1482 }
1483
1484 dpde = p[k][j+1][i] - p[k][j][i];
1485
1486 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1487 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1488 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1489 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1490 }
1491 }
1492 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1493 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1494 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1495 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1496 }
1497 }
1498 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1499 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1500 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1501 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1502 }
1503 }
1504 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1505 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1506 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1507 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1508 }
1509 }
1510 else {
1511 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1512 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1513 }
1514
1515 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1516
1517
1518 rhs[k][j][i].y -=
1519 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1520 jcsi[k][j][i].y * jeta[k][j][i].y +
1521 jcsi[k][j][i].z * jeta[k][j][i].z) +
1522 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1523 jeta[k][j][i].y * jeta[k][j][i].y +
1524 jeta[k][j][i].z * jeta[k][j][i].z) +
1525 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1526 jzet[k][j][i].y * jeta[k][j][i].y +
1527 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1528
1529 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1530 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1531 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1532 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1533 }
1534 }
1535 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1536 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1537 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1538 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1539 }
1540 }
1541 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1542 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1543 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1544 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1545 }
1546 }
1547 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1548 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1549 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1550 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1551 }
1552 }
1553 else {
1554 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1555 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1556 }
1557
1558 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1559 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1560 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1561 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1562 }
1563 }
1564 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1565 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1566 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1567 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1568 }
1569 }
1570 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1571 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1572 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1573 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1574 }
1575 }
1576 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1577 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1578 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1579 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1580 }
1581 }
1582 else {
1583 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1584 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1585 }
1586
1587 dpdz = (p[k+1][j][i] - p[k][j][i]);
1588
1589 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1590
1591 rhs[k][j][i].z -=
1592 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1593 kcsi[k][j][i].y * kzet[k][j][i].y +
1594 kcsi[k][j][i].z * kzet[k][j][i].z) +
1595 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1596 keta[k][j][i].y * kzet[k][j][i].y +
1597 keta[k][j][i].z * kzet[k][j][i].z) +
1598 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1599 kzet[k][j][i].y * kzet[k][j][i].y +
1600 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1601
1602 }
1603 }
1604 }
1605
1606
1607 //Mohsen March 2012//
1608
1609 // rhs.x at boundaries for periodic bc at i direction//
1610 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs==0){
1611 for (k=lzs; k<lze; k++) {
1612 for (j=lys; j<lye; j++) {
1613 i=xs;
1614 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1615
1616 dpdc = p[k][j][i+1] - p[k][j][i];
1617
1618 if ((j==my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1619 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1620 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1621 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1622 }
1623 }
1624 else if ((j==my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1625 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1626 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1627 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1628 }
1629 }
1630 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1631 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1632 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1633 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1634 }
1635 }
1636 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1637 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1638 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1639 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1640 }
1641 }
1642 else {
1643 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1644 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1645 }
1646
1647 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1648 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1649 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1650 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1651 }
1652 }
1653 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1654 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1655 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1656 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1657 }
1658 }
1659 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1660 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1661 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1662 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1663 }
1664 }
1665 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1666 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1667 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1668 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1669 }
1670 }
1671 else {
1672 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1673 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1674 }
1675
1676 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1677 rhs[k][j][i].x -=
1678 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1679 icsi[k][j][i].y * icsi[k][j][i].y +
1680 icsi[k][j][i].z * icsi[k][j][i].z)+
1681 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1682 ieta[k][j][i].y * icsi[k][j][i].y +
1683 ieta[k][j][i].z * icsi[k][j][i].z)+
1684 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1685 izet[k][j][i].y * icsi[k][j][i].y +
1686 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1687 }
1688 }
1689 }
1690
1691// rhs.y at boundaries for periodic bc at j direction//
1692 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys==0){
1693 for (k=lzs; k<lze; k++) {
1694 for (i=lxs; i<lxe; i++) {
1695
1696 j=ys;
1697 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1698
1699 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1700 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1701 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1702 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1703 }
1704 }
1705 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1706 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1707 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1708 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1709 }
1710 }
1711 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1712 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1713 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1714 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1715 }
1716 }
1717 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1718 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1719 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1720 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1721 }
1722 }
1723 else {
1724 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1725 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1726 }
1727
1728 dpde = p[k][j+1][i] - p[k][j][i];
1729
1730 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1731 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1732 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1733 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1734 }
1735 }
1736 else if ((k == mz-2 || k==1 ) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1737 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1738 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1739 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1740 }
1741 }
1742 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1743 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1744 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1745 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1746 }
1747 }
1748 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1749 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1750 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1751 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1752 }
1753 }
1754 else {
1755 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1756 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1757 }
1758
1759 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1760
1761 rhs[k][j][i].y -=
1762 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1763 jcsi[k][j][i].y * jeta[k][j][i].y +
1764 jcsi[k][j][i].z * jeta[k][j][i].z)+
1765 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1766 jeta[k][j][i].y * jeta[k][j][i].y +
1767 jeta[k][j][i].z * jeta[k][j][i].z)+
1768 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1769 jzet[k][j][i].y * jeta[k][j][i].y +
1770 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1771
1772 }
1773 }
1774 }
1775
1776 // rhs.z at boundaries for periodic bc at k direction//
1777 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC&& zs==0){
1778 for (j=lys; j<lye; j++) {
1779 for (i=lxs; i<lxe; i++) {
1780
1781 k=zs;
1782 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1783
1784 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1785 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1786 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1787 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1788 }
1789 }
1790 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1791 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1792 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1793 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1794 }
1795 }
1796 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1797 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1798 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1799 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1800 }
1801 }
1802 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1803 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1804 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1805 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1806 }
1807 }
1808 else {
1809 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1810 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1811 }
1812
1813 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1814 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1815 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1816 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1817 }
1818 }
1819 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1820 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1821 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1822 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1823 }
1824 }
1825 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1826 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1827 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1828 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1829 }
1830 }
1831 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1832 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1833 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1834 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1835 }
1836 }
1837 else {
1838 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1839 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1840 }
1841
1842 dpdz = (p[k+1][j][i] - p[k][j][i]);
1843
1844 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1845
1846 rhs[k][j][i].z -=
1847 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1848 kcsi[k][j][i].y * kzet[k][j][i].y +
1849 kcsi[k][j][i].z * kzet[k][j][i].z)+
1850 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1851 keta[k][j][i].y * kzet[k][j][i].y +
1852 keta[k][j][i].z * kzet[k][j][i].z)+
1853 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1854 kzet[k][j][i].y * kzet[k][j][i].y +
1855 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1856
1857 }
1858 }
1859 }
1860
1861 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1862
1863 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure Gradient added to RHS .\n");
1864 PetscInt TwoD = simCtx->TwoD;
1865
1866 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases initiated .\n");
1867
1868 // 7. Final clean-up for immersed boundaries and 2D cases
1869 for (k=lzs; k<lze; k++) {
1870 for (j=lys; j<lye; j++) {
1871 for (i=lxs; i<lxe; i++) {
1872 if (TwoD==1)
1873 rhs[k][j][i].x =0.;
1874 else if (TwoD==2)
1875 rhs[k][j][i].y =0.;
1876 else if (TwoD==3)
1877 rhs[k][j][i].z =0.;
1878
1879 if (nvert[k][j][i]>0.1) {
1880 rhs[k][j][i].x = 0;
1881 rhs[k][j][i].y = 0;
1882 rhs[k][j][i].z = 0;
1883 }
1884 if (nvert[k][j][i+1]>0.1) {
1885 rhs[k][j][i].x=0;
1886 }
1887 if (nvert[k][j+1][i]>0.1) {
1888 rhs[k][j][i].y=0;
1889 }
1890 if (nvert[k+1][j][i]>0.1) {
1891 rhs[k][j][i].z=0;
1892 }
1893 }
1894 }
1895 }
1896 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases complete .\n");
1897
1898 // ========================================================================
1899
1900 // --- Restore all PETSc array pointers ---
1901 // DMDAVecRestoreArray(fda, user->lUcont, &ucont);
1902 DMDAVecRestoreArray(fda, Rhs, &rhs);
1903 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Rhs restored successfully! .\n");
1904
1905 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1906 DMDAVecRestoreArray(fda, user->lEta, &eta);
1907 DMDAVecRestoreArray(fda, user->lZet, &zet);
1908 DMDAVecRestoreArray(da, user->lAj, &aj);
1909 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Face metrics restored successfully! .\n");
1910
1911 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1912 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1913 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1914 DMDAVecRestoreArray(da, user->lIAj, &iaj);
1915 LOG_ALLOW(GLOBAL,LOG_DEBUG,"I Face metrics restored successfully! .\n");
1916
1917 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1918 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1919 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1920 DMDAVecRestoreArray(da, user->lJAj, &jaj);
1921 LOG_ALLOW(GLOBAL,LOG_DEBUG,"J Face metrics restored successfully! .\n");
1922
1923 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1924 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1925 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1926 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1927 LOG_ALLOW(GLOBAL,LOG_DEBUG,"K Face metrics restored successfully! .\n");
1928
1929 DMDAVecRestoreArray(da, user->lP, &p);
1930 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure restored successfully! .\n");
1931
1932 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1933 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Nvert restored successfully! .\n");
1934
1935 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Cell Centered scalars restored successfully! .\n");
1936
1937 // --- Destroy temporary work vectors ---
1938 ierr = VecDestroy(&Conv); CHKERRQ(ierr);
1939 ierr = VecDestroy(&Visc); CHKERRQ(ierr);
1940 ierr = VecDestroy(&Rc); CHKERRQ(ierr);
1941 ierr = VecDestroy(&Rct); CHKERRQ(ierr);
1942 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Temporary work vectors destroyed successfully! .\n");
1943
1944 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: RHS computation complete.\n",
1945 simCtx->rank, user->_this);
1946
1947 PROFILE_FUNCTION_END;
1948 PetscFunctionReturn(0);
1949}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
ComputeBodyForces
PetscErrorCode ComputeBodyForces(UserCtx *user, Vec Rct)
Internal helper implementation: ComputeBodyForces().
Definition rhs.c:1162
Viscous
PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
Implementation of Viscous().
Definition rhs.c:519
Convection
PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
Implementation of Convection().
Definition rhs.c:13
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2247
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1361
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:674
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:674
SimCtx::rank
PetscMPIInt rank
Definition variables.h:646
UserCtx::_this
PetscInt _this
Definition variables.h:824
SimCtx::invicid
PetscInt invicid
Definition variables.h:674
UserCtx::lUcont
Vec lUcont
Definition variables.h:837
UserCtx::info
DMDALocalInfo info
Definition variables.h:818
UserCtx::lUcat
Vec lUcat
Definition variables.h:837
SimCtx::rotateframe
PetscInt rotateframe
Definition variables.h:674
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeEulerianDiffusivity()

PetscErrorCode ComputeEulerianDiffusivity ( UserCtx *  user)

Computes the effective diffusivity scalar field (Gamma_eff) on the Eulerian grid.

This function calculates the total diffusivity used to drive the stochastic motion of particles (Scalar FDF). It combines molecular diffusion and turbulent diffusion.

Formula: Gamma_eff = nu/Sc + nu_t/Sc_t

Where:

  • nu = 1/Re (kinematic viscosity)
  • nu_t (eddy viscosity from LES/RANS model)
  • Sc (molecular Schmidt number)
  • Sc_t (turbulent Schmidt number)
Note
If turbulence models are disabled, nu_t is assumed to be 0.
This function updates the local ghost values of lDiffusivity at the end to ensure gradients can be computed correctly at subdomain boundaries.
Parameters
[in,out]userPointer to the user context containing grid data and simulation parameters.
Returns
PetscErrorCode 0 on success.

Computes the effective diffusivity scalar field (Gamma_eff) on the Eulerian grid.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/rhs.h.

See also
ComputeEulerianDiffusivity()

Definition at line 1959 of file rhs.c.

1960{
1961 PetscErrorCode ierr;
1962 DM da = user->da;
1963 PetscInt i, j, k, xs, ys, zs, xm, ym, zm, xe, ye, ze;
1964 PetscInt lxs, lys, lzs, lxe, lye, lze;
1965
1966 // Pointers for 3D grid access
1967 PetscReal ***diff_arr; // Output: Diffusivity field
1968 PetscReal ***nut_arr; // Input: Eddy Viscosity field (optional)
1969
1970 // Physics parameters
1971 PetscReal nu_molecular, gamma_molecular;
1972 PetscReal nu_turbulent, gamma_turbulent;
1973 PetscReal Sc, Sct;
1974 PetscBool use_turbulence_model;
1975
1976 PetscFunctionBeginUser;
1977 PROFILE_FUNCTION_BEGIN;
1978
1979 if (VerificationDiffusivityOverrideActive(user->simCtx)) {
1980 ierr = ApplyVerificationDiffusivityOverride(user); CHKERRQ(ierr);
1981 PROFILE_FUNCTION_END;
1982 PetscFunctionReturn(0);
1983 }
1984
1985 // ------------------------------------------------------------------------
1986 // 1. Parameter Setup & Safety Checks
1987 // ------------------------------------------------------------------------
1988
1989 // Determine Molecular Viscosity (nu = 1/Re)
1990 // Guard against division by zero if Re is not set or infinite (inviscid)
1991 if (user->simCtx->ren > 1.0e-12) {
1992 nu_molecular = 1.0 / user->simCtx->ren;
1993 } else {
1994 nu_molecular = 0.0;
1995 }
1996
1997 // Set Schmidt Numbers (Default to 1.0 if not provided to prevent NaN)
1998 Sc = (user->simCtx->schmidt_number > 1.0e-6) ? user->simCtx->schmidt_number : 1.0;
1999 Sct = (user->simCtx->Turbulent_schmidt_number > 1.0e-6) ? user->simCtx->Turbulent_schmidt_number : 0.7;
2000
2001 // Pre-calculate molecular component
2002 gamma_molecular = nu_molecular / Sc;
2003
2004 // Check if a turbulence model is active (LES or RANS)
2005 use_turbulence_model = (user->simCtx->les || user->simCtx->rans) ? PETSC_TRUE : PETSC_FALSE;
2006
2007 // ------------------------------------------------------------------------
2008 // 2. Data Access
2009 // ------------------------------------------------------------------------
2010
2011 // Get local grid boundaries
2012 DMDALocalInfo info;
2013 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
2014
2015 xs = info.xs; ys = info.ys; zs = info.zs;
2016 xm = info.xm; ym = info.ym; zm = info.zm;
2017 xe = xs + xm; ye = ys + ym; ze = zs + zm;
2018
2019 lxs = (xs == 0)? xs + 1 : xs;
2020 lys = (ys == 0)? ys + 1 : ys;
2021 lzs = (zs == 0)? zs + 1 : zs;
2022 lxe = (xe == info.mx)? xe - 1 : xe;
2023 lye = (ye == info.my)? ye - 1 : ye;
2024 lze = (ze == info.mz)? ze - 1 : ze;
2025
2026 // Get write access to the output Diffusivity array
2027 ierr = DMDAVecGetArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
2028
2029 // Get read access to Eddy Viscosity only if turbulence is active
2030 if (use_turbulence_model) {
2031 ierr = DMDAVecGetArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
2032 }
2033
2034 // ------------------------------------------------------------------------
2035 // 3. Calculation Loop
2036 // ------------------------------------------------------------------------
2037
2038 for (k = lzs; k < lze; k++) {
2039 for (j = lys; j < lye; j++) {
2040 for (i = lxs; i < lxe; i++) {
2041
2042 gamma_turbulent = 0.0;
2043
2044 if (use_turbulence_model) {
2045 // Fetch local eddy viscosity
2046 nu_turbulent = nut_arr[k][j][i];
2047
2048 // NUMERICAL SAFETY:
2049 // Some turbulence models (dynamic SGS) can locally produce
2050 // slightly negative viscosity. We clamp this to 0 to prevent
2051 // negative diffusivity, which crashes the Langevin sqrt().
2052 if (nu_turbulent < 0.0) {
2053 nu_turbulent = 0.0;
2054 }
2055
2056 gamma_turbulent = nu_turbulent / Sct;
2057 }
2058
2059 // Sum components
2060 diff_arr[k][j][i] = gamma_molecular + gamma_turbulent;
2061 }
2062 }
2063 }
2064
2065 // ------------------------------------------------------------------------
2066 // 4. Cleanup & Synchronization
2067 // ------------------------------------------------------------------------
2068
2069 // Restore arrays
2070 if (use_turbulence_model) {
2071 ierr = DMDAVecRestoreArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
2072 }
2073 ierr = DMDAVecRestoreArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
2074
2075 // Update Ghost Points
2076 // This is required because downstream operations (Drift Gradient Calculation
2077 // and Particle Interpolation) will need access to the halo regions of this field.
2078 ierr = UpdateLocalGhosts(user,"Diffusivity"); CHKERRQ(ierr);
2079 PROFILE_FUNCTION_END;
2080 PetscFunctionReturn(0);
2081}
SimCtx::schmidt_number
PetscReal schmidt_number
Definition variables.h:709
SimCtx::Turbulent_schmidt_number
PetscReal Turbulent_schmidt_number
Definition variables.h:709
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:865
UserCtx::Diffusivity
Vec Diffusivity
Definition variables.h:840
ApplyVerificationDiffusivityOverride
PetscErrorCode ApplyVerificationDiffusivityOverride(UserCtx *user)
Populates the Eulerian diffusivity field from a verification-only source override.
Definition verification_sources.c:31
VerificationDiffusivityOverrideActive
PetscBool VerificationDiffusivityOverrideActive(const SimCtx *simCtx)
Reports whether a verification-only diffusivity override is active.
Definition verification_sources.c:11
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeEulerianDiffusivityGradient()

PetscErrorCode ComputeEulerianDiffusivityGradient ( UserCtx *  user)

Computes the Eulerian gradient of the effective diffusivity field.

Reads the scalar diffusivity field and writes a vector gradient field used by particle stochastic transport updates.

Parameters
[in,out]userPointer to user context containing diffusivity vectors.
Returns
PetscErrorCode 0 on success.

Computes the Eulerian gradient of the effective diffusivity field.

Local to this translation unit.

Definition at line 2089 of file rhs.c.

2090{
2091 PetscErrorCode ierr;
2092 DM da = user->da, fda = user->fda;
2093 DMDALocalInfo info = user->info;
2094
2095 // 1. Determine Global Dimensions
2096 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2097
2098 // 2. Determine Local Loop Bounds (Skip Ghosts/Unused Indices)
2099 // Grid uses indices 1 to M-2 for physical cells.
2100 // Index 0 and M-1 are ghost/boundary holders.
2101
2102 // Start: If we own the global start (0), skip it and start at 1.
2103 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2104 // End: If we own the global end (mx), stop before it (mx-1), so loop covers mx-2.
2105 PetscInt lxe = (info.xs + info.xm == mx) ? mx - 1 : info.xs + info.xm;
2106
2107 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2108 PetscInt lye = (info.ys + info.ym == my) ? my - 1 : info.ys + info.ym;
2109
2110 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2111 PetscInt lze = (info.zs + info.zm == mz) ? mz - 1 : info.zs + info.zm;
2112
2113 PetscInt i, j, k;
2114
2115 // Pointers
2116 PetscReal ***diff; // Input (Scalar)
2117 Cmpnts ***grad_diff; // Output (Vector)
2118 Cmpnts ***csi, ***eta, ***zet;
2119 PetscReal ***aj;
2120
2121 // Boundary Flags (Check if Periodic)
2122 PetscBool p_x = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC);
2123 PetscBool p_y = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC);
2124 PetscBool p_z = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC);
2125
2126 PetscFunctionBeginUser;
2127 PROFILE_FUNCTION_BEGIN;
2128
2129 // 3. Update Ghosts for Diffusivity
2130 // Required so that Central Differences at i=2 can read i=1,
2131 // and Central Differences at Periodic Boundaries work correctly.
2132 ierr = UpdateLocalGhosts(user, "Diffusivity"); CHKERRQ(ierr);
2133
2134 // 4. Get Arrays (Read Only)
2135 ierr = DMDAVecGetArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2136 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2137 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2138 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2139 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2140
2141 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity and Metrics arrays accessed successfully! .\n");
2142
2143 // 5. Get Output Array (Read/Write)
2144 ierr = DMDAVecGetArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2145
2146 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity Gradient array accessed successfully! .\n");
2147
2148 // 6. Loop over Physical Domain
2149 for (k = lzs; k < lze; k++) {
2150 for (j = lys; j < lye; j++) {
2151 for (i = lxs; i < lxe; i++) {
2152
2153 PetscReal dGdCsi, dGdEta, dGdZet;
2154
2155 // ---------------------------------------------------------
2156 // I-Direction (Csi)
2157 // ---------------------------------------------------------
2158 if (!p_x && i == 1) {
2159 // Physical Start: 2nd Order Forward Difference
2160 // Stencil: [-3, 4, -1] / 2
2161 dGdCsi = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j][i+1] - diff[k][j][i+2]) * 0.5;
2162 }
2163 else if (!p_x && i == mx - 2) {
2164 // Physical End: 2nd Order Backward Difference
2165 // Stencil: [1, -4, 3] / 2
2166 dGdCsi = (3.0 * diff[k][j][i] - 4.0 * diff[k][j][i-1] + diff[k][j][i-2]) * 0.5;
2167 }
2168 else {
2169 // Interior / Periodic: Central Difference
2170 // Stencil: [-1, 0, 1] / 2
2171 dGdCsi = (diff[k][j][i+1] - diff[k][j][i-1]) * 0.5;
2172 }
2173
2174 // ---------------------------------------------------------
2175 // J-Direction (Eta)
2176 // ---------------------------------------------------------
2177 if (!p_y && j == 1) {
2178 // Physical Start: Forward
2179 dGdEta = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j+1][i] - diff[k][j+2][i]) * 0.5;
2180 }
2181 else if (!p_y && j == my - 2) {
2182 // Physical End: Backward
2183 dGdEta = (3.0 * diff[k][j][i] - 4.0 * diff[k][j-1][i] + diff[k][j-2][i]) * 0.5;
2184 }
2185 else {
2186 // Interior: Central
2187 dGdEta = (diff[k][j+1][i] - diff[k][j-1][i]) * 0.5;
2188 }
2189
2190 // ---------------------------------------------------------
2191 // K-Direction (Zet)
2192 // ---------------------------------------------------------
2193 if (!p_z && k == 1) {
2194 // Physical Start: Forward
2195 dGdZet = (-3.0 * diff[k][j][i] + 4.0 * diff[k+1][j][i] - diff[k+2][j][i]) * 0.5;
2196 }
2197 else if (!p_z && k == mz - 2) {
2198 // Physical End: Backward
2199 dGdZet = (3.0 * diff[k][j][i] - 4.0 * diff[k-1][j][i] + diff[k-2][j][i]) * 0.5;
2200 }
2201 else {
2202 // Interior: Central
2203 dGdZet = (diff[k+1][j][i] - diff[k-1][j][i]) * 0.5;
2204 }
2205
2206 // ---------------------------------------------------------
2207 // Transform to Physical Space (Cartesian Gradient)
2208 // ---------------------------------------------------------
2209 TransformScalarDerivativesToPhysical(aj[k][j][i],
2210 csi[k][j][i], eta[k][j][i], zet[k][j][i],
2211 dGdCsi, dGdEta, dGdZet,
2212 &grad_diff[k][j][i]);
2213 }
2214 }
2215 }
2216
2217 // 7. Restore Arrays
2218 ierr = DMDAVecRestoreArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2219 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2220 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2221 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2222 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2223 ierr = DMDAVecRestoreArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2224
2225 // 8. Update Ghosts for the Result
2226 // Important: Particles near subdomain boundaries will need to interpolate
2227 // this gradient vector, so the ghosts of the vector field must be filled.
2228 ierr = UpdateLocalGhosts(user, "DiffusivityGradient"); CHKERRQ(ierr);
2229
2230 PROFILE_FUNCTION_END;
2231 PetscFunctionReturn(0);
2232}
TransformScalarDerivativesToPhysical
void TransformScalarDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, PetscReal dPhi_dcsi, PetscReal dPhi_deta, PetscReal dPhi_dzet, Cmpnts *gradPhi)
Transforms scalar derivatives from computational space to physical space.
Definition setup.c:2808
UserCtx::DiffusivityGradient
Vec DiffusivityGradient
Definition variables.h:841
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:840
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.8