PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Macros | Functions
poisson.c File Reference
#include "poisson.h"
#include "logging.h"
Include dependency graph for poisson.c:

Go to the source code of this file.

Macros

#define __FUNCT__   "GhostNodeVelocity"
 
#define GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user)
 
#define __FUNCT__   "GridRestriction"
 
#define CP   0
 
#define EP   1
 
#define WP   2
 
#define NP   3
 
#define SP   4
 
#define TP   5
 
#define BP   6
 
#define NE   7
 
#define SE   8
 
#define NW   9
 
#define SW   10
 
#define TN   11
 
#define BN   12
 
#define TS   13
 
#define BS   14
 
#define TE   15
 
#define BE   16
 
#define TW   17
 
#define BW   18
 
#define __FUNCT__   "Projection"
 
#define __FUNCT__   "UpdatePressure"
 
#define __FUNCT__   "PoissonNullSpaceFunction"
 
#define __FUNCT__   "PoissonLHSNew"
 
#define __FUNCT__   "PoissonRHS"
 
#define __FUNCT__   "PoissonSolver_MG"
 

Functions

static PetscErrorCode GhostNodeVelocity (UserCtx *user)
 
static PetscInt Gidx (PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
 
static PetscErrorCode GridRestriction (PetscInt i, PetscInt j, PetscInt k, PetscInt *ih, PetscInt *jh, PetscInt *kh, UserCtx *user)
 
PetscErrorCode Projection (UserCtx *user)
 Performs the projection step to enforce an incompressible velocity field.
 
PetscErrorCode UpdatePressure (UserCtx *user)
 Updates the pressure field and handles periodic boundary conditions.
 
static PetscErrorCode mymatmultadd (Mat mat, Vec v1, Vec v2)
 
PetscErrorCode PoissonNullSpaceFunction (MatNullSpace nullsp, Vec X, void *ctx)
 The callback function for PETSc's MatNullSpace object.
 
PetscErrorCode MyInterpolation (Mat A, Vec X, Vec F)
 The callback function for the multigrid interpolation operator (MatShell).
 
static PetscErrorCode RestrictResidual_SolidAware (Mat A, Vec X, Vec F)
 
PetscErrorCode MyRestriction (Mat A, Vec X, Vec F)
 The callback function for the multigrid restriction operator (MatShell).
 
PetscErrorCode PoissonLHSNew (UserCtx *user)
 Assembles the Left-Hand Side (LHS) matrix for the Pressure Poisson Equation.
 
PetscErrorCode PoissonRHS (UserCtx *user, Vec B)
 Computes the Right-Hand-Side (RHS) of the Poisson equation, which is the divergence of the intermediate velocity field.
 
PetscErrorCode VolumeFlux_rev (UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
 A specialized version of VolumeFlux, likely for reversed normals.
 
PetscErrorCode VolumeFlux (UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
 Calculates the net flux across the immersed boundary surface.
 
PetscErrorCode FullyBlocked (UserCtx *user)
 
PetscErrorCode MyNvertRestriction (UserCtx *user_h, UserCtx *user_c)
 
PetscErrorCode PoissonSolver_MG (UserMG *usermg)
 Solves the pressure-Poisson equation using a geometric multigrid method.
 

Macro Definition Documentation

◆ __FUNCT__ [1/8]

#define __FUNCT__   "GhostNodeVelocity"

Definition at line 6 of file poisson.c.

◆ GridInterpolation

#define GridInterpolation (   i,
  j,
  k,
  ic,
  jc,
  kc,
  ia,
  ja,
  ka,
  user 
)

Definition at line 798 of file poisson.c.

799 { \
800 ic = i; \
801 ia = 0; \
802 } \
803 else { \
804 ic = (i+1) / 2; \
805 ia = (i - 2 * (ic)) == 0 ? 1 : -1; \
806 if (i==1 || i==mx-2) ia = 0; \
807 }\
808 if ((user->jsc)) { \
809 jc = j; \
810 ja = 0; \
811 } \
812 else { \
813 jc = (j+1) / 2; \
814 ja = (j - 2 * (jc)) == 0 ? 1 : -1; \
815 if (j==1 || j==my-2) ja = 0; \
816 } \
817 if ((user->ksc)) { \
818 kc = k; \
819 ka = 0; \
820 } \
821 else { \
822 kc = (k+1) / 2; \
823 ka = (k - 2 * (kc)) == 0 ? 1 : -1; \
824 if (k==1 || k==mz-2) ka = 0; \
825 } \
826 if (ka==-1 && nvert_c[kc-1][jc][ic] > 0.1) ka = 0; \
827 else if (ka==1 && nvert_c[kc+1][jc][ic] > 0.1) ka = 0; \
828 if (ja==-1 && nvert_c[kc][jc-1][ic] > 0.1) ja = 0; \
829 else if (ja==1 && nvert_c[kc][jc+1][ic] > 0.1) ja = 0; \
830 if (ia==-1 && nvert_c[kc][jc][ic-1] > 0.1) ia = 0; \
831 else if (ia==1 && nvert_c[kc][jc][ic+1] > 0.1) ia = 0;

◆ __FUNCT__ [2/8]

#define __FUNCT__   "GridRestriction"

Definition at line 6 of file poisson.c.

◆ CP

#define CP   0

Definition at line 886 of file poisson.c.

◆ EP

#define EP   1

Definition at line 888 of file poisson.c.

◆ WP

#define WP   2

Definition at line 889 of file poisson.c.

◆ NP

#define NP   3

Definition at line 890 of file poisson.c.

◆ SP

#define SP   4

Definition at line 891 of file poisson.c.

◆ TP

#define TP   5

Definition at line 892 of file poisson.c.

◆ BP

#define BP   6

Definition at line 893 of file poisson.c.

◆ NE

#define NE   7

Definition at line 896 of file poisson.c.

◆ SE

#define SE   8

Definition at line 897 of file poisson.c.

◆ NW

#define NW   9

Definition at line 898 of file poisson.c.

◆ SW

#define SW   10

Definition at line 899 of file poisson.c.

◆ TN

#define TN   11

Definition at line 900 of file poisson.c.

◆ BN

#define BN   12

Definition at line 901 of file poisson.c.

◆ TS

#define TS   13

Definition at line 902 of file poisson.c.

◆ BS

#define BS   14

Definition at line 903 of file poisson.c.

◆ TE

#define TE   15

Definition at line 904 of file poisson.c.

◆ BE

#define BE   16

Definition at line 905 of file poisson.c.

◆ TW

#define TW   17

Definition at line 906 of file poisson.c.

◆ BW

#define BW   18

Definition at line 907 of file poisson.c.

◆ __FUNCT__ [3/8]

#define __FUNCT__   "Projection"

Definition at line 6 of file poisson.c.

◆ __FUNCT__ [4/8]

#define __FUNCT__   "UpdatePressure"

Definition at line 6 of file poisson.c.

◆ __FUNCT__ [5/8]

#define __FUNCT__   "PoissonNullSpaceFunction"

Definition at line 6 of file poisson.c.

◆ __FUNCT__ [6/8]

#define __FUNCT__   "PoissonLHSNew"

Definition at line 6 of file poisson.c.

◆ __FUNCT__ [7/8]

#define __FUNCT__   "PoissonRHS"

Definition at line 6 of file poisson.c.

◆ __FUNCT__ [8/8]

#define __FUNCT__   "PoissonSolver_MG"

Definition at line 6 of file poisson.c.

Function Documentation

◆ GhostNodeVelocity()

static PetscErrorCode GhostNodeVelocity ( UserCtx *  user)
static

Definition at line 8 of file poisson.c.

9{
10
11 // --- CONTEXT ACQUISITION BLOCK ---
12 // Get the master simulation context from the UserCtx.
13 SimCtx *simCtx = user->simCtx;
14
15 // Create local variables to mirror the legacy globals for minimal code changes.
16 PetscInt wallfunction = simCtx->wallfunction;
17 const PetscInt les = simCtx->les;
18 const PetscInt rans = simCtx->rans;
19 const PetscInt ti = simCtx->step;
20 const PetscReal U_bc = simCtx->U_bc;
21 // --- END CONTEXT ACQUISITION BLOCK ---
22
23 DM da = user->da, fda = user->fda;
24 DMDALocalInfo info = user->info;
25 PetscInt xs = info.xs, xe = info.xs + info.xm;
26 PetscInt ys = info.ys, ye = info.ys + info.ym;
27 PetscInt zs = info.zs, ze = info.zs + info.zm;
28 PetscInt mx = info.mx, my = info.my, mz = info.mz;
29
30 PetscInt i, j, k;
31 PetscInt lxs, lxe, lys, lye, lzs, lze;
32
33 Cmpnts ***ubcs, ***ucat,***lucat, ***csi, ***eta, ***zet,***cent;
34
35
36 PetscReal Un, nx,ny,nz,A;
37
38 lxs = xs; lxe = xe;
39 lys = ys; lye = ye;
40 lzs = zs; lze = ze;
41
42 PetscInt gxs, gxe, gys, gye, gzs, gze;
43
44 gxs = info.gxs; gxe = gxs + info.gxm;
45 gys = info.gys; gye = gys + info.gym;
46 gzs = info.gzs; gze = gzs + info.gzm;
47
48 if (xs==0) lxs = xs+1;
49 if (ys==0) lys = ys+1;
50 if (zs==0) lzs = zs+1;
51
52 if (xe==mx) lxe = xe-1;
53 if (ye==my) lye = ye-1;
54 if (ze==mz) lze = ze-1;
55
56 PetscFunctionBeginUser;
57
58 PROFILE_FUNCTION_BEGIN;
59
60 DMDAVecGetArray(fda, user->lCsi, &csi);
61 DMDAVecGetArray(fda, user->lEta, &eta);
62 DMDAVecGetArray(fda, user->lZet, &zet);
63
64
65/* ================================================================================== */
66/* WALL FUNCTION */
67/* ================================================================================== */
68 if (wallfunction) {
69 PetscReal ***ustar, ***aj, ***nvert;
70 DMDAVecGetArray(fda, user->Ucat, &ucat);
71 DMDAVecGetArray(da, user->lUstar, &ustar);
72 DMDAVecGetArray(da, user->lAj, &aj);
73 DMDAVecGetArray(da, user->lNvert, &nvert);
74
75 // wall function for boundary
76 //Dec 2015
77 for (k=lzs; k<lze; k++)
78 for (j=lys; j<lye; j++)
79 for (i=lxs; i<lxe; i++) {
80 if( nvert[k][j][i]<1.1 && ( (user->bctype[0]==-1 && i==1) || (user->bctype[1]==-1 && i==mx-2) ) ) {
81 double area = sqrt( csi[k][j][i].x*csi[k][j][i].x + csi[k][j][i].y*csi[k][j][i].y + csi[k][j][i].z*csi[k][j][i].z );
82 double sb, sc;
83 double ni[3], nj[3], nk[3];
84 Cmpnts Uc, Ua;
85
86 Ua.x = Ua.y = Ua.z = 0;
87 sb = 0.5/aj[k][j][i]/area;
88
89 if(i==1) {
90 sc = 2*sb + 0.5/aj[k][j][i+1]/area;
91 Uc = ucat[k][j][i+1];
92 }
93 else {
94 sc = 2*sb + 0.5/aj[k][j][i-1]/area;
95 Uc = ucat[k][j][i-1];
96 }
97
98 CalculateFaceNormalAndArea(csi[k][j][i], eta[k][j][i], zet[k][j][i], ni, nj, nk,NULL,NULL,NULL); // Passing Null pointers for area calculation.
99 if(i==mx-2) ni[0]*=-1, ni[1]*=-1, ni[2]*=-1;
100
101 //if(i==1) printf("%f %f, %f %f %f, %e %e %e, %e\n", sc, sb, Ua.x, Ua.y, Ua.z, Uc.x, Uc.y, Uc.z, ucat[k][j][i+2].z);
102 wall_function (user, sc, sb, Ua, Uc, &ucat[k][j][i], &ustar[k][j][i], ni[0], ni[1], ni[2]);
103 nvert[k][j][i]=1; /* set nvert to 1 to exclude from rhs */
104 // ucat[k][j][i] = lucat[k][j][i];
105 }
106
107 if( nvert[k][j][i]<1.1 && ( (user->bctype[2]==-1 && j==1) || (user->bctype[3]==-1 && j==my-2) ) ) {
108 double area = sqrt( eta[k][j][i].x*eta[k][j][i].x + eta[k][j][i].y*eta[k][j][i].y + eta[k][j][i].z*eta[k][j][i].z );
109 double sb, sc;
110 double ni[3], nj[3], nk[3];
111 Cmpnts Uc, Ua;
112
113 Ua.x = Ua.y = Ua.z = 0;
114 sb = 0.5/aj[k][j][i]/area;
115
116 if(j==1) {
117 sc = 2*sb + 0.5/aj[k][j+1][i]/area;
118 Uc = ucat[k][j+1][i];
119 }
120 else {
121 sc = 2*sb + 0.5/aj[k][j-1][i]/area;
122 Uc = ucat[k][j-1][i];
123 }
124
125 CalculateFaceNormalAndArea(csi[k][j][i], eta[k][j][i], zet[k][j][i], ni, nj, nk,NULL,NULL,NULL); // Passing Null pointers for Area.
126 //if(j==my-2) nj[0]*=-1, nj[1]*=-1, nj[2]*=-1;
127
128 wall_function (user, sc, sb, Ua, Uc, &ucat[k][j][i], &ustar[k][j][i], nj[0], nj[1], nj[2]);
129 nvert[k][j][i]=1; /* set nvert to 1 to exclude from rhs */
130 // ucat[k][j][i] = lucat[k][j][i];
131 }
132
133 }
134 DMDAVecRestoreArray(da, user->lAj, &aj);
135 DMDAVecRestoreArray(da, user->lNvert, &nvert);
136 DMDAVecRestoreArray(da, user->lUstar, &ustar);
137 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
138
139 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
140 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
141 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
142/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
143
144 } // if wallfunction
145
146/* ================================================================================== */
147/* Driven Cavity */
148/* ================================================================================== */
149
150 DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs);
151 DMDAVecGetArray(fda, user->Ucat, &ucat);
152
153 /* Designed for driven cavity problem (top(k=kmax) wall moving)
154 u_x = 1 at k==kmax */
155 if (user->bctype[5]==2) {
156 if (ze==mz) {
157 k = ze-1;
158 for (j=lys; j<lye; j++) {
159 for (i=lxs; i<lxe; i++) {
160 ubcs[k][j][i].x = 0.;// - ucat[k-1][j][i].x;
161 ubcs[k][j][i].y = 1.;//- ucat[k-1][j][i].y;
162 ubcs[k][j][i].z = 0.;//- ucat[k-1][j][i].z;
163 }
164 }
165 }
166 }
167 /* ==================================================================================== */
168 /* Cylinder O-grid */
169 /* ==================================================================================== */
170 if (user->bctype[3]==12) {
171 /* Designed to test O-grid for flow over a cylinder at jmax velocity is 1 (horizontal)
172 u_x = 1 at k==kmax */
173 if (ye==my) {
174 j = ye-1;
175 for (k=lzs; k<lze; k++) {
176 for (i=lxs; i<lxe; i++) {
177 ubcs[k][j][i].x = 0.;
178 ubcs[k][j][i].y = 0.;
179 ubcs[k][j][i].z = 1.;
180 }
181 }
182 }
183 }
184/* ==================================================================================== */
185/* Annulus */
186/* ==================================================================================== */
187/* designed to test periodic boundary condition for c grid j=0 rotates */
188
189 if (user->bctype[2]==11) {
190 DMDAVecGetArray(fda, user->Cent, &cent);
191 if (ys==0){
192 j=0;
193 for (k=lzs; k<lze; k++) {
194 for (i=lxs; i<lxe; i++) {
195 ubcs[k][j][i].x =cent[k][j+1][i].y/sqrt(cent[k][j+1][i].x*cent[k][j+1][i].x+cent[k][j+1][i].y*cent[k][j+1][i].y);;
196 // ubcs[k][j][i].x=0.0;
197 ubcs[k][j][i].y =-cent[k][j+1][i].x/sqrt(cent[k][j+1][i].x*cent[k][j+1][i].x+cent[k][j+1][i].y*cent[k][j+1][i].y);
198 // ubcs[k][j][i].y = -cent[k][j+1][i].z/sqrt(cent[k][j+1][i].z*cent[k][j+1][i].z+cent[k][j+1][i].y*cent[k][j+1][i].y);
199 ubcs[k][j][i].z = 0.0;
200 // ubcs[k][j][i].z =cent[k][j+1][i].y/sqrt(cent[k][j+1][i].z*cent[k][j+1][i].z+cent[k][j+1][i].y*cent[k][j+1][i].y);
201 }
202 }
203 }
204 DMDAVecRestoreArray(fda, user->Cent, &cent);
205 }
206
207/* ================================================================================== */
208/* SYMMETRY BC */
209/* ================================================================================== */
210
211 if (user->bctype[0]==3) {
212 if (xs==0) {
213 i= xs;
214
215 for (k=lzs; k<lze; k++) {
216 for (j=lys; j<lye; j++) {
217 A=sqrt(csi[k][j][i].z*csi[k][j][i].z +
218 csi[k][j][i].y*csi[k][j][i].y +
219 csi[k][j][i].x*csi[k][j][i].x);
220 nx=csi[k][j][i].x/A;
221 ny=csi[k][j][i].y/A;
222 nz=csi[k][j][i].z/A;
223 Un=ucat[k][j][i+1].x*nx+ucat[k][j][i+1].y*ny+ucat[k][j][i+1].z*nz;
224 ubcs[k][j][i].x = ucat[k][j][i+1].x-Un*nx;//-V_frame.x;
225 ubcs[k][j][i].y = ucat[k][j][i+1].y-Un*ny;
226 ubcs[k][j][i].z = ucat[k][j][i+1].z-Un*nz;
227 }
228 }
229 }
230 }
231
232 if (user->bctype[1]==3) {
233 if (xe==mx) {
234 i= xe-1;
235
236 for (k=lzs; k<lze; k++) {
237 for (j=lys; j<lye; j++) {
238 A=sqrt(csi[k][j][i-1].z*csi[k][j][i-1].z +
239 csi[k][j][i-1].y*csi[k][j][i-1].y +
240 csi[k][j][i-1].x*csi[k][j][i-1].x);
241 nx=csi[k][j][i-1].x/A;
242 ny=csi[k][j][i-1].y/A;
243 nz=csi[k][j][i-1].z/A;
244 Un=ucat[k][j][i-1].x*nx+ucat[k][j][i-1].y*ny+ucat[k][j][i-1].z*nz;
245 ubcs[k][j][i].x = ucat[k][j][i-1].x-Un*nx;
246 ubcs[k][j][i].y = ucat[k][j][i-1].y-Un*ny;
247 ubcs[k][j][i].z = ucat[k][j][i-1].z-Un*nz;
248 }
249 }
250 }
251 }
252
253 if (user->bctype[2]==3) {
254 if (ys==0) {
255 j= ys;
256
257 for (k=lzs; k<lze; k++) {
258 for (i=lxs; i<lxe; i++) {
259 A=sqrt(eta[k][j][i].z*eta[k][j][i].z +
260 eta[k][j][i].y*eta[k][j][i].y +
261 eta[k][j][i].x*eta[k][j][i].x);
262 nx=eta[k][j][i].x/A;
263 ny=eta[k][j][i].y/A;
264 nz=eta[k][j][i].z/A;
265 Un=ucat[k][j+1][i].x*nx+ucat[k][j+1][i].y*ny+ucat[k][j+1][i].z*nz;
266 ubcs[k][j][i].x = ucat[k][j+1][i].x-Un*nx;
267 ubcs[k][j][i].y = ucat[k][j+1][i].y-Un*ny;
268 ubcs[k][j][i].z = ucat[k][j+1][i].z-Un*nz;
269 }
270 }
271 }
272 }
273
274 if (user->bctype[3]==3) {
275 if (ye==my) {
276 j=ye-1;
277
278 for (k=lzs; k<lze; k++) {
279 for (i=lxs; i<lxe; i++) {
280 A=sqrt(eta[k][j-1][i].z*eta[k][j-1][i].z +
281 eta[k][j-1][i].y*eta[k][j-1][i].y +
282 eta[k][j-1][i].x*eta[k][j-1][i].x);
283 nx=eta[k][j-1][i].x/A;
284 ny=eta[k][j-1][i].y/A;
285 nz=eta[k][j-1][i].z/A;
286 Un=ucat[k][j-1][i].x*nx+ucat[k][j-1][i].y*ny+ucat[k][j-1][i].z*nz;
287 ubcs[k][j][i].x = ucat[k][j-1][i].x-Un*nx;
288 ubcs[k][j][i].y = ucat[k][j-1][i].y-Un*ny;
289 ubcs[k][j][i].z = ucat[k][j-1][i].z-Un*nz;
290 }
291 }
292 }
293 }
294
295
296 if (user->bctype[4]==3) {
297 if (zs==0) {
298 k = 0;
299 for (j=lys; j<lye; j++) {
300 for (i=lxs; i<lxe; i++) {
301 A=sqrt(zet[k][j][i].z*zet[k][j][i].z +
302 zet[k][j][i].y*zet[k][j][i].y +
303 zet[k][j][i].x*zet[k][j][i].x);
304 nx=zet[k][j][i].x/A;
305 ny=zet[k][j][i].y/A;
306 nz=zet[k][j][i].z/A;
307 Un=ucat[k][j][i].x*nx+ucat[k][j][i].y*ny+ucat[k][j][i].z*nz;
308 ubcs[k][j][i].x = ucat[k][j][i].x-Un*nx;
309 ubcs[k][j][i].y = ucat[k][j][i].y-Un*ny;
310 ubcs[k][j][i].z = ucat[k][j][i].z-Un*nz;
311 }
312 }
313 }
314 }
315
316 if (user->bctype[5]==3) {
317 if (ze==mz) {
318 k = ze-1;
319 for (j=lys; j<lye; j++) {
320 for (i=lxs; i<lxe; i++) {
321 A=sqrt(zet[k-1][j][i].z*zet[k-1][j][i].z +
322 zet[k-1][j][i].y*zet[k-1][j][i].y +
323 zet[k-1][j][i].x*zet[k-1][j][i].x);
324 nx=zet[k-1][j][i].x/A;
325 ny=zet[k-1][j][i].y/A;
326 nz=zet[k-1][j][i].z/A;
327 Un=ucat[k-1][j][i].x*nx+ucat[k-1][j][i].y*ny+ucat[k-1][j][i].z*nz;
328 ubcs[k][j][i].x = ucat[k-1][j][i].x-Un*nx;
329 ubcs[k][j][i].y = ucat[k-1][j][i].y-Un*ny;
330 ubcs[k][j][i].z = ucat[k-1][j][i].z-Un*nz;
331 }
332 }
333 }
334 }
335/* ==================================================================================== */
336 /* Rheology */
337 /* ==================================================================================== */
338
339 // PetscPrintf(PETSC_COMM_WORLD, "moving plate velocity for rheology setup is %le \n",U_bc);
340
341 if (user->bctype[2]==13){
342 if (ys==0){
343 j=0;
344 for (k=lzs; k<lze; k++) {
345 for (i=lxs; i<lxe; i++) {
346 ubcs[k][j][i].x = 0.;
347 //ubcs[k][j][i].x = -U_bc;
348 ubcs[k][j][i].y = 0.;
349 ubcs[k][j][i].z =-U_bc;
350 //ubcs[k][j][i].z =0.0;
351 }
352 }
353 }
354 }
355 if (user->bctype[3]==13){
356 if (ye==my){
357 j=ye-1;
358 for (k=lzs; k<lze; k++) {
359 for (i=lxs; i<lxe; i++) {
360 ubcs[k][j][i].x = 0.;
361 //ubcs[k][j][i].x = U_bc;
362 ubcs[k][j][i].y = 0.;
363 ubcs[k][j][i].z = U_bc;
364 //ubcs[k][j][i].z =0.0;
365 }
366 }
367 }
368 }
369 if (user->bctype[4]==13){
370 if (zs==0){
371 k=0;
372 for (j=lys; j<lye; j++) {
373 for (i=lxs; i<lxe; i++) {
374 ubcs[k][j][i].x =-U_bc;
375 ubcs[k][j][i].y = 0.;
376 ubcs[k][j][i].z = 0.;
377 }
378 }
379 }
380 }
381 if (user->bctype[5]==13){
382 if (ze==mz){
383 k=ze-1;
384 for (j=lys; j<lye; j++) {
385 for (i=lxs; i<lxe; i++) {
386 ubcs[k][j][i].x = U_bc;
387 ubcs[k][j][i].y = 0.;
388 ubcs[k][j][i].z = 0.;
389 }
390 }
391 }
392 }
393/* ================================================================================== */
394 // boundary conditions on ghost nodes
395/* ================================================================================== */
396//Mohsen Aug 2012
397 if (xs==0 && user->bctype[0]!=7) {
398 i = xs;
399 for (k=zs; k<ze; k++) {
400 for (j=ys; j<ye; j++) {
401 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j][i+1].x;
402 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j][i+1].y;
403 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j][i+1].z;
404 }
405 }
406 }
407
408 if (xe==mx && user->bctype[0]!=7) {
409 i = xe-1;
410 for (k=zs; k<ze; k++) {
411 for (j=ys; j<ye; j++) {
412 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j][i-1].x;
413 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j][i-1].y;
414 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j][i-1].z;
415 }
416 }
417 }
418
419 if (ys==0 && user->bctype[2]!=7) {
420 j = ys;
421 for (k=zs; k<ze; k++) {
422 for (i=xs; i<xe; i++) {
423 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j+1][i].x;
424 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j+1][i].y;
425 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j+1][i].z;
426 }
427 }
428 }
429
430 if (ye==my && user->bctype[2]!=7) {
431 j = ye-1;
432 for (k=zs; k<ze; k++) {
433 for (i=xs; i<xe; i++) {
434 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j-1][i].x;
435 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j-1][i].y;
436 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j-1][i].z;
437 }
438 }
439 }
440
441 if (zs==0 && user->bctype[4]!=7) {
442 k = zs;
443 for (j=ys; j<ye; j++) {
444 for (i=xs; i<xe; i++) {
445 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k+1][j][i].x;
446 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k+1][j][i].y;
447 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k+1][j][i].z;
448 }
449 }
450 }
451
452 if (ze==mz && user->bctype[4]!=7) {
453 k = ze-1;
454 for (j=ys; j<ye; j++) {
455 for (i=xs; i<xe; i++) {
456 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k-1][j][i].x;
457 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k-1][j][i].y;
458 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k-1][j][i].z;
459 }
460 }
461 }
462 DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs);
463 DMDAVecRestoreArray(fda, user->Ucat,&ucat);
464/* ================================================================================== */
465/* Periodic BC Mohsen Aug 2012
466/* ================================================================================== */
467
468 if (user->bctype[0]==7 || user->bctype[2]==7 || user->bctype[4]==7 ){
469
470 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
471 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
472
473 DMDAVecGetArray(fda, user->lUcat, &lucat);
474 DMDAVecGetArray(fda, user->Ucat, &ucat);
475
476 if (user->bctype[0]==7 || user->bctype[1]==7){
477 if (xs==0){
478 i=xs;
479 for (k=zs; k<ze; k++) {
480 for (j=ys; j<ye; j++) {
481 if(k>0 && k<user->KM && j>0 && j<user->JM){
482 ucat[k][j][i]=lucat[k][j][i-2];
483 }
484 }
485 }
486 }
487 }
488 if (user->bctype[2]==7 || user->bctype[3]==7){
489 if (ys==0){
490 j=ys;
491 for (k=zs; k<ze; k++) {
492 for (i=xs; i<xe; i++) {
493 if(k>0 && k<user->KM && i>0 && i<user->IM){
494 ucat[k][j][i]=lucat[k][j-2][i];
495 }
496 }
497 }
498 }
499 }
500 if (user->bctype[4]==7 || user->bctype[5]==7){
501 if (zs==0){
502 k=zs;
503 for (j=ys; j<ye; j++) {
504 for (i=xs; i<xe; i++) {
505 if(j>0 && j<user->JM && i>0 && i<user->IM){
506 ucat[k][j][i]=lucat[k-2][j][i];
507 }
508 }
509 }
510 }
511 }
512 if (user->bctype[0]==7 || user->bctype[1]==7){
513 if (xe==mx){
514 i=mx-1;
515 for (k=zs; k<ze; k++) {
516 for (j=ys; j<ye; j++) {
517 if(k>0 && k<user->KM && j>0 && j<user->JM){
518 ucat[k][j][i]=lucat[k][j][i+2];
519 }
520 }
521 }
522 }
523 }
524 if (user->bctype[2]==7 || user->bctype[3]==7){
525 if (ye==my){
526 j=my-1;
527 for (k=zs; k<ze; k++) {
528 for (i=xs; i<xe; i++) {
529 if(k>0 && k<user->KM && i>0 && i<user->IM){
530 ucat[k][j][i]=lucat[k][j+2][i];
531 }
532 }
533 }
534 }
535 }
536 if (user->bctype[4]==7 || user->bctype[5]==7){
537 if (ze==mz){
538 k=mz-1;
539 for (j=ys; j<ye; j++) {
540 for (i=xs; i<xe; i++) {
541 if(j>0 && j<user->JM && i>0 && i<user->IM){
542 ucat[k][j][i].x=lucat[k+2][j][i].x;
543 }
544 }
545 }
546 }
547 }
548
549 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
550 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
551
552 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
553 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
554 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
555/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
556 }
557
558 // velocity on the corner point
559 DMDAVecGetArray(fda, user->Ucat, &ucat);
560
561 if (zs==0 ) {
562 k=0;
563 if (xs==0) {
564 i=0;
565 for (j=ys; j<ye; j++) {
566 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j][i+1].x);
567 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j][i+1].y);
568 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j][i+1].z);
569 }
570 }
571 if (xe == mx) {
572 i=mx-1;
573 for (j=ys; j<ye; j++) {
574 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j][i-1].x);
575 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j][i-1].y);
576 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j][i-1].z);
577 }
578 }
579
580 if (ys==0) {
581 j=0;
582 for (i=xs; i<xe; i++) {
583 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j+1][i].x);
584 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j+1][i].y);
585 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j+1][i].z);
586 }
587 }
588
589 if (ye==my) {
590 j=my-1;
591 for (i=xs; i<xe; i++) {
592 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j-1][i].x);
593 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j-1][i].y);
594 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j-1][i].z);
595 }
596 }
597 }
598
599 if (ze==mz) {
600 k=mz-1;
601 if (xs==0) {
602 i=0;
603 for (j=ys; j<ye; j++) {
604 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j][i+1].x);
605 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j][i+1].y);
606 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j][i+1].z);
607 }
608 }
609 if (xe == mx) {
610 i=mx-1;
611 for (j=ys; j<ye; j++) {
612 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j][i-1].x);
613 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j][i-1].y);
614 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j][i-1].z);
615 }
616 }
617 if (ys==0) {
618 j=0;
619 for (i=xs; i<xe; i++) {
620 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j+1][i].x);
621 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j+1][i].y);
622 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j+1][i].z);
623 }
624 }
625
626 if (ye==my) {
627 j=my-1;
628 for (i=xs; i<xe; i++) {
629 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j-1][i].x);
630 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j-1][i].y);
631 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j-1][i].z);
632 }
633 }
634 }
635
636 if (ys==0) {
637 j=0;
638 if (xs==0) {
639 i=0;
640 for (k=zs; k<ze; k++) {
641 ucat[k][j][i].x = 0.5*(ucat[k][j+1][i].x+ucat[k][j][i+1].x);
642 ucat[k][j][i].y = 0.5*(ucat[k][j+1][i].y+ucat[k][j][i+1].y);
643 ucat[k][j][i].z = 0.5*(ucat[k][j+1][i].z+ucat[k][j][i+1].z);
644 }
645 }
646
647 if (xe==mx) {
648 i=mx-1;
649 for (k=zs; k<ze; k++) {
650 ucat[k][j][i].x = 0.5*(ucat[k][j+1][i].x+ucat[k][j][i-1].x);
651 ucat[k][j][i].y = 0.5*(ucat[k][j+1][i].y+ucat[k][j][i-1].y);
652 ucat[k][j][i].z = 0.5*(ucat[k][j+1][i].z+ucat[k][j][i-1].z);
653 }
654 }
655 }
656
657 if (ye==my) {
658 j=my-1;
659 if (xs==0) {
660 i=0;
661 for (k=zs; k<ze; k++) {
662 ucat[k][j][i].x = 0.5*(ucat[k][j-1][i].x+ucat[k][j][i+1].x);
663 ucat[k][j][i].y = 0.5*(ucat[k][j-1][i].y+ucat[k][j][i+1].y);
664 ucat[k][j][i].z = 0.5*(ucat[k][j-1][i].z+ucat[k][j][i+1].z);
665 }
666 }
667
668 if (xe==mx) {
669 i=mx-1;
670 for (k=zs; k<ze; k++) {
671 ucat[k][j][i].x = 0.5*(ucat[k][j-1][i].x+ucat[k][j][i-1].x);
672 ucat[k][j][i].y = 0.5*(ucat[k][j-1][i].y+ucat[k][j][i-1].y);
673 ucat[k][j][i].z = 0.5*(ucat[k][j-1][i].z+ucat[k][j][i-1].z);
674 }
675 }
676 }
677
678 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
679
680 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
681 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
682
683 if (user->bctype[0]==7 || user->bctype[2]==7 || user->bctype[4]==7 ){
684 //i-direction
685 DMDAVecGetArray(fda, user->lUcat, &lucat);
686 DMDAVecGetArray(fda, user->Ucat, &ucat);
687
688 if (user->bctype[0]==7){
689 if (xs==0){
690 i=xs;
691 for (k=zs; k<ze; k++) {
692 for (j=ys; j<ye; j++) {
693 ucat[k][j][i]=lucat[k][j][i-2];
694 }
695 }
696 }
697 }
698 if (user->bctype[1]==7){
699 if (xe==mx){
700 i=xe-1;
701 for (k=zs; k<ze; k++) {
702 for (j=ys; j<ye; j++) {
703 ucat[k][j][i]=lucat[k][j][i+2];
704 }
705 }
706 }
707 }
708 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
709 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
710
711 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
712/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
713
714 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
715 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
716
717 //j-direction
718
719 DMDAVecGetArray(fda, user->lUcat, &lucat);
720 DMDAVecGetArray(fda, user->Ucat, &ucat);
721
722 if (user->bctype[2]==7){
723 if (ys==0){
724 j=ys;
725 for (k=zs; k<ze; k++) {
726 for (i=xs; i<xe; i++) {
727 ucat[k][j][i]=lucat[k][j-2][i];
728 }
729 }
730 }
731 }
732
733 if (user->bctype[3]==7){
734 if (ye==my){
735 j=my-1;
736 for (k=zs; k<ze; k++) {
737 for (i=xs; i<xe; i++) {
738 ucat[k][j][i]=lucat[k][j+2][i];
739 }
740 }
741 }
742 }
743
744 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
745 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
746
747 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
748/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
749
750 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
751 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
752
753 //k-direction
754
755 DMDAVecGetArray(fda, user->lUcat, &lucat);
756 DMDAVecGetArray(fda, user->Ucat, &ucat);
757
758 if (user->bctype[4]==7){
759 if (zs==0){
760 k=zs;
761 for (j=ys; j<ye; j++) {
762 for (i=xs; i<xe; i++) {
763 ucat[k][j][i]=lucat[k-2][j][i];
764 }
765 }
766 }
767 }
768 if (user->bctype[5]==7){
769 if (ze==mz){
770 k=mz-1;
771 for (j=ys; j<ye; j++) {
772 for (i=xs; i<xe; i++) {
773 ucat[k][j][i].x=lucat[k+2][j][i].x;
774 }
775 }
776 }
777 }
778
779 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
780 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
781
782 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
783/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
784
785 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
786 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
787 }
788
789 DMDAVecRestoreArray(fda, user->lCsi, &csi);
790 DMDAVecRestoreArray(fda, user->lEta, &eta);
791 DMDAVecRestoreArray(fda, user->lZet, &zet);
792
793 PROFILE_FUNCTION_END;
794 PetscFunctionReturn(0);
795
796}
CalculateFaceNormalAndArea
PetscErrorCode CalculateFaceNormalAndArea(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
Computes the unit normal vectors and areas of the three faces of a computational cell.
Definition Metric.c:257
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:731
UserCtx::lNvert
Vec lNvert
Definition variables.h:688
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:664
SimCtx::rans
PetscInt rans
Definition variables.h:609
UserCtx::lZet
Vec lZet
Definition variables.h:705
UserCtx::lUstar
Vec lUstar
Definition variables.h:683
BCS::Ubcs
Vec Ubcs
Boundary condition velocity values. (Comment: "An ugly hack, waste of memory")
Definition variables.h:121
Cmpnts::x
PetscScalar x
Definition variables.h:101
UserCtx::Bcs
BCS Bcs
Definition variables.h:682
UserCtx::lCsi
Vec lCsi
Definition variables.h:705
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::Ucat
Vec Ucat
Definition variables.h:688
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:610
UserCtx::bctype
PetscInt bctype[6]
Definition variables.h:684
SimCtx::U_bc
PetscReal U_bc
Definition variables.h:604
SimCtx::step
PetscInt step
Definition variables.h:546
UserCtx::lAj
Vec lAj
Definition variables.h:705
UserCtx::info
DMDALocalInfo info
Definition variables.h:668
UserCtx::lUcat
Vec lUcat
Definition variables.h:688
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:705
UserCtx::Cent
Vec Cent
Definition variables.h:705
SimCtx::les
PetscInt les
Definition variables.h:609
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:538
wall_function
void wall_function(UserCtx *user, double sc, double sb, Cmpnts Ua, Cmpnts Uc, Cmpnts *Ub, PetscReal *ustar, double nx, double ny, double nz)
Definition wallfunction.c:245
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Gidx()

static PetscInt Gidx ( PetscInt  i,
PetscInt  j,
PetscInt  k,
UserCtx *  user 
)
static

Definition at line 833 of file poisson.c.

835{
836 PetscInt nidx;
837 DMDALocalInfo info = user->info;
838
839 PetscInt mx = info.mx, my = info.my;
840
841 AO ao;
842 DMDAGetAO(user->da, &ao);
843 nidx=i+j*mx+k*mx*my;
844
845 AOApplicationToPetsc(ao,1,&nidx);
846
847 return (nidx);
848}
Here is the caller graph for this function:

◆ GridRestriction()

static PetscErrorCode GridRestriction ( PetscInt  i,
PetscInt  j,
PetscInt  k,
PetscInt *  ih,
PetscInt *  jh,
PetscInt *  kh,
UserCtx *  user 
)
static

Definition at line 853 of file poisson.c.

856{
857 PetscFunctionBeginUser;
858 PROFILE_FUNCTION_BEGIN;
859 if ((user->isc)) {
860 *ih = i;
861 }
862 else {
863 *ih = 2 * i;
864 }
865
866 if ((user->jsc)) {
867 *jh = j;
868 }
869 else {
870 *jh = 2 * j;
871 }
872
873 if ((user->ksc)) {
874 *kh = k;
875 }
876 else {
877 *kh = 2 * k;
878 }
879
880 PROFILE_FUNCTION_END;
881 PetscFunctionReturn(0);
882}
UserCtx::isc
PetscInt isc
Definition variables.h:674
UserCtx::ksc
PetscInt ksc
Definition variables.h:674
UserCtx::jsc
PetscInt jsc
Definition variables.h:674
Here is the caller graph for this function:

◆ Projection()

PetscErrorCode Projection ( UserCtx *  user)

Performs the projection step to enforce an incompressible velocity field.

Corrects the contravariant velocity field Ucont to be divergence-free using the gradient of the pressure correction field Phi.

This function executes the final "correction" stage of a fractional-step (projection) method for solving the incompressible Navier-Stokes equations. After solving the Pressure Poisson Equation to get a pressure correction Phi, this function uses Phi to correct the intermediate velocity field, making it divergence-free.

The main steps are:

  1. Calculate Pressure Gradient: At each velocity point (on the cell faces), it computes the gradient of the pressure correction field (∇Φ). This is done using finite differences on a generalized curvilinear grid.
  2. Correct Velocity Field: It updates the contravariant velocity components Ucont by subtracting the scaled pressure gradient. The update rule is: U_new = U_intermediate - Δt * G * ∇Φ, where G is the metric tensor g_ij that transforms the gradient into contravariant components.
  3. Boundary-Aware Gradients: The gradient calculation is highly sensitive to boundaries. The code uses complex, dynamic stencils that switch from centered differences to one-sided differences if a standard stencil would use a point inside an immersed solid (nvert > 0.1). This preserves accuracy at the fluid-solid interface.
  4. Periodic Boundary Updates: Includes dedicated loops to correctly calculate the pressure gradient at the domain edges for periodic boundary conditions.
  5. Finalization: After correcting Ucont, it calls helper functions to convert the velocity back to Cartesian coordinates (Contra2Cart) and update all ghost cell values.
Parameters
[in,out]userPointer to the UserCtx struct, containing simulation context, grid metrics, and the velocity/pressure fields.
Returns
PetscErrorCode 0 on success, or a non-zero error code on failure.

Definition at line 948 of file poisson.c.

949{
950 PetscFunctionBeginUser;
951 PROFILE_FUNCTION_BEGIN;
952 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering Projection step to correct velocity field.\n");
953
954 //================================================================================
955 // Section 1: Initialization and Data Acquisition
956 //================================================================================
957
958 // --- Get simulation and grid context ---
959 SimCtx *simCtx = user->simCtx;
960 DM da = user->da, fda = user->fda;
961 DMDALocalInfo info = user->info;
962
963 // --- Grid dimensions ---
964 PetscInt mx = info.mx, my = info.my, mz = info.mz;
965 PetscInt xs = info.xs, xe = info.xs + info.xm;
966 PetscInt ys = info.ys, ye = info.ys + info.ym;
967 PetscInt zs = info.zs, ze = info.zs + info.zm;
968
969 // --- Loop bounds (excluding outer ghost layers) ---
970 PetscInt lxs = (xs == 0) ? xs + 1 : xs;
971 PetscInt lxe = (xe == mx) ? xe - 1 : xe;
972 PetscInt lys = (ys == 0) ? ys + 1 : ys;
973 PetscInt lye = (ye == my) ? ye - 1 : ye;
974 PetscInt lzs = (zs == 0) ? zs + 1 : zs;
975 PetscInt lze = (ze == mz) ? ze - 1 : ze;
976
977 // --- Get direct pointer access to grid metric and field data ---
978 Cmpnts ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
979 PetscReal ***iaj, ***jaj, ***kaj, ***p, ***nvert;
980 Cmpnts ***ucont;
981 DMDAVecGetArray(fda, user->lICsi, &icsi); DMDAVecGetArray(fda, user->lIEta, &ieta); DMDAVecGetArray(fda, user->lIZet, &izet);
982 DMDAVecGetArray(fda, user->lJCsi, &jcsi); DMDAVecGetArray(fda, user->lJEta, &jeta); DMDAVecGetArray(fda, user->lJZet, &jzet);
983 DMDAVecGetArray(fda, user->lKCsi, &kcsi); DMDAVecGetArray(fda, user->lKEta, &keta); DMDAVecGetArray(fda, user->lKZet, &kzet);
984 DMDAVecGetArray(da, user->lIAj, &iaj); DMDAVecGetArray(da, user->lJAj, &jaj); DMDAVecGetArray(da, user->lKAj, &kaj);
985 DMDAVecGetArray(da, user->lNvert, &nvert);
986 DMDAVecGetArray(da, user->lPhi, &p); // Note: using lPhi, which is the pressure correction
987 //DMDAVecGetArray(da,user->lP,&p);
988 DMDAVecGetArray(fda, user->Ucont, &ucont);
989
990 // --- Constants for clarity ---
991 const PetscReal IBM_FLUID_THRESHOLD = 0.1;
992 const PetscReal scale = simCtx->dt * simCtx->st / COEF_TIME_ACCURACY;
993
994 LOG_ALLOW(GLOBAL,LOG_DEBUG," Starting velocity correction: Scale = %le .\n",scale);
995
996 //================================================================================
997 // Section 2: Correct Velocity Components
998 //================================================================================
999 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Calculating pressure gradients and correcting velocity components.\n");
1000
1001 // --- Main loop over interior domain points ---
1002 for (PetscInt k = lzs; k < lze; k++) {
1003 for (PetscInt j = lys; j < lye; j++) {
1004 for (PetscInt i = lxs; i < lxe; i++) {
1005
1006 // --- Correct U_contravariant (x-component of velocity) ---
1007 PetscInt i_end = (user->bctype[0] == 7) ? mx - 1 : mx - 2;
1008 if (i < i_end) {
1009
1010 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k][j][i + 1] > IBM_FLUID_THRESHOLD)) {
1011 // Compute pressure derivatives (dp/d_csi, dp/d_eta, dp/d_zet) at the i-face
1012
1013 PetscReal dpdc = p[k][j][i + 1] - p[k][j][i];
1014 PetscReal dpde = 0.0, dpdz = 0.0;
1015
1016 // Boundary-aware stencil for dp/d_eta
1017 if ((j==my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1018 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1019 dpde = (p[k][j][i] + p[k][j][i+1] -
1020 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1021 }
1022 }
1023
1024 else if ((j==my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1025 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) { dpde = (p[k][j][i] + p[k][j][i+1] - p[k][j-1][i] - p[k][j-1][i+1]) * 0.5; }
1026 }
1027
1028 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1029 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
1030 }
1031
1032 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1033 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
1034 }
1035
1036 else { dpde = (p[k][j+1][i] + p[k][j+1][i+1] - p[k][j-1][i] - p[k][j-1][i+1]) * 0.25; }
1037
1038 // Boundary-aware stencil for dp/d_zet
1039 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1040 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) { dpdz = (p[k][j][i] + p[k][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.5; }
1041 }
1042
1043 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1044 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) { dpdz = (p[k][j][i] + p[k][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.5; }
1045 }
1046
1047 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1048 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
1049 }
1050
1051 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1052 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k][j][i] - p[k][j][i+1]) * 0.5; }
1053 }
1054
1055 else { dpdz = (p[k+1][j][i] + p[k+1][j][i+1] - p[k-1][j][i] - p[k-1][j][i+1]) * 0.25; }
1056
1057 // Apply the correction: U_new = U_old - dt * (g11*dpdc + g12*dpde + g13*dpdz)
1058
1059
1060
1061 PetscReal grad_p_x = (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x + icsi[k][j][i].y * icsi[k][j][i].y
1062 + icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
1063 dpde * (ieta[k][j][i].x * icsi[k][j][i].x + ieta[k][j][i].y * icsi[k][j][i].y
1064 + ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
1065 dpdz * (izet[k][j][i].x * icsi[k][j][i].x + izet[k][j][i].y * icsi[k][j][i].y
1066 + izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i]);
1067
1068 PetscReal correction = grad_p_x*scale;
1069 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Csi Direction: %le.\n",correction);
1070 ucont[k][j][i].x -= correction;
1071
1072 }
1073 }
1074
1075 // --- Correct V_contravariant (y-component of velocity) ---
1076 PetscInt j_end = (user->bctype[2] == 7) ? my - 1 : my - 2;
1077 if (j < j_end) {
1078 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k][j + 1][i] > IBM_FLUID_THRESHOLD)) {
1079 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1080 dpde = p[k][j + 1][i] - p[k][j][i];
1081
1082 // Boundary-aware stencil for dp/d_csi
1083 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1084 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) { dpdc = (p[k][j][i] + p[k][j+1][i] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.5; }
1085 } else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1086 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) { dpdc = (p[k][j][i] + p[k][j+1][i] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.5; }
1087 } else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1088 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
1089 } else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1090 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
1091 } else { dpdc = (p[k][j][i+1] + p[k][j+1][i+1] - p[k][j][i-1] - p[k][j+1][i-1]) * 0.25; }
1092
1093 // Boundary-aware stencil for dp/d_zet
1094 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1095 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) { dpdz = (p[k][j][i] + p[k][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.5; }
1096 } else if ((k == mz-2 || k==1 ) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1097 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) { dpdz = (p[k][j][i] + p[k][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.5; }
1098 } else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1099 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
1100 } else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1101 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k][j][i] - p[k][j+1][i]) * 0.5; }
1102 } else { dpdz = (p[k+1][j][i] + p[k+1][j+1][i] - p[k-1][j][i] - p[k-1][j+1][i]) * 0.25; }
1103
1104 PetscReal grad_p_y = (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x + jcsi[k][j][i].y * jeta[k][j][i].y + jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
1105 dpde * (jeta[k][j][i].x * jeta[k][j][i].x + jeta[k][j][i].y * jeta[k][j][i].y + jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
1106 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x + jzet[k][j][i].y * jeta[k][j][i].y + jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i]);
1107
1108 PetscReal correction = grad_p_y*scale;
1109 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Eta Direction: %le.\n",correction);
1110 ucont[k][j][i].y -= correction;
1111 }
1112 }
1113
1114 // --- Correct W_contravariant (z-component of velocity) ---
1115 PetscInt k_end = (user->bctype[4] == 7) ? mz - 1 : mz - 2;
1116 if (k < k_end) {
1117 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k + 1][j][i] > IBM_FLUID_THRESHOLD)) {
1118 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1119 dpdz = p[k + 1][j][i] - p[k][j][i];
1120
1121 // Boundary-aware stencil for dp/d_csi
1122 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1123 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) { dpdc = (p[k][j][i] + p[k+1][j][i] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.5; }
1124 } else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1125 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) { dpdc = (p[k][j][i] + p[k+1][j][i] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.5; }
1126 } else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1127 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
1128 } else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1129 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
1130 } else { dpdc = (p[k][j][i+1] + p[k+1][j][i+1] - p[k][j][i-1] - p[k+1][j][i-1]) * 0.25; }
1131
1132 // Boundary-aware stencil for dp/d_eta
1133 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1134 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) { dpde = (p[k][j][i] + p[k+1][j][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.5; }
1135 } else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1136 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) { dpde = (p[k][j][i] + p[k+1][j][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.5; }
1137 } else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1138 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
1139 } else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1140 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j][i] - p[k+1][j][i]) * 0.5; }
1141 } else { dpde = (p[k][j+1][i] + p[k+1][j+1][i] - p[k][j-1][i] - p[k+1][j-1][i]) * 0.25; }
1142
1143 PetscReal grad_p_z = (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x + kcsi[k][j][i].y * kzet[k][j][i].y + kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
1144 dpde * (keta[k][j][i].x * kzet[k][j][i].x + keta[k][j][i].y * kzet[k][j][i].y + keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
1145 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x + kzet[k][j][i].y * kzet[k][j][i].y + kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i]);
1146
1147 // ========================= DEBUG PRINT =========================
1148 LOG_LOOP_ALLOW_EXACT(GLOBAL, LOG_DEBUG, k, 5,
1149 "[k=%d, j=%d, i=%d] ---- Neighbor Pressures ----\n"
1150 " Central Z-Neighbors: p[k+1][j][i] = %g | p[k][j][i] = %g\n"
1151 " Eta-Stencil (Y-dir): p[k][j-1][i] = %g, p[k+1][j-1][i] = %g | p[k][j+1][i] = %g, p[k+1][j+1][i] = %g\n"
1152 " Csi-Stencil (X-dir): p[k][j][i-1] = %g, p[k+1][j][i-1] = %g | p[k][j][i+1] = %g, p[k+1][j][i+1] = %g\n",
1153 k, j, i,
1154 p[k + 1][j][i], p[k][j][i],
1155 p[k][j - 1][i], p[k + 1][j - 1][i], p[k][j + 1][i], p[k + 1][j + 1][i],
1156 p[k][j][i - 1], p[k + 1][j][i - 1], p[k][j][i + 1], p[k + 1][j][i + 1]);
1157 // ======================= END DEBUG PRINT =======================
1158
1159 LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," dpdc: %le | dpde: %le | dpdz: %le.\n",dpdc,dpde,dpdz);
1160 PetscReal correction = grad_p_z*scale;
1161 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Zet Direction: %le.\n",correction);
1162 ucont[k][j][i].z -= correction;
1163 }
1164 }
1165 }
1166 }
1167 }
1168
1169 // --- Explicit correction for periodic boundaries (if necessary) ---
1170 // The main loop handles the interior, but this handles the first physical layer at periodic boundaries.
1171 // Note: This logic is largely duplicated from the main loop and could be merged, but is preserved for fidelity.
1172 if (user->bctype[0] == 7 && xs == 0) {
1173 for (PetscInt k=lzs; k<lze; k++) {
1174 for (PetscInt j=lys; j<lye; j++) {
1175 PetscInt i=xs;
1176
1177 PetscReal dpdc = p[k][j][i+1] - p[k][j][i];
1178
1179 PetscReal dpde = 0.;
1180 PetscReal dpdz = 0.;
1181
1182 if ((j==my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1183 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1184 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1185 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1186 }
1187 }
1188 else if ((j==my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1189 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1190 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1191 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1192 }
1193 }
1194 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1195 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1196 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1197 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1198 }
1199 }
1200 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1201 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1202 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1203 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1204 }
1205 }
1206 else {
1207 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1208 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1209 }
1210
1211 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1212 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1213 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1214 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1215 }
1216 }
1217 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1218 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1219 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1220 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1221 }
1222 }
1223 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1224 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1225 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1226 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1227 }
1228 }
1229 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1230 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1231 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1232 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1233 }
1234 }
1235 else {
1236 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1237 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1238 }
1239
1240
1241
1242 if (!(nvert[k][j][i] + nvert[k][j][i+1])) {
1243 ucont[k][j][i].x -=
1244 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1245 icsi[k][j][i].y * icsi[k][j][i].y +
1246 icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
1247 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1248 ieta[k][j][i].y * icsi[k][j][i].y +
1249 ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
1250 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1251 izet[k][j][i].y * icsi[k][j][i].y +
1252 izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i])
1253 * scale;
1254
1255 }
1256 }
1257 }
1258 }
1259 if (user->bctype[2] == 7 && ys == 0) {
1260
1261 for (PetscInt k=lzs; k<lze; k++) {
1262 for (PetscInt i=lxs; i<lxe; i++) {
1263 PetscInt j=ys;
1264
1265 PetscReal dpdc = 0.;
1266 PetscReal dpdz = 0.;
1267 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1268 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1269 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1270 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1271 }
1272 }
1273 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1274 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1275 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1276 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1277 }
1278 }
1279 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1280 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1281 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1282 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1283 }
1284 }
1285 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1286 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1287 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1288 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1289 }
1290 }
1291 else {
1292 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1293 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1294 }
1295
1296 PetscReal dpde = p[k][j+1][i] - p[k][j][i];
1297
1298 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1299 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1300 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1301 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1302 }
1303 }
1304 else if ((k == mz-2 || k==1 ) && user->bctype[0]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1305 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1306 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1307 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1308 }
1309 }
1310 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1311 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1312 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1313 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1314 }
1315 }
1316 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1317 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1318 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1319 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1320 }
1321 }
1322 else {
1323 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1324 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1325 }
1326
1327 if (!(nvert[k][j][i] + nvert[k][j+1][i])) {
1328 ucont[k][j][i].y -=
1329 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1330 jcsi[k][j][i].y * jeta[k][j][i].y +
1331 jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
1332 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1333 jeta[k][j][i].y * jeta[k][j][i].y +
1334 jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
1335 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1336 jzet[k][j][i].y * jeta[k][j][i].y +
1337 jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i])
1338 * scale;
1339 }
1340 }
1341 }
1342 }
1343
1344 if (user->bctype[4] == 7 && zs == 0) {
1345 for (PetscInt j=lys; j<lye; j++) {
1346 for (PetscInt i=lxs; i<lxe; i++) {
1347
1348 PetscInt k=zs;
1349 PetscReal dpdc = 0.;
1350 PetscReal dpde = 0.;
1351
1352 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1353 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1354 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1355 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1356 }
1357 }
1358 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1359 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1360 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1361 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1362 }
1363 }
1364 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1365 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1366 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1367 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1368 }
1369 }
1370 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1371 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1372 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1373 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1374 }
1375 }
1376 else {
1377 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1378 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1379 }
1380
1381 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1382 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1383 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1384 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1385 }
1386 }
1387 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1388 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1389 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1390 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1391 }
1392 }
1393 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1394 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1395 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1396 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1397 }
1398 }
1399 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1400 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1401 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1402 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1403 }
1404 }
1405 else {
1406 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1407 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1408 }
1409
1410 PetscReal dpdz = p[k+1][j][i] - p[k][j][i];
1411
1412 if (!(nvert[k][j][i] + nvert[k+1][j][i])) {
1413
1414 ucont[k][j][i].z -=
1415 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1416 kcsi[k][j][i].y * kzet[k][j][i].y +
1417 kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
1418 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1419 keta[k][j][i].y * kzet[k][j][i].y +
1420 keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
1421 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1422 kzet[k][j][i].y * kzet[k][j][i].y +
1423 kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i])
1424 * scale;
1425
1426 }
1427 }
1428 }
1429 }
1430
1431 // --- Optional: Enforce specific mass flow rate for channel flow simulations ---
1432 if (simCtx->channelz) {
1433 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Applying channel flow mass flux correction.\n");
1434 // DMDAVecGetArray(fda, user->lCsi, &csi); DMDAVecGetArray(fda, user->lEta, &eta); DMDAVecGetArray(fda, user->lZet, &zet);
1435 //DMDAVecGetArray(da, user->lAj, &aj);
1436 // This entire block was commented out in the original code. It calculates the current
1437 // total flux in the z-direction, compares it to a target flux (user->simCtx->FluxInSum),
1438 // and computes a uniform velocity correction `ratio` to enforce the target.
1439 // The implementation for applying the `ratio` was also commented out.
1440 // Preserving this as-is.
1441 }
1442
1443 //================================================================================
1444 // Section 3: Finalization and Cleanup
1445 //================================================================================
1446
1447 // --- Restore access to all PETSc vector arrays ---
1448 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
1449 // DMDAVecRestoreArray(fda, user->lCsi, &csi); DMDAVecRestoreArray(fda, user->lEta, &eta); DMDAVecRestoreArray(fda, user->lZet, &zet);
1450 //DMDAVecRestoreArray(da, user->lAj, &aj);
1451 DMDAVecRestoreArray(fda, user->lICsi, &icsi); DMDAVecRestoreArray(fda, user->lIEta, &ieta); DMDAVecRestoreArray(fda, user->lIZet, &izet);
1452 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi); DMDAVecRestoreArray(fda, user->lJEta, &jeta); DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1453 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi); DMDAVecRestoreArray(fda, user->lKEta, &keta); DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1454 DMDAVecRestoreArray(da, user->lIAj, &iaj); DMDAVecRestoreArray(da, user->lJAj, &jaj); DMDAVecRestoreArray(da, user->lKAj, &kaj);
1455 DMDAVecRestoreArray(da, user->lP, &p);
1456 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1457
1458 // --- Update ghost cells for the newly corrected velocity field ---
1459 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating ghost cells for corrected velocity.\n");
1460 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
1461 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
1462
1463 // --- Convert velocity to Cartesian and update ghost nodes ---
1464 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Converting velocity to Cartesian and finalizing ghost nodes.\n");
1465 Contra2Cart(user);
1466 GhostNodeVelocity(user);
1467
1468 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting Projection step.\n");
1469 PROFILE_FUNCTION_END;
1470 PetscFunctionReturn(0);
1471}
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
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:201
LOG_LOOP_ALLOW_EXACT
#define LOG_LOOP_ALLOW_EXACT(scope, level, var, val, fmt,...)
Logs a custom message if a variable equals a specific value.
Definition logging.h:349
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
GhostNodeVelocity
static PetscErrorCode GhostNodeVelocity(UserCtx *user)
Definition poisson.c:8
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2084
UserCtx::lIEta
Vec lIEta
Definition variables.h:707
UserCtx::lIZet
Vec lIZet
Definition variables.h:707
SimCtx::channelz
PetscInt channelz
Definition variables.h:593
SimCtx::st
PetscReal st
Definition variables.h:581
UserCtx::lIAj
Vec lIAj
Definition variables.h:707
UserCtx::lKEta
Vec lKEta
Definition variables.h:709
SimCtx::dt
PetscReal dt
Definition variables.h:552
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:708
UserCtx::Ucont
Vec Ucont
Definition variables.h:688
UserCtx::lPhi
Vec lPhi
Definition variables.h:688
UserCtx::lKZet
Vec lKZet
Definition variables.h:709
UserCtx::lJEta
Vec lJEta
Definition variables.h:708
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:709
UserCtx::lJZet
Vec lJZet
Definition variables.h:708
UserCtx::lUcont
Vec lUcont
Definition variables.h:688
UserCtx::lICsi
Vec lICsi
Definition variables.h:707
UserCtx::lJAj
Vec lJAj
Definition variables.h:708
UserCtx::lKAj
Vec lKAj
Definition variables.h:709
COEF_TIME_ACCURACY
#define COEF_TIME_ACCURACY
Coefficient controlling the temporal accuracy scheme (e.g., 1.5 for BDF2).
Definition variables.h:57
Here is the call graph for this function:
Here is the caller graph for this function:

◆ UpdatePressure()

PetscErrorCode UpdatePressure ( UserCtx *  user)

Updates the pressure field and handles periodic boundary conditions.

Updates the pressure field P with the pressure correction Phi computed by the Poisson solver.

This function is a core part of the projection method in a CFD solver. It is called after the Pressure Poisson Equation has been solved to obtain a pressure correction field, Phi.

The function performs two main operations:

  1. Core Pressure Update: It updates the main pressure field P by adding the pressure correction Phi at every grid point in the fluid domain. The operation is P_new = P_old + Phi.
  2. Periodic Boundary Synchronization: If any of the domain boundaries are periodic (bctype == 7), this function manually updates the pressure values in the ghost cells. It copies values from the physical cells on one side of the domain to the corresponding ghost cells on the opposite side. This ensures that subsequent calculations involving pressure gradients are correct across periodic boundaries. The refactored version consolidates the original code's redundant updates into a single, efficient pass.
Parameters
[in,out]userPointer to the UserCtx struct, containing simulation context and the pressure vectors P and Phi.
Returns
PetscErrorCode 0 on success, or a non-zero error code on failure.

Definition at line 1502 of file poisson.c.

1503{
1504 PetscFunctionBeginUser;
1505 PROFILE_FUNCTION_BEGIN;
1506 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering UpdatePressure.\n");
1507
1508 //================================================================================
1509 // Section 1: Initialization and Data Acquisition
1510 //================================================================================
1511 DM da = user->da;
1512 DMDALocalInfo info = user->info;
1513
1514 // Local grid extents for the main update loop
1515 PetscInt xs = info.xs, xe = info.xs + info.xm;
1516 PetscInt ys = info.ys, ye = info.ys + info.ym;
1517 PetscInt zs = info.zs, ze = info.zs + info.zm;
1518 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1519
1520 PetscInt i,j,k;
1521
1522 // --- Get direct pointer access to PETSc vector data for performance ---
1523 PetscReal ***p, ***phi;
1524 DMDAVecGetArray(da, user->P, &p);
1525 DMDAVecGetArray(da, user->Phi, &phi);
1526
1527 //================================================================================
1528 // Section 2: Core Pressure Update
1529 //================================================================================
1530 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Performing core pressure update (P_new = P_old + Phi).\n");
1531 for (PetscInt k = zs; k < ze; k++) {
1532 for (PetscInt j = ys; j < ye; j++) {
1533 for (PetscInt i = xs; i < xe; i++) {
1534 // This is the fundamental pressure update in a projection method.
1535 p[k][j][i] += phi[k][j][i];
1536 }
1537 }
1538 }
1539
1540 // Restore arrays now that the core computation is done.
1541 DMDAVecRestoreArray(da, user->Phi, &phi);
1542 DMDAVecRestoreArray(da, user->P, &p);
1543
1544
1545 //================================================================================
1546 // Section 3: Handle Periodic Boundary Condition Synchronization
1547 //================================================================================
1548 // This block is executed only if at least one boundary is periodic.
1549 // The original code contained many redundant Get/Restore and update calls.
1550 // This refactored version performs the same effective logic but in a single,
1551 // efficient pass, which is numerically equivalent and much cleaner.
1552 if (user->bctype[0] == 7 || user->bctype[1] == 7 ||
1553 user->bctype[2] == 7 || user->bctype[3] == 7 ||
1554 user->bctype[4] == 7 || user->bctype[5] == 7)
1555 {
1556 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Synchronizing ghost cells for periodic boundaries.\n");
1557
1558 // First, update the local vectors (including ghost regions) with the new global data.
1559 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
1560 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
1561 DMGlobalToLocalBegin(da, user->Phi, INSERT_VALUES, user->lPhi);
1562 DMGlobalToLocalEnd(da, user->Phi, INSERT_VALUES, user->lPhi);
1563
1564 // Get pointers to all necessary local and global arrays ONCE.
1565 PetscReal ***lp, ***lphi;
1566 DMDAVecGetArray(da, user->lP, &lp);
1567 DMDAVecGetArray(da, user->lPhi, &lphi);
1568 DMDAVecGetArray(da, user->P, &p);
1569 DMDAVecGetArray(da, user->Phi, &phi);
1570
1571 // --- X-Direction Periodic Update ---
1572 if (user->bctype[0] == 7 || user->bctype[1] == 7) {
1573 // Update left boundary physical cells from right boundary ghost cells
1574 if (xs == 0) {
1575 PetscInt i = 0;
1576 for (k = zs; k < ze; k++) {
1577 for (j = ys; j < ye; j++) {
1578 // Note: Accessing lp[...][i-2] reads from the ghost cell region.
1579 p[k][j][i] = lp[k][j][i - 2];
1580 phi[k][j][i] = lphi[k][j][i - 2];
1581 }
1582 }
1583 }
1584 // Update right boundary physical cells from left boundary ghost cells
1585 if (xe == mx) {
1586 PetscInt i = mx - 1;
1587 for (k = zs; k < ze; k++) {
1588 for (j = ys; j < ye; j++) {
1589 p[k][j][i] = lp[k][j][i + 2];
1590 phi[k][j][i] = lphi[k][j][i + 2];
1591 }
1592 }
1593 }
1594 }
1595
1596 // --- Y-Direction Periodic Update ---
1597 if (user->bctype[2] == 7 || user->bctype[3] == 7) {
1598 // Update bottom boundary
1599 if (ys == 0) {
1600 PetscInt j = 0;
1601 for (k = zs; k < ze; k++) {
1602 for (i = xs; i < xe; i++) {
1603 p[k][j][i] = lp[k][j - 2][i];
1604 phi[k][j][i] = lphi[k][j - 2][i];
1605 }
1606 }
1607 }
1608 // Update top boundary
1609 if (ye == my) {
1610 PetscInt j = my - 1;
1611 for (k = zs; k < ze; k++) {
1612 for (i = xs; i < xe; i++) {
1613 p[k][j][i] = lp[k][j + 2][i];
1614 phi[k][j][i] = lphi[k][j + 2][i];
1615 }
1616 }
1617 }
1618 }
1619
1620 // --- Z-Direction Periodic Update ---
1621 if (user->bctype[4] == 7 || user->bctype[5] == 7) {
1622 // Update front boundary
1623 if (zs == 0) {
1624 PetscInt k = 0;
1625 for (j = ys; j < ye; j++) {
1626 for (i = xs; i < xe; i++) {
1627 p[k][j][i] = lp[k - 2][j][i];
1628 phi[k][j][i] = lphi[k - 2][j][i];
1629 }
1630 }
1631 }
1632 // Update back boundary
1633 if (ze == mz) {
1634 PetscInt k = mz - 1;
1635 for (j = ys; j < ye; j++) {
1636 for (i = xs; i < xe; i++) {
1637 p[k][j][i] = lp[k + 2][j][i];
1638 phi[k][j][i] = lphi[k + 2][j][i];
1639 }
1640 }
1641 }
1642 }
1643
1644 // Restore all arrays ONCE.
1645 DMDAVecRestoreArray(da, user->lP, &lp);
1646 DMDAVecRestoreArray(da, user->lPhi, &lphi);
1647 DMDAVecRestoreArray(da, user->P, &p);
1648 DMDAVecRestoreArray(da, user->Phi, &phi);
1649
1650 // After manually updating the physical boundary cells, we must call
1651 // DMGlobalToLocal again to ensure all processes have the updated ghost
1652 // values for the *next* function that needs them.
1653 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
1654 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
1655 DMGlobalToLocalBegin(da, user->Phi, INSERT_VALUES, user->lPhi);
1656 DMGlobalToLocalEnd(da, user->Phi, INSERT_VALUES, user->lPhi);
1657 }
1658
1659 //================================================================================
1660 // Section 4: Final Cleanup (pointers already restored)
1661 //================================================================================
1662
1663 UpdateLocalGhosts(user,"P");
1664 UpdateLocalGhosts(user,"Phi");
1665
1666 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting UpdatePressure.\n");
1667 PROFILE_FUNCTION_END;
1668 PetscFunctionReturn(0);
1669}
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1091
UserCtx::Phi
Vec Phi
Definition variables.h:688
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mymatmultadd()

static PetscErrorCode mymatmultadd ( Mat  mat,
Vec  v1,
Vec  v2 
)
static

Definition at line 1673 of file poisson.c.

1674{
1675 Vec vt;
1676 VecDuplicate(v2, &vt);
1677 MatMult(mat, v1, vt);
1678 VecAYPX(v2, 1., vt);
1679 VecDestroy(&vt);
1680 return(0);
1681}

◆ PoissonNullSpaceFunction()

PetscErrorCode PoissonNullSpaceFunction ( MatNullSpace  nullsp,
Vec  X,
void *  ctx 
)

The callback function for PETSc's MatNullSpace object.

This function removes the null space from the Poisson solution vector by ensuring the average pressure is zero, which is necessary for problems with pure Neumann boundary conditions.

Parameters
nullspThe MatNullSpace context.
XThe vector to be corrected.
ctxA void pointer to the UserCtx.
Returns
PetscErrorCode 0 on success.

Definition at line 1686 of file poisson.c.

1687{
1688 UserCtx *user = (UserCtx*)ctx;
1689
1690 DM da = user->da;
1691
1692 DMDALocalInfo info = user->info;
1693 PetscInt xs = info.xs, xe = info.xs + info.xm;
1694 PetscInt ys = info.ys, ye = info.ys + info.ym;
1695 PetscInt zs = info.zs, ze = info.zs + info.zm;
1696 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1697 PetscInt lxs, lxe, lys, lye, lzs, lze;
1698
1699 PetscReal ***x, ***nvert;
1700 PetscInt i, j, k;
1701
1702/* /\* First remove a constant from the Vec field X *\/ */
1703
1704
1705 /* Then apply boundary conditions */
1706 DMDAVecGetArray(da, X, &x);
1707 DMDAVecGetArray(da, user->lNvert, &nvert);
1708
1709 lxs = xs; lxe = xe;
1710 lys = ys; lye = ye;
1711 lzs = zs; lze = ze;
1712
1713 if (xs==0) lxs = xs+1;
1714 if (ys==0) lys = ys+1;
1715 if (zs==0) lzs = zs+1;
1716
1717 if (xe==mx) lxe = xe-1;
1718 if (ye==my) lye = ye-1;
1719 if (ze==mz) lze = ze-1;
1720
1721 PetscReal lsum, sum;
1722 PetscReal lnum, num;
1723
1724 if (user->multinullspace) PetscPrintf(PETSC_COMM_WORLD, "MultiNullSpace!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1725 if (!user->multinullspace) {
1726 lsum = 0;
1727 lnum = 0;
1728 for (k=lzs; k<lze; k++) {
1729 for (j=lys; j<lye; j++) {
1730 for (i=lxs; i<lxe; i++) {
1731 if (nvert[k][j][i] < 0.1) {
1732 lsum += x[k][j][i];
1733 lnum ++;
1734 }
1735 }
1736 }
1737 }
1738
1739 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1740 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1741 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1742/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1743 sum = sum / (-1.0 * num);
1744
1745 for (k=lzs; k<lze; k++) {
1746 for (j=lys; j<lye; j++) {
1747 for (i=lxs; i<lxe; i++) {
1748 if (nvert[k][j][i] < 0.1) {
1749 x[k][j][i] +=sum;
1750 }
1751 }
1752 }
1753 }
1754 }
1755 else {
1756 lsum = 0;
1757 lnum = 0;
1758 for (j=lys; j<lye; j++) {
1759 for (i=lxs; i<lxe; i++) {
1760 for (k=lzs; k<lze; k++) {
1761 if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1762 lsum += x[k][j][i];
1763 lnum ++;
1764 }
1765 }
1766 }
1767 }
1768 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1769 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1770 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1771/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1772 sum /= -num;
1773 for (j=lys; j<lye; j++) {
1774 for (i=lxs; i<lxe; i++) {
1775 for (k=lzs; k<lze; k++) {
1776 if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1777 x[k][j][i] += sum;
1778 }
1779 }
1780 }
1781 }
1782
1783 lsum = 0;
1784 lnum = 0;
1785 for (j=lys; j<lye; j++) {
1786 for (i=lxs; i<lxe; i++) {
1787 for (k=lzs; k<lze; k++) {
1788 if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1789 lsum += x[k][j][i];
1790 lnum ++;
1791 }
1792 }
1793 }
1794 }
1795 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1796 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1797 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1798/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1799 sum /= -num;
1800 for (j=lys; j<lye; j++) {
1801 for (i=lxs; i<lxe; i++) {
1802 for (k=lzs; k<lze; k++) {
1803 if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1804 x[k][j][i] += sum;
1805 }
1806 }
1807 }
1808 }
1809
1810 } //if multinullspace
1811 if (zs == 0) {
1812 k = 0;
1813 for (j=ys; j<ye; j++) {
1814 for (i=xs; i<xe; i++) {
1815 x[k][j][i] = 0.;
1816 }
1817 }
1818 }
1819
1820 if (ze == mz) {
1821 k = mz-1;
1822 for (j=ys; j<ye; j++) {
1823 for (i=xs; i<xe; i++) {
1824 x[k][j][i] = 0.;
1825 }
1826 }
1827 }
1828
1829 if (ys == 0) {
1830 j = 0;
1831 for (k=zs; k<ze; k++) {
1832 for (i=xs; i<xe; i++) {
1833 x[k][j][i] = 0.;
1834 }
1835 }
1836 }
1837
1838 if (ye == my) {
1839 j = my-1;
1840 for (k=zs; k<ze; k++) {
1841 for (i=xs; i<xe; i++) {
1842 x[k][j][i] = 0.;
1843 }
1844 }
1845 }
1846
1847 if (xs == 0) {
1848 i = 0;
1849 for (k=zs; k<ze; k++) {
1850 for (j=ys; j<ye; j++) {
1851 x[k][j][i] = 0.;
1852 }
1853 }
1854 }
1855
1856 if (xe == mx) {
1857 i = mx-1;
1858 for (k=zs; k<ze; k++) {
1859 for (j=ys; j<ye; j++) {
1860 x[k][j][i] = 0.;
1861 }
1862 }
1863 }
1864
1865 for (k=zs; k<ze; k++) {
1866 for (j=ys; j<ye; j++) {
1867 for (i=xs; i<xe; i++) {
1868 if (nvert[k][j][i] > 0.1)
1869 x[k][j][i] = 0.;
1870 }
1871 }
1872 }
1873 DMDAVecRestoreArray(da, X, &x);
1874 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1875
1876 return 0;
1877}
UserCtx::KSKE
PetscInt * KSKE
Definition variables.h:697
UserCtx::multinullspace
PetscBool multinullspace
Definition variables.h:698
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:661
Here is the caller graph for this function:

◆ MyInterpolation()

PetscErrorCode MyInterpolation ( Mat  A,
Vec  X,
Vec  F 
)

The callback function for the multigrid interpolation operator (MatShell).

Defines the coarse-to-fine grid transfer for the pressure correction.

Parameters
AThe shell matrix context.
XThe coarse-grid source vector.
FThe fine-grid destination vector.
Returns
PetscErrorCode 0 on success.

Definition at line 1879 of file poisson.c.

1880{
1881 UserCtx *user;
1882
1883 MatShellGetContext(A, (void**)&user);
1884
1885
1886
1887 DM da = user->da;
1888
1889 DM da_c = *user->da_c;
1890
1891 DMDALocalInfo info = user->info;
1892 PetscInt xs = info.xs, xe = info.xs + info.xm;
1893 PetscInt ys = info.ys, ye = info.ys + info.ym;
1894 PetscInt zs = info.zs, ze = info.zs + info.zm;
1895 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1896 PetscInt lxs, lxe, lys, lye, lzs, lze;
1897
1898 PetscReal ***f, ***x, ***nvert, ***nvert_c;
1899 PetscInt i, j, k, ic, jc, kc, ia, ja, ka;
1900
1901 lxs = xs; lxe = xe;
1902 lys = ys; lye = ye;
1903 lzs = zs; lze = ze;
1904
1905 if (xs==0) lxs = xs+1;
1906 if (ys==0) lys = ys+1;
1907 if (zs==0) lzs = zs+1;
1908
1909 if (xe==mx) lxe = xe-1;
1910 if (ye==my) lye = ye-1;
1911 if (ze==mz) lze = ze-1;
1912
1913
1914 DMDAVecGetArray(da, F, &f);
1915
1916
1917 Vec lX;
1918 DMCreateLocalVector(da_c, &lX);
1919
1920 DMGlobalToLocalBegin(da_c, X, INSERT_VALUES, lX);
1921 DMGlobalToLocalEnd(da_c, X, INSERT_VALUES, lX);
1922 DMDAVecGetArray(da_c, lX, &x);
1923
1924 DMDAVecGetArray(da, user->lNvert, &nvert);
1925 DMDAVecGetArray(da_c, *(user->lNvert_c), &nvert_c);
1926 for (k=lzs; k<lze; k++) {
1927 for (j=lys; j<lye; j++) {
1928 for (i=lxs; i<lxe; i++) {
1929
1930 GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user);
1931
1932 f[k][j][i] = (x[kc ][jc ][ic ] * 9 +
1933 x[kc ][jc+ja][ic ] * 3 +
1934 x[kc ][jc ][ic+ia] * 3 +
1935 x[kc ][jc+ja][ic+ia]) * 3./64. +
1936 (x[kc+ka][jc ][ic ] * 9 +
1937 x[kc+ka][jc+ja][ic ] * 3 +
1938 x[kc+ka][jc ][ic+ia] * 3 +
1939 x[kc+ka][jc+ja][ic+ia]) /64.;
1940 }
1941 }
1942 }
1943
1944 for (k=zs; k<ze; k++) {
1945 for (j=ys; j<ye; j++) {
1946 for (i=xs; i<xe; i++) {
1947
1948 if (i==0) {
1949 f[k][j][i] = 0.;//-f[k][j][i+1];
1950 }
1951 else if (i==mx-1) {
1952 f[k][j][i] = 0.;//-f[k][j][i-1];
1953 }
1954 else if (j==0) {
1955 f[k][j][i] = 0.;//-f[k][j+1][i];
1956 }
1957 else if (j==my-1) {
1958 f[k][j][i] = 0.;//-f[k][j-1][i];
1959 }
1960 else if (k==0) {
1961 f[k][j][i] = 0.;//-f[k+1][j][i];
1962 }
1963 else if (k==mz-1) {
1964 f[k][j][i] = 0.;//-f[k-1][j][i];
1965 }
1966 if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
1967
1968 }
1969 }
1970 }
1971
1972 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1973 DMDAVecRestoreArray(da_c, *(user->lNvert_c), &nvert_c);
1974
1975 DMDAVecRestoreArray(da_c, lX, &x);
1976
1977 VecDestroy(&lX);
1978 DMDAVecRestoreArray(da, F, &f);
1979
1980
1981
1982 return 0;
1983
1984}
GridInterpolation
#define GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user)
Definition poisson.c:798
UserCtx::da_c
DM * da_c
Definition variables.h:723
UserCtx::lNvert_c
Vec * lNvert_c
Definition variables.h:724
Here is the caller graph for this function:

◆ RestrictResidual_SolidAware()

static PetscErrorCode RestrictResidual_SolidAware ( Mat  A,
Vec  X,
Vec  F 
)
static

Definition at line 1986 of file poisson.c.

1987{
1988 UserCtx *user;
1989 MatShellGetContext(A, (void**)&user);
1990
1991 DM da = user->da;
1992 DM da_f = *user->da_f;
1993
1994 DMDALocalInfo info;
1995 DMDAGetLocalInfo(da, &info);
1996 PetscInt xs = info.xs, xe = info.xs + info.xm;
1997 PetscInt ys = info.ys, ye = info.ys + info.ym;
1998 PetscInt zs = info.zs, ze = info.zs + info.zm;
1999 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2000
2001 PetscReal ***f, ***x, ***nvert;
2002 PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
2003
2004 DMDAVecGetArray(da, F, &f);
2005
2006 Vec lX;
2007 DMCreateLocalVector(da_f, &lX);
2008 DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
2009 DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);
2010 DMDAVecGetArray(da_f, lX, &x);
2011
2012 DMDAVecGetArray(da, user->lNvert, &nvert);
2013
2014 PetscReal ***nvert_f;
2015 DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
2016
2017 if ((user->isc)) ia = 0;
2018 else ia = 1;
2019
2020 if ((user->jsc)) ja = 0;
2021 else ja = 1;
2022
2023 if ((user->ksc)) ka = 0;
2024 else ka = 1;
2025
2026 for (k=zs; k<ze; k++) {
2027 for (j=ys; j<ye; j++) {
2028 for (i=xs; i<xe; i++) {
2029 // --- CORRECTED LOGIC ---
2030 // First, check if the current point is a boundary point.
2031 // If it is, it does not contribute to the coarse grid residual.
2032 if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1 || nvert[k][j][i] > 0.1) {
2033 f[k][j][i] = 0.0;
2034 }
2035 // Only if it's a true interior fluid point, perform the restriction.
2036 else {
2037 GridRestriction(i, j, k, &ih, &jh, &kh, user);
2038 f[k][j][i] = 0.125 *
2039 (x[kh ][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih ]) +
2040 x[kh ][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih-ia]) +
2041 x[kh ][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih ]) +
2042 x[kh-ka][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih ]) +
2043 x[kh ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih-ia]) +
2044 x[kh-ka][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih ]) +
2045 x[kh-ka][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih-ia]) +
2046 x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
2047 }
2048 }
2049 }
2050 }
2051
2052 DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
2053 DMDAVecRestoreArray(da_f, lX, &x);
2054 VecDestroy(&lX);
2055 DMDAVecRestoreArray(da, F, &f);
2056 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2057
2058 return 0;
2059}
GridRestriction
static PetscErrorCode GridRestriction(PetscInt i, PetscInt j, PetscInt k, PetscInt *ih, PetscInt *jh, PetscInt *kh, UserCtx *user)
Definition poisson.c:853
UserCtx::user_f
UserCtx * user_f
Definition variables.h:722
UserCtx::da_f
DM * da_f
Definition variables.h:723
Here is the call graph for this function:
Here is the caller graph for this function:

◆ MyRestriction()

PetscErrorCode MyRestriction ( Mat  A,
Vec  X,
Vec  F 
)

The callback function for the multigrid restriction operator (MatShell).

Defines the fine-to-coarse grid transfer for the Poisson residual.

Parameters
AThe shell matrix context.
XThe fine-grid source vector.
FThe coarse-grid destination vector.
Returns
PetscErrorCode 0 on success.

Definition at line 2061 of file poisson.c.

2062{
2063 UserCtx *user;
2064
2065 MatShellGetContext(A, (void**)&user);
2066
2067
2068 DM da = user->da;
2069
2070 DM da_f = *user->da_f;
2071
2072 DMDALocalInfo info;
2073 DMDAGetLocalInfo(da, &info);
2074 PetscInt xs = info.xs, xe = info.xs + info.xm;
2075 PetscInt ys = info.ys, ye = info.ys + info.ym;
2076 PetscInt zs = info.zs, ze = info.zs + info.zm;
2077 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2078 // PetscInt lxs, lxe, lys, lye, lzs, lze;
2079
2080 PetscReal ***f, ***x, ***nvert;
2081 PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
2082
2083 DMDAVecGetArray(da, F, &f);
2084
2085 Vec lX;
2086
2087 DMCreateLocalVector(da_f, &lX);
2088 DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
2089 DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);
2090 DMDAVecGetArray(da_f, lX, &x);
2091
2092 DMDAVecGetArray(da, user->lNvert, &nvert);
2093
2094 PetscReal ***nvert_f;
2095 DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
2096
2097 if ((user->isc)) ia = 0;
2098 else ia = 1;
2099
2100 if ((user->jsc)) ja = 0;
2101 else ja = 1;
2102
2103 if ((user->ksc)) ka = 0;
2104 else ka = 1;
2105
2106 for (k=zs; k<ze; k++) {
2107 for (j=ys; j<ye; j++) {
2108 for (i=xs; i<xe; i++) {
2109 if (k==0) {
2110 f[k][j][i] = 0.;
2111 }
2112 else if (k==mz-1) {
2113 f[k][j][i] = 0.;
2114 }
2115 else if (j==0) {
2116 f[k][j][i] = 0.;
2117 }
2118 else if (j==my-1) {
2119 f[k][j][i] = 0.;
2120 }
2121 else if (i==0) {
2122 f[k][j][i] = 0.;
2123 }
2124 else if (i==mx-1) {
2125 f[k][j][i] = 0.;
2126 }
2127 else {
2128 GridRestriction(i, j, k, &ih, &jh, &kh, user);
2129 f[k][j][i] = 0.125 *
2130 (x[kh ][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih ]) +
2131 x[kh ][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih-ia]) +
2132 x[kh ][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih ]) +
2133 x[kh-ka][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih ]) +
2134 x[kh ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih-ia]) +
2135 x[kh-ka][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih ]) +
2136 x[kh-ka][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih-ia]) +
2137 x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
2138
2139
2140
2141 if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
2142 }
2143 }
2144 }
2145 }
2146
2147
2148 DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
2149
2150 DMDAVecRestoreArray(da_f, lX, &x);
2151 VecDestroy(&lX);
2152
2153 DMDAVecRestoreArray(da, F, &f);
2154 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2155
2156
2157 return 0;
2158}
Here is the call graph for this function:

◆ PoissonLHSNew()

PetscErrorCode PoissonLHSNew ( UserCtx *  user)

Assembles the Left-Hand Side (LHS) matrix for the Pressure Poisson Equation.

Assembles the Left-Hand-Side (LHS) matrix (Laplacian operator) for the Poisson equation on a single grid level.

This function constructs the sparse matrix A (the LHS) for the linear system Ax = B, which is the Pressure Poisson Equation (PPE). The matrix A represents a discrete version of the negative Laplacian operator (-∇²), tailored for a general curvilinear, staggered grid.

The assembly process is highly complex and follows these main steps:

  1. Matrix Initialization: On the first call, it allocates a sparse PETSc matrix A pre-allocating space for a 19-point stencil per row. On subsequent calls, it simply zeroes out the existing matrix.
  2. Metric Tensor Calculation: It first computes the 9 components of the metric tensor (g11, g12, ..., g33) at the cell faces. These g_ij coefficients are essential for defining the Laplacian operator in generalized curvilinear coordinates and account for grid stretching and non-orthogonality.
  3. Stencil Assembly Loop: The function iterates through every grid point (i,j,k). For each point, it determines the matrix row entries based on its status:
    • Boundary/Solid Point: If the point is on a domain boundary or inside an immersed solid (nvert > 0.1), it sets an identity row (A(row,row) = 1), effectively removing it from the linear solve.
    • Fluid Point: For a fluid point, it computes the 19 coefficients of the finite volume stencil. This involves summing contributions from each of the six faces of the control volume around the point.
  4. Boundary-Aware Stencils: The stencil calculation is critically dependent on the state of neighboring cells. The code contains intricate logic to check if neighbors are fluid or solid. If a standard centered-difference stencil would cross into a solid, the scheme is automatically adapted to a one-sided difference to maintain accuracy at the fluid-solid interface.
  5. Matrix Finalization: After all rows corresponding to the local processor's grid points have been set, MatAssemblyBegin/End is called to finalize the global matrix, making it ready for use by a linear solver.
Parameters
[in,out]userPointer to the UserCtx struct, containing all simulation context, grid data, and the matrix A.
Returns
PetscErrorCode 0 on success, or a non-zero error code on failure.

Definition at line 2207 of file poisson.c.

2208{
2209 PetscFunctionBeginUser;
2210 PROFILE_FUNCTION_BEGIN;
2211 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonLHSNew to assemble Laplacian matrix.\n");
2212 PetscErrorCode ierr;
2213 //================================================================================
2214 // Section 1: Initialization and Data Acquisition
2215 //================================================================================
2216
2217
2218 // --- Get simulation and grid context ---
2219 SimCtx *simCtx = user->simCtx;
2220 DM da = user->da, fda = user->fda;
2221 DMDALocalInfo info = user->info;
2222 PetscInt IM = user->IM, JM = user->JM, KM = user->KM;
2223 PetscInt i,j,k;
2224
2225 // --- Grid dimensions ---
2226 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2227 PetscInt xs = info.xs, xe = info.xs + info.xm;
2228 PetscInt ys = info.ys, ye = info.ys + info.ym;
2229 PetscInt zs = info.zs, ze = info.zs + info.zm;
2230 PetscInt gxs = info.gxs, gxe = gxs + info.gxm;
2231 PetscInt gys = info.gys, gye = gys + info.gym;
2232 PetscInt gzs = info.gzs, gze = gzs + info.gzm;
2233
2234 // --- Define constants for clarity ---
2235 const PetscReal IBM_FLUID_THRESHOLD = 0.1;
2236
2237 // --- Allocate the LHS matrix A on the first call ---
2238 if (!user->assignedA) {
2239 LOG_ALLOW(GLOBAL, LOG_INFO, "First call: Creating LHS matrix 'A' with 19-point stencil preallocation.\n");
2240 PetscInt N = mx * my * mz; // Total size
2241 PetscInt M; // Local size
2242 VecGetLocalSize(user->Phi, &M);
2243 // Create a sparse AIJ matrix, preallocating for 19 non-zeros per row (d=diagonal, o=off-diagonal)
2244 MatCreateAIJ(PETSC_COMM_WORLD, M, M, N, N, 19, PETSC_NULL, 19, PETSC_NULL, &(user->A));
2245 user->assignedA = PETSC_TRUE;
2246 }
2247
2248 // Zero out matrix entries from the previous solve
2249 MatZeroEntries(user->A);
2250
2251 // --- Get direct pointer access to grid metric data ---
2252 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
2253 PetscReal ***aj, ***iaj, ***jaj, ***kaj, ***nvert;
2254 DMDAVecGetArray(fda, user->lCsi, &csi); DMDAVecGetArray(fda, user->lEta, &eta); DMDAVecGetArray(fda, user->lZet, &zet);
2255 DMDAVecGetArray(fda, user->lICsi, &icsi); DMDAVecGetArray(fda, user->lIEta, &ieta); DMDAVecGetArray(fda, user->lIZet, &izet);
2256 DMDAVecGetArray(fda, user->lJCsi, &jcsi); DMDAVecGetArray(fda, user->lJEta, &jeta); DMDAVecGetArray(fda, user->lJZet, &jzet);
2257 DMDAVecGetArray(fda, user->lKCsi, &kcsi); DMDAVecGetArray(fda, user->lKEta, &keta); DMDAVecGetArray(fda, user->lKZet, &kzet);
2258 DMDAVecGetArray(da, user->lAj, &aj); DMDAVecGetArray(da, user->lIAj, &iaj); DMDAVecGetArray(da, user->lJAj, &jaj); DMDAVecGetArray(da, user->lKAj, &kaj);
2259 DMDAVecGetArray(da, user->lNvert, &nvert);
2260
2261 // --- Create temporary vectors for the metric tensor components G_ij ---
2262 Vec G11, G12, G13, G21, G22, G23, G31, G32, G33;
2263 PetscReal ***g11, ***g12, ***g13, ***g21, ***g22, ***g23, ***g31, ***g32, ***g33;
2264 VecDuplicate(user->lAj, &G11); VecDuplicate(user->lAj, &G12); VecDuplicate(user->lAj, &G13);
2265 VecDuplicate(user->lAj, &G21); VecDuplicate(user->lAj, &G22); VecDuplicate(user->lAj, &G23);
2266 VecDuplicate(user->lAj, &G31); VecDuplicate(user->lAj, &G32); VecDuplicate(user->lAj, &G33);
2267 DMDAVecGetArray(da, G11, &g11); DMDAVecGetArray(da, G12, &g12); DMDAVecGetArray(da, G13, &g13);
2268 DMDAVecGetArray(da, G21, &g21); DMDAVecGetArray(da, G22, &g22); DMDAVecGetArray(da, G23, &g23);
2269 DMDAVecGetArray(da, G31, &g31); DMDAVecGetArray(da, G32, &g32); DMDAVecGetArray(da, G33, &g33);
2270
2271 //================================================================================
2272 // Section 2: Pre-compute Metric Tensor Coefficients (g_ij)
2273 //================================================================================
2274 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pre-computing metric tensor components (g_ij).\n");
2275 for (k = gzs; k < gze; k++) {
2276 for (j = gys; j < gye; j++) {
2277 for (i = gxs; i < gxe; i++) {
2278 // These coefficients represent the dot products of the grid's contravariant base vectors,
2279 // scaled by face area. They are the core of the Laplacian operator on a curvilinear grid.
2280 if(i>-1 && j>-1 && k>-1 && i<IM+1 && j<JM+1 && k<KM+1){
2281 g11[k][j][i] = (icsi[k][j][i].x * icsi[k][j][i].x + icsi[k][j][i].y * icsi[k][j][i].y + icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
2282 g12[k][j][i] = (ieta[k][j][i].x * icsi[k][j][i].x + ieta[k][j][i].y * icsi[k][j][i].y + ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
2283 g13[k][j][i] = (izet[k][j][i].x * icsi[k][j][i].x + izet[k][j][i].y * icsi[k][j][i].y + izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i];
2284 g21[k][j][i] = (jcsi[k][j][i].x * jeta[k][j][i].x + jcsi[k][j][i].y * jeta[k][j][i].y + jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
2285 g22[k][j][i] = (jeta[k][j][i].x * jeta[k][j][i].x + jeta[k][j][i].y * jeta[k][j][i].y + jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
2286 g23[k][j][i] = (jzet[k][j][i].x * jeta[k][j][i].x + jzet[k][j][i].y * jeta[k][j][i].y + jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i];
2287 g31[k][j][i] = (kcsi[k][j][i].x * kzet[k][j][i].x + kcsi[k][j][i].y * kzet[k][j][i].y + kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
2288 g32[k][j][i] = (keta[k][j][i].x * kzet[k][j][i].x + keta[k][j][i].y * kzet[k][j][i].y + keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
2289 g33[k][j][i] = (kzet[k][j][i].x * kzet[k][j][i].x + kzet[k][j][i].y * kzet[k][j][i].y + kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i];
2290 }
2291 }
2292 }
2293 }
2294
2295 //================================================================================
2296 // Section 3: Assemble the LHS Matrix A
2297 //================================================================================
2298 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling the LHS matrix A using a 19-point stencil.\n");
2299
2300 // --- Define domain boundaries for stencil logic, accounting for periodic BCs ---
2301 PetscInt x_str, x_end, y_str, y_end, z_str, z_end;
2302 if (user->bctype[0] == 7) { x_end = mx - 1; x_str = 0; }
2303 else { x_end = mx - 2; x_str = 1; }
2304 if (user->bctype[2] == 7) { y_end = my - 1; y_str = 0; }
2305 else { y_end = my - 2; y_str = 1; }
2306 if (user->bctype[4] == 7) { z_end = mz - 1; z_str = 0; }
2307 else { z_end = mz - 2; z_str = 1; }
2308
2309 // --- Main assembly loop over all local grid points ---
2310 for (k = zs; k < ze; k++) {
2311 for (j = ys; j < ye; j++) {
2312 for (i = xs; i < xe; i++) {
2313 PetscScalar vol[19]; // Holds the 19 stencil coefficient values for the current row
2314 PetscInt idx[19]; // Holds the 19 global column indices for the current row
2315 PetscInt row = Gidx(i, j, k, user); // Global index for the current row
2316
2317 // --- Handle Domain Boundary and Immersed Solid Points ---
2318 // For these points, we don't solve the Poisson equation. We set an identity
2319 // row (A_ii = 1) to effectively fix the pressure value (usually to 0).
2320 if (i == 0 || i == mx - 1 || j == 0 || j == my - 1 || k == 0 || k == mz - 1 || nvert[k][j][i] > IBM_FLUID_THRESHOLD) {
2321 vol[CP] = 1.0;
2322 idx[CP] = row;
2323 MatSetValues(user->A, 1, &row, 1, &idx[CP], &vol[CP], INSERT_VALUES);
2324 }
2325 // --- Handle Fluid Points ---
2326 else {
2327 for (PetscInt m = 0; m < 19; m++) {
2328 vol[m] = 0.0;
2329 }
2330
2331 /************************************************************************
2332 * EAST FACE CONTRIBUTION (between i and i+1)
2333 ************************************************************************/
2334 if (nvert[k][j][i + 1] < IBM_FLUID_THRESHOLD && i != x_end) { // East neighbor is fluid
2335 // Primary derivative term: d/d_csi (g11 * dP/d_csi)
2336 vol[CP] -= g11[k][j][i];
2337 vol[EP] += g11[k][j][i];
2338
2339 // Cross-derivative term: d/d_csi (g12 * dP/d_eta).
2340 // This requires an average of dP/d_eta. If a neighbor is solid, the stencil
2341 // dynamically switches to a one-sided difference to avoid using solid points.
2342 if ((j == my-2 && user->bctype[2]!=7) || nvert[k][j+1][i] + nvert[k][j+1][i+1] > 0.1) {
2343 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
2344 vol[CP] += g12[k][j][i] * 0.5; vol[EP] += g12[k][j][i] * 0.5;
2345 vol[SP] -= g12[k][j][i] * 0.5; vol[SE] -= g12[k][j][i] * 0.5;
2346 }
2347 }
2348 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k][j+1][i+1] > 0.1) {
2349 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
2350 vol[CP] += g12[k][j][i] * 0.5; vol[EP] += g12[k][j][i] * 0.5;
2351 vol[SP] -= g12[k][j][i] * 0.5; vol[SE] -= g12[k][j][i] * 0.5;
2352 }
2353 }
2354 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
2355 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
2356 vol[NP] += g12[k][j][i] * 0.5; vol[NE] += g12[k][j][i] * 0.5;
2357 vol[CP] -= g12[k][j][i] * 0.5; vol[EP] -= g12[k][j][i] * 0.5;
2358 }
2359 }
2360 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
2361 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
2362 vol[NP] += g12[k][j][i] * 0.5; vol[NE] += g12[k][j][i] * 0.5;
2363 vol[CP] -= g12[k][j][i] * 0.5; vol[EP] -= g12[k][j][i] * 0.5;
2364 }
2365 }
2366 else { // Centered difference
2367 vol[NP] += g12[k][j][i] * 0.25; vol[NE] += g12[k][j][i] * 0.25;
2368 vol[SP] -= g12[k][j][i] * 0.25; vol[SE] -= g12[k][j][i] * 0.25;
2369 }
2370
2371 // Cross-derivative term: d/d_csi (g13 * dP/d_zet)
2372 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
2373 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
2374 vol[CP] += g13[k][j][i] * 0.5; vol[EP] += g13[k][j][i] * 0.5;
2375 vol[BP] -= g13[k][j][i] * 0.5; vol[BE] -= g13[k][j][i] * 0.5;
2376 }
2377 }
2378 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
2379 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
2380 vol[CP] += g13[k][j][i] * 0.5; vol[EP] += g13[k][j][i] * 0.5;
2381 vol[BP] -= g13[k][j][i] * 0.5; vol[BE] -= g13[k][j][i] * 0.5;
2382 }
2383 }
2384 else if ((k == 1 && user->bctype[4]!=7) || nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
2385 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
2386 vol[TP] += g13[k][j][i] * 0.5; vol[TE] += g13[k][j][i] * 0.5;
2387 vol[CP] -= g13[k][j][i] * 0.5; vol[EP] -= g13[k][j][i] * 0.5;
2388 }
2389 }
2390 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
2391 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
2392 vol[TP] += g13[k][j][i] * 0.5; vol[TE] += g13[k][j][i] * 0.5;
2393 vol[CP] -= g13[k][j][i] * 0.5; vol[EP] -= g13[k][j][i] * 0.5;
2394 }
2395 }
2396 else { // Centered difference
2397 vol[TP] += g13[k][j][i] * 0.25; vol[TE] += g13[k][j][i] * 0.25;
2398 vol[BP] -= g13[k][j][i] * 0.25; vol[BE] -= g13[k][j][i] * 0.25;
2399 }
2400 }
2401
2402 /************************************************************************
2403 * WEST FACE CONTRIBUTION (between i-1 and i)
2404 ************************************************************************/
2405 if (nvert[k][j][i-1] < IBM_FLUID_THRESHOLD && i != x_str) {
2406 vol[CP] -= g11[k][j][i-1];
2407 vol[WP] += g11[k][j][i-1];
2408
2409 if ((j == my-2 && user->bctype[2]!=7) || nvert[k][j+1][i] + nvert[k][j+1][i-1] > 0.1) {
2410 if (nvert[k][j-1][i] + nvert[k][j-1][i-1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
2411 vol[CP] -= g12[k][j][i-1] * 0.5; vol[WP] -= g12[k][j][i-1] * 0.5;
2412 vol[SP] += g12[k][j][i-1] * 0.5; vol[SW] += g12[k][j][i-1] * 0.5;
2413 }
2414 }
2415 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k][j+1][i-1] > 0.1) {
2416 if (nvert[k][j-1][i] + nvert[k][j-1][i-1] < 0.1) {
2417 vol[CP] -= g12[k][j][i-1] * 0.5; vol[WP] -= g12[k][j][i-1] * 0.5;
2418 vol[SP] += g12[k][j][i-1] * 0.5; vol[SW] += g12[k][j][i-1] * 0.5;
2419 }
2420 }
2421 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k][j-1][i-1] > 0.1) {
2422 if (nvert[k][j+1][i] + nvert[k][j+1][i-1] < 0.1) {
2423 vol[NP] -= g12[k][j][i-1] * 0.5; vol[NW] -= g12[k][j][i-1] * 0.5;
2424 vol[CP] += g12[k][j][i-1] * 0.5; vol[WP] += g12[k][j][i-1] * 0.5;
2425 }
2426 }
2427 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i-1] > 0.1) {
2428 if (nvert[k][j+1][i] + nvert[k][j+1][i-1] < 0.1) {
2429 vol[NP] -= g12[k][j][i-1] * 0.5; vol[NW] -= g12[k][j][i-1] * 0.5;
2430 vol[CP] += g12[k][j][i-1] * 0.5; vol[WP] += g12[k][j][i-1] * 0.5;
2431 }
2432 }
2433 else {
2434 vol[NP] -= g12[k][j][i-1] * 0.25; vol[NW] -= g12[k][j][i-1] * 0.25;
2435 vol[SP] += g12[k][j][i-1] * 0.25; vol[SW] += g12[k][j][i-1] * 0.25;
2436 }
2437
2438 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i-1] > 0.1) {
2439 if (nvert[k-1][j][i] + nvert[k-1][j][i-1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
2440 vol[CP] -= g13[k][j][i-1] * 0.5; vol[WP] -= g13[k][j][i-1] * 0.5;
2441 vol[BP] += g13[k][j][i-1] * 0.5; vol[BW] += g13[k][j][i-1] * 0.5;
2442 }
2443 }
2444 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i-1] > 0.1) {
2445 if (nvert[k-1][j][i] + nvert[k-1][j][i-1] < 0.1) {
2446 vol[CP] -= g13[k][j][i-1] * 0.5; vol[WP] -= g13[k][j][i-1] * 0.5;
2447 vol[BP] += g13[k][j][i-1] * 0.5; vol[BW] += g13[k][j][i-1] * 0.5;
2448 }
2449 }
2450 else if ((k == 1 && user->bctype[4]!=7) || nvert[k-1][j][i] + nvert[k-1][j][i-1] > 0.1) {
2451 if (nvert[k+1][j][i] + nvert[k+1][j][i-1] < 0.1) {
2452 vol[TP] -= g13[k][j][i-1] * 0.5; vol[TW] -= g13[k][j][i-1] * 0.5;
2453 vol[CP] += g13[k][j][i-1] * 0.5; vol[WP] += g13[k][j][i-1] * 0.5;
2454 }
2455 }
2456 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i-1] > 0.1) {
2457 if (nvert[k+1][j][i] + nvert[k+1][j][i-1] < 0.1) {
2458 vol[TP] -= g13[k][j][i-1] * 0.5; vol[TW] -= g13[k][j][i-1] * 0.5;
2459 vol[CP] += g13[k][j][i-1] * 0.5; vol[WP] += g13[k][j][i-1] * 0.5;
2460 }
2461 }
2462 else {
2463 vol[TP] -= g13[k][j][i-1] * 0.25; vol[TW] -= g13[k][j][i-1] * 0.25;
2464 vol[BP] += g13[k][j][i-1] * 0.25; vol[BW] += g13[k][j][i-1] * 0.25;
2465 }
2466 }
2467
2468 /************************************************************************
2469 * NORTH FACE CONTRIBUTION (between j and j+1)
2470 ************************************************************************/
2471 if (nvert[k][j+1][i] < IBM_FLUID_THRESHOLD && j != y_end) {
2472 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
2473 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
2474 vol[CP] += g21[k][j][i] * 0.5; vol[NP] += g21[k][j][i] * 0.5;
2475 vol[WP] -= g21[k][j][i] * 0.5; vol[NW] -= g21[k][j][i] * 0.5;
2476 }
2477 }
2478 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
2479 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
2480 vol[CP] += g21[k][j][i] * 0.5; vol[NP] += g21[k][j][i] * 0.5;
2481 vol[WP] -= g21[k][j][i] * 0.5; vol[NW] -= g21[k][j][i] * 0.5;
2482 }
2483 }
2484 else if ((i == 1 && user->bctype[0]!=7) || nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
2485 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
2486 vol[EP] += g21[k][j][i] * 0.5; vol[NE] += g21[k][j][i] * 0.5;
2487 vol[CP] -= g21[k][j][i] * 0.5; vol[NP] -= g21[k][j][i] * 0.5;
2488 }
2489 }
2490 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
2491 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
2492 vol[EP] += g21[k][j][i] * 0.5; vol[NE] += g21[k][j][i] * 0.5;
2493 vol[CP] -= g21[k][j][i] * 0.5; vol[NP] -= g21[k][j][i] * 0.5;
2494 }
2495 }
2496 else {
2497 vol[EP] += g21[k][j][i] * 0.25; vol[NE] += g21[k][j][i] * 0.25;
2498 vol[WP] -= g21[k][j][i] * 0.25; vol[NW] -= g21[k][j][i] * 0.25;
2499 }
2500
2501 vol[CP] -= g22[k][j][i];
2502 vol[NP] += g22[k][j][i];
2503
2504 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
2505 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
2506 vol[CP] += g23[k][j][i] * 0.5; vol[NP] += g23[k][j][i] * 0.5;
2507 vol[BP] -= g23[k][j][i] * 0.5; vol[BN] -= g23[k][j][i] * 0.5;
2508 }
2509 }
2510 else if ((k == mz-2 || k==1 ) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
2511 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
2512 vol[CP] += g23[k][j][i] * 0.5; vol[NP] += g23[k][j][i] * 0.5;
2513 vol[BP] -= g23[k][j][i] * 0.5; vol[BN] -= g23[k][j][i] * 0.5;
2514 }
2515 }
2516 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
2517 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
2518 vol[TP] += g23[k][j][i] * 0.5; vol[TN] += g23[k][j][i] * 0.5;
2519 vol[CP] -= g23[k][j][i] * 0.5; vol[NP] -= g23[k][j][i] * 0.5;
2520 }
2521 }
2522 else if ((k == 1 || k==mz-2 ) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
2523 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
2524 vol[TP] += g23[k][j][i] * 0.5; vol[TN] += g23[k][j][i] * 0.5;
2525 vol[CP] -= g23[k][j][i] * 0.5; vol[NP] -= g23[k][j][i] * 0.5;
2526 }
2527 }
2528 else {
2529 vol[TP] += g23[k][j][i] * 0.25; vol[TN] += g23[k][j][i] * 0.25;
2530 vol[BP] -= g23[k][j][i] * 0.25; vol[BN] -= g23[k][j][i] * 0.25;
2531 }
2532 }
2533
2534 /************************************************************************
2535 * SOUTH FACE CONTRIBUTION (between j-1 and j)
2536 ************************************************************************/
2537 if (nvert[k][j-1][i] < IBM_FLUID_THRESHOLD && j != y_str) {
2538 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j-1][i+1] > 0.1) {
2539 if (nvert[k][j][i-1] + nvert[k][j-1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
2540 vol[CP] -= g21[k][j-1][i] * 0.5; vol[SP] -= g21[k][j-1][i] * 0.5;
2541 vol[WP] += g21[k][j-1][i] * 0.5; vol[SW] += g21[k][j-1][i] * 0.5;
2542 }
2543 }
2544 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j-1][i+1] > 0.1) {
2545 if (nvert[k][j][i-1] + nvert[k][j-1][i-1] < 0.1) {
2546 vol[CP] -= g21[k][j-1][i] * 0.5; vol[SP] -= g21[k][j-1][i] * 0.5;
2547 vol[WP] += g21[k][j-1][i] * 0.5; vol[SW] += g21[k][j-1][i] * 0.5;
2548 }
2549 }
2550 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j-1][i-1] > 0.1) {
2551 if (nvert[k][j][i+1] + nvert[k][j-1][i+1] < 0.1) {
2552 vol[EP] -= g21[k][j-1][i] * 0.5; vol[SE] -= g21[k][j-1][i] * 0.5;
2553 vol[CP] += g21[k][j-1][i] * 0.5; vol[SP] += g21[k][j-1][i] * 0.5;
2554 }
2555 }
2556 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j-1][i-1] > 0.1) {
2557 if (nvert[k][j][i+1] + nvert[k][j-1][i+1] < 0.1) {
2558 vol[EP] -= g21[k][j-1][i] * 0.5; vol[SE] -= g21[k][j-1][i] * 0.5;
2559 vol[CP] += g21[k][j-1][i] * 0.5; vol[SP] += g21[k][j-1][i] * 0.5;
2560 }
2561 }
2562 else {
2563 vol[EP] -= g21[k][j-1][i] * 0.25; vol[SE] -= g21[k][j-1][i] * 0.25;
2564 vol[WP] += g21[k][j-1][i] * 0.25; vol[SW] += g21[k][j-1][i] * 0.25;
2565 }
2566
2567 vol[CP] -= g22[k][j-1][i];
2568 vol[SP] += g22[k][j-1][i];
2569
2570 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j-1][i] > 0.1) {
2571 if (nvert[k-1][j][i] + nvert[k-1][j-1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
2572 vol[CP] -= g23[k][j-1][i] * 0.5; vol[SP] -= g23[k][j-1][i] * 0.5;
2573 vol[BP] += g23[k][j-1][i] * 0.5; vol[BS] += g23[k][j-1][i] * 0.5;
2574 }
2575 }
2576 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j-1][i] > 0.1) {
2577 if (nvert[k-1][j][i] + nvert[k-1][j-1][i] < 0.1 ) {
2578 vol[CP] -= g23[k][j-1][i] * 0.5; vol[SP] -= g23[k][j-1][i] * 0.5;
2579 vol[BP] += g23[k][j-1][i] * 0.5; vol[BS] += g23[k][j-1][i] * 0.5;
2580 }
2581 }
2582 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j-1][i] > 0.1) {
2583 if (nvert[k+1][j][i] + nvert[k+1][j-1][i] < 0.1) {
2584 vol[TP] -= g23[k][j-1][i] * 0.5; vol[TS] -= g23[k][j-1][i] * 0.5;
2585 vol[CP] += g23[k][j-1][i] * 0.5; vol[SP] += g23[k][j-1][i] * 0.5;
2586 }
2587 }
2588 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j-1][i] > 0.1) {
2589 if (nvert[k+1][j][i] + nvert[k+1][j-1][i] < 0.1) {
2590 vol[TP] -= g23[k][j-1][i] * 0.5; vol[TS] -= g23[k][j-1][i] * 0.5;
2591 vol[CP] += g23[k][j-1][i] * 0.5; vol[SP] += g23[k][j-1][i] * 0.5;
2592 }
2593 }
2594 else {
2595 vol[TP] -= g23[k][j-1][i] * 0.25; vol[TS] -= g23[k][j-1][i] * 0.25;
2596 vol[BP] += g23[k][j-1][i] * 0.25; vol[BS] += g23[k][j-1][i] * 0.25;
2597 }
2598 }
2599
2600 /************************************************************************
2601 * TOP FACE CONTRIBUTION (between k and k+1)
2602 ************************************************************************/
2603 if (nvert[k+1][j][i] < IBM_FLUID_THRESHOLD && k != z_end) {
2604 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
2605 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
2606 vol[CP] += g31[k][j][i] * 0.5; vol[TP] += g31[k][j][i] * 0.5;
2607 vol[WP] -= g31[k][j][i] * 0.5; vol[TW] -= g31[k][j][i] * 0.5;
2608 }
2609 }
2610 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
2611 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
2612 vol[CP] += g31[k][j][i] * 0.5; vol[TP] += g31[k][j][i] * 0.5;
2613 vol[WP] -= g31[k][j][i] * 0.5; vol[TW] -= g31[k][j][i] * 0.5;
2614 }
2615 }
2616 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
2617 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
2618 vol[EP] += g31[k][j][i] * 0.5; vol[TE] += g31[k][j][i] * 0.5;
2619 vol[CP] -= g31[k][j][i] * 0.5; vol[TP] -= g31[k][j][i] * 0.5;
2620 }
2621 }
2622 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
2623 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
2624 vol[EP] += g31[k][j][i] * 0.5; vol[TE] += g31[k][j][i] * 0.5;
2625 vol[CP] -= g31[k][j][i] * 0.5; vol[TP] -= g31[k][j][i] * 0.5;
2626 }
2627 }
2628 else {
2629 vol[EP] += g31[k][j][i] * 0.25; vol[TE] += g31[k][j][i] * 0.25;
2630 vol[WP] -= g31[k][j][i] * 0.25; vol[TW] -= g31[k][j][i] * 0.25;
2631 }
2632
2633 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
2634 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
2635 vol[CP] += g32[k][j][i] * 0.5; vol[TP] += g32[k][j][i] * 0.5;
2636 vol[SP] -= g32[k][j][i] * 0.5; vol[TS] -= g32[k][j][i] * 0.5;
2637 }
2638 }
2639 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
2640 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
2641 vol[CP] += g32[k][j][i] * 0.5; vol[TP] += g32[k][j][i] * 0.5;
2642 vol[SP] -= g32[k][j][i] * 0.5; vol[TS] -= g32[k][j][i] * 0.5;
2643 }
2644 }
2645 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
2646 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
2647 vol[NP] += g32[k][j][i] * 0.5; vol[TN] += g32[k][j][i] * 0.5;
2648 vol[CP] -= g32[k][j][i] * 0.5; vol[TP] -= g32[k][j][i] * 0.5;
2649 }
2650 }
2651 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
2652 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
2653 vol[NP] += g32[k][j][i] * 0.5; vol[TN] += g32[k][j][i] * 0.5;
2654 vol[CP] -= g32[k][j][i] * 0.5; vol[TP] -= g32[k][j][i] * 0.5;
2655 }
2656 }
2657 else {
2658 vol[NP] += g32[k][j][i] * 0.25; vol[TN] += g32[k][j][i] * 0.25;
2659 vol[SP] -= g32[k][j][i] * 0.25; vol[TS] -= g32[k][j][i] * 0.25;
2660 }
2661
2662 vol[CP] -= g33[k][j][i];
2663 vol[TP] += g33[k][j][i];
2664 }
2665
2666 /************************************************************************
2667 * BOTTOM FACE CONTRIBUTION (between k-1 and k)
2668 ************************************************************************/
2669 if (nvert[k-1][j][i] < IBM_FLUID_THRESHOLD && k != z_str) {
2670 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k-1][j][i+1] > 0.1) {
2671 if (nvert[k][j][i-1] + nvert[k-1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
2672 vol[CP] -= g31[k-1][j][i] * 0.5; vol[BP] -= g31[k-1][j][i] * 0.5;
2673 vol[WP] += g31[k-1][j][i] * 0.5; vol[BW] += g31[k-1][j][i] * 0.5;
2674 }
2675 }
2676 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k-1][j][i+1] > 0.1) {
2677 if (nvert[k][j][i-1] + nvert[k-1][j][i-1] < 0.1) {
2678 vol[CP] -= g31[k-1][j][i] * 0.5; vol[BP] -= g31[k-1][j][i] * 0.5;
2679 vol[WP] += g31[k-1][j][i] * 0.5; vol[BW] += g31[k-1][j][i] * 0.5;
2680 }
2681 }
2682 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k-1][j][i-1] > 0.1) {
2683 if (nvert[k][j][i+1] + nvert[k-1][j][i+1] < 0.1) {
2684 vol[EP] -= g31[k-1][j][i] * 0.5; vol[BE] -= g31[k-1][j][i] * 0.5;
2685 vol[CP] += g31[k-1][j][i] * 0.5; vol[BP] += g31[k-1][j][i] * 0.5;
2686 }
2687 }
2688 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k-1][j][i-1] > 0.1) {
2689 if (nvert[k][j][i+1] + nvert[k-1][j][i+1] < 0.1) {
2690 vol[EP] -= g31[k-1][j][i] * 0.5; vol[BE] -= g31[k-1][j][i] * 0.5;
2691 vol[CP] += g31[k-1][j][i] * 0.5; vol[BP] += g31[k-1][j][i] * 0.5;
2692 }
2693 }
2694 else {
2695 vol[EP] -= g31[k-1][j][i] * 0.25; vol[BE] -= g31[k-1][j][i] * 0.25;
2696 vol[WP] += g31[k-1][j][i] * 0.25; vol[BW] += g31[k-1][j][i] * 0.25;
2697 }
2698
2699 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k-1][j+1][i] > 0.1) {
2700 if (nvert[k][j-1][i] + nvert[k-1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
2701 vol[CP] -= g32[k-1][j][i] * 0.5; vol[BP] -= g32[k-1][j][i] * 0.5;
2702 vol[SP] += g32[k-1][j][i] * 0.5; vol[BS] += g32[k-1][j][i] * 0.5;
2703 }
2704 }
2705 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k-1][j+1][i] > 0.1) {
2706 if (nvert[k][j-1][i] + nvert[k-1][j-1][i] < 0.1) {
2707 vol[CP] -= g32[k-1][j][i] * 0.5; vol[BP] -= g32[k-1][j][i] * 0.5;
2708 vol[SP] += g32[k-1][j][i] * 0.5; vol[BS] += g32[k-1][j][i] * 0.5;
2709 }
2710 }
2711 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k-1][j-1][i] > 0.1) {
2712 if (nvert[k][j+1][i] + nvert[k-1][j+1][i] < 0.1) {
2713 vol[NP] -= g32[k-1][j][i] * 0.5; vol[BN] -= g32[k-1][j][i] * 0.5;
2714 vol[CP] += g32[k-1][j][i] * 0.5; vol[BP] += g32[k-1][j][i] * 0.5;
2715 }
2716 }
2717 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k-1][j-1][i] > 0.1) {
2718 if (nvert[k][j+1][i] + nvert[k-1][j+1][i] < 0.1) {
2719 vol[NP] -= g32[k-1][j][i] * 0.5; vol[BN] -= g32[k-1][j][i] * 0.5;
2720 vol[CP] += g32[k-1][j][i] * 0.5; vol[BP] += g32[k-1][j][i] * 0.5;
2721 }
2722 }
2723 else {
2724 vol[NP] -= g32[k-1][j][i] * 0.25; vol[BN] -= g32[k-1][j][i] * 0.25;
2725 vol[SP] += g32[k-1][j][i] * 0.25; vol[BS] += g32[k-1][j][i] * 0.25;
2726 }
2727
2728 vol[CP] -= g33[k-1][j][i];
2729 vol[BP] += g33[k-1][j][i];
2730 }
2731
2732 // --- Final scaling and insertion into the matrix ---
2733
2734 // Scale all stencil coefficients by the negative cell volume (-aj).
2735 for (PetscInt m = 0; m < 19; m++) {
2736 vol[m] *= -aj[k][j][i];
2737 }
2738
2739 // Set the global column indices for the 19 stencil points, handling periodic BCs.
2740 idx[CP] = Gidx(i, j, k, user);
2741 if (user->bctype[0]==7 && i==mx-2) idx[EP] = Gidx(1, j, k, user); else idx[EP] = Gidx(i+1, j, k, user);
2742 if (user->bctype[0]==7 && i==1) idx[WP] = Gidx(mx-2, j, k, user); else idx[WP] = Gidx(i-1, j, k, user);
2743 if (user->bctype[2]==7 && j==my-2) idx[NP] = Gidx(i, 1, k, user); else idx[NP] = Gidx(i, j+1, k, user);
2744 if (user->bctype[2]==7 && j==1) idx[SP] = Gidx(i, my-2, k, user); else idx[SP] = Gidx(i, j-1, k, user);
2745 if (user->bctype[4]==7 && k==mz-2) idx[TP] = Gidx(i, j, 1, user); else idx[TP] = Gidx(i, j, k+1, user);
2746 if (user->bctype[4]==7 && k==1) idx[BP] = Gidx(i, j, mz-2, user); else idx[BP] = Gidx(i, j, k-1, user);
2747 if (user->bctype[0]==7 && user->bctype[2]==7 && i==mx-2 && j==my-2) idx[NE] = Gidx(1, 1, k, user); else if (user->bctype[0]==7 && i==mx-2) idx[NE] = Gidx(1, j+1, k, user); else if (user->bctype[2]==7 && j==my-2) idx[NE] = Gidx(i+1, 1, k, user); else idx[NE] = Gidx(i+1, j+1, k, user);
2748 if (user->bctype[0]==7 && user->bctype[2]==7 && i==mx-2 && j==1) idx[SE] = Gidx(1, my-2, k, user); else if (user->bctype[0]==7 && i==mx-2) idx[SE] = Gidx(1, j-1, k, user); else if (user->bctype[2]==7 && j==1) idx[SE] = Gidx(i+1, my-2, k, user); else idx[SE] = Gidx(i+1, j-1, k, user);
2749 if (user->bctype[0]==7 && user->bctype[2]==7 && i==1 && j==my-2) idx[NW] = Gidx(mx-2, 1, k, user); else if (user->bctype[0]==7 && i==1) idx[NW] = Gidx(mx-2, j+1, k, user); else if (user->bctype[2]==7 && j==my-2) idx[NW] = Gidx(i-1, 1, k, user); else idx[NW] = Gidx(i-1, j+1, k, user);
2750 if (user->bctype[0]==7 && user->bctype[2]==7 && i==1 && j==1) idx[SW] = Gidx(mx-2, my-2, k, user); else if (user->bctype[0]==7 && i==1) idx[SW] = Gidx(mx-2, j-1, k, user); else if (user->bctype[2]==7 && j==1) idx[SW] = Gidx(i-1, my-2, k, user); else idx[SW] = Gidx(i-1, j-1, k, user);
2751 if (user->bctype[2]==7 && user->bctype[4]==7 && j==my-2 && k==mz-2) idx[TN] = Gidx(i, 1, 1, user); else if (user->bctype[2]==7 && j==my-2) idx[TN] = Gidx(i, 1, k+1, user); else if (user->bctype[4]==7 && k==mz-2) idx[TN] = Gidx(i, j+1, 1, user); else idx[TN] = Gidx(i, j+1, k+1, user);
2752 if (user->bctype[2]==7 && user->bctype[4]==7 && j==my-2 && k==1) idx[BN] = Gidx(i, 1, mz-2, user); else if(user->bctype[2]==7 && j==my-2) idx[BN] = Gidx(i, 1, k-1, user); else if (user->bctype[4]==7 && k==1) idx[BN] = Gidx(i, j+1, mz-2, user); else idx[BN] = Gidx(i, j+1, k-1, user);
2753 if (user->bctype[2]==7 && user->bctype[4]==7 && j==1 && k==mz-2) idx[TS] = Gidx(i, my-2, 1, user); else if (user->bctype[2]==7 && j==1) idx[TS] = Gidx(i, my-2, k+1, user); else if (user->bctype[4]==7 && k==mz-2) idx[TS] = Gidx(i, j-1, 1, user); else idx[TS] = Gidx(i, j-1, k+1, user);
2754 if (user->bctype[2]==7 && user->bctype[4]==7 && j==1 && k==1) idx[BS] = Gidx(i, my-2, mz-2, user); else if (user->bctype[2]==7 && j==1) idx[BS] = Gidx(i, my-2, k-1, user); else if (user->bctype[4]==7 && k==1) idx[BS] = Gidx(i, j-1, mz-2, user); else idx[BS] = Gidx(i, j-1, k-1, user);
2755 if (user->bctype[0]==7 && user->bctype[4]==7 && i==mx-2 && k==mz-2) idx[TE] = Gidx(1, j, 1, user); else if(user->bctype[0]==7 && i==mx-2) idx[TE] = Gidx(1, j, k+1, user); else if(user->bctype[4]==7 && k==mz-2) idx[TE] = Gidx(i+1, j, 1, user); else idx[TE] = Gidx(i+1, j, k+1, user);
2756 if (user->bctype[0]==7 && user->bctype[4]==7 && i==mx-2 && k==1) idx[BE] = Gidx(1, j, mz-2, user); else if(user->bctype[0]==7 && i==mx-2) idx[BE] = Gidx(1, j, k-1, user); else if(user->bctype[4]==7 && k==1) idx[BE] = Gidx(i+1, j, mz-2, user); else idx[BE] = Gidx(i+1, j, k-1, user);
2757 if (user->bctype[0]==7 && user->bctype[4]==7 && i==1 && k==mz-2) idx[TW] = Gidx(mx-2, j, 1, user); else if(user->bctype[0]==7 && i==1) idx[TW] = Gidx(mx-2, j, k+1, user); else if (user->bctype[4]==7 && k==mz-2) idx[TW] = Gidx(i-1, j, 1, user); else idx[TW] = Gidx(i-1, j, k+1, user);
2758 if (user->bctype[0]==7 && user->bctype[4]==7 && i==1 && k==1) idx[BW] = Gidx(mx-2, j, mz-2, user); else if (user->bctype[0]==7 && i==1) idx[BW] = Gidx(mx-2, j, k-1, user); else if (user->bctype[4]==7 && k==1) idx[BW] = Gidx(i-1, j, mz-2, user); else idx[BW] = Gidx(i-1, j, k-1, user);
2759
2760 // Insert the computed row into the matrix A.
2761 MatSetValues(user->A, 1, &row, 19, idx, vol, INSERT_VALUES);
2762 }
2763 }
2764 }
2765 }
2766
2767 //================================================================================
2768 // Section 4: Finalize Matrix and Cleanup
2769 //================================================================================
2770
2771 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finalizing matrix assembly.\n");
2772 MatAssemblyBegin(user->A, MAT_FINAL_ASSEMBLY);
2773 MatAssemblyEnd(user->A, MAT_FINAL_ASSEMBLY);
2774
2775 PetscReal max_A;
2776
2777 ierr = MatNorm(user->A,NORM_INFINITY,&max_A);CHKERRQ(ierr);
2778
2779 LOG_ALLOW(GLOBAL,LOG_DEBUG," Max value in A matrix for level %d = %le.\n",user->thislevel,max_A);
2780
2781 // if (get_log_level() >= LOG_DEBUG) {
2782 // ierr = MatView(user->A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
2783 // }
2784
2785 // --- Restore access to all PETSc vectors and destroy temporary ones ---
2786 DMDAVecRestoreArray(da, G11, &g11); DMDAVecRestoreArray(da, G12, &g12); DMDAVecRestoreArray(da, G13, &g13);
2787 DMDAVecRestoreArray(da, G21, &g21); DMDAVecRestoreArray(da, G22, &g22); DMDAVecRestoreArray(da, G23, &g23);
2788 DMDAVecRestoreArray(da, G31, &g31); DMDAVecRestoreArray(da, G32, &g32); DMDAVecRestoreArray(da, G33, &g33);
2789
2790 VecDestroy(&G11); VecDestroy(&G12); VecDestroy(&G13);
2791 VecDestroy(&G21); VecDestroy(&G22); VecDestroy(&G23);
2792 VecDestroy(&G31); VecDestroy(&G32); VecDestroy(&G33);
2793
2794 DMDAVecRestoreArray(fda, user->lCsi, &csi); DMDAVecRestoreArray(fda, user->lEta, &eta); DMDAVecRestoreArray(fda, user->lZet, &zet);
2795 DMDAVecRestoreArray(fda, user->lICsi, &icsi); DMDAVecRestoreArray(fda, user->lIEta, &ieta); DMDAVecRestoreArray(fda, user->lIZet, &izet);
2796 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi); DMDAVecRestoreArray(fda, user->lJEta, &jeta); DMDAVecRestoreArray(fda, user->lJZet, &jzet);
2797 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi); DMDAVecRestoreArray(fda, user->lKEta, &keta); DMDAVecRestoreArray(fda, user->lKZet, &kzet);
2798 DMDAVecRestoreArray(da, user->lAj, &aj); DMDAVecRestoreArray(da, user->lIAj, &iaj); DMDAVecRestoreArray(da, user->lJAj, &jaj); DMDAVecRestoreArray(da, user->lKAj, &kaj);
2799 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2800
2801 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting PoissonLHSNew.\n");
2802 PROFILE_FUNCTION_END;
2803 PetscFunctionReturn(0);
2804}
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:32
TW
#define TW
Definition poisson.c:906
SE
#define SE
Definition poisson.c:897
BN
#define BN
Definition poisson.c:901
WP
#define WP
Definition poisson.c:889
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Definition poisson.c:833
SW
#define SW
Definition poisson.c:899
BS
#define BS
Definition poisson.c:903
NE
#define NE
Definition poisson.c:896
CP
#define CP
Definition poisson.c:886
BE
#define BE
Definition poisson.c:905
BP
#define BP
Definition poisson.c:893
BW
#define BW
Definition poisson.c:907
TE
#define TE
Definition poisson.c:904
TS
#define TS
Definition poisson.c:902
NP
#define NP
Definition poisson.c:890
EP
#define EP
Definition poisson.c:888
TN
#define TN
Definition poisson.c:900
SP
#define SP
Definition poisson.c:891
TP
#define TP
Definition poisson.c:892
NW
#define NW
Definition poisson.c:898
UserCtx::KM
PetscInt KM
Definition variables.h:670
UserCtx::assignedA
PetscBool assignedA
Definition variables.h:701
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:721
UserCtx::JM
PetscInt JM
Definition variables.h:670
UserCtx::IM
PetscInt IM
Definition variables.h:670
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PoissonRHS()

PetscErrorCode PoissonRHS ( UserCtx *  user,
Vec  B 
)

Computes the Right-Hand-Side (RHS) of the Poisson equation, which is the divergence of the intermediate velocity field.

Parameters
userThe UserCtx for the grid level.
BThe PETSc Vec where the RHS result will be stored.
Returns
PetscErrorCode 0 on success.

Definition at line 2810 of file poisson.c.

2811{
2812 DMDALocalInfo info = user->info;
2813 PetscInt xs = info.xs, xe = info.xs + info.xm;
2814 PetscInt ys = info.ys, ye = info.ys + info.ym;
2815 PetscInt zs = info.zs, ze = info.zs + info.zm;
2816 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2817
2818 PetscInt i, j, k;
2819 PetscReal ***nvert, ***aj, ***rb, dt = user->simCtx->dt;
2820 struct Components{
2821 PetscReal x;
2822 PetscReal y;
2823 PetscReal z;
2824 } *** ucont;
2825
2826 PROFILE_FUNCTION_BEGIN;
2827
2828 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonRHS to compute pressure equation RHS.\n");
2829
2830 DMDAVecGetArray(user->da, B, &rb);
2831 DMDAVecGetArray(user->fda, user->lUcont, &ucont);
2832 DMDAVecGetArray(user->da, user->lNvert, &nvert);
2833 DMDAVecGetArray(user->da, user->lAj, &aj);
2834
2835
2836 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing RHS values for each cell.\n");
2837
2838 for (k=zs; k<ze; k++) {
2839 for (j=ys; j<ye; j++) {
2840 for (i=xs; i<xe; i++) {
2841
2842 if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1) {
2843 rb[k][j][i] = 0.;
2844 }
2845 else if (nvert[k][j][i] > 0.1) {
2846 rb[k][j][i] = 0;
2847 }
2848 else {
2849 rb[k][j][i] = -(ucont[k][j][i].x - ucont[k][j][i-1].x +
2850 ucont[k][j][i].y - ucont[k][j-1][i].y +
2851 ucont[k][j][i].z - ucont[k-1][j][i].z) / dt
2852 * aj[k][j][i] / user->simCtx->st * COEF_TIME_ACCURACY;
2853
2854 }
2855 }
2856 }
2857 }
2858
2859
2860 // --- Check the solvability condition for the Poisson equation ---
2861 // The global sum of the RHS (proportional to the total divergence) must be zero.
2862 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying solvability condition (sum of RHS terms).\n");
2863 PetscReal lsum=0., sum=0.;
2864
2865 for (k=zs; k<ze; k++) {
2866 for (j=ys; j<ye; j++) {
2867 for (i=xs; i<xe; i++) {
2868
2869 lsum += rb[k][j][i] / aj[k][j][i]* dt/COEF_TIME_ACCURACY;
2870
2871 }
2872 }
2873 }
2874
2875 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2876
2877 LOG_ALLOW(GLOBAL, LOG_INFO, "Global Sum of RHS (Divergence Check): %le\n", sum);
2878
2879 user->simCtx->summationRHS = sum;
2880
2881 DMDAVecRestoreArray(user->fda, user->lUcont, &ucont);
2882 DMDAVecRestoreArray(user->da, user->lNvert, &nvert);
2883 DMDAVecRestoreArray(user->da, user->lAj, &aj);
2884 DMDAVecRestoreArray(user->da, B, &rb);
2885
2886
2887 PROFILE_FUNCTION_END;
2888 return 0;
2889}
SimCtx::summationRHS
PetscReal summationRHS
Definition variables.h:637
Here is the caller graph for this function:

◆ VolumeFlux_rev()

PetscErrorCode VolumeFlux_rev ( UserCtx *  user,
PetscReal *  ibm_Flux,
PetscReal *  ibm_Area,
PetscInt  flg 
)

A specialized version of VolumeFlux, likely for reversed normals.

Parameters
userThe UserCtx for the grid level.
ibm_Flux(Output) The calculated net flux.
ibm_Area(Output) The total surface area of the IB.
flgA flag controlling the correction behavior.
Returns
PetscErrorCode 0 on success.

Definition at line 2891 of file poisson.c.

2893{
2894
2895 // --- CONTEXT ACQUISITION BLOCK ---
2896 // Get the master simulation context from the UserCtx.
2897 SimCtx *simCtx = user->simCtx;
2898
2899 // Create a local variable to mirror the legacy global for minimal code changes.
2900 const PetscInt NumberOfBodies = simCtx->NumberOfBodies;
2901 // --- END CONTEXT ACQUISITION BLOCK ---
2902
2903 DM da = user->da, fda = user->fda;
2904
2905 DMDALocalInfo info = user->info;
2906
2907 PetscInt xs = info.xs, xe = info.xs + info.xm;
2908 PetscInt ys = info.ys, ye = info.ys + info.ym;
2909 PetscInt zs = info.zs, ze = info.zs + info.zm;
2910 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2911
2912 PetscInt i, j, k;
2913 PetscInt lxs, lys, lzs, lxe, lye, lze;
2914
2915 lxs = xs; lxe = xe;
2916 lys = ys; lye = ye;
2917 lzs = zs; lze = ze;
2918
2919 if (xs==0) lxs = xs+1;
2920 if (ys==0) lys = ys+1;
2921 if (zs==0) lzs = zs+1;
2922
2923 if (xe==mx) lxe = xe-1;
2924 if (ye==my) lye = ye-1;
2925 if (ze==mz) lze = ze-1;
2926
2927 PetscReal ***nvert, ibmval=1.5;
2928 Cmpnts ***ucor, ***csi, ***eta, ***zet;
2929 DMDAVecGetArray(fda, user->Ucont, &ucor);
2930 DMDAVecGetArray(fda, user->lCsi, &csi);
2931 DMDAVecGetArray(fda, user->lEta, &eta);
2932 DMDAVecGetArray(fda, user->lZet, &zet);
2933 DMDAVecGetArray(da, user->lNvert, &nvert);
2934
2935 PetscReal libm_Flux, libm_area;
2936 libm_Flux = 0;
2937 libm_area = 0;
2938 for (k=lzs; k<lze; k++) {
2939 for (j=lys; j<lye; j++) {
2940 for (i=lxs; i<lxe; i++) {
2941 if (nvert[k][j][i] < 0.1) {
2942 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
2943 libm_Flux += ucor[k][j][i].x;
2944 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2945 csi[k][j][i].y * csi[k][j][i].y +
2946 csi[k][j][i].z * csi[k][j][i].z);
2947
2948 }
2949 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
2950 libm_Flux += ucor[k][j][i].y;
2951 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2952 eta[k][j][i].y * eta[k][j][i].y +
2953 eta[k][j][i].z * eta[k][j][i].z);
2954 }
2955 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
2956 libm_Flux += ucor[k][j][i].z;
2957 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2958 zet[k][j][i].y * zet[k][j][i].y +
2959 zet[k][j][i].z * zet[k][j][i].z);
2960 }
2961 }
2962
2963 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
2964 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
2965 libm_Flux -= ucor[k][j][i].x;
2966 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2967 csi[k][j][i].y * csi[k][j][i].y +
2968 csi[k][j][i].z * csi[k][j][i].z);
2969
2970 }
2971 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
2972 libm_Flux -= ucor[k][j][i].y;
2973 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2974 eta[k][j][i].y * eta[k][j][i].y +
2975 eta[k][j][i].z * eta[k][j][i].z);
2976 }
2977 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
2978 libm_Flux -= ucor[k][j][i].z;
2979 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2980 zet[k][j][i].y * zet[k][j][i].y +
2981 zet[k][j][i].z * zet[k][j][i].z);
2982 }
2983 }
2984
2985 }
2986 }
2987 }
2988
2989 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2990 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2991
2992 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
2993/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
2994 PetscPrintf(PETSC_COMM_WORLD, "IBMFlux REV %le %le\n", *ibm_Flux, *ibm_Area);
2995
2996 PetscReal correction;
2997
2998 if (*ibm_Area > 1.e-15) {
2999 if (flg)
3000 correction = (*ibm_Flux + user->FluxIntpSum) / *ibm_Area;
3001 else
3002 correction = *ibm_Flux / *ibm_Area;
3003 }
3004 else {
3005 correction = 0;
3006 }
3007
3008 for (k=lzs; k<lze; k++) {
3009 for (j=lys; j<lye; j++) {
3010 for (i=lxs; i<lxe; i++) {
3011 if (nvert[k][j][i] < 0.1) {
3012 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
3013 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
3014 csi[k][j][i].y * csi[k][j][i].y +
3015 csi[k][j][i].z * csi[k][j][i].z) *
3016 correction;
3017
3018 }
3019 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
3020 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
3021 eta[k][j][i].y * eta[k][j][i].y +
3022 eta[k][j][i].z * eta[k][j][i].z) *
3023 correction;
3024 }
3025 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
3026 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
3027 zet[k][j][i].y * zet[k][j][i].y +
3028 zet[k][j][i].z * zet[k][j][i].z) *
3029 correction;
3030 }
3031 }
3032
3033 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
3034 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
3035 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3036 csi[k][j][i].y * csi[k][j][i].y +
3037 csi[k][j][i].z * csi[k][j][i].z) *
3038 correction;
3039
3040 }
3041 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
3042 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3043 eta[k][j][i].y * eta[k][j][i].y +
3044 eta[k][j][i].z * eta[k][j][i].z) *
3045 correction;
3046 }
3047 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
3048 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3049 zet[k][j][i].y * zet[k][j][i].y +
3050 zet[k][j][i].z * zet[k][j][i].z) *
3051 correction;
3052 }
3053 }
3054
3055 }
3056 }
3057 }
3058
3059
3060
3061 libm_Flux = 0;
3062 libm_area = 0;
3063 for (k=lzs; k<lze; k++) {
3064 for (j=lys; j<lye; j++) {
3065 for (i=lxs; i<lxe; i++) {
3066 if (nvert[k][j][i] < 0.1) {
3067 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
3068 libm_Flux += ucor[k][j][i].x;
3069 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3070 csi[k][j][i].y * csi[k][j][i].y +
3071 csi[k][j][i].z * csi[k][j][i].z);
3072
3073 }
3074 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
3075 libm_Flux += ucor[k][j][i].y;
3076 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3077 eta[k][j][i].y * eta[k][j][i].y +
3078 eta[k][j][i].z * eta[k][j][i].z);
3079 }
3080 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
3081 libm_Flux += ucor[k][j][i].z;
3082 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3083 zet[k][j][i].y * zet[k][j][i].y +
3084 zet[k][j][i].z * zet[k][j][i].z);
3085 }
3086 }
3087
3088 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
3089 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
3090 libm_Flux -= ucor[k][j][i].x;
3091 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3092 csi[k][j][i].y * csi[k][j][i].y +
3093 csi[k][j][i].z * csi[k][j][i].z);
3094
3095 }
3096 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
3097 libm_Flux -= ucor[k][j][i].y;
3098 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3099 eta[k][j][i].y * eta[k][j][i].y +
3100 eta[k][j][i].z * eta[k][j][i].z);
3101 }
3102 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
3103 libm_Flux -= ucor[k][j][i].z;
3104 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3105 zet[k][j][i].y * zet[k][j][i].y +
3106 zet[k][j][i].z * zet[k][j][i].z);
3107 }
3108 }
3109
3110 }
3111 }
3112 }
3113
3114 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3115 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3116
3117 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
3118/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
3119 PetscPrintf(PETSC_COMM_WORLD, "IBMFlux22 REV %le %le\n", *ibm_Flux, *ibm_Area);
3120
3121 DMDAVecRestoreArray(da, user->lNvert, &nvert);
3122 DMDAVecRestoreArray(fda, user->lCsi, &csi);
3123 DMDAVecRestoreArray(fda, user->lEta, &eta);
3124 DMDAVecRestoreArray(fda, user->lZet, &zet);
3125 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
3126
3127 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3128 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3129 return 0;
3130}
SimCtx::NumberOfBodies
PetscInt NumberOfBodies
Definition variables.h:587
UserCtx::FluxIntpSum
PetscReal FluxIntpSum
Definition variables.h:685
Here is the caller graph for this function:

◆ VolumeFlux()

PetscErrorCode VolumeFlux ( UserCtx *  user,
PetscReal *  ibm_Flux,
PetscReal *  ibm_Area,
PetscInt  flg 
)

Calculates the net flux across the immersed boundary surface.

Parameters
userThe UserCtx for the grid level.
ibm_Flux(Output) The calculated net flux.
ibm_Area(Output) The total surface area of the IB.
flgA flag controlling the correction behavior.
Returns
PetscErrorCode 0 on success.

Definition at line 3133 of file poisson.c.

3134{
3135 // --- CONTEXT ACQUISITION BLOCK ---
3136 // Get the master simulation context from the UserCtx.
3137 SimCtx *simCtx = user->simCtx;
3138
3139 // Create local variables to mirror the legacy globals for minimal code changes.
3140 const PetscInt NumberOfBodies = simCtx->NumberOfBodies;
3141 const PetscInt channelz = simCtx->channelz;
3142 const PetscInt fish = simCtx->fish;
3143 // --- END CONTEXT ACQUISITION BLOCK ---
3144
3145 DM da = user->da, fda = user->fda;
3146
3147 DMDALocalInfo info = user->info;
3148
3149 PetscInt xs = info.xs, xe = info.xs + info.xm;
3150 PetscInt ys = info.ys, ye = info.ys + info.ym;
3151 PetscInt zs = info.zs, ze = info.zs + info.zm;
3152 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3153
3154 PetscInt i, j, k,ibi;
3155 PetscInt lxs, lys, lzs, lxe, lye, lze;
3156
3157 PetscInt gxs, gxe, gys, gye, gzs, gze;
3158
3159 gxs = info.gxs; gxe = gxs + info.gxm;
3160 gys = info.gys; gye = gys + info.gym;
3161 gzs = info.gzs; gze = gzs + info.gzm;
3162
3163 lxs = xs; lxe = xe;
3164 lys = ys; lye = ye;
3165 lzs = zs; lze = ze;
3166
3167 if (xs==0) lxs = xs+1;
3168 if (ys==0) lys = ys+1;
3169 if (zs==0) lzs = zs+1;
3170
3171 if (xe==mx) lxe = xe-1;
3172 if (ye==my) lye = ye-1;
3173 if (ze==mz) lze = ze-1;
3174
3175 PetscReal epsilon=1.e-8;
3176 PetscReal ***nvert, ibmval=1.9999;
3177
3178 struct Components {
3179 PetscReal x;
3180 PetscReal y;
3181 PetscReal z;
3182 }***ucor, ***lucor,***csi, ***eta, ***zet;
3183
3184
3185 PetscInt xend=mx-2 ,yend=my-2,zend=mz-2;
3186
3187 if (user->bctype[0]==7) xend=mx-1;
3188 if (user->bctype[2]==7) yend=my-1;
3189 if (user->bctype[4]==7) zend=mz-1;
3190
3191 DMDAVecGetArray(fda, user->Ucont, &ucor);
3192 DMDAVecGetArray(fda, user->lCsi, &csi);
3193 DMDAVecGetArray(fda, user->lEta, &eta);
3194 DMDAVecGetArray(fda, user->lZet, &zet);
3195 DMDAVecGetArray(da, user->lNvert, &nvert);
3196
3197 PetscReal libm_Flux, libm_area, libm_Flux_abs=0., ibm_Flux_abs;
3198 libm_Flux = 0;
3199 libm_area = 0;
3200
3201 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering VolumeFlux to enforce no-penetration condition.\n");
3202
3203 //Mohsen March 2017
3204 PetscReal *lIB_Flux, *lIB_area,*IB_Flux,*IB_Area;
3205 if (NumberOfBodies > 1) {
3206
3207 lIB_Flux=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3208 lIB_area=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3209 IB_Flux=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3210 IB_Area=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3211
3212
3213 for (ibi=0; ibi<NumberOfBodies; ibi++) {
3214 lIB_Flux[ibi]=0.0;
3215 lIB_area[ibi]=0.0;
3216 IB_Flux[ibi]=0.0;
3217 IB_Area[ibi]=0.0;
3218 }
3219 }
3220
3221
3222 //================================================================================
3223 // PASS 1: Calculate Uncorrected Flux and Area
3224 // This pass measures the total fluid "leakage" across the immersed boundary.
3225 //================================================================================
3226 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 1: Measuring uncorrected flux and area.\n");
3227
3228 for (k=lzs; k<lze; k++) {
3229 for (j=lys; j<lye; j++) {
3230 for (i=lxs; i<lxe; i++) {
3231 if (nvert[k][j][i] < 0.1) {
3232 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] < ibmval && i < xend) {
3233
3234 if (fabs(ucor[k][j][i].x)>epsilon) {
3235 libm_Flux += ucor[k][j][i].x;
3236 if (flg==3)
3237 libm_Flux_abs += fabs(ucor[k][j][i].x)/sqrt(csi[k][j][i].x * csi[k][j][i].x +
3238 csi[k][j][i].y * csi[k][j][i].y +
3239 csi[k][j][i].z * csi[k][j][i].z);
3240 else
3241 libm_Flux_abs += fabs(ucor[k][j][i].x);
3242
3243 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3244 csi[k][j][i].y * csi[k][j][i].y +
3245 csi[k][j][i].z * csi[k][j][i].z);
3246
3247 if (NumberOfBodies > 1) {
3248
3249 ibi=(int)((nvert[k][j][i+1]-1.0)*1001);
3250 lIB_Flux[ibi] += ucor[k][j][i].x;
3251 lIB_area[ibi] += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3252 csi[k][j][i].y * csi[k][j][i].y +
3253 csi[k][j][i].z * csi[k][j][i].z);
3254 }
3255 } else
3256 ucor[k][j][i].x=0.;
3257
3258 }
3259 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
3260
3261 if (fabs(ucor[k][j][i].y)>epsilon) {
3262 libm_Flux += ucor[k][j][i].y;
3263 if (flg==3)
3264 libm_Flux_abs += fabs(ucor[k][j][i].y)/sqrt(eta[k][j][i].x * eta[k][j][i].x +
3265 eta[k][j][i].y * eta[k][j][i].y +
3266 eta[k][j][i].z * eta[k][j][i].z);
3267 else
3268 libm_Flux_abs += fabs(ucor[k][j][i].y);
3269 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3270 eta[k][j][i].y * eta[k][j][i].y +
3271 eta[k][j][i].z * eta[k][j][i].z);
3272 if (NumberOfBodies > 1) {
3273
3274 ibi=(int)((nvert[k][j+1][i]-1.0)*1001);
3275
3276 lIB_Flux[ibi] += ucor[k][j][i].y;
3277 lIB_area[ibi] += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3278 eta[k][j][i].y * eta[k][j][i].y +
3279 eta[k][j][i].z * eta[k][j][i].z);
3280 }
3281 } else
3282 ucor[k][j][i].y=0.;
3283 }
3284 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
3285
3286 if (fabs(ucor[k][j][i].z)>epsilon) {
3287 libm_Flux += ucor[k][j][i].z;
3288 if (flg==3)
3289 libm_Flux_abs += fabs(ucor[k][j][i].z)/sqrt(zet[k][j][i].x * zet[k][j][i].x +
3290 zet[k][j][i].y * zet[k][j][i].y +
3291 zet[k][j][i].z * zet[k][j][i].z);
3292 else
3293 libm_Flux_abs += fabs(ucor[k][j][i].z);
3294 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3295 zet[k][j][i].y * zet[k][j][i].y +
3296 zet[k][j][i].z * zet[k][j][i].z);
3297
3298 if (NumberOfBodies > 1) {
3299
3300 ibi=(int)((nvert[k+1][j][i]-1.0)*1001);
3301 lIB_Flux[ibi] += ucor[k][j][i].z;
3302 lIB_area[ibi] += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3303 zet[k][j][i].y * zet[k][j][i].y +
3304 zet[k][j][i].z * zet[k][j][i].z);
3305 }
3306 }else
3307 ucor[k][j][i].z=0.;
3308 }
3309 }
3310
3311 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
3312
3313 if (nvert[k][j][i+1] < 0.1 && i < xend) {
3314 if (fabs(ucor[k][j][i].x)>epsilon) {
3315 libm_Flux -= ucor[k][j][i].x;
3316 if (flg==3)
3317 libm_Flux_abs += fabs(ucor[k][j][i].x)/sqrt(csi[k][j][i].x * csi[k][j][i].x +
3318 csi[k][j][i].y * csi[k][j][i].y +
3319 csi[k][j][i].z * csi[k][j][i].z);
3320 else
3321 libm_Flux_abs += fabs(ucor[k][j][i].x);
3322 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3323 csi[k][j][i].y * csi[k][j][i].y +
3324 csi[k][j][i].z * csi[k][j][i].z);
3325 if (NumberOfBodies > 1) {
3326
3327 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3328 lIB_Flux[ibi] -= ucor[k][j][i].x;
3329 lIB_area[ibi] += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3330 csi[k][j][i].y * csi[k][j][i].y +
3331 csi[k][j][i].z * csi[k][j][i].z);
3332 }
3333
3334 }else
3335 ucor[k][j][i].x=0.;
3336 }
3337 if (nvert[k][j+1][i] < 0.1 && j < yend) {
3338 if (fabs(ucor[k][j][i].y)>epsilon) {
3339 libm_Flux -= ucor[k][j][i].y;
3340 if (flg==3)
3341 libm_Flux_abs += fabs(ucor[k][j][i].y)/ sqrt(eta[k][j][i].x * eta[k][j][i].x +
3342 eta[k][j][i].y * eta[k][j][i].y +
3343 eta[k][j][i].z * eta[k][j][i].z);
3344 else
3345 libm_Flux_abs += fabs(ucor[k][j][i].y);
3346 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3347 eta[k][j][i].y * eta[k][j][i].y +
3348 eta[k][j][i].z * eta[k][j][i].z);
3349 if (NumberOfBodies > 1) {
3350
3351 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3352 lIB_Flux[ibi] -= ucor[k][j][i].y;
3353 lIB_area[ibi] += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3354 eta[k][j][i].y * eta[k][j][i].y +
3355 eta[k][j][i].z * eta[k][j][i].z);
3356 }
3357 }else
3358 ucor[k][j][i].y=0.;
3359 }
3360 if (nvert[k+1][j][i] < 0.1 && k < zend) {
3361 if (fabs(ucor[k][j][i].z)>epsilon) {
3362 libm_Flux -= ucor[k][j][i].z;
3363 if (flg==3)
3364 libm_Flux_abs += fabs(ucor[k][j][i].z)/sqrt(zet[k][j][i].x * zet[k][j][i].x +
3365 zet[k][j][i].y * zet[k][j][i].y +
3366 zet[k][j][i].z * zet[k][j][i].z);
3367 else
3368 libm_Flux_abs += fabs(ucor[k][j][i].z);
3369 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3370 zet[k][j][i].y * zet[k][j][i].y +
3371 zet[k][j][i].z * zet[k][j][i].z);
3372 if (NumberOfBodies > 1) {
3373
3374 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3375 lIB_Flux[ibi] -= ucor[k][j][i].z;
3376 lIB_area[ibi] += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3377 zet[k][j][i].y * zet[k][j][i].y +
3378 zet[k][j][i].z * zet[k][j][i].z);
3379 }
3380 }else
3381 ucor[k][j][i].z=0.;
3382 }
3383 }
3384
3385 }
3386 }
3387 }
3388
3389 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3390 MPI_Allreduce(&libm_Flux_abs, &ibm_Flux_abs,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3391 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3392
3393 if (NumberOfBodies > 1) {
3394 MPI_Allreduce(lIB_Flux,IB_Flux,NumberOfBodies,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
3395 MPI_Allreduce(lIB_area,IB_Area,NumberOfBodies,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
3396 }
3397
3398 PetscReal correction;
3399
3400 PetscReal *Correction;
3401 if (NumberOfBodies > 1) {
3402 Correction=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3403 for (ibi=0; ibi<NumberOfBodies; ibi++) Correction[ibi]=0.0;
3404 }
3405
3406 if (*ibm_Area > 1.e-15) {
3407 if (flg>1)
3408 correction = (*ibm_Flux + user->FluxIntpSum)/ ibm_Flux_abs;
3409 else if (flg)
3410 correction = (*ibm_Flux + user->FluxIntpSum) / *ibm_Area;
3411 else
3412 correction = *ibm_Flux / *ibm_Area;
3413 if (NumberOfBodies > 1)
3414 for (ibi=0; ibi<NumberOfBodies; ibi++) if (IB_Area[ibi]>1.e-15) Correction[ibi] = IB_Flux[ibi] / IB_Area[ibi];
3415 }
3416 else {
3417 correction = 0;
3418 }
3419 // --- Log the uncorrected results and calculated correction ---
3420 LOG_ALLOW(GLOBAL, LOG_INFO, "IBM Uncorrected Flux: %g, Area: %g, Correction: %g\n", *ibm_Flux, *ibm_Area, correction);
3421 if (NumberOfBodies>1){
3422 for (ibi=0; ibi<NumberOfBodies; ibi++) LOG_ALLOW(GLOBAL, LOG_INFO, " [Body %d] Uncorrected Flux: %g, Area: %g, Correction: %g\n", ibi, IB_Flux[ibi], IB_Area[ibi], Correction[ibi]);
3423 }
3424
3425 //================================================================================
3426 // PASS 2: Apply Correction to Velocity Field
3427 // This pass modifies the velocity at the boundary to enforce zero net flux.
3428 //================================================================================
3429 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 2: Applying velocity corrections at the boundary.\n");
3430
3431 for (k=lzs; k<lze; k++) {
3432 for (j=lys; j<lye; j++) {
3433 for (i=lxs; i<lxe; i++) {
3434 if (nvert[k][j][i] < 0.1) {
3435 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] <ibmval && i < xend) {
3436 if (fabs(ucor[k][j][i].x)>epsilon){
3437 if (flg==3)
3438 ucor[k][j][i].x -=correction*fabs(ucor[k][j][i].x)/
3439 sqrt(csi[k][j][i].x * csi[k][j][i].x +
3440 csi[k][j][i].y * csi[k][j][i].y +
3441 csi[k][j][i].z * csi[k][j][i].z);
3442 else if (flg==2)
3443 ucor[k][j][i].x -=correction*fabs(ucor[k][j][i].x);
3444 else if (NumberOfBodies > 1) {
3445 ibi=(int)((nvert[k][j][i+1]-1.0)*1001);
3446 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
3447 csi[k][j][i].y * csi[k][j][i].y +
3448 csi[k][j][i].z * csi[k][j][i].z) *
3449 Correction[ibi];
3450 }
3451 else
3452 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
3453 csi[k][j][i].y * csi[k][j][i].y +
3454 csi[k][j][i].z * csi[k][j][i].z) *
3455 correction;
3456 }
3457 }
3458 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
3459 if (fabs(ucor[k][j][i].y)>epsilon) {
3460 if (flg==3)
3461 ucor[k][j][i].y -=correction*fabs(ucor[k][j][i].y)/
3462 sqrt(eta[k][j][i].x * eta[k][j][i].x +
3463 eta[k][j][i].y * eta[k][j][i].y +
3464 eta[k][j][i].z * eta[k][j][i].z);
3465 else if (flg==2)
3466 ucor[k][j][i].y -=correction*fabs(ucor[k][j][i].y);
3467 else if (NumberOfBodies > 1) {
3468 ibi=(int)((nvert[k][j+1][i]-1.0)*1001);
3469 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
3470 eta[k][j][i].y * eta[k][j][i].y +
3471 eta[k][j][i].z * eta[k][j][i].z) *
3472 Correction[ibi];
3473 }
3474 else
3475 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
3476 eta[k][j][i].y * eta[k][j][i].y +
3477 eta[k][j][i].z * eta[k][j][i].z) *
3478 correction;
3479 }
3480 }
3481 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
3482 if (fabs(ucor[k][j][i].z)>epsilon) {
3483 if (flg==3)
3484 ucor[k][j][i].z -= correction*fabs(ucor[k][j][i].z)/
3485 sqrt(zet[k][j][i].x * zet[k][j][i].x +
3486 zet[k][j][i].y * zet[k][j][i].y +
3487 zet[k][j][i].z * zet[k][j][i].z);
3488 else if (flg==2)
3489 ucor[k][j][i].z -= correction*fabs(ucor[k][j][i].z);
3490 else if (NumberOfBodies > 1) {
3491 ibi=(int)((nvert[k+1][j][i]-1.0)*1001);
3492 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
3493 zet[k][j][i].y * zet[k][j][i].y +
3494 zet[k][j][i].z * zet[k][j][i].z) *
3495 Correction[ibi];
3496 }
3497 else
3498 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
3499 zet[k][j][i].y * zet[k][j][i].y +
3500 zet[k][j][i].z * zet[k][j][i].z) *
3501 correction;
3502 }
3503 }
3504 }
3505
3506 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
3507 if (nvert[k][j][i+1] < 0.1 && i < xend) {
3508 if (fabs(ucor[k][j][i].x)>epsilon) {
3509 if (flg==3)
3510 ucor[k][j][i].x += correction*fabs(ucor[k][j][i].x)/
3511 sqrt(csi[k][j][i].x * csi[k][j][i].x +
3512 csi[k][j][i].y * csi[k][j][i].y +
3513 csi[k][j][i].z * csi[k][j][i].z);
3514 else if (flg==2)
3515 ucor[k][j][i].x += correction*fabs(ucor[k][j][i].x);
3516 else if (NumberOfBodies > 1) {
3517 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3518 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3519 csi[k][j][i].y * csi[k][j][i].y +
3520 csi[k][j][i].z * csi[k][j][i].z) *
3521 Correction[ibi];
3522 }
3523 else
3524 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3525 csi[k][j][i].y * csi[k][j][i].y +
3526 csi[k][j][i].z * csi[k][j][i].z) *
3527 correction;
3528 }
3529 }
3530 if (nvert[k][j+1][i] < 0.1 && j < yend) {
3531 if (fabs(ucor[k][j][i].y)>epsilon) {
3532 if (flg==3)
3533 ucor[k][j][i].y +=correction*fabs(ucor[k][j][i].y)/
3534 sqrt(eta[k][j][i].x * eta[k][j][i].x +
3535 eta[k][j][i].y * eta[k][j][i].y +
3536 eta[k][j][i].z * eta[k][j][i].z);
3537 else if (flg==2)
3538 ucor[k][j][i].y +=correction*fabs(ucor[k][j][i].y);
3539 else if (NumberOfBodies > 1) {
3540 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3541 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3542 eta[k][j][i].y * eta[k][j][i].y +
3543 eta[k][j][i].z * eta[k][j][i].z) *
3544 Correction[ibi];
3545 }
3546 else
3547 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3548 eta[k][j][i].y * eta[k][j][i].y +
3549 eta[k][j][i].z * eta[k][j][i].z) *
3550 correction;
3551 }
3552 }
3553 if (nvert[k+1][j][i] < 0.1 && k < zend) {
3554 if (fabs(ucor[k][j][i].z)>epsilon) {
3555 if (flg==3)
3556 ucor[k][j][i].z += correction*fabs(ucor[k][j][i].z)/
3557 sqrt(zet[k][j][i].x * zet[k][j][i].x +
3558 zet[k][j][i].y * zet[k][j][i].y +
3559 zet[k][j][i].z * zet[k][j][i].z);
3560 else if (flg==2)
3561 ucor[k][j][i].z += correction*fabs(ucor[k][j][i].z);
3562 else if (NumberOfBodies > 1) {
3563 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3564 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3565 zet[k][j][i].y * zet[k][j][i].y +
3566 zet[k][j][i].z * zet[k][j][i].z) *
3567 Correction[ibi];
3568 }
3569 else
3570 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3571 zet[k][j][i].y * zet[k][j][i].y +
3572 zet[k][j][i].z * zet[k][j][i].z) *
3573 correction;
3574 }
3575 }
3576 }
3577
3578 }
3579 }
3580 }
3581
3582 //================================================================================
3583 // PASS 3: Verification
3584 // This optional pass recalculates the flux to confirm the correction was successful.
3585 //================================================================================
3586 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 3: Verifying corrected flux.\n");
3587
3588 libm_Flux = 0;
3589 libm_area = 0;
3590 for (k=lzs; k<lze; k++) {
3591 for (j=lys; j<lye; j++) {
3592 for (i=lxs; i<lxe; i++) {
3593 if (nvert[k][j][i] < 0.1) {
3594 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] < ibmval && i < xend) {
3595 libm_Flux += ucor[k][j][i].x;
3596 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3597 csi[k][j][i].y * csi[k][j][i].y +
3598 csi[k][j][i].z * csi[k][j][i].z);
3599
3600 }
3601 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
3602 libm_Flux += ucor[k][j][i].y;
3603 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3604 eta[k][j][i].y * eta[k][j][i].y +
3605 eta[k][j][i].z * eta[k][j][i].z);
3606 }
3607 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
3608 libm_Flux += ucor[k][j][i].z;
3609 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3610 zet[k][j][i].y * zet[k][j][i].y +
3611 zet[k][j][i].z * zet[k][j][i].z);
3612 }
3613 }
3614
3615 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
3616 if (nvert[k][j][i+1] < 0.1 && i < xend) {
3617 libm_Flux -= ucor[k][j][i].x;
3618 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3619 csi[k][j][i].y * csi[k][j][i].y +
3620 csi[k][j][i].z * csi[k][j][i].z);
3621
3622 }
3623 if (nvert[k][j+1][i] < 0.1 && j < yend) {
3624 libm_Flux -= ucor[k][j][i].y;
3625 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3626 eta[k][j][i].y * eta[k][j][i].y +
3627 eta[k][j][i].z * eta[k][j][i].z);
3628 }
3629 if (nvert[k+1][j][i] < 0.1 && k < zend) {
3630 libm_Flux -= ucor[k][j][i].z;
3631 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3632 zet[k][j][i].y * zet[k][j][i].y +
3633 zet[k][j][i].z * zet[k][j][i].z);
3634 }
3635 }
3636
3637 }
3638 }
3639 }
3640
3641 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3642 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3643
3644 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
3645/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
3646 LOG_ALLOW(GLOBAL, LOG_INFO, "IBM Corrected (Verified) Flux: %g, Area: %g\n", *ibm_Flux, *ibm_Area);
3647
3648
3649 if (user->bctype[0]==7 || user->bctype[1]==7){
3650 if (xe==mx){
3651 i=mx-2;
3652 for (k=lzs; k<lze; k++) {
3653 for (j=lys; j<lye; j++) {
3654 // if(j>0 && k>0 && j<user->JM && k<user->KM){
3655 if ((nvert[k][j][i]>ibmval && nvert[k][j][i+1]<0.1) || (nvert[k][j][i]<0.1 && nvert[k][j][i+1]>ibmval)) ucor[k][j][i].x=0.0;
3656
3657 // }
3658 }
3659 }
3660 }
3661 }
3662
3663 if (user->bctype[2]==7 || user->bctype[3]==7){
3664 if (ye==my){
3665 j=my-2;
3666 for (k=lzs; k<lze; k++) {
3667 for (i=lxs; i<lxe; i++) {
3668 // if(i>0 && k>0 && i<user->IM && k<user->KM){
3669 if ((nvert[k][j][i]>ibmval && nvert[k][j+1][i]<0.1) || (nvert[k][j][i]<0.1 && nvert[k][j+1][i]>ibmval)) ucor[k][j][i].y=0.0;
3670 // }
3671 }
3672 }
3673 }
3674 }
3675
3676 if (user->bctype[4]==7 || user->bctype[5]==7){
3677 if (ze==mz){
3678 k=mz-2;
3679 for (j=lys; j<lye; j++) {
3680 for (i=lxs; i<lxe; i++) {
3681 // if(i>0 && j>0 && i<user->IM && j<user->JM){
3682 if ((nvert[k][j][i]>ibmval && nvert[k+1][j][i]<0.1) || (nvert[k][j][i]<0.1 && nvert[k+1][j][i]>ibmval)) ucor[k][j][i].z=0.0;
3683 // }
3684 }
3685 }
3686 }
3687 }
3688
3689
3690 DMDAVecRestoreArray(da, user->lNvert, &nvert);
3691 DMDAVecRestoreArray(fda, user->lCsi, &csi);
3692 DMDAVecRestoreArray(fda, user->lEta, &eta);
3693 DMDAVecRestoreArray(fda, user->lZet, &zet);
3694 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
3695
3696 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3697 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3698
3699 /* periodci boundary condition update corrected flux */
3700 //Mohsen Dec 2015
3701 DMDAVecGetArray(fda, user->lUcont, &lucor);
3702 DMDAVecGetArray(fda, user->Ucont, &ucor);
3703
3704 if (user->bctype[0]==7 || user->bctype[1]==7){
3705 if (xs==0){
3706 i=xs;
3707 for (k=zs; k<ze; k++) {
3708 for (j=ys; j<ye; j++) {
3709 if(j>0 && k>0 && j<user->JM && k<user->KM){
3710 ucor[k][j][i].x=lucor[k][j][i-2].x;
3711 }
3712 }
3713 }
3714 }
3715 }
3716 if (user->bctype[2]==7 || user->bctype[3]==7){
3717 if (ys==0){
3718 j=ys;
3719 for (k=zs; k<ze; k++) {
3720 for (i=xs; i<xe; i++) {
3721 if(i>0 && k>0 && i<user->IM && k<user->KM){
3722 ucor[k][j][i].y=lucor[k][j-2][i].y;
3723 }
3724 }
3725 }
3726 }
3727 }
3728 if (user->bctype[4]==7 || user->bctype[5]==7){
3729 if (zs==0){
3730 k=zs;
3731 for (j=ys; j<ye; j++) {
3732 for (i=xs; i<xe; i++) {
3733 if(i>0 && j>0 && i<user->IM && j<user->JM){
3734 ucor[k][j][i].z=lucor[k-2][j][i].z;
3735 }
3736 }
3737 }
3738 }
3739 }
3740 DMDAVecRestoreArray(fda, user->lUcont, &lucor);
3741 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
3742
3743 /* DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); */
3744/* DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); */
3745 DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont);
3746 DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont);
3747
3748 if (NumberOfBodies > 1) {
3749 free(lIB_Flux);
3750 free(lIB_area);
3751 free(IB_Flux);
3752 free(IB_Area);
3753 free(Correction);
3754 }
3755
3756 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting VolumeFlux.\n");
3757
3758 return 0;
3759}
SimCtx::fish
PetscInt fish
Definition variables.h:572
Here is the caller graph for this function:

◆ FullyBlocked()

PetscErrorCode FullyBlocked ( UserCtx *  user)

Definition at line 3761 of file poisson.c.

3762{
3763 DM da = user->da;
3764 Vec nNvert;
3765 DMDALocalInfo info = user->info;
3766
3767 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3768
3769 PetscInt i, j, k;
3770
3771 PetscInt *KSKE = user->KSKE;
3772 PetscReal ***nvert;
3773 PetscBool *Blocked;
3774
3775 DMDACreateNaturalVector(da, &nNvert);
3776 DMDAGlobalToNaturalBegin(da, user->Nvert, INSERT_VALUES, nNvert);
3777 DMDAGlobalToNaturalEnd(da, user->Nvert, INSERT_VALUES, nNvert);
3778
3779 VecScatter ctx;
3780 Vec Zvert;
3781 VecScatterCreateToZero(nNvert, &ctx, &Zvert);
3782
3783 VecScatterBegin(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
3784 VecScatterEnd(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
3785
3786 VecScatterDestroy(&ctx);
3787 VecDestroy(&nNvert);
3788
3789 PetscInt rank;
3790 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
3791
3792 if (!rank) {
3793
3794 VecGetArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
3795 PetscMalloc(mx*my*sizeof(PetscBool), &Blocked);
3796 for (j=1; j<my-1; j++) {
3797 for (i=1; i<mx-1; i++) {
3798 Blocked[j*mx+i] = PETSC_FALSE;
3799 for (k=0; k<mz; k++) {
3800 if (nvert[k][j][i] > 0.1) {
3801 if (!Blocked[j*mx+i]) {
3802 KSKE[2*(j*mx+i)] = k;
3803 Blocked[j*mx+i] = PETSC_TRUE;
3804 }
3805 else {
3806 KSKE[2*(j*mx+i)] = PetscMin(KSKE[2*(j*mx+i)], k);
3807 }
3808 }
3809 }
3810 }
3811 }
3812
3813
3814 user->multinullspace = PETSC_TRUE;
3815 for (j=1; j<my-1; j++) {
3816 for (i=1; i<mx-1; i++) {
3817 if (!Blocked[j*mx+i]) {
3818 user->multinullspace = PETSC_FALSE;
3819 break;
3820 }
3821 }
3822 }
3823 PetscFree(Blocked);
3824 VecRestoreArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
3825 MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD);
3826 if (user->multinullspace) {
3827 MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD);
3828
3829 }
3830 }
3831 else {
3832 MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD);
3833 if (user->multinullspace) {
3834 MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD);
3835 }
3836 }
3837
3838
3839
3840 VecDestroy(&Zvert);
3841 return 0;
3842}
UserCtx::Nvert
Vec Nvert
Definition variables.h:688
Here is the caller graph for this function:

◆ MyNvertRestriction()

PetscErrorCode MyNvertRestriction ( UserCtx *  user_h,
UserCtx *  user_c 
)

Definition at line 3844 of file poisson.c.

3845{
3846 // DA da = user_c->da, fda = user_c->fda;
3847
3848
3849
3850 DMDALocalInfo info = user_c->info;
3851 PetscInt xs = info.xs, xe = info.xs + info.xm;
3852 PetscInt ys = info.ys, ye = info.ys + info.ym;
3853 PetscInt zs = info.zs, ze = info.zs + info.zm;
3854 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3855
3856 PetscInt i,j,k;
3857 PetscInt ih, jh, kh, ia, ja, ka;
3858 PetscInt lxs, lxe, lys, lye, lzs, lze;
3859
3860 PetscReal ***nvert, ***nvert_h;
3861
3862 DMDAVecGetArray(user_h->da, user_h->lNvert, &nvert_h);
3863 DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert);
3864
3865 lxs = xs; lxe = xe;
3866 lys = ys; lye = ye;
3867 lzs = zs; lze = ze;
3868
3869 if (xs==0) lxs = xs+1;
3870 if (ys==0) lys = ys+1;
3871 if (zs==0) lzs = zs+1;
3872
3873 if (xe==mx) lxe = xe-1;
3874 if (ye==my) lye = ye-1;
3875 if (ze==mz) lze = ze-1;
3876
3877 if ((user_c->isc)) ia = 0;
3878 else ia = 1;
3879
3880 if ((user_c->jsc)) ja = 0;
3881 else ja = 1;
3882
3883 if ((user_c->ksc)) ka = 0;
3884 else ka = 1;
3885
3886 VecSet(user_c->Nvert, 0.);
3887 if (user_c->thislevel > 0) {
3888 for (k=lzs; k<lze; k++) {
3889 for (j=lys; j<lye; j++) {
3890 for (i=lxs; i<lxe; i++) {
3891 GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
3892 if (nvert_h[kh ][jh ][ih ] *
3893 nvert_h[kh ][jh ][ih-ia] *
3894 nvert_h[kh ][jh-ja][ih ] *
3895 nvert_h[kh-ka][jh ][ih ] *
3896 nvert_h[kh ][jh-ja][ih-ia] *
3897 nvert_h[kh-ka][jh ][ih-ia] *
3898 nvert_h[kh-ka][jh-ja][ih ] *
3899 nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
3900 nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
3901 }
3902 }
3903 }
3904 }
3905 }
3906 else {
3907 for (k=lzs; k<lze; k++) {
3908 for (j=lys; j<lye; j++) {
3909 for (i=lxs; i<lxe; i++) {
3910 GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
3911 if (nvert_h[kh ][jh ][ih ] *
3912 nvert_h[kh ][jh ][ih-ia] *
3913 nvert_h[kh ][jh-ja][ih ] *
3914 nvert_h[kh-ka][jh ][ih ] *
3915 nvert_h[kh ][jh-ja][ih-ia] *
3916 nvert_h[kh-ka][jh ][ih-ia] *
3917 nvert_h[kh-ka][jh-ja][ih ] *
3918 nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
3919 nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
3920 }
3921 }
3922 }
3923 }
3924 }
3925 DMDAVecRestoreArray(user_h->da, user_h->lNvert, &nvert_h);
3926 DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert);
3927
3928 DMGlobalToLocalBegin(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3929 DMGlobalToLocalEnd(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3930 //Mohsen Dec 2015
3931 DMDAVecGetArray(user_c->da, user_c->lNvert, &nvert);
3932 DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert_h);
3933
3934 for (k=lzs; k<lze; k++) {
3935 for (j=lys; j<lye; j++) {
3936 for (i=lxs; i<lxe; i++) {
3937 if (nvert_h[k][j][i] < 0.1) {
3938 if (nvert[k][j][i+1] + nvert[k][j][i-1] > 1.1 &&
3939 nvert[k][j+1][i] + nvert[k][j-1][i] > 1.1 &&
3940 nvert[k+1][j][i] + nvert[k-1][j][i] > 1.1) {
3941 nvert_h[k][j][i] = 1.;
3942 }
3943 }
3944 }
3945 }
3946 }
3947
3948 DMDAVecRestoreArray(user_c->da, user_c->lNvert, &nvert);
3949 DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert_h);
3950 DMGlobalToLocalBegin(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3951 DMGlobalToLocalEnd(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3952 /* DMLocalToGlobalBegin(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
3953/* DMLocalToGlobalEnd(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
3954 return 0;
3955}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PoissonSolver_MG()

PetscErrorCode PoissonSolver_MG ( UserMG *  usermg)

Solves the pressure-Poisson equation using a geometric multigrid method.

This function orchestrates the entire multigrid V-cycle for the pressure correction equation. It assembles the Laplacian matrix on all grid levels, sets up the KSP solvers, smoothers, restriction/interpolation operators, and executes the solve.

Parameters
usermgThe UserMG context containing the entire multigrid hierarchy.
Returns
PetscErrorCode 0 on success.

Definition at line 3959 of file poisson.c.

3960{
3961 // --- CONTEXT ACQUISITION BLOCK ---
3962 // Get the master simulation context from the first block's UserCtx on the finest level.
3963 // This provides access to all former global variables.
3964 SimCtx *simCtx = usermg->mgctx[0].user[0].simCtx;
3965
3966 // Create local variables to mirror the legacy globals for minimal code changes.
3967 const PetscInt block_number = simCtx->block_number;
3968 const PetscInt immersed = simCtx->immersed;
3969 const PetscInt MHV = simCtx->MHV;
3970 const PetscInt LV = simCtx->LV;
3971 PetscMPIInt rank = simCtx->rank;
3972 // --- END CONTEXT ACQUISITION BLOCK ---
3973
3974 PetscErrorCode ierr;
3975 PetscInt l, bi;
3976 MGCtx *mgctx = usermg->mgctx;
3977 KSP mgksp, subksp;
3978 PC mgpc, subpc;
3979 UserCtx *user;
3980
3981 PetscFunctionBeginUser; // Moved to after variable declarations
3982 PROFILE_FUNCTION_BEGIN;
3983 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting Multigrid Poisson Solve...\n");
3984
3985 for (bi = 0; bi < block_number; bi++) {
3986
3987 // ====================================================================
3988 // SECTION: Immersed Boundary Specific Setup (Conditional)
3989 // ====================================================================
3990 if (immersed) {
3991 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Performing IBM pre-solve setup (Nvert restriction, etc.).\n", bi);
3992 for (l = usermg->mglevels - 1; l > 0; l--) {
3993 mgctx[l].user[bi].multinullspace = PETSC_FALSE;
3994 MyNvertRestriction(&mgctx[l].user[bi], &mgctx[l-1].user[bi]);
3995 }
3996 // Coarsest level check for disconnected domains due to IBM
3997 l = 0;
3998 user = mgctx[l].user;
3999 ierr = PetscMalloc1(user[bi].info.mx * user[bi].info.my * 2, &user[bi].KSKE); CHKERRQ(ierr);
4000 FullyBlocked(&user[bi]);
4001 }
4002
4003
4004 l = usermg->mglevels - 1;
4005 user = mgctx[l].user;
4006
4007 // --- 1. Compute RHS of the Poisson Equation ---
4008 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Computing Poisson RHS...\n", bi);
4009 ierr = VecDuplicate(user[bi].P, &user[bi].B); CHKERRQ(ierr);
4010
4011 PetscReal ibm_Flux, ibm_Area;
4012 PetscInt flg = immersed - 1;
4013
4014 // Calculate volume flux source terms (often from IBM)
4015 VolumeFlux(&user[bi], &ibm_Flux, &ibm_Area, flg);
4016 if (MHV || LV) {
4017 flg = ((MHV > 1 || LV) && bi == 0) ? 1 : 0;
4018 VolumeFlux_rev(&user[bi], &ibm_Flux, &ibm_Area, flg);
4019 }
4020 // Calculate the main divergence term
4021 PoissonRHS(&user[bi], user[bi].B);
4022
4023 // --- 2. Assemble LHS Matrix (Laplacian) on all MG levels ---
4024 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Assembling Poisson LHS on all levels...\n", bi);
4025 for (l = usermg->mglevels - 1; l >= 0; l--) {
4026 user = mgctx[l].user;
4027 LOG_ALLOW(GLOBAL,LOG_DEBUG," Calculating LHS for Level %d.\n",l);
4028 PoissonLHSNew(&user[bi]);
4029 }
4030
4031 // --- 3. Setup PETSc KSP and PCMG (Multigrid Preconditioner) ---
4032 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Configuring KSP and PCMG...\n", bi);
4033
4034 ierr = KSPCreate(PETSC_COMM_WORLD, &mgksp); CHKERRQ(ierr);
4035 ierr = KSPAppendOptionsPrefix(mgksp, "ps_"); CHKERRQ(ierr);
4036
4037 // =======================================================================
4038 DualMonitorCtx *monctx;
4039 char filen[128];
4040 PetscBool has_monitor_option;
4041
4042 // 1. Allocate the context and set it up.
4043 ierr = PetscNew(&monctx); CHKERRQ(ierr);
4044
4045 monctx->step = simCtx->step;
4046 monctx->block_id = bi;
4047 monctx->file_handle = NULL;
4048
4049 // Only rank 0 handles the file.
4050 if (!rank) {
4051 sprintf(filen, "%s/Poisson_Solver_Convergence_History_Block_%d.log", simCtx->log_dir,bi);
4052 // On the very first step of the entire simulation, TRUNCATE the file.
4053 if (simCtx->step == simCtx->StartStep) {
4054 monctx->file_handle = fopen(filen, "w");
4055 } else { // For all subsequent steps, APPEND to the file.
4056 monctx->file_handle = fopen(filen, "a");
4057 }
4058
4059 if (monctx->file_handle) {
4060 PetscFPrintf(PETSC_COMM_SELF, monctx->file_handle, "--- Convergence for Timestep %d, Block %d ---\n", (int)simCtx->step, bi);
4061 } else {
4062 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open KSP monitor log file: %s", filen);
4063 }
4064 }
4065
4066 ierr = PetscOptionsHasName(NULL, NULL, "-ps_ksp_pic_monitor_true_residual", &has_monitor_option); CHKERRQ(ierr);
4067 monctx->log_to_console = has_monitor_option;
4068
4069 ierr = KSPMonitorSet(mgksp, DualKSPMonitor, monctx, DualMonitorDestroy); CHKERRQ(ierr);
4070 // =======================================================================
4071
4072 ierr = KSPGetPC(mgksp, &mgpc); CHKERRQ(ierr);
4073 ierr = PCSetType(mgpc, PCMG); CHKERRQ(ierr);
4074
4075 ierr = PCMGSetLevels(mgpc, usermg->mglevels, PETSC_NULL); CHKERRQ(ierr);
4076 ierr = PCMGSetCycleType(mgpc, PC_MG_CYCLE_V); CHKERRQ(ierr);
4077 ierr = PCMGSetType(mgpc, PC_MG_MULTIPLICATIVE); CHKERRQ(ierr);
4078
4079 // --- 4. Define Restriction and Interpolation Operators for MG ---
4080 for (l = usermg->mglevels - 1; l > 0; l--) {
4081
4082 // Get stable pointers directly from the main mgctx array.
4083 // These pointers point to memory that will persist.
4084 UserCtx *fine_user_ctx = &mgctx[l].user[bi];
4085 UserCtx *coarse_user_ctx = &mgctx[l-1].user[bi];
4086
4087 // --- Configure the context pointers ---
4088 // The coarse UserCtx needs to know about the fine grid for restriction.
4089 coarse_user_ctx->da_f = &(fine_user_ctx->da);
4090 coarse_user_ctx->user_f = fine_user_ctx;
4091
4092 // The fine UserCtx needs to know about the coarse grid for interpolation.
4093 fine_user_ctx->da_c = &(coarse_user_ctx->da);
4094 fine_user_ctx->user_c = coarse_user_ctx;
4095 fine_user_ctx->lNvert_c = &(coarse_user_ctx->lNvert);
4096
4097 // --- Get matrix dimensions ---
4098 PetscInt m_c = (coarse_user_ctx->info.xm * coarse_user_ctx->info.ym * coarse_user_ctx->info.zm);
4099 PetscInt m_f = (fine_user_ctx->info.xm * fine_user_ctx->info.ym * fine_user_ctx->info.zm);
4100 PetscInt M_c = (coarse_user_ctx->info.mx * coarse_user_ctx->info.my * coarse_user_ctx->info.mz);
4101 PetscInt M_f = (fine_user_ctx->info.mx * fine_user_ctx->info.my * fine_user_ctx->info.mz);
4102
4103 LOG_ALLOW(GLOBAL,LOG_DEBUG,"level = %d; m_c = %d; m_f = %d; M_c = %d; M_f = %d.\n",l,m_c,m_f,M_c,M_f);
4104 // --- Create the MatShell objects ---
4105 // Pass the STABLE pointer coarse_user_ctx as the context for restriction.
4106 ierr = MatCreateShell(PETSC_COMM_WORLD, m_c, m_f, M_c, M_f, (void*)coarse_user_ctx, &fine_user_ctx->MR); CHKERRQ(ierr);
4107
4108 // Pass the STABLE pointer fine_user_ctx as the context for interpolation.
4109 ierr = MatCreateShell(PETSC_COMM_WORLD, m_f, m_c, M_f, M_c, (void*)fine_user_ctx, &fine_user_ctx->MP); CHKERRQ(ierr);
4110
4111 // --- Set the operations for the MatShells ---
4112 ierr = MatShellSetOperation(fine_user_ctx->MR, MATOP_MULT, (void(*)(void))RestrictResidual_SolidAware); CHKERRQ(ierr);
4113 ierr = MatShellSetOperation(fine_user_ctx->MP, MATOP_MULT, (void(*)(void))MyInterpolation); CHKERRQ(ierr);
4114
4115 // --- Register the operators with PCMG ---
4116 ierr = PCMGSetRestriction(mgpc, l, fine_user_ctx->MR); CHKERRQ(ierr);
4117 ierr = PCMGSetInterpolation(mgpc, l, fine_user_ctx->MP); CHKERRQ(ierr);
4118
4119 }
4120
4121 // --- 5. Configure Solvers on Each MG Level ---
4122 for (l = usermg->mglevels - 1; l >= 0; l--) {
4123 user = mgctx[l].user;
4124 if (l > 0) { // Smoother for fine levels
4125 ierr = PCMGGetSmoother(mgpc, l, &subksp); CHKERRQ(ierr);
4126 } else { // Direct or iterative solver for the coarsest level
4127 ierr = PCMGGetCoarseSolve(mgpc, &subksp); CHKERRQ(ierr);
4128 ierr = KSPSetTolerances(subksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, 40); CHKERRQ(ierr);
4129 }
4130
4131 ierr = KSPSetOperators(subksp, user[bi].A, user[bi].A); CHKERRQ(ierr);
4132 ierr = KSPSetFromOptions(subksp); CHKERRQ(ierr);
4133 ierr = KSPGetPC(subksp, &subpc); CHKERRQ(ierr);
4134 ierr = PCSetType(subpc, PCBJACOBI); CHKERRQ(ierr);
4135
4136 KSP *subsubksp;
4137 PC subsubpc;
4138 PetscInt nlocal;
4139
4140 // This logic is required for both the smoother and the coarse solve
4141 // since both use PCBJACOBI.
4142 ierr = KSPSetUp(subksp); CHKERRQ(ierr); // Set up KSP to allow access to sub-KSPs
4143 ierr = PCBJacobiGetSubKSP(subpc, &nlocal, NULL, &subsubksp); CHKERRQ(ierr);
4144
4145 for (PetscInt abi = 0; abi < nlocal; abi++) {
4146 ierr = KSPGetPC(subsubksp[abi], &subsubpc); CHKERRQ(ierr);
4147 // Add the critical shift amount
4148 ierr = PCFactorSetShiftAmount(subsubpc, 1.e-10); CHKERRQ(ierr);
4149 }
4150
4151 ierr = MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &user[bi].nullsp); CHKERRQ(ierr);
4152 ierr = MatNullSpaceSetFunction(user[bi].nullsp, PoissonNullSpaceFunction, &user[bi]); CHKERRQ(ierr);
4153 ierr = MatSetNullSpace(user[bi].A, user[bi].nullsp); CHKERRQ(ierr);
4154
4155 ierr = PCMGSetResidual(mgpc, l, PCMGResidualDefault, user[bi].A); CHKERRQ(ierr);
4156 ierr = KSPSetUp(subksp); CHKERRQ(ierr);
4157
4158 if (l < usermg->mglevels - 1) {
4159 ierr = MatCreateVecs(user[bi].A, &user[bi].R, PETSC_NULL); CHKERRQ(ierr);
4160 ierr = PCMGSetRhs(mgpc, l, user[bi].R); CHKERRQ(ierr);
4161 }
4162 }
4163
4164 // --- 6. Set Final KSP Operators and Solve ---
4165 l = usermg->mglevels - 1;
4166 user = mgctx[l].user;
4167
4168 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Setting KSP operators and solving...\n", bi);
4169 ierr = KSPSetOperators(mgksp, user[bi].A, user[bi].A); CHKERRQ(ierr);
4170 ierr = MatSetNullSpace(user[bi].A, user[bi].nullsp); CHKERRQ(ierr);
4171 ierr = KSPSetFromOptions(mgksp); CHKERRQ(ierr);
4172 ierr = KSPSetUp(mgksp); CHKERRQ(ierr);
4173 ierr = KSPSolve(mgksp, user[bi].B, user[bi].Phi); CHKERRQ(ierr);
4174
4175 // --- 7. Cleanup for this block ---
4176 for (l = usermg->mglevels - 1; l >= 0; l--) {
4177 user = mgctx[l].user;
4178 MatNullSpaceDestroy(&user[bi].nullsp);
4179 MatDestroy(&user[bi].A);
4180 user[bi].assignedA = PETSC_FALSE;
4181 if (l > 0) {
4182 MatDestroy(&user[bi].MR);
4183 MatDestroy(&user[bi].MP);
4184 } else if (l==0 && immersed) {
4185 PetscFree(user[bi].KSKE);
4186 }
4187 if (l < usermg->mglevels - 1) {
4188 VecDestroy(&user[bi].R);
4189 }
4190 }
4191
4192 KSPDestroy(&mgksp);
4193 VecDestroy(&mgctx[usermg->mglevels-1].user[bi].B);
4194
4195 } // End of loop over blocks
4196
4197 LOG_ALLOW(GLOBAL, LOG_INFO, "Multigrid Poisson Solve complete.\n");
4198 PROFILE_FUNCTION_END;
4199 PetscFunctionReturn(0);
4200}
DualMonitorDestroy
PetscErrorCode DualMonitorDestroy(void **ctx)
Destroys the DualMonitorCtx.
Definition logging.c:770
DualMonitorCtx::log_to_console
PetscBool log_to_console
Definition logging.h:59
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:46
DualMonitorCtx::step
PetscInt step
Definition logging.h:61
DualKSPMonitor
PetscErrorCode DualKSPMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
A custom KSP monitor that logs to a file and optionally to the console.
Definition logging.c:802
DualMonitorCtx::file_handle
FILE * file_handle
Definition logging.h:58
DualMonitorCtx::block_id
PetscInt block_id
Definition logging.h:62
DualMonitorCtx
Context for a dual-purpose KSP monitor.
Definition logging.h:57
MyNvertRestriction
PetscErrorCode MyNvertRestriction(UserCtx *user_h, UserCtx *user_c)
Definition poisson.c:3844
PoissonNullSpaceFunction
PetscErrorCode PoissonNullSpaceFunction(MatNullSpace nullsp, Vec X, void *ctx)
The callback function for PETSc's MatNullSpace object.
Definition poisson.c:1686
PoissonLHSNew
PetscErrorCode PoissonLHSNew(UserCtx *user)
Assembles the Left-Hand Side (LHS) matrix for the Pressure Poisson Equation.
Definition poisson.c:2207
VolumeFlux_rev
PetscErrorCode VolumeFlux_rev(UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
A specialized version of VolumeFlux, likely for reversed normals.
Definition poisson.c:2891
RestrictResidual_SolidAware
static PetscErrorCode RestrictResidual_SolidAware(Mat A, Vec X, Vec F)
Definition poisson.c:1986
VolumeFlux
PetscErrorCode VolumeFlux(UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
Calculates the net flux across the immersed boundary surface.
Definition poisson.c:3133
PoissonRHS
PetscErrorCode PoissonRHS(UserCtx *user, Vec B)
Computes the Right-Hand-Side (RHS) of the Poisson equation, which is the divergence of the intermedia...
Definition poisson.c:2810
MyInterpolation
PetscErrorCode MyInterpolation(Mat A, Vec X, Vec F)
The callback function for the multigrid interpolation operator (MatShell).
Definition poisson.c:1879
FullyBlocked
PetscErrorCode FullyBlocked(UserCtx *user)
Definition poisson.c:3761
SimCtx::MHV
PetscInt MHV
Definition variables.h:572
MGCtx::user
UserCtx * user
Definition variables.h:427
SimCtx::rank
PetscMPIInt rank
Definition variables.h:541
SimCtx::LV
PetscInt LV
Definition variables.h:572
SimCtx::block_number
PetscInt block_number
Definition variables.h:593
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:548
UserCtx::user_c
UserCtx * user_c
Definition variables.h:722
SimCtx::log_dir
char log_dir[PETSC_MAX_PATH_LEN]
Definition variables.h:562
UserMG::mglevels
PetscInt mglevels
Definition variables.h:434
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:437
SimCtx::immersed
PetscInt immersed
Definition variables.h:567
MGCtx
Context for Multigrid operations.
Definition variables.h:426
Here is the call graph for this function:
Here is the caller graph for this function: