PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
poisson.c
Go to the documentation of this file.
1// In src/poisson.c (or wherever PoissonSolver_MG is)
2#include "poisson.h" // The new header for this file
3#include "logging.h"
4
5static void Calculate_Covariant_metrics(double g[3][3], double G[3][3])
6{
7 /*
8 | csi.x csi.y csi.z |-1 | x.csi x.eta x.zet |
9 | eta.x eta.y eta.z | = | y.csi y.eta y.zet |
10 | zet.x zet.y zet.z | | z.csi z.eta z.zet |
11
12 */
13 const double a11=g[0][0], a12=g[0][1], a13=g[0][2];
14 const double a21=g[1][0], a22=g[1][1], a23=g[1][2];
15 const double a31=g[2][0], a32=g[2][1], a33=g[2][2];
16
17 double det= a11*(a33*a22-a32*a23)-a21*(a33*a12-a32*a13)+a31*(a23*a12-a22*a13);
18
19 G[0][0] = (a33*a22-a32*a23)/det, G[0][1] = - (a33*a12-a32*a13)/det, G[0][2] = (a23*a12-a22*a13)/det;
20 G[1][0] = -(a33*a21-a31*a23)/det, G[1][1] = (a33*a11-a31*a13)/det, G[1][2] = - (a23*a11-a21*a13)/det;
21 G[2][0] = (a32*a21-a31*a22)/det, G[2][1] = - (a32*a11-a31*a12)/det, G[2][2] = (a22*a11-a21*a12)/det;
22};
23
24static void Calculate_normal(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3])
25{
26 double g[3][3];
27 double G[3][3];
28
29 g[0][0]=csi.x, g[0][1]=csi.y, g[0][2]=csi.z;
30 g[1][0]=eta.x, g[1][1]=eta.y, g[1][2]=eta.z;
31 g[2][0]=zet.x, g[2][1]=zet.y, g[2][2]=zet.z;
32
33 Calculate_Covariant_metrics(g, G);
34 double xcsi=G[0][0], ycsi=G[1][0], zcsi=G[2][0];
35 double xeta=G[0][1], yeta=G[1][1], zeta=G[2][1];
36 double xzet=G[0][2], yzet=G[1][2], zzet=G[2][2];
37
38 double nx_i = xcsi, ny_i = ycsi, nz_i = zcsi;
39 double nx_j = xeta, ny_j = yeta, nz_j = zeta;
40 double nx_k = xzet, ny_k = yzet, nz_k = zzet;
41
42 double sum_i=sqrt(nx_i*nx_i+ny_i*ny_i+nz_i*nz_i);
43 double sum_j=sqrt(nx_j*nx_j+ny_j*ny_j+nz_j*nz_j);
44 double sum_k=sqrt(nx_k*nx_k+ny_k*ny_k+nz_k*nz_k);
45
46// *Ai = sqrt( g[0][0]*g[0][0] + g[0][1]*g[0][1] + g[0][2]*g[0][2] ); // area
47// *Aj = sqrt( g[1][0]*g[1][0] + g[1][1]*g[1][1] + g[1][2]*g[1][2] );
48// *Ak =sqrt( g[2][0]*g[2][0] + g[2][1]*g[2][1] + g[2][2]*g[2][2] );
49
50 nx_i /= sum_i, ny_i /= sum_i, nz_i /= sum_i;
51 nx_j /= sum_j, ny_j /= sum_j, nz_j /= sum_j;
52 nx_k /= sum_k, ny_k /= sum_k, nz_k /= sum_k;
53
54 ni[0] = nx_i, ni[1] = ny_i, ni[2] = nz_i;
55 nj[0] = nx_j, nj[1] = ny_j, nj[2] = nz_j;
56 nk[0] = nx_k, nk[1] = ny_k, nk[2] = nz_k;
57}
58
59static PetscErrorCode GhostNodeVelocity(UserCtx *user)
60{
61
62 // --- CONTEXT ACQUISITION BLOCK ---
63 // Get the master simulation context from the UserCtx.
64 SimCtx *simCtx = user->simCtx;
65
66 // Create local variables to mirror the legacy globals for minimal code changes.
67 PetscInt wallfunction = simCtx->wallfunction;
68 const PetscInt les = simCtx->les;
69 const PetscInt rans = simCtx->rans;
70 const PetscInt ti = simCtx->step;
71 const PetscReal U_bc = simCtx->U_bc;
72 // --- END CONTEXT ACQUISITION BLOCK ---
73
74 DM da = user->da, fda = user->fda;
75 DMDALocalInfo info = user->info;
76 PetscInt xs = info.xs, xe = info.xs + info.xm;
77 PetscInt ys = info.ys, ye = info.ys + info.ym;
78 PetscInt zs = info.zs, ze = info.zs + info.zm;
79 PetscInt mx = info.mx, my = info.my, mz = info.mz;
80
81 PetscInt i, j, k;
82 PetscInt lxs, lxe, lys, lye, lzs, lze;
83
84 Cmpnts ***ubcs, ***ucat,***lucat, ***csi, ***eta, ***zet,***cent;
85
86
87 PetscReal Un, nx,ny,nz,A;
88
89 lxs = xs; lxe = xe;
90 lys = ys; lye = ye;
91 lzs = zs; lze = ze;
92
93 PetscInt gxs, gxe, gys, gye, gzs, gze;
94
95 gxs = info.gxs; gxe = gxs + info.gxm;
96 gys = info.gys; gye = gys + info.gym;
97 gzs = info.gzs; gze = gzs + info.gzm;
98
99 if (xs==0) lxs = xs+1;
100 if (ys==0) lys = ys+1;
101 if (zs==0) lzs = zs+1;
102
103 if (xe==mx) lxe = xe-1;
104 if (ye==my) lye = ye-1;
105 if (ze==mz) lze = ze-1;
106
107 DMDAVecGetArray(fda, user->lCsi, &csi);
108 DMDAVecGetArray(fda, user->lEta, &eta);
109 DMDAVecGetArray(fda, user->lZet, &zet);
110
111
112/* ================================================================================== */
113/* WALL FUNCTION */
114/* ================================================================================== */
115 if (wallfunction) {
116 PetscReal ***ustar, ***aj, ***nvert;
117 DMDAVecGetArray(fda, user->Ucat, &ucat);
118 DMDAVecGetArray(da, user->lUstar, &ustar);
119 DMDAVecGetArray(da, user->lAj, &aj);
120 DMDAVecGetArray(da, user->lNvert, &nvert);
121
122 // wall function for boundary
123 //Dec 2015
124 for (k=lzs; k<lze; k++)
125 for (j=lys; j<lye; j++)
126 for (i=lxs; i<lxe; i++) {
127 if( nvert[k][j][i]<1.1 && ( (user->bctype[0]==-1 && i==1) || (user->bctype[1]==-1 && i==mx-2) ) ) {
128 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 );
129 double sb, sc;
130 double ni[3], nj[3], nk[3];
131 Cmpnts Uc, Ua;
132
133 Ua.x = Ua.y = Ua.z = 0;
134 sb = 0.5/aj[k][j][i]/area;
135
136 if(i==1) {
137 sc = 2*sb + 0.5/aj[k][j][i+1]/area;
138 Uc = ucat[k][j][i+1];
139 }
140 else {
141 sc = 2*sb + 0.5/aj[k][j][i-1]/area;
142 Uc = ucat[k][j][i-1];
143 }
144
145 Calculate_normal(csi[k][j][i], eta[k][j][i], zet[k][j][i], ni, nj, nk);
146 if(i==mx-2) ni[0]*=-1, ni[1]*=-1, ni[2]*=-1;
147
148 //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);
149 wall_function (user, sc, sb, Ua, Uc, &ucat[k][j][i], &ustar[k][j][i], ni[0], ni[1], ni[2]);
150 nvert[k][j][i]=1; /* set nvert to 1 to exclude from rhs */
151 // ucat[k][j][i] = lucat[k][j][i];
152 }
153
154 if( nvert[k][j][i]<1.1 && ( (user->bctype[2]==-1 && j==1) || (user->bctype[3]==-1 && j==my-2) ) ) {
155 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 );
156 double sb, sc;
157 double ni[3], nj[3], nk[3];
158 Cmpnts Uc, Ua;
159
160 Ua.x = Ua.y = Ua.z = 0;
161 sb = 0.5/aj[k][j][i]/area;
162
163 if(j==1) {
164 sc = 2*sb + 0.5/aj[k][j+1][i]/area;
165 Uc = ucat[k][j+1][i];
166 }
167 else {
168 sc = 2*sb + 0.5/aj[k][j-1][i]/area;
169 Uc = ucat[k][j-1][i];
170 }
171
172 Calculate_normal(csi[k][j][i], eta[k][j][i], zet[k][j][i], ni, nj, nk);
173 //if(j==my-2) nj[0]*=-1, nj[1]*=-1, nj[2]*=-1;
174
175 wall_function (user, sc, sb, Ua, Uc, &ucat[k][j][i], &ustar[k][j][i], nj[0], nj[1], nj[2]);
176 nvert[k][j][i]=1; /* set nvert to 1 to exclude from rhs */
177 // ucat[k][j][i] = lucat[k][j][i];
178 }
179
180 }
181 DMDAVecRestoreArray(da, user->lAj, &aj);
182 DMDAVecRestoreArray(da, user->lNvert, &nvert);
183 DMDAVecRestoreArray(da, user->lUstar, &ustar);
184 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
185
186 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
187 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
188 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
189/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
190
191 } // if wallfunction
192
193/* ================================================================================== */
194/* Driven Cavity */
195/* ================================================================================== */
196
197 DMDAVecGetArray(fda, user->Bcs.Ubcs, &ubcs);
198 DMDAVecGetArray(fda, user->Ucat, &ucat);
199
200 /* Designed for driven cavity problem (top(k=kmax) wall moving)
201 u_x = 1 at k==kmax */
202 if (user->bctype[5]==2) {
203 if (ze==mz) {
204 k = ze-1;
205 for (j=lys; j<lye; j++) {
206 for (i=lxs; i<lxe; i++) {
207 ubcs[k][j][i].x = 0.;// - ucat[k-1][j][i].x;
208 ubcs[k][j][i].y = 1.;//- ucat[k-1][j][i].y;
209 ubcs[k][j][i].z = 0.;//- ucat[k-1][j][i].z;
210 }
211 }
212 }
213 }
214 /* ==================================================================================== */
215 /* Cylinder O-grid */
216 /* ==================================================================================== */
217 if (user->bctype[3]==12) {
218 /* Designed to test O-grid for flow over a cylinder at jmax velocity is 1 (horizontal)
219 u_x = 1 at k==kmax */
220 if (ye==my) {
221 j = ye-1;
222 for (k=lzs; k<lze; k++) {
223 for (i=lxs; i<lxe; i++) {
224 ubcs[k][j][i].x = 0.;
225 ubcs[k][j][i].y = 0.;
226 ubcs[k][j][i].z = 1.;
227 }
228 }
229 }
230 }
231/* ==================================================================================== */
232/* Annulus */
233/* ==================================================================================== */
234/* designed to test periodic boundary condition for c grid j=0 rotates */
235
236 if (user->bctype[2]==11) {
237 DMDAVecGetArray(fda, user->Cent, &cent);
238 if (ys==0){
239 j=0;
240 for (k=lzs; k<lze; k++) {
241 for (i=lxs; i<lxe; i++) {
242 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);;
243 // ubcs[k][j][i].x=0.0;
244 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);
245 // 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);
246 ubcs[k][j][i].z = 0.0;
247 // 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);
248 }
249 }
250 }
251 DMDAVecRestoreArray(fda, user->Cent, &cent);
252 }
253
254/* ================================================================================== */
255/* SYMMETRY BC */
256/* ================================================================================== */
257
258 if (user->bctype[0]==3) {
259 if (xs==0) {
260 i= xs;
261
262 for (k=lzs; k<lze; k++) {
263 for (j=lys; j<lye; j++) {
264 A=sqrt(csi[k][j][i].z*csi[k][j][i].z +
265 csi[k][j][i].y*csi[k][j][i].y +
266 csi[k][j][i].x*csi[k][j][i].x);
267 nx=csi[k][j][i].x/A;
268 ny=csi[k][j][i].y/A;
269 nz=csi[k][j][i].z/A;
270 Un=ucat[k][j][i+1].x*nx+ucat[k][j][i+1].y*ny+ucat[k][j][i+1].z*nz;
271 ubcs[k][j][i].x = ucat[k][j][i+1].x-Un*nx;//-V_frame.x;
272 ubcs[k][j][i].y = ucat[k][j][i+1].y-Un*ny;
273 ubcs[k][j][i].z = ucat[k][j][i+1].z-Un*nz;
274 }
275 }
276 }
277 }
278
279 if (user->bctype[1]==3) {
280 if (xe==mx) {
281 i= xe-1;
282
283 for (k=lzs; k<lze; k++) {
284 for (j=lys; j<lye; j++) {
285 A=sqrt(csi[k][j][i-1].z*csi[k][j][i-1].z +
286 csi[k][j][i-1].y*csi[k][j][i-1].y +
287 csi[k][j][i-1].x*csi[k][j][i-1].x);
288 nx=csi[k][j][i-1].x/A;
289 ny=csi[k][j][i-1].y/A;
290 nz=csi[k][j][i-1].z/A;
291 Un=ucat[k][j][i-1].x*nx+ucat[k][j][i-1].y*ny+ucat[k][j][i-1].z*nz;
292 ubcs[k][j][i].x = ucat[k][j][i-1].x-Un*nx;
293 ubcs[k][j][i].y = ucat[k][j][i-1].y-Un*ny;
294 ubcs[k][j][i].z = ucat[k][j][i-1].z-Un*nz;
295 }
296 }
297 }
298 }
299
300 if (user->bctype[2]==3) {
301 if (ys==0) {
302 j= ys;
303
304 for (k=lzs; k<lze; k++) {
305 for (i=lxs; i<lxe; i++) {
306 A=sqrt(eta[k][j][i].z*eta[k][j][i].z +
307 eta[k][j][i].y*eta[k][j][i].y +
308 eta[k][j][i].x*eta[k][j][i].x);
309 nx=eta[k][j][i].x/A;
310 ny=eta[k][j][i].y/A;
311 nz=eta[k][j][i].z/A;
312 Un=ucat[k][j+1][i].x*nx+ucat[k][j+1][i].y*ny+ucat[k][j+1][i].z*nz;
313 ubcs[k][j][i].x = ucat[k][j+1][i].x-Un*nx;
314 ubcs[k][j][i].y = ucat[k][j+1][i].y-Un*ny;
315 ubcs[k][j][i].z = ucat[k][j+1][i].z-Un*nz;
316 }
317 }
318 }
319 }
320
321 if (user->bctype[3]==3) {
322 if (ye==my) {
323 j=ye-1;
324
325 for (k=lzs; k<lze; k++) {
326 for (i=lxs; i<lxe; i++) {
327 A=sqrt(eta[k][j-1][i].z*eta[k][j-1][i].z +
328 eta[k][j-1][i].y*eta[k][j-1][i].y +
329 eta[k][j-1][i].x*eta[k][j-1][i].x);
330 nx=eta[k][j-1][i].x/A;
331 ny=eta[k][j-1][i].y/A;
332 nz=eta[k][j-1][i].z/A;
333 Un=ucat[k][j-1][i].x*nx+ucat[k][j-1][i].y*ny+ucat[k][j-1][i].z*nz;
334 ubcs[k][j][i].x = ucat[k][j-1][i].x-Un*nx;
335 ubcs[k][j][i].y = ucat[k][j-1][i].y-Un*ny;
336 ubcs[k][j][i].z = ucat[k][j-1][i].z-Un*nz;
337 }
338 }
339 }
340 }
341
342
343 if (user->bctype[4]==3) {
344 if (zs==0) {
345 k = 0;
346 for (j=lys; j<lye; j++) {
347 for (i=lxs; i<lxe; i++) {
348 A=sqrt(zet[k][j][i].z*zet[k][j][i].z +
349 zet[k][j][i].y*zet[k][j][i].y +
350 zet[k][j][i].x*zet[k][j][i].x);
351 nx=zet[k][j][i].x/A;
352 ny=zet[k][j][i].y/A;
353 nz=zet[k][j][i].z/A;
354 Un=ucat[k][j][i].x*nx+ucat[k][j][i].y*ny+ucat[k][j][i].z*nz;
355 ubcs[k][j][i].x = ucat[k][j][i].x-Un*nx;
356 ubcs[k][j][i].y = ucat[k][j][i].y-Un*ny;
357 ubcs[k][j][i].z = ucat[k][j][i].z-Un*nz;
358 }
359 }
360 }
361 }
362
363 if (user->bctype[5]==3) {
364 if (ze==mz) {
365 k = ze-1;
366 for (j=lys; j<lye; j++) {
367 for (i=lxs; i<lxe; i++) {
368 A=sqrt(zet[k-1][j][i].z*zet[k-1][j][i].z +
369 zet[k-1][j][i].y*zet[k-1][j][i].y +
370 zet[k-1][j][i].x*zet[k-1][j][i].x);
371 nx=zet[k-1][j][i].x/A;
372 ny=zet[k-1][j][i].y/A;
373 nz=zet[k-1][j][i].z/A;
374 Un=ucat[k-1][j][i].x*nx+ucat[k-1][j][i].y*ny+ucat[k-1][j][i].z*nz;
375 ubcs[k][j][i].x = ucat[k-1][j][i].x-Un*nx;
376 ubcs[k][j][i].y = ucat[k-1][j][i].y-Un*ny;
377 ubcs[k][j][i].z = ucat[k-1][j][i].z-Un*nz;
378 }
379 }
380 }
381 }
382/* ==================================================================================== */
383 /* Rheology */
384 /* ==================================================================================== */
385
386 // PetscPrintf(PETSC_COMM_WORLD, "moving plate velocity for rheology setup is %le \n",U_bc);
387
388 if (user->bctype[2]==13){
389 if (ys==0){
390 j=0;
391 for (k=lzs; k<lze; k++) {
392 for (i=lxs; i<lxe; i++) {
393 ubcs[k][j][i].x = 0.;
394 //ubcs[k][j][i].x = -U_bc;
395 ubcs[k][j][i].y = 0.;
396 ubcs[k][j][i].z =-U_bc;
397 //ubcs[k][j][i].z =0.0;
398 }
399 }
400 }
401 }
402 if (user->bctype[3]==13){
403 if (ye==my){
404 j=ye-1;
405 for (k=lzs; k<lze; k++) {
406 for (i=lxs; i<lxe; i++) {
407 ubcs[k][j][i].x = 0.;
408 //ubcs[k][j][i].x = U_bc;
409 ubcs[k][j][i].y = 0.;
410 ubcs[k][j][i].z = U_bc;
411 //ubcs[k][j][i].z =0.0;
412 }
413 }
414 }
415 }
416 if (user->bctype[4]==13){
417 if (zs==0){
418 k=0;
419 for (j=lys; j<lye; j++) {
420 for (i=lxs; i<lxe; i++) {
421 ubcs[k][j][i].x =-U_bc;
422 ubcs[k][j][i].y = 0.;
423 ubcs[k][j][i].z = 0.;
424 }
425 }
426 }
427 }
428 if (user->bctype[5]==13){
429 if (ze==mz){
430 k=ze-1;
431 for (j=lys; j<lye; j++) {
432 for (i=lxs; i<lxe; i++) {
433 ubcs[k][j][i].x = U_bc;
434 ubcs[k][j][i].y = 0.;
435 ubcs[k][j][i].z = 0.;
436 }
437 }
438 }
439 }
440/* ================================================================================== */
441 // boundary conditions on ghost nodes
442/* ================================================================================== */
443//Mohsen Aug 2012
444 if (xs==0 && user->bctype[0]!=7) {
445 i = xs;
446 for (k=zs; k<ze; k++) {
447 for (j=ys; j<ye; j++) {
448 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j][i+1].x;
449 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j][i+1].y;
450 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j][i+1].z;
451 }
452 }
453 }
454
455 if (xe==mx && user->bctype[0]!=7) {
456 i = xe-1;
457 for (k=zs; k<ze; k++) {
458 for (j=ys; j<ye; j++) {
459 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j][i-1].x;
460 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j][i-1].y;
461 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j][i-1].z;
462 }
463 }
464 }
465
466 if (ys==0 && user->bctype[2]!=7) {
467 j = ys;
468 for (k=zs; k<ze; k++) {
469 for (i=xs; i<xe; i++) {
470 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j+1][i].x;
471 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j+1][i].y;
472 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j+1][i].z;
473 }
474 }
475 }
476
477 if (ye==my && user->bctype[2]!=7) {
478 j = ye-1;
479 for (k=zs; k<ze; k++) {
480 for (i=xs; i<xe; i++) {
481 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k][j-1][i].x;
482 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k][j-1][i].y;
483 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k][j-1][i].z;
484 }
485 }
486 }
487
488 if (zs==0 && user->bctype[4]!=7) {
489 k = zs;
490 for (j=ys; j<ye; j++) {
491 for (i=xs; i<xe; i++) {
492 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k+1][j][i].x;
493 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k+1][j][i].y;
494 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k+1][j][i].z;
495 }
496 }
497 }
498
499 if (ze==mz && user->bctype[4]!=7) {
500 k = ze-1;
501 for (j=ys; j<ye; j++) {
502 for (i=xs; i<xe; i++) {
503 ucat[k][j][i].x = 2 * ubcs[k][j][i].x - ucat[k-1][j][i].x;
504 ucat[k][j][i].y = 2 * ubcs[k][j][i].y - ucat[k-1][j][i].y;
505 ucat[k][j][i].z = 2 * ubcs[k][j][i].z - ucat[k-1][j][i].z;
506 }
507 }
508 }
509 DMDAVecRestoreArray(fda, user->Bcs.Ubcs, &ubcs);
510 DMDAVecRestoreArray(fda, user->Ucat,&ucat);
511/* ================================================================================== */
512/* Periodic BC Mohsen Aug 2012
513/* ================================================================================== */
514
515 if (user->bctype[0]==7 || user->bctype[2]==7 || user->bctype[4]==7 ){
516
517 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
518 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
519
520 DMDAVecGetArray(fda, user->lUcat, &lucat);
521 DMDAVecGetArray(fda, user->Ucat, &ucat);
522
523 if (user->bctype[0]==7 || user->bctype[1]==7){
524 if (xs==0){
525 i=xs;
526 for (k=zs; k<ze; k++) {
527 for (j=ys; j<ye; j++) {
528 if(k>0 && k<user->KM && j>0 && j<user->JM){
529 ucat[k][j][i]=lucat[k][j][i-2];
530 }
531 }
532 }
533 }
534 }
535 if (user->bctype[2]==7 || user->bctype[3]==7){
536 if (ys==0){
537 j=ys;
538 for (k=zs; k<ze; k++) {
539 for (i=xs; i<xe; i++) {
540 if(k>0 && k<user->KM && i>0 && i<user->IM){
541 ucat[k][j][i]=lucat[k][j-2][i];
542 }
543 }
544 }
545 }
546 }
547 if (user->bctype[4]==7 || user->bctype[5]==7){
548 if (zs==0){
549 k=zs;
550 for (j=ys; j<ye; j++) {
551 for (i=xs; i<xe; i++) {
552 if(j>0 && j<user->JM && i>0 && i<user->IM){
553 ucat[k][j][i]=lucat[k-2][j][i];
554 }
555 }
556 }
557 }
558 }
559 if (user->bctype[0]==7 || user->bctype[1]==7){
560 if (xe==mx){
561 i=mx-1;
562 for (k=zs; k<ze; k++) {
563 for (j=ys; j<ye; j++) {
564 if(k>0 && k<user->KM && j>0 && j<user->JM){
565 ucat[k][j][i]=lucat[k][j][i+2];
566 }
567 }
568 }
569 }
570 }
571 if (user->bctype[2]==7 || user->bctype[3]==7){
572 if (ye==my){
573 j=my-1;
574 for (k=zs; k<ze; k++) {
575 for (i=xs; i<xe; i++) {
576 if(k>0 && k<user->KM && i>0 && i<user->IM){
577 ucat[k][j][i]=lucat[k][j+2][i];
578 }
579 }
580 }
581 }
582 }
583 if (user->bctype[4]==7 || user->bctype[5]==7){
584 if (ze==mz){
585 k=mz-1;
586 for (j=ys; j<ye; j++) {
587 for (i=xs; i<xe; i++) {
588 if(j>0 && j<user->JM && i>0 && i<user->IM){
589 ucat[k][j][i].x=lucat[k+2][j][i].x;
590 }
591 }
592 }
593 }
594 }
595
596 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
597 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
598
599 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
600 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
601 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
602/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
603 }
604
605 // velocity on the corner point
606 DMDAVecGetArray(fda, user->Ucat, &ucat);
607
608 if (zs==0 ) {
609 k=0;
610 if (xs==0) {
611 i=0;
612 for (j=ys; j<ye; j++) {
613 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j][i+1].x);
614 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j][i+1].y);
615 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j][i+1].z);
616 }
617 }
618 if (xe == mx) {
619 i=mx-1;
620 for (j=ys; j<ye; j++) {
621 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j][i-1].x);
622 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j][i-1].y);
623 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j][i-1].z);
624 }
625 }
626
627 if (ys==0) {
628 j=0;
629 for (i=xs; i<xe; i++) {
630 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j+1][i].x);
631 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j+1][i].y);
632 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j+1][i].z);
633 }
634 }
635
636 if (ye==my) {
637 j=my-1;
638 for (i=xs; i<xe; i++) {
639 ucat[k][j][i].x = 0.5*(ucat[k+1][j][i].x+ucat[k][j-1][i].x);
640 ucat[k][j][i].y = 0.5*(ucat[k+1][j][i].y+ucat[k][j-1][i].y);
641 ucat[k][j][i].z = 0.5*(ucat[k+1][j][i].z+ucat[k][j-1][i].z);
642 }
643 }
644 }
645
646 if (ze==mz) {
647 k=mz-1;
648 if (xs==0) {
649 i=0;
650 for (j=ys; j<ye; j++) {
651 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j][i+1].x);
652 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j][i+1].y);
653 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j][i+1].z);
654 }
655 }
656 if (xe == mx) {
657 i=mx-1;
658 for (j=ys; j<ye; j++) {
659 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j][i-1].x);
660 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j][i-1].y);
661 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j][i-1].z);
662 }
663 }
664 if (ys==0) {
665 j=0;
666 for (i=xs; i<xe; i++) {
667 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j+1][i].x);
668 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j+1][i].y);
669 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j+1][i].z);
670 }
671 }
672
673 if (ye==my) {
674 j=my-1;
675 for (i=xs; i<xe; i++) {
676 ucat[k][j][i].x = 0.5*(ucat[k-1][j][i].x+ucat[k][j-1][i].x);
677 ucat[k][j][i].y = 0.5*(ucat[k-1][j][i].y+ucat[k][j-1][i].y);
678 ucat[k][j][i].z = 0.5*(ucat[k-1][j][i].z+ucat[k][j-1][i].z);
679 }
680 }
681 }
682
683 if (ys==0) {
684 j=0;
685 if (xs==0) {
686 i=0;
687 for (k=zs; k<ze; k++) {
688 ucat[k][j][i].x = 0.5*(ucat[k][j+1][i].x+ucat[k][j][i+1].x);
689 ucat[k][j][i].y = 0.5*(ucat[k][j+1][i].y+ucat[k][j][i+1].y);
690 ucat[k][j][i].z = 0.5*(ucat[k][j+1][i].z+ucat[k][j][i+1].z);
691 }
692 }
693
694 if (xe==mx) {
695 i=mx-1;
696 for (k=zs; k<ze; k++) {
697 ucat[k][j][i].x = 0.5*(ucat[k][j+1][i].x+ucat[k][j][i-1].x);
698 ucat[k][j][i].y = 0.5*(ucat[k][j+1][i].y+ucat[k][j][i-1].y);
699 ucat[k][j][i].z = 0.5*(ucat[k][j+1][i].z+ucat[k][j][i-1].z);
700 }
701 }
702 }
703
704 if (ye==my) {
705 j=my-1;
706 if (xs==0) {
707 i=0;
708 for (k=zs; k<ze; k++) {
709 ucat[k][j][i].x = 0.5*(ucat[k][j-1][i].x+ucat[k][j][i+1].x);
710 ucat[k][j][i].y = 0.5*(ucat[k][j-1][i].y+ucat[k][j][i+1].y);
711 ucat[k][j][i].z = 0.5*(ucat[k][j-1][i].z+ucat[k][j][i+1].z);
712 }
713 }
714
715 if (xe==mx) {
716 i=mx-1;
717 for (k=zs; k<ze; k++) {
718 ucat[k][j][i].x = 0.5*(ucat[k][j-1][i].x+ucat[k][j][i-1].x);
719 ucat[k][j][i].y = 0.5*(ucat[k][j-1][i].y+ucat[k][j][i-1].y);
720 ucat[k][j][i].z = 0.5*(ucat[k][j-1][i].z+ucat[k][j][i-1].z);
721 }
722 }
723 }
724
725 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
726
727 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
728 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
729
730 if (user->bctype[0]==7 || user->bctype[2]==7 || user->bctype[4]==7 ){
731 //i-direction
732 DMDAVecGetArray(fda, user->lUcat, &lucat);
733 DMDAVecGetArray(fda, user->Ucat, &ucat);
734
735 if (user->bctype[0]==7){
736 if (xs==0){
737 i=xs;
738 for (k=zs; k<ze; k++) {
739 for (j=ys; j<ye; j++) {
740 ucat[k][j][i]=lucat[k][j][i-2];
741 }
742 }
743 }
744 }
745 if (user->bctype[1]==7){
746 if (xe==mx){
747 i=xe-1;
748 for (k=zs; k<ze; k++) {
749 for (j=ys; j<ye; j++) {
750 ucat[k][j][i]=lucat[k][j][i+2];
751 }
752 }
753 }
754 }
755 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
756 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
757
758 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
759/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
760
761 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
762 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
763
764 //j-direction
765
766 DMDAVecGetArray(fda, user->lUcat, &lucat);
767 DMDAVecGetArray(fda, user->Ucat, &ucat);
768
769 if (user->bctype[2]==7){
770 if (ys==0){
771 j=ys;
772 for (k=zs; k<ze; k++) {
773 for (i=xs; i<xe; i++) {
774 ucat[k][j][i]=lucat[k][j-2][i];
775 }
776 }
777 }
778 }
779
780 if (user->bctype[3]==7){
781 if (ye==my){
782 j=my-1;
783 for (k=zs; k<ze; k++) {
784 for (i=xs; i<xe; i++) {
785 ucat[k][j][i]=lucat[k][j+2][i];
786 }
787 }
788 }
789 }
790
791 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
792 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
793
794 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
795/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
796
797 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
798 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
799
800 //k-direction
801
802 DMDAVecGetArray(fda, user->lUcat, &lucat);
803 DMDAVecGetArray(fda, user->Ucat, &ucat);
804
805 if (user->bctype[4]==7){
806 if (zs==0){
807 k=zs;
808 for (j=ys; j<ye; j++) {
809 for (i=xs; i<xe; i++) {
810 ucat[k][j][i]=lucat[k-2][j][i];
811 }
812 }
813 }
814 }
815 if (user->bctype[5]==7){
816 if (ze==mz){
817 k=mz-1;
818 for (j=ys; j<ye; j++) {
819 for (i=xs; i<xe; i++) {
820 ucat[k][j][i].x=lucat[k+2][j][i].x;
821 }
822 }
823 }
824 }
825
826 DMDAVecRestoreArray(fda, user->lUcat, &lucat);
827 DMDAVecRestoreArray(fda, user->Ucat, &ucat);
828
829 /* DMLocalToGlobalBegin(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
830/* DMLocalToGlobalEnd(fda, user->lUcat, INSERT_VALUES, user->Ucat); */
831
832 DMGlobalToLocalBegin(fda, user->Ucat, INSERT_VALUES, user->lUcat);
833 DMGlobalToLocalEnd(fda, user->Ucat, INSERT_VALUES, user->lUcat);
834 }
835
836 DMDAVecRestoreArray(fda, user->lCsi, &csi);
837 DMDAVecRestoreArray(fda, user->lEta, &eta);
838 DMDAVecRestoreArray(fda, user->lZet, &zet);
839
840 return(0);
841
842}
843
844#define GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user) \
845 if ((user->isc)) { \
846 ic = i; \
847 ia = 0; \
848 } \
849 else { \
850 ic = (i+1) / 2; \
851 ia = (i - 2 * (ic)) == 0 ? 1 : -1; \
852 if (i==1 || i==mx-2) ia = 0; \
853 }\
854 if ((user->jsc)) { \
855 jc = j; \
856 ja = 0; \
857 } \
858 else { \
859 jc = (j+1) / 2; \
860 ja = (j - 2 * (jc)) == 0 ? 1 : -1; \
861 if (j==1 || j==my-2) ja = 0; \
862 } \
863 if ((user->ksc)) { \
864 kc = k; \
865 ka = 0; \
866 } \
867 else { \
868 kc = (k+1) / 2; \
869 ka = (k - 2 * (kc)) == 0 ? 1 : -1; \
870 if (k==1 || k==mz-2) ka = 0; \
871 } \
872 if (ka==-1 && nvert_c[kc-1][jc][ic] > 0.1) ka = 0; \
873 else if (ka==1 && nvert_c[kc+1][jc][ic] > 0.1) ka = 0; \
874 if (ja==-1 && nvert_c[kc][jc-1][ic] > 0.1) ja = 0; \
875 else if (ja==1 && nvert_c[kc][jc+1][ic] > 0.1) ja = 0; \
876 if (ia==-1 && nvert_c[kc][jc][ic-1] > 0.1) ia = 0; \
877 else if (ia==1 && nvert_c[kc][jc][ic+1] > 0.1) ia = 0;
878
879static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
880
881{
882 PetscInt nidx;
883 DMDALocalInfo info = user->info;
884
885 PetscInt mx = info.mx, my = info.my;
886
887 AO ao;
888 DMDAGetAO(user->da, &ao);
889 nidx=i+j*mx+k*mx*my;
890
891 AOApplicationToPetsc(ao,1,&nidx);
892
893 return (nidx);
894}
895
896static PetscErrorCode GridRestriction(PetscInt i, PetscInt j, PetscInt k,
897 PetscInt *ih, PetscInt *jh, PetscInt *kh,
898 UserCtx *user)
899{
900 if ((user->isc)) {
901 *ih = i;
902 }
903 else {
904 *ih = 2 * i;
905 }
906
907 if ((user->jsc)) {
908 *jh = j;
909 }
910 else {
911 *jh = 2 * j;
912 }
913
914 if ((user->ksc)) {
915 *kh = k;
916 }
917 else {
918 *kh = 2 * k;
919 }
920
921 return 0;
922}
923
924
925//Cell Center
926#define CP 0
927// Face Centers
928#define EP 1
929#define WP 2
930#define NP 3
931#define SP 4
932#define TP 5
933#define BP 6
934
935// Edge Centers
936#define NE 7
937#define SE 8
938#define NW 9
939#define SW 10
940#define TN 11
941#define BN 12
942#define TS 13
943#define BS 14
944#define TE 15
945#define BE 16
946#define TW 17
947#define BW 18
948
949#undef __FUNCT__
950#define __FUNCT__ "Projection"
951/**
952 * @brief Performs the projection step to enforce an incompressible velocity field.
953 *
954 * @details
955 * This function executes the final "correction" stage of a fractional-step (projection)
956 * method for solving the incompressible Navier-Stokes equations. After solving the
957 * Pressure Poisson Equation to get a pressure correction `Phi`, this function uses
958 * `Phi` to correct the intermediate velocity field, making it divergence-free.
959 *
960 * The main steps are:
961 * 1. **Calculate Pressure Gradient**: At each velocity point (on the cell faces), it
962 * computes the gradient of the pressure correction field (`∇Φ`). This is done using
963 * finite differences on a generalized curvilinear grid.
964 *
965 * 2. **Correct Velocity Field**: It updates the contravariant velocity components `Ucont`
966 * by subtracting the scaled pressure gradient. The update rule is:
967 * `U_new = U_intermediate - Δt * G * ∇Φ`, where `G` is the metric tensor `g_ij`
968 * that transforms the gradient into contravariant components.
969 *
970 * 3. **Boundary-Aware Gradients**: The gradient calculation is highly sensitive to
971 * boundaries. The code uses complex, dynamic stencils that switch from centered
972 * differences to one-sided differences if a standard stencil would use a point
973 * inside an immersed solid (`nvert > 0.1`). This preserves accuracy at the
974 * fluid-solid interface.
975 *
976 * 4. **Periodic Boundary Updates**: Includes dedicated loops to correctly calculate
977 * the pressure gradient at the domain edges for periodic boundary conditions.
978 *
979 * 5. **Finalization**: After correcting `Ucont`, it calls helper functions to convert
980 * the velocity back to Cartesian coordinates (`Contra2Cart`) and update all
981 * ghost cell values.
982 *
983 * @param[in, out] user Pointer to the UserCtx struct, containing simulation context,
984 * grid metrics, and the velocity/pressure fields.
985 *
986 * @return PetscErrorCode 0 on success, or a non-zero error code on failure.
987 */
988PetscErrorCode Projection(UserCtx *user)
989{
990 PetscFunctionBeginUser;
991 PROFILE_FUNCTION_BEGIN;
992 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering Projection step to correct velocity field.\n");
993
994 //================================================================================
995 // Section 1: Initialization and Data Acquisition
996 //================================================================================
997
998 // --- Get simulation and grid context ---
999 SimCtx *simCtx = user->simCtx;
1000 DM da = user->da, fda = user->fda;
1001 DMDALocalInfo info = user->info;
1002
1003 // --- Grid dimensions ---
1004 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1005 PetscInt xs = info.xs, xe = info.xs + info.xm;
1006 PetscInt ys = info.ys, ye = info.ys + info.ym;
1007 PetscInt zs = info.zs, ze = info.zs + info.zm;
1008
1009 // --- Loop bounds (excluding outer ghost layers) ---
1010 PetscInt lxs = (xs == 0) ? xs + 1 : xs;
1011 PetscInt lxe = (xe == mx) ? xe - 1 : xe;
1012 PetscInt lys = (ys == 0) ? ys + 1 : ys;
1013 PetscInt lye = (ye == my) ? ye - 1 : ye;
1014 PetscInt lzs = (zs == 0) ? zs + 1 : zs;
1015 PetscInt lze = (ze == mz) ? ze - 1 : ze;
1016
1017 // --- Get direct pointer access to grid metric and field data ---
1018 Cmpnts ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1019 PetscReal ***iaj, ***jaj, ***kaj, ***p, ***nvert;
1020 Cmpnts ***ucont;
1021 DMDAVecGetArray(fda, user->lICsi, &icsi); DMDAVecGetArray(fda, user->lIEta, &ieta); DMDAVecGetArray(fda, user->lIZet, &izet);
1022 DMDAVecGetArray(fda, user->lJCsi, &jcsi); DMDAVecGetArray(fda, user->lJEta, &jeta); DMDAVecGetArray(fda, user->lJZet, &jzet);
1023 DMDAVecGetArray(fda, user->lKCsi, &kcsi); DMDAVecGetArray(fda, user->lKEta, &keta); DMDAVecGetArray(fda, user->lKZet, &kzet);
1024 DMDAVecGetArray(da, user->lIAj, &iaj); DMDAVecGetArray(da, user->lJAj, &jaj); DMDAVecGetArray(da, user->lKAj, &kaj);
1025 DMDAVecGetArray(da, user->lNvert, &nvert);
1026 DMDAVecGetArray(da, user->lPhi, &p); // Note: using lPhi, which is the pressure correction
1027 //DMDAVecGetArray(da,user->lP,&p);
1028 DMDAVecGetArray(fda, user->Ucont, &ucont);
1029
1030 // --- Constants for clarity ---
1031 const PetscReal IBM_FLUID_THRESHOLD = 0.1;
1032 const PetscReal scale = simCtx->dt * simCtx->st / COEF_TIME_ACCURACY;
1033
1034 LOG_ALLOW(GLOBAL,LOG_DEBUG," Starting velocity correction: Scale = %le .\n",scale);
1035
1036 //================================================================================
1037 // Section 2: Correct Velocity Components
1038 //================================================================================
1039 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Calculating pressure gradients and correcting velocity components.\n");
1040
1041 // --- Main loop over interior domain points ---
1042 for (PetscInt k = lzs; k < lze; k++) {
1043 for (PetscInt j = lys; j < lye; j++) {
1044 for (PetscInt i = lxs; i < lxe; i++) {
1045
1046 // --- Correct U_contravariant (x-component of velocity) ---
1047 PetscInt i_end = (user->bctype[0] == 7) ? mx - 1 : mx - 2;
1048 if (i < i_end) {
1049
1050 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k][j][i + 1] > IBM_FLUID_THRESHOLD)) {
1051 // Compute pressure derivatives (dp/d_csi, dp/d_eta, dp/d_zet) at the i-face
1052
1053 PetscReal dpdc = p[k][j][i + 1] - p[k][j][i];
1054 PetscReal dpde = 0.0, dpdz = 0.0;
1055
1056 // Boundary-aware stencil for dp/d_eta
1057 if ((j==my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1058 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1059 dpde = (p[k][j][i] + p[k][j][i+1] -
1060 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1061 }
1062 }
1063
1064 else if ((j==my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1065 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; }
1066 }
1067
1068 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1069 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; }
1070 }
1071
1072 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1073 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; }
1074 }
1075
1076 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; }
1077
1078 // Boundary-aware stencil for dp/d_zet
1079 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1080 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; }
1081 }
1082
1083 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1084 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; }
1085 }
1086
1087 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1088 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; }
1089 }
1090
1091 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1092 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; }
1093 }
1094
1095 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; }
1096
1097 // Apply the correction: U_new = U_old - dt * (g11*dpdc + g12*dpde + g13*dpdz)
1098
1099
1100
1101 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
1102 + icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
1103 dpde * (ieta[k][j][i].x * icsi[k][j][i].x + ieta[k][j][i].y * icsi[k][j][i].y
1104 + ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
1105 dpdz * (izet[k][j][i].x * icsi[k][j][i].x + izet[k][j][i].y * icsi[k][j][i].y
1106 + izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i]);
1107
1108 PetscReal correction = grad_p_x*scale;
1109 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Csi Direction: %le.\n",correction);
1110 ucont[k][j][i].x -= correction;
1111
1112 }
1113 }
1114
1115 // --- Correct V_contravariant (y-component of velocity) ---
1116 PetscInt j_end = (user->bctype[2] == 7) ? my - 1 : my - 2;
1117 if (j < j_end) {
1118 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k][j + 1][i] > IBM_FLUID_THRESHOLD)) {
1119 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1120 dpde = p[k][j + 1][i] - p[k][j][i];
1121
1122 // Boundary-aware stencil for dp/d_csi
1123 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1124 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; }
1125 } else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1126 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; }
1127 } else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1128 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; }
1129 } else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1130 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; }
1131 } 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; }
1132
1133 // Boundary-aware stencil for dp/d_zet
1134 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1135 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; }
1136 } else if ((k == mz-2 || k==1 ) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1137 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; }
1138 } else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1139 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; }
1140 } else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1141 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; }
1142 } 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; }
1143
1144 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] +
1145 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] +
1146 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]);
1147
1148 PetscReal correction = grad_p_y*scale;
1149 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Eta Direction: %le.\n",correction);
1150 ucont[k][j][i].y -= correction;
1151 }
1152 }
1153
1154 // --- Correct W_contravariant (z-component of velocity) ---
1155 PetscInt k_end = (user->bctype[4] == 7) ? mz - 1 : mz - 2;
1156 if (k < k_end) {
1157 if (!(nvert[k][j][i] > IBM_FLUID_THRESHOLD || nvert[k + 1][j][i] > IBM_FLUID_THRESHOLD)) {
1158 PetscReal dpdc = 0.0, dpde = 0.0, dpdz = 0.0;
1159 dpdz = p[k + 1][j][i] - p[k][j][i];
1160
1161 // Boundary-aware stencil for dp/d_csi
1162 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1163 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; }
1164 } else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1165 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; }
1166 } else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1167 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; }
1168 } else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1169 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; }
1170 } 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; }
1171
1172 // Boundary-aware stencil for dp/d_eta
1173 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1174 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; }
1175 } else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1176 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; }
1177 } else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1178 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; }
1179 } else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1180 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; }
1181 } 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; }
1182
1183 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] +
1184 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] +
1185 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]);
1186
1187 // ========================= DEBUG PRINT =========================
1188 LOG_LOOP_ALLOW_EXACT(GLOBAL, LOG_DEBUG, k, 5,
1189 "[k=%d, j=%d, i=%d] ---- Neighbor Pressures ----\n"
1190 " Central Z-Neighbors: p[k+1][j][i] = %g | p[k][j][i] = %g\n"
1191 " 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"
1192 " 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",
1193 k, j, i,
1194 p[k + 1][j][i], p[k][j][i],
1195 p[k][j - 1][i], p[k + 1][j - 1][i], p[k][j + 1][i], p[k + 1][j + 1][i],
1196 p[k][j][i - 1], p[k + 1][j][i - 1], p[k][j][i + 1], p[k + 1][j][i + 1]);
1197 // ======================= END DEBUG PRINT =======================
1198
1199 LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," dpdc: %le | dpde: %le | dpdz: %le.\n",dpdc,dpde,dpdz);
1200 PetscReal correction = grad_p_z*scale;
1201 //LOG_LOOP_ALLOW_EXACT(GLOBAL,LOG_DEBUG,k,5," Flux correction in Zet Direction: %le.\n",correction);
1202 ucont[k][j][i].z -= correction;
1203 }
1204 }
1205 }
1206 }
1207 }
1208
1209 // --- Explicit correction for periodic boundaries (if necessary) ---
1210 // The main loop handles the interior, but this handles the first physical layer at periodic boundaries.
1211 // Note: This logic is largely duplicated from the main loop and could be merged, but is preserved for fidelity.
1212 if (user->bctype[0] == 7 && xs == 0) {
1213 for (PetscInt k=lzs; k<lze; k++) {
1214 for (PetscInt j=lys; j<lye; j++) {
1215 PetscInt i=xs;
1216
1217 PetscReal dpdc = p[k][j][i+1] - p[k][j][i];
1218
1219 PetscReal dpde = 0.;
1220 PetscReal dpdz = 0.;
1221
1222 if ((j==my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1223 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1224 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1225 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1226 }
1227 }
1228 else if ((j==my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1229 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1230 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1231 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1232 }
1233 }
1234 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1235 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1236 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1237 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1238 }
1239 }
1240 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1241 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1242 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1243 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1244 }
1245 }
1246 else {
1247 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1248 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1249 }
1250
1251 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1252 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1253 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1254 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1255 }
1256 }
1257 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1258 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1259 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1260 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1261 }
1262 }
1263 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1264 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1265 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1266 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1267 }
1268 }
1269 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1270 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1271 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1272 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1273 }
1274 }
1275 else {
1276 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1277 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1278 }
1279
1280
1281
1282 if (!(nvert[k][j][i] + nvert[k][j][i+1])) {
1283 ucont[k][j][i].x -=
1284 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1285 icsi[k][j][i].y * icsi[k][j][i].y +
1286 icsi[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
1287 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1288 ieta[k][j][i].y * icsi[k][j][i].y +
1289 ieta[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i] +
1290 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1291 izet[k][j][i].y * icsi[k][j][i].y +
1292 izet[k][j][i].z * icsi[k][j][i].z) * iaj[k][j][i])
1293 * scale;
1294
1295 }
1296 }
1297 }
1298 }
1299 if (user->bctype[2] == 7 && ys == 0) {
1300
1301 for (PetscInt k=lzs; k<lze; k++) {
1302 for (PetscInt i=lxs; i<lxe; i++) {
1303 PetscInt j=ys;
1304
1305 PetscReal dpdc = 0.;
1306 PetscReal dpdz = 0.;
1307 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1308 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1309 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1310 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1311 }
1312 }
1313 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1314 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1315 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1316 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1317 }
1318 }
1319 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1320 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1321 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1322 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1323 }
1324 }
1325 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1326 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1327 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1328 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1329 }
1330 }
1331 else {
1332 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1333 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1334 }
1335
1336 PetscReal dpde = p[k][j+1][i] - p[k][j][i];
1337
1338 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1339 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1340 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1341 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1342 }
1343 }
1344 else if ((k == mz-2 || k==1 ) && user->bctype[0]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1345 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1346 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1347 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1348 }
1349 }
1350 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1351 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1352 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1353 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1354 }
1355 }
1356 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1357 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1358 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1359 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1360 }
1361 }
1362 else {
1363 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1364 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1365 }
1366
1367 if (!(nvert[k][j][i] + nvert[k][j+1][i])) {
1368 ucont[k][j][i].y -=
1369 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1370 jcsi[k][j][i].y * jeta[k][j][i].y +
1371 jcsi[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
1372 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1373 jeta[k][j][i].y * jeta[k][j][i].y +
1374 jeta[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i] +
1375 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1376 jzet[k][j][i].y * jeta[k][j][i].y +
1377 jzet[k][j][i].z * jeta[k][j][i].z) * jaj[k][j][i])
1378 * scale;
1379 }
1380 }
1381 }
1382 }
1383
1384 if (user->bctype[4] == 7 && zs == 0) {
1385 for (PetscInt j=lys; j<lye; j++) {
1386 for (PetscInt i=lxs; i<lxe; i++) {
1387
1388 PetscInt k=zs;
1389 PetscReal dpdc = 0.;
1390 PetscReal dpde = 0.;
1391
1392 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1393 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1394 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1395 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1396 }
1397 }
1398 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1399 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1400 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1401 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1402 }
1403 }
1404 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1405 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1406 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1407 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1408 }
1409 }
1410 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1411 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1412 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1413 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1414 }
1415 }
1416 else {
1417 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1418 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1419 }
1420
1421 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1422 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1423 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1424 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1425 }
1426 }
1427 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1428 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1429 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1430 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1431 }
1432 }
1433 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1434 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1435 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1436 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1437 }
1438 }
1439 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1440 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1441 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1442 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1443 }
1444 }
1445 else {
1446 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1447 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1448 }
1449
1450 PetscReal dpdz = p[k+1][j][i] - p[k][j][i];
1451
1452 if (!(nvert[k][j][i] + nvert[k+1][j][i])) {
1453
1454 ucont[k][j][i].z -=
1455 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1456 kcsi[k][j][i].y * kzet[k][j][i].y +
1457 kcsi[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
1458 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1459 keta[k][j][i].y * kzet[k][j][i].y +
1460 keta[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i] +
1461 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1462 kzet[k][j][i].y * kzet[k][j][i].y +
1463 kzet[k][j][i].z * kzet[k][j][i].z) * kaj[k][j][i])
1464 * scale;
1465
1466 }
1467 }
1468 }
1469 }
1470
1471 // --- Optional: Enforce specific mass flow rate for channel flow simulations ---
1472 if (simCtx->channelz) {
1473 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Applying channel flow mass flux correction.\n");
1474 // DMDAVecGetArray(fda, user->lCsi, &csi); DMDAVecGetArray(fda, user->lEta, &eta); DMDAVecGetArray(fda, user->lZet, &zet);
1475 //DMDAVecGetArray(da, user->lAj, &aj);
1476 // This entire block was commented out in the original code. It calculates the current
1477 // total flux in the z-direction, compares it to a target flux (user->simCtx->FluxInSum),
1478 // and computes a uniform velocity correction `ratio` to enforce the target.
1479 // The implementation for applying the `ratio` was also commented out.
1480 // Preserving this as-is.
1481 }
1482
1483 //================================================================================
1484 // Section 3: Finalization and Cleanup
1485 //================================================================================
1486
1487 // --- Restore access to all PETSc vector arrays ---
1488 DMDAVecRestoreArray(fda, user->Ucont, &ucont);
1489 // DMDAVecRestoreArray(fda, user->lCsi, &csi); DMDAVecRestoreArray(fda, user->lEta, &eta); DMDAVecRestoreArray(fda, user->lZet, &zet);
1490 //DMDAVecRestoreArray(da, user->lAj, &aj);
1491 DMDAVecRestoreArray(fda, user->lICsi, &icsi); DMDAVecRestoreArray(fda, user->lIEta, &ieta); DMDAVecRestoreArray(fda, user->lIZet, &izet);
1492 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi); DMDAVecRestoreArray(fda, user->lJEta, &jeta); DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1493 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi); DMDAVecRestoreArray(fda, user->lKEta, &keta); DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1494 DMDAVecRestoreArray(da, user->lIAj, &iaj); DMDAVecRestoreArray(da, user->lJAj, &jaj); DMDAVecRestoreArray(da, user->lKAj, &kaj);
1495 DMDAVecRestoreArray(da, user->lP, &p);
1496 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1497
1498 // --- Update ghost cells for the newly corrected velocity field ---
1499 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Updating ghost cells for corrected velocity.\n");
1500 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
1501 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
1502
1503 // --- Convert velocity to Cartesian and update ghost nodes ---
1504 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Converting velocity to Cartesian and finalizing ghost nodes.\n");
1505 Contra2Cart(user);
1506 GhostNodeVelocity(user);
1507
1508 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting Projection step.\n");
1509 PROFILE_FUNCTION_END;
1510 PetscFunctionReturn(0);
1511}
1512
1513#undef __FUNCT__
1514#define __FUNCT__ "UpdatePressure"
1515/**
1516 * @brief Updates the pressure field and handles periodic boundary conditions.
1517 *
1518 * @details
1519 * This function is a core part of the projection method in a CFD solver. It is
1520 * called after the Pressure Poisson Equation has been solved to obtain a pressure
1521 * correction field, `Phi`.
1522 *
1523 * The function performs two main operations:
1524 *
1525 * 1. **Core Pressure Update**: It updates the main pressure field `P` by adding the
1526 * pressure correction `Phi` at every grid point in the fluid domain. The
1527 * operation is `P_new = P_old + Phi`.
1528 *
1529 * 2. **Periodic Boundary Synchronization**: If any of the domain boundaries are
1530 * periodic (`bctype == 7`), this function manually updates the pressure values
1531 * in the ghost cells. It copies values from the physical cells on one side of
1532 * the domain to the corresponding ghost cells on the opposite side. This ensures
1533 * that subsequent calculations involving pressure gradients are correct across
1534 * periodic boundaries. The refactored version consolidates the original code's
1535 * redundant updates into a single, efficient pass.
1536 *
1537 * @param[in, out] user Pointer to the UserCtx struct, containing simulation context and
1538 * the pressure vectors `P` and `Phi`.
1539 *
1540 * @return PetscErrorCode 0 on success, or a non-zero error code on failure.
1541 */
1542PetscErrorCode UpdatePressure(UserCtx *user)
1543{
1544 PetscFunctionBeginUser;
1545 PROFILE_FUNCTION_BEGIN;
1546 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering UpdatePressure.\n");
1547
1548 //================================================================================
1549 // Section 1: Initialization and Data Acquisition
1550 //================================================================================
1551 DM da = user->da;
1552 DMDALocalInfo info = user->info;
1553
1554 // Local grid extents for the main update loop
1555 PetscInt xs = info.xs, xe = info.xs + info.xm;
1556 PetscInt ys = info.ys, ye = info.ys + info.ym;
1557 PetscInt zs = info.zs, ze = info.zs + info.zm;
1558 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1559
1560 PetscInt i,j,k;
1561
1562 // --- Get direct pointer access to PETSc vector data for performance ---
1563 PetscReal ***p, ***phi;
1564 DMDAVecGetArray(da, user->P, &p);
1565 DMDAVecGetArray(da, user->Phi, &phi);
1566
1567 //================================================================================
1568 // Section 2: Core Pressure Update
1569 //================================================================================
1570 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Performing core pressure update (P_new = P_old + Phi).\n");
1571 for (PetscInt k = zs; k < ze; k++) {
1572 for (PetscInt j = ys; j < ye; j++) {
1573 for (PetscInt i = xs; i < xe; i++) {
1574 // This is the fundamental pressure update in a projection method.
1575 p[k][j][i] += phi[k][j][i];
1576 }
1577 }
1578 }
1579
1580 // Restore arrays now that the core computation is done.
1581 DMDAVecRestoreArray(da, user->Phi, &phi);
1582 DMDAVecRestoreArray(da, user->P, &p);
1583
1584
1585 //================================================================================
1586 // Section 3: Handle Periodic Boundary Condition Synchronization
1587 //================================================================================
1588 // This block is executed only if at least one boundary is periodic.
1589 // The original code contained many redundant Get/Restore and update calls.
1590 // This refactored version performs the same effective logic but in a single,
1591 // efficient pass, which is numerically equivalent and much cleaner.
1592 if (user->bctype[0] == 7 || user->bctype[1] == 7 ||
1593 user->bctype[2] == 7 || user->bctype[3] == 7 ||
1594 user->bctype[4] == 7 || user->bctype[5] == 7)
1595 {
1596 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Synchronizing ghost cells for periodic boundaries.\n");
1597
1598 // First, update the local vectors (including ghost regions) with the new global data.
1599 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
1600 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
1601 DMGlobalToLocalBegin(da, user->Phi, INSERT_VALUES, user->lPhi);
1602 DMGlobalToLocalEnd(da, user->Phi, INSERT_VALUES, user->lPhi);
1603
1604 // Get pointers to all necessary local and global arrays ONCE.
1605 PetscReal ***lp, ***lphi;
1606 DMDAVecGetArray(da, user->lP, &lp);
1607 DMDAVecGetArray(da, user->lPhi, &lphi);
1608 DMDAVecGetArray(da, user->P, &p);
1609 DMDAVecGetArray(da, user->Phi, &phi);
1610
1611 // --- X-Direction Periodic Update ---
1612 if (user->bctype[0] == 7 || user->bctype[1] == 7) {
1613 // Update left boundary physical cells from right boundary ghost cells
1614 if (xs == 0) {
1615 PetscInt i = 0;
1616 for (k = zs; k < ze; k++) {
1617 for (j = ys; j < ye; j++) {
1618 // Note: Accessing lp[...][i-2] reads from the ghost cell region.
1619 p[k][j][i] = lp[k][j][i - 2];
1620 phi[k][j][i] = lphi[k][j][i - 2];
1621 }
1622 }
1623 }
1624 // Update right boundary physical cells from left boundary ghost cells
1625 if (xe == mx) {
1626 PetscInt i = mx - 1;
1627 for (k = zs; k < ze; k++) {
1628 for (j = ys; j < ye; j++) {
1629 p[k][j][i] = lp[k][j][i + 2];
1630 phi[k][j][i] = lphi[k][j][i + 2];
1631 }
1632 }
1633 }
1634 }
1635
1636 // --- Y-Direction Periodic Update ---
1637 if (user->bctype[2] == 7 || user->bctype[3] == 7) {
1638 // Update bottom boundary
1639 if (ys == 0) {
1640 PetscInt j = 0;
1641 for (k = zs; k < ze; k++) {
1642 for (i = xs; i < xe; i++) {
1643 p[k][j][i] = lp[k][j - 2][i];
1644 phi[k][j][i] = lphi[k][j - 2][i];
1645 }
1646 }
1647 }
1648 // Update top boundary
1649 if (ye == my) {
1650 PetscInt j = my - 1;
1651 for (k = zs; k < ze; k++) {
1652 for (i = xs; i < xe; i++) {
1653 p[k][j][i] = lp[k][j + 2][i];
1654 phi[k][j][i] = lphi[k][j + 2][i];
1655 }
1656 }
1657 }
1658 }
1659
1660 // --- Z-Direction Periodic Update ---
1661 if (user->bctype[4] == 7 || user->bctype[5] == 7) {
1662 // Update front boundary
1663 if (zs == 0) {
1664 PetscInt k = 0;
1665 for (j = ys; j < ye; j++) {
1666 for (i = xs; i < xe; i++) {
1667 p[k][j][i] = lp[k - 2][j][i];
1668 phi[k][j][i] = lphi[k - 2][j][i];
1669 }
1670 }
1671 }
1672 // Update back boundary
1673 if (ze == mz) {
1674 PetscInt k = mz - 1;
1675 for (j = ys; j < ye; j++) {
1676 for (i = xs; i < xe; i++) {
1677 p[k][j][i] = lp[k + 2][j][i];
1678 phi[k][j][i] = lphi[k + 2][j][i];
1679 }
1680 }
1681 }
1682 }
1683
1684 // Restore all arrays ONCE.
1685 DMDAVecRestoreArray(da, user->lP, &lp);
1686 DMDAVecRestoreArray(da, user->lPhi, &lphi);
1687 DMDAVecRestoreArray(da, user->P, &p);
1688 DMDAVecRestoreArray(da, user->Phi, &phi);
1689
1690 // After manually updating the physical boundary cells, we must call
1691 // DMGlobalToLocal again to ensure all processes have the updated ghost
1692 // values for the *next* function that needs them.
1693 DMGlobalToLocalBegin(da, user->P, INSERT_VALUES, user->lP);
1694 DMGlobalToLocalEnd(da, user->P, INSERT_VALUES, user->lP);
1695 DMGlobalToLocalBegin(da, user->Phi, INSERT_VALUES, user->lPhi);
1696 DMGlobalToLocalEnd(da, user->Phi, INSERT_VALUES, user->lPhi);
1697 }
1698
1699 //================================================================================
1700 // Section 4: Final Cleanup (pointers already restored)
1701 //================================================================================
1702
1703 UpdateLocalGhosts(user,"P");
1704 UpdateLocalGhosts(user,"Phi");
1705
1706 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting UpdatePressure.\n");
1707 PROFILE_FUNCTION_END;
1708 PetscFunctionReturn(0);
1709}
1710
1711static PetscErrorCode mymatmultadd(Mat mat, Vec v1, Vec v2)
1712{
1713
1714 Vec vt;
1715 VecDuplicate(v2, &vt);
1716 MatMult(mat, v1, vt);
1717 VecAYPX(v2, 1., vt);
1718 VecDestroy(&vt);
1719 return(0);
1720}
1721
1722PetscErrorCode PoissonNullSpaceFunction(MatNullSpace nullsp,Vec X, void *ctx)
1723{
1724 UserCtx *user = (UserCtx*)ctx;
1725
1726 DM da = user->da;
1727
1728 DMDALocalInfo info = user->info;
1729 PetscInt xs = info.xs, xe = info.xs + info.xm;
1730 PetscInt ys = info.ys, ye = info.ys + info.ym;
1731 PetscInt zs = info.zs, ze = info.zs + info.zm;
1732 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1733 PetscInt lxs, lxe, lys, lye, lzs, lze;
1734
1735 PetscReal ***x, ***nvert;
1736 PetscInt i, j, k;
1737
1738/* /\* First remove a constant from the Vec field X *\/ */
1739
1740
1741 /* Then apply boundary conditions */
1742 DMDAVecGetArray(da, X, &x);
1743 DMDAVecGetArray(da, user->lNvert, &nvert);
1744
1745 lxs = xs; lxe = xe;
1746 lys = ys; lye = ye;
1747 lzs = zs; lze = ze;
1748
1749 if (xs==0) lxs = xs+1;
1750 if (ys==0) lys = ys+1;
1751 if (zs==0) lzs = zs+1;
1752
1753 if (xe==mx) lxe = xe-1;
1754 if (ye==my) lye = ye-1;
1755 if (ze==mz) lze = ze-1;
1756
1757 PetscReal lsum, sum;
1758 PetscReal lnum, num;
1759
1760 if (user->multinullspace) PetscPrintf(PETSC_COMM_WORLD, "MultiNullSpace!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1761 if (!user->multinullspace) {
1762 lsum = 0;
1763 lnum = 0;
1764 for (k=lzs; k<lze; k++) {
1765 for (j=lys; j<lye; j++) {
1766 for (i=lxs; i<lxe; i++) {
1767 if (nvert[k][j][i] < 0.1) {
1768 lsum += x[k][j][i];
1769 lnum ++;
1770 }
1771 }
1772 }
1773 }
1774
1775 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1776 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1777 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1778/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1779 sum = sum / (-1.0 * num);
1780
1781 for (k=lzs; k<lze; k++) {
1782 for (j=lys; j<lye; j++) {
1783 for (i=lxs; i<lxe; i++) {
1784 if (nvert[k][j][i] < 0.1) {
1785 x[k][j][i] +=sum;
1786 }
1787 }
1788 }
1789 }
1790 }
1791 else {
1792 lsum = 0;
1793 lnum = 0;
1794 for (j=lys; j<lye; j++) {
1795 for (i=lxs; i<lxe; i++) {
1796 for (k=lzs; k<lze; k++) {
1797 if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1798 lsum += x[k][j][i];
1799 lnum ++;
1800 }
1801 }
1802 }
1803 }
1804 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1805 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1806 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1807/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1808 sum /= -num;
1809 for (j=lys; j<lye; j++) {
1810 for (i=lxs; i<lxe; i++) {
1811 for (k=lzs; k<lze; k++) {
1812 if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1813 x[k][j][i] += sum;
1814 }
1815 }
1816 }
1817 }
1818
1819 lsum = 0;
1820 lnum = 0;
1821 for (j=lys; j<lye; j++) {
1822 for (i=lxs; i<lxe; i++) {
1823 for (k=lzs; k<lze; k++) {
1824 if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1825 lsum += x[k][j][i];
1826 lnum ++;
1827 }
1828 }
1829 }
1830 }
1831 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1832 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1833 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1834/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1835 sum /= -num;
1836 for (j=lys; j<lye; j++) {
1837 for (i=lxs; i<lxe; i++) {
1838 for (k=lzs; k<lze; k++) {
1839 if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1840 x[k][j][i] += sum;
1841 }
1842 }
1843 }
1844 }
1845
1846 } //if multinullspace
1847 if (zs == 0) {
1848 k = 0;
1849 for (j=ys; j<ye; j++) {
1850 for (i=xs; i<xe; i++) {
1851 x[k][j][i] = 0.;
1852 }
1853 }
1854 }
1855
1856 if (ze == mz) {
1857 k = mz-1;
1858 for (j=ys; j<ye; j++) {
1859 for (i=xs; i<xe; i++) {
1860 x[k][j][i] = 0.;
1861 }
1862 }
1863 }
1864
1865 if (ys == 0) {
1866 j = 0;
1867 for (k=zs; k<ze; k++) {
1868 for (i=xs; i<xe; i++) {
1869 x[k][j][i] = 0.;
1870 }
1871 }
1872 }
1873
1874 if (ye == my) {
1875 j = my-1;
1876 for (k=zs; k<ze; k++) {
1877 for (i=xs; i<xe; i++) {
1878 x[k][j][i] = 0.;
1879 }
1880 }
1881 }
1882
1883 if (xs == 0) {
1884 i = 0;
1885 for (k=zs; k<ze; k++) {
1886 for (j=ys; j<ye; j++) {
1887 x[k][j][i] = 0.;
1888 }
1889 }
1890 }
1891
1892 if (xe == mx) {
1893 i = mx-1;
1894 for (k=zs; k<ze; k++) {
1895 for (j=ys; j<ye; j++) {
1896 x[k][j][i] = 0.;
1897 }
1898 }
1899 }
1900
1901 for (k=zs; k<ze; k++) {
1902 for (j=ys; j<ye; j++) {
1903 for (i=xs; i<xe; i++) {
1904 if (nvert[k][j][i] > 0.1)
1905 x[k][j][i] = 0.;
1906 }
1907 }
1908 }
1909 DMDAVecRestoreArray(da, X, &x);
1910 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1911
1912 return 0;
1913}
1914
1915PetscErrorCode MyInterpolation(Mat A, Vec X, Vec F)
1916{
1917 UserCtx *user;
1918
1919 MatShellGetContext(A, (void**)&user);
1920
1921
1922
1923 DM da = user->da;
1924
1925 DM da_c = *user->da_c;
1926
1927 DMDALocalInfo info = user->info;
1928 PetscInt xs = info.xs, xe = info.xs + info.xm;
1929 PetscInt ys = info.ys, ye = info.ys + info.ym;
1930 PetscInt zs = info.zs, ze = info.zs + info.zm;
1931 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1932 PetscInt lxs, lxe, lys, lye, lzs, lze;
1933
1934 PetscReal ***f, ***x, ***nvert, ***nvert_c;
1935 PetscInt i, j, k, ic, jc, kc, ia, ja, ka;
1936
1937 lxs = xs; lxe = xe;
1938 lys = ys; lye = ye;
1939 lzs = zs; lze = ze;
1940
1941 if (xs==0) lxs = xs+1;
1942 if (ys==0) lys = ys+1;
1943 if (zs==0) lzs = zs+1;
1944
1945 if (xe==mx) lxe = xe-1;
1946 if (ye==my) lye = ye-1;
1947 if (ze==mz) lze = ze-1;
1948
1949
1950 DMDAVecGetArray(da, F, &f);
1951
1952
1953 Vec lX;
1954 DMCreateLocalVector(da_c, &lX);
1955
1956 DMGlobalToLocalBegin(da_c, X, INSERT_VALUES, lX);
1957 DMGlobalToLocalEnd(da_c, X, INSERT_VALUES, lX);
1958 DMDAVecGetArray(da_c, lX, &x);
1959
1960 DMDAVecGetArray(da, user->lNvert, &nvert);
1961 DMDAVecGetArray(da_c, *(user->lNvert_c), &nvert_c);
1962 for (k=lzs; k<lze; k++) {
1963 for (j=lys; j<lye; j++) {
1964 for (i=lxs; i<lxe; i++) {
1965
1966 GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user);
1967
1968 f[k][j][i] = (x[kc ][jc ][ic ] * 9 +
1969 x[kc ][jc+ja][ic ] * 3 +
1970 x[kc ][jc ][ic+ia] * 3 +
1971 x[kc ][jc+ja][ic+ia]) * 3./64. +
1972 (x[kc+ka][jc ][ic ] * 9 +
1973 x[kc+ka][jc+ja][ic ] * 3 +
1974 x[kc+ka][jc ][ic+ia] * 3 +
1975 x[kc+ka][jc+ja][ic+ia]) /64.;
1976 }
1977 }
1978 }
1979
1980 for (k=zs; k<ze; k++) {
1981 for (j=ys; j<ye; j++) {
1982 for (i=xs; i<xe; i++) {
1983
1984 if (i==0) {
1985 f[k][j][i] = 0.;//-f[k][j][i+1];
1986 }
1987 else if (i==mx-1) {
1988 f[k][j][i] = 0.;//-f[k][j][i-1];
1989 }
1990 else if (j==0) {
1991 f[k][j][i] = 0.;//-f[k][j+1][i];
1992 }
1993 else if (j==my-1) {
1994 f[k][j][i] = 0.;//-f[k][j-1][i];
1995 }
1996 else if (k==0) {
1997 f[k][j][i] = 0.;//-f[k+1][j][i];
1998 }
1999 else if (k==mz-1) {
2000 f[k][j][i] = 0.;//-f[k-1][j][i];
2001 }
2002 if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
2003
2004 }
2005 }
2006 }
2007
2008 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2009 DMDAVecRestoreArray(da_c, *(user->lNvert_c), &nvert_c);
2010
2011 DMDAVecRestoreArray(da_c, lX, &x);
2012
2013 VecDestroy(&lX);
2014 DMDAVecRestoreArray(da, F, &f);
2015
2016
2017
2018 return 0;
2019
2020}
2021
2022static PetscErrorCode RestrictResidual_SolidAware(Mat A, Vec X, Vec F)
2023{
2024 UserCtx *user;
2025 MatShellGetContext(A, (void**)&user);
2026
2027 DM da = user->da;
2028 DM da_f = *user->da_f;
2029
2030 DMDALocalInfo info;
2031 DMDAGetLocalInfo(da, &info);
2032 PetscInt xs = info.xs, xe = info.xs + info.xm;
2033 PetscInt ys = info.ys, ye = info.ys + info.ym;
2034 PetscInt zs = info.zs, ze = info.zs + info.zm;
2035 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2036
2037 PetscReal ***f, ***x, ***nvert;
2038 PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
2039
2040 DMDAVecGetArray(da, F, &f);
2041
2042 Vec lX;
2043 DMCreateLocalVector(da_f, &lX);
2044 DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
2045 DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);
2046 DMDAVecGetArray(da_f, lX, &x);
2047
2048 DMDAVecGetArray(da, user->lNvert, &nvert);
2049
2050 PetscReal ***nvert_f;
2051 DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
2052
2053 if ((user->isc)) ia = 0;
2054 else ia = 1;
2055
2056 if ((user->jsc)) ja = 0;
2057 else ja = 1;
2058
2059 if ((user->ksc)) ka = 0;
2060 else ka = 1;
2061
2062 for (k=zs; k<ze; k++) {
2063 for (j=ys; j<ye; j++) {
2064 for (i=xs; i<xe; i++) {
2065 // --- CORRECTED LOGIC ---
2066 // First, check if the current point is a boundary point.
2067 // If it is, it does not contribute to the coarse grid residual.
2068 if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1 || nvert[k][j][i] > 0.1) {
2069 f[k][j][i] = 0.0;
2070 }
2071 // Only if it's a true interior fluid point, perform the restriction.
2072 else {
2073 GridRestriction(i, j, k, &ih, &jh, &kh, user);
2074 f[k][j][i] = 0.125 *
2075 (x[kh ][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih ]) +
2076 x[kh ][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih-ia]) +
2077 x[kh ][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih ]) +
2078 x[kh-ka][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih ]) +
2079 x[kh ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih-ia]) +
2080 x[kh-ka][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih ]) +
2081 x[kh-ka][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih-ia]) +
2082 x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
2083 }
2084 }
2085 }
2086 }
2087
2088 DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
2089 DMDAVecRestoreArray(da_f, lX, &x);
2090 VecDestroy(&lX);
2091 DMDAVecRestoreArray(da, F, &f);
2092 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2093
2094 return 0;
2095}
2096
2097PetscErrorCode MyRestriction(Mat A, Vec X, Vec F)
2098{
2099 UserCtx *user;
2100
2101 MatShellGetContext(A, (void**)&user);
2102
2103
2104 DM da = user->da;
2105
2106 DM da_f = *user->da_f;
2107
2108 DMDALocalInfo info;
2109 DMDAGetLocalInfo(da, &info);
2110 PetscInt xs = info.xs, xe = info.xs + info.xm;
2111 PetscInt ys = info.ys, ye = info.ys + info.ym;
2112 PetscInt zs = info.zs, ze = info.zs + info.zm;
2113 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2114 // PetscInt lxs, lxe, lys, lye, lzs, lze;
2115
2116 PetscReal ***f, ***x, ***nvert;
2117 PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
2118
2119 DMDAVecGetArray(da, F, &f);
2120
2121 Vec lX;
2122
2123 DMCreateLocalVector(da_f, &lX);
2124 DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
2125 DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);
2126 DMDAVecGetArray(da_f, lX, &x);
2127
2128 DMDAVecGetArray(da, user->lNvert, &nvert);
2129
2130 PetscReal ***nvert_f;
2131 DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
2132
2133 if ((user->isc)) ia = 0;
2134 else ia = 1;
2135
2136 if ((user->jsc)) ja = 0;
2137 else ja = 1;
2138
2139 if ((user->ksc)) ka = 0;
2140 else ka = 1;
2141
2142 for (k=zs; k<ze; k++) {
2143 for (j=ys; j<ye; j++) {
2144 for (i=xs; i<xe; i++) {
2145 if (k==0) {
2146 f[k][j][i] = 0.;
2147 }
2148 else if (k==mz-1) {
2149 f[k][j][i] = 0.;
2150 }
2151 else if (j==0) {
2152 f[k][j][i] = 0.;
2153 }
2154 else if (j==my-1) {
2155 f[k][j][i] = 0.;
2156 }
2157 else if (i==0) {
2158 f[k][j][i] = 0.;
2159 }
2160 else if (i==mx-1) {
2161 f[k][j][i] = 0.;
2162 }
2163 else {
2164 GridRestriction(i, j, k, &ih, &jh, &kh, user);
2165 f[k][j][i] = 0.125 *
2166 (x[kh ][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih ]) +
2167 x[kh ][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih-ia]) +
2168 x[kh ][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih ]) +
2169 x[kh-ka][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih ]) +
2170 x[kh ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih-ia]) +
2171 x[kh-ka][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih ]) +
2172 x[kh-ka][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih-ia]) +
2173 x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
2174
2175
2176
2177 if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
2178 }
2179 }
2180 }
2181 }
2182
2183
2184 DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
2185
2186 DMDAVecRestoreArray(da_f, lX, &x);
2187 VecDestroy(&lX);
2188
2189 DMDAVecRestoreArray(da, F, &f);
2190 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2191
2192
2193 return 0;
2194}
2195
2196/**
2197 * @brief Assembles the Left-Hand Side (LHS) matrix for the Pressure Poisson Equation.
2198 *
2199 * @details
2200 * This function constructs the sparse matrix `A` (the LHS) for the linear system `Ax = B`,
2201 * which is the Pressure Poisson Equation (PPE). The matrix `A` represents a discrete
2202 * version of the negative Laplacian operator (-∇²), tailored for a general curvilinear,
2203 * staggered grid.
2204 *
2205 * The assembly process is highly complex and follows these main steps:
2206 *
2207 * 1. **Matrix Initialization**: On the first call, it allocates a sparse PETSc matrix `A`
2208 * pre-allocating space for a 19-point stencil per row. On subsequent calls, it
2209 * simply zeroes out the existing matrix.
2210 *
2211 * 2. **Metric Tensor Calculation**: It first computes the 9 components of the metric
2212 * tensor (`g11`, `g12`, ..., `g33`) at the cell faces. These `g_ij` coefficients
2213 * are essential for defining the Laplacian operator in generalized curvilinear
2214 * coordinates and account for grid stretching and non-orthogonality.
2215 *
2216 * 3. **Stencil Assembly Loop**: The function iterates through every grid point `(i,j,k)`.
2217 * For each point, it determines the matrix row entries based on its status:
2218 * - **Boundary/Solid Point**: If the point is on a domain boundary or inside an
2219 * immersed solid (`nvert > 0.1`), it sets an identity row (`A(row,row) = 1`),
2220 * effectively removing it from the linear solve.
2221 * - **Fluid Point**: For a fluid point, it computes the 19 coefficients of the
2222 * finite volume stencil. This involves summing contributions from each of the
2223 * six faces of the control volume around the point.
2224 *
2225 * 4. **Boundary-Aware Stencils**: The stencil calculation is critically dependent on the
2226 * state of neighboring cells. The code contains intricate logic to check if neighbors
2227 * are fluid or solid. If a standard centered-difference stencil would cross into a
2228 * solid, the scheme is automatically adapted to a one-sided difference to maintain
2229 * accuracy at the fluid-solid interface.
2230 *
2231 * 5. **Matrix Finalization**: After all rows corresponding to the local processor's
2232 * grid points have been set, `MatAssemblyBegin/End` is called to finalize the
2233 * global matrix, making it ready for use by a linear solver.
2234 *
2235 * @param[in, out] user Pointer to the UserCtx struct, containing all simulation context,
2236 * grid data, and the matrix `A`.
2237 *
2238 * @return PetscErrorCode 0 on success, or a non-zero error code on failure.
2239 */
2240PetscErrorCode PoissonLHSNew(UserCtx *user)
2241{
2242 PetscFunctionBeginUser;
2243 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonLHSNew to assemble Laplacian matrix.\n");
2244 PetscErrorCode ierr;
2245 //================================================================================
2246 // Section 1: Initialization and Data Acquisition
2247 //================================================================================
2248
2249
2250 // --- Get simulation and grid context ---
2251 SimCtx *simCtx = user->simCtx;
2252 DM da = user->da, fda = user->fda;
2253 DMDALocalInfo info = user->info;
2254 PetscInt IM = user->IM, JM = user->JM, KM = user->KM;
2255 PetscInt i,j,k;
2256
2257 // --- Grid dimensions ---
2258 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2259 PetscInt xs = info.xs, xe = info.xs + info.xm;
2260 PetscInt ys = info.ys, ye = info.ys + info.ym;
2261 PetscInt zs = info.zs, ze = info.zs + info.zm;
2262 PetscInt gxs = info.gxs, gxe = gxs + info.gxm;
2263 PetscInt gys = info.gys, gye = gys + info.gym;
2264 PetscInt gzs = info.gzs, gze = gzs + info.gzm;
2265
2266 // --- Define constants for clarity ---
2267 const PetscReal IBM_FLUID_THRESHOLD = 0.1;
2268
2269 // --- Allocate the LHS matrix A on the first call ---
2270 if (!user->assignedA) {
2271 LOG_ALLOW(GLOBAL, LOG_INFO, "First call: Creating LHS matrix 'A' with 19-point stencil preallocation.\n");
2272 PetscInt N = mx * my * mz; // Total size
2273 PetscInt M; // Local size
2274 VecGetLocalSize(user->Phi, &M);
2275 // Create a sparse AIJ matrix, preallocating for 19 non-zeros per row (d=diagonal, o=off-diagonal)
2276 MatCreateAIJ(PETSC_COMM_WORLD, M, M, N, N, 19, PETSC_NULL, 19, PETSC_NULL, &(user->A));
2277 user->assignedA = PETSC_TRUE;
2278 }
2279
2280 // Zero out matrix entries from the previous solve
2281 MatZeroEntries(user->A);
2282
2283 // --- Get direct pointer access to grid metric data ---
2284 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
2285 PetscReal ***aj, ***iaj, ***jaj, ***kaj, ***nvert;
2286 DMDAVecGetArray(fda, user->lCsi, &csi); DMDAVecGetArray(fda, user->lEta, &eta); DMDAVecGetArray(fda, user->lZet, &zet);
2287 DMDAVecGetArray(fda, user->lICsi, &icsi); DMDAVecGetArray(fda, user->lIEta, &ieta); DMDAVecGetArray(fda, user->lIZet, &izet);
2288 DMDAVecGetArray(fda, user->lJCsi, &jcsi); DMDAVecGetArray(fda, user->lJEta, &jeta); DMDAVecGetArray(fda, user->lJZet, &jzet);
2289 DMDAVecGetArray(fda, user->lKCsi, &kcsi); DMDAVecGetArray(fda, user->lKEta, &keta); DMDAVecGetArray(fda, user->lKZet, &kzet);
2290 DMDAVecGetArray(da, user->lAj, &aj); DMDAVecGetArray(da, user->lIAj, &iaj); DMDAVecGetArray(da, user->lJAj, &jaj); DMDAVecGetArray(da, user->lKAj, &kaj);
2291 DMDAVecGetArray(da, user->lNvert, &nvert);
2292
2293 // --- Create temporary vectors for the metric tensor components G_ij ---
2294 Vec G11, G12, G13, G21, G22, G23, G31, G32, G33;
2295 PetscReal ***g11, ***g12, ***g13, ***g21, ***g22, ***g23, ***g31, ***g32, ***g33;
2296 VecDuplicate(user->lAj, &G11); VecDuplicate(user->lAj, &G12); VecDuplicate(user->lAj, &G13);
2297 VecDuplicate(user->lAj, &G21); VecDuplicate(user->lAj, &G22); VecDuplicate(user->lAj, &G23);
2298 VecDuplicate(user->lAj, &G31); VecDuplicate(user->lAj, &G32); VecDuplicate(user->lAj, &G33);
2299 DMDAVecGetArray(da, G11, &g11); DMDAVecGetArray(da, G12, &g12); DMDAVecGetArray(da, G13, &g13);
2300 DMDAVecGetArray(da, G21, &g21); DMDAVecGetArray(da, G22, &g22); DMDAVecGetArray(da, G23, &g23);
2301 DMDAVecGetArray(da, G31, &g31); DMDAVecGetArray(da, G32, &g32); DMDAVecGetArray(da, G33, &g33);
2302
2303 //================================================================================
2304 // Section 2: Pre-compute Metric Tensor Coefficients (g_ij)
2305 //================================================================================
2306 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pre-computing metric tensor components (g_ij).\n");
2307 for (k = gzs; k < gze; k++) {
2308 for (j = gys; j < gye; j++) {
2309 for (i = gxs; i < gxe; i++) {
2310 // These coefficients represent the dot products of the grid's contravariant base vectors,
2311 // scaled by face area. They are the core of the Laplacian operator on a curvilinear grid.
2312 if(i>-1 && j>-1 && k>-1 && i<IM+1 && j<JM+1 && k<KM+1){
2313 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];
2314 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];
2315 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];
2316 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];
2317 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];
2318 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];
2319 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];
2320 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];
2321 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];
2322 }
2323 }
2324 }
2325 }
2326
2327 //================================================================================
2328 // Section 3: Assemble the LHS Matrix A
2329 //================================================================================
2330 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Assembling the LHS matrix A using a 19-point stencil.\n");
2331
2332 // --- Define domain boundaries for stencil logic, accounting for periodic BCs ---
2333 PetscInt x_str, x_end, y_str, y_end, z_str, z_end;
2334 if (user->bctype[0] == 7) { x_end = mx - 1; x_str = 0; }
2335 else { x_end = mx - 2; x_str = 1; }
2336 if (user->bctype[2] == 7) { y_end = my - 1; y_str = 0; }
2337 else { y_end = my - 2; y_str = 1; }
2338 if (user->bctype[4] == 7) { z_end = mz - 1; z_str = 0; }
2339 else { z_end = mz - 2; z_str = 1; }
2340
2341 // --- Main assembly loop over all local grid points ---
2342 for (k = zs; k < ze; k++) {
2343 for (j = ys; j < ye; j++) {
2344 for (i = xs; i < xe; i++) {
2345 PetscScalar vol[19]; // Holds the 19 stencil coefficient values for the current row
2346 PetscInt idx[19]; // Holds the 19 global column indices for the current row
2347 PetscInt row = Gidx(i, j, k, user); // Global index for the current row
2348
2349 // --- Handle Domain Boundary and Immersed Solid Points ---
2350 // For these points, we don't solve the Poisson equation. We set an identity
2351 // row (A_ii = 1) to effectively fix the pressure value (usually to 0).
2352 if (i == 0 || i == mx - 1 || j == 0 || j == my - 1 || k == 0 || k == mz - 1 || nvert[k][j][i] > IBM_FLUID_THRESHOLD) {
2353 vol[CP] = 1.0;
2354 idx[CP] = row;
2355 MatSetValues(user->A, 1, &row, 1, &idx[CP], &vol[CP], INSERT_VALUES);
2356 }
2357 // --- Handle Fluid Points ---
2358 else {
2359 for (PetscInt m = 0; m < 19; m++) {
2360 vol[m] = 0.0;
2361 }
2362
2363 /************************************************************************
2364 * EAST FACE CONTRIBUTION (between i and i+1)
2365 ************************************************************************/
2366 if (nvert[k][j][i + 1] < IBM_FLUID_THRESHOLD && i != x_end) { // East neighbor is fluid
2367 // Primary derivative term: d/d_csi (g11 * dP/d_csi)
2368 vol[CP] -= g11[k][j][i];
2369 vol[EP] += g11[k][j][i];
2370
2371 // Cross-derivative term: d/d_csi (g12 * dP/d_eta).
2372 // This requires an average of dP/d_eta. If a neighbor is solid, the stencil
2373 // dynamically switches to a one-sided difference to avoid using solid points.
2374 if ((j == my-2 && user->bctype[2]!=7) || nvert[k][j+1][i] + nvert[k][j+1][i+1] > 0.1) {
2375 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
2376 vol[CP] += g12[k][j][i] * 0.5; vol[EP] += g12[k][j][i] * 0.5;
2377 vol[SP] -= g12[k][j][i] * 0.5; vol[SE] -= g12[k][j][i] * 0.5;
2378 }
2379 }
2380 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k][j+1][i+1] > 0.1) {
2381 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
2382 vol[CP] += g12[k][j][i] * 0.5; vol[EP] += g12[k][j][i] * 0.5;
2383 vol[SP] -= g12[k][j][i] * 0.5; vol[SE] -= g12[k][j][i] * 0.5;
2384 }
2385 }
2386 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
2387 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
2388 vol[NP] += g12[k][j][i] * 0.5; vol[NE] += g12[k][j][i] * 0.5;
2389 vol[CP] -= g12[k][j][i] * 0.5; vol[EP] -= g12[k][j][i] * 0.5;
2390 }
2391 }
2392 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
2393 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
2394 vol[NP] += g12[k][j][i] * 0.5; vol[NE] += g12[k][j][i] * 0.5;
2395 vol[CP] -= g12[k][j][i] * 0.5; vol[EP] -= g12[k][j][i] * 0.5;
2396 }
2397 }
2398 else { // Centered difference
2399 vol[NP] += g12[k][j][i] * 0.25; vol[NE] += g12[k][j][i] * 0.25;
2400 vol[SP] -= g12[k][j][i] * 0.25; vol[SE] -= g12[k][j][i] * 0.25;
2401 }
2402
2403 // Cross-derivative term: d/d_csi (g13 * dP/d_zet)
2404 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
2405 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
2406 vol[CP] += g13[k][j][i] * 0.5; vol[EP] += g13[k][j][i] * 0.5;
2407 vol[BP] -= g13[k][j][i] * 0.5; vol[BE] -= g13[k][j][i] * 0.5;
2408 }
2409 }
2410 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
2411 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
2412 vol[CP] += g13[k][j][i] * 0.5; vol[EP] += g13[k][j][i] * 0.5;
2413 vol[BP] -= g13[k][j][i] * 0.5; vol[BE] -= g13[k][j][i] * 0.5;
2414 }
2415 }
2416 else if ((k == 1 && user->bctype[4]!=7) || nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
2417 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
2418 vol[TP] += g13[k][j][i] * 0.5; vol[TE] += g13[k][j][i] * 0.5;
2419 vol[CP] -= g13[k][j][i] * 0.5; vol[EP] -= g13[k][j][i] * 0.5;
2420 }
2421 }
2422 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
2423 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
2424 vol[TP] += g13[k][j][i] * 0.5; vol[TE] += g13[k][j][i] * 0.5;
2425 vol[CP] -= g13[k][j][i] * 0.5; vol[EP] -= g13[k][j][i] * 0.5;
2426 }
2427 }
2428 else { // Centered difference
2429 vol[TP] += g13[k][j][i] * 0.25; vol[TE] += g13[k][j][i] * 0.25;
2430 vol[BP] -= g13[k][j][i] * 0.25; vol[BE] -= g13[k][j][i] * 0.25;
2431 }
2432 }
2433
2434 /************************************************************************
2435 * WEST FACE CONTRIBUTION (between i-1 and i)
2436 ************************************************************************/
2437 if (nvert[k][j][i-1] < IBM_FLUID_THRESHOLD && i != x_str) {
2438 vol[CP] -= g11[k][j][i-1];
2439 vol[WP] += g11[k][j][i-1];
2440
2441 if ((j == my-2 && user->bctype[2]!=7) || nvert[k][j+1][i] + nvert[k][j+1][i-1] > 0.1) {
2442 if (nvert[k][j-1][i] + nvert[k][j-1][i-1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
2443 vol[CP] -= g12[k][j][i-1] * 0.5; vol[WP] -= g12[k][j][i-1] * 0.5;
2444 vol[SP] += g12[k][j][i-1] * 0.5; vol[SW] += g12[k][j][i-1] * 0.5;
2445 }
2446 }
2447 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k][j+1][i-1] > 0.1) {
2448 if (nvert[k][j-1][i] + nvert[k][j-1][i-1] < 0.1) {
2449 vol[CP] -= g12[k][j][i-1] * 0.5; vol[WP] -= g12[k][j][i-1] * 0.5;
2450 vol[SP] += g12[k][j][i-1] * 0.5; vol[SW] += g12[k][j][i-1] * 0.5;
2451 }
2452 }
2453 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k][j-1][i-1] > 0.1) {
2454 if (nvert[k][j+1][i] + nvert[k][j+1][i-1] < 0.1) {
2455 vol[NP] -= g12[k][j][i-1] * 0.5; vol[NW] -= g12[k][j][i-1] * 0.5;
2456 vol[CP] += g12[k][j][i-1] * 0.5; vol[WP] += g12[k][j][i-1] * 0.5;
2457 }
2458 }
2459 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i-1] > 0.1) {
2460 if (nvert[k][j+1][i] + nvert[k][j+1][i-1] < 0.1) {
2461 vol[NP] -= g12[k][j][i-1] * 0.5; vol[NW] -= g12[k][j][i-1] * 0.5;
2462 vol[CP] += g12[k][j][i-1] * 0.5; vol[WP] += g12[k][j][i-1] * 0.5;
2463 }
2464 }
2465 else {
2466 vol[NP] -= g12[k][j][i-1] * 0.25; vol[NW] -= g12[k][j][i-1] * 0.25;
2467 vol[SP] += g12[k][j][i-1] * 0.25; vol[SW] += g12[k][j][i-1] * 0.25;
2468 }
2469
2470 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i-1] > 0.1) {
2471 if (nvert[k-1][j][i] + nvert[k-1][j][i-1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
2472 vol[CP] -= g13[k][j][i-1] * 0.5; vol[WP] -= g13[k][j][i-1] * 0.5;
2473 vol[BP] += g13[k][j][i-1] * 0.5; vol[BW] += g13[k][j][i-1] * 0.5;
2474 }
2475 }
2476 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i-1] > 0.1) {
2477 if (nvert[k-1][j][i] + nvert[k-1][j][i-1] < 0.1) {
2478 vol[CP] -= g13[k][j][i-1] * 0.5; vol[WP] -= g13[k][j][i-1] * 0.5;
2479 vol[BP] += g13[k][j][i-1] * 0.5; vol[BW] += g13[k][j][i-1] * 0.5;
2480 }
2481 }
2482 else if ((k == 1 && user->bctype[4]!=7) || nvert[k-1][j][i] + nvert[k-1][j][i-1] > 0.1) {
2483 if (nvert[k+1][j][i] + nvert[k+1][j][i-1] < 0.1) {
2484 vol[TP] -= g13[k][j][i-1] * 0.5; vol[TW] -= g13[k][j][i-1] * 0.5;
2485 vol[CP] += g13[k][j][i-1] * 0.5; vol[WP] += g13[k][j][i-1] * 0.5;
2486 }
2487 }
2488 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i-1] > 0.1) {
2489 if (nvert[k+1][j][i] + nvert[k+1][j][i-1] < 0.1) {
2490 vol[TP] -= g13[k][j][i-1] * 0.5; vol[TW] -= g13[k][j][i-1] * 0.5;
2491 vol[CP] += g13[k][j][i-1] * 0.5; vol[WP] += g13[k][j][i-1] * 0.5;
2492 }
2493 }
2494 else {
2495 vol[TP] -= g13[k][j][i-1] * 0.25; vol[TW] -= g13[k][j][i-1] * 0.25;
2496 vol[BP] += g13[k][j][i-1] * 0.25; vol[BW] += g13[k][j][i-1] * 0.25;
2497 }
2498 }
2499
2500 /************************************************************************
2501 * NORTH FACE CONTRIBUTION (between j and j+1)
2502 ************************************************************************/
2503 if (nvert[k][j+1][i] < IBM_FLUID_THRESHOLD && j != y_end) {
2504 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
2505 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
2506 vol[CP] += g21[k][j][i] * 0.5; vol[NP] += g21[k][j][i] * 0.5;
2507 vol[WP] -= g21[k][j][i] * 0.5; vol[NW] -= g21[k][j][i] * 0.5;
2508 }
2509 }
2510 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
2511 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
2512 vol[CP] += g21[k][j][i] * 0.5; vol[NP] += g21[k][j][i] * 0.5;
2513 vol[WP] -= g21[k][j][i] * 0.5; vol[NW] -= g21[k][j][i] * 0.5;
2514 }
2515 }
2516 else if ((i == 1 && user->bctype[0]!=7) || nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
2517 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
2518 vol[EP] += g21[k][j][i] * 0.5; vol[NE] += g21[k][j][i] * 0.5;
2519 vol[CP] -= g21[k][j][i] * 0.5; vol[NP] -= g21[k][j][i] * 0.5;
2520 }
2521 }
2522 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
2523 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
2524 vol[EP] += g21[k][j][i] * 0.5; vol[NE] += g21[k][j][i] * 0.5;
2525 vol[CP] -= g21[k][j][i] * 0.5; vol[NP] -= g21[k][j][i] * 0.5;
2526 }
2527 }
2528 else {
2529 vol[EP] += g21[k][j][i] * 0.25; vol[NE] += g21[k][j][i] * 0.25;
2530 vol[WP] -= g21[k][j][i] * 0.25; vol[NW] -= g21[k][j][i] * 0.25;
2531 }
2532
2533 vol[CP] -= g22[k][j][i];
2534 vol[NP] += g22[k][j][i];
2535
2536 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
2537 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
2538 vol[CP] += g23[k][j][i] * 0.5; vol[NP] += g23[k][j][i] * 0.5;
2539 vol[BP] -= g23[k][j][i] * 0.5; vol[BN] -= g23[k][j][i] * 0.5;
2540 }
2541 }
2542 else if ((k == mz-2 || k==1 ) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
2543 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
2544 vol[CP] += g23[k][j][i] * 0.5; vol[NP] += g23[k][j][i] * 0.5;
2545 vol[BP] -= g23[k][j][i] * 0.5; vol[BN] -= g23[k][j][i] * 0.5;
2546 }
2547 }
2548 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
2549 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
2550 vol[TP] += g23[k][j][i] * 0.5; vol[TN] += g23[k][j][i] * 0.5;
2551 vol[CP] -= g23[k][j][i] * 0.5; vol[NP] -= g23[k][j][i] * 0.5;
2552 }
2553 }
2554 else if ((k == 1 || k==mz-2 ) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
2555 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
2556 vol[TP] += g23[k][j][i] * 0.5; vol[TN] += g23[k][j][i] * 0.5;
2557 vol[CP] -= g23[k][j][i] * 0.5; vol[NP] -= g23[k][j][i] * 0.5;
2558 }
2559 }
2560 else {
2561 vol[TP] += g23[k][j][i] * 0.25; vol[TN] += g23[k][j][i] * 0.25;
2562 vol[BP] -= g23[k][j][i] * 0.25; vol[BN] -= g23[k][j][i] * 0.25;
2563 }
2564 }
2565
2566 /************************************************************************
2567 * SOUTH FACE CONTRIBUTION (between j-1 and j)
2568 ************************************************************************/
2569 if (nvert[k][j-1][i] < IBM_FLUID_THRESHOLD && j != y_str) {
2570 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j-1][i+1] > 0.1) {
2571 if (nvert[k][j][i-1] + nvert[k][j-1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
2572 vol[CP] -= g21[k][j-1][i] * 0.5; vol[SP] -= g21[k][j-1][i] * 0.5;
2573 vol[WP] += g21[k][j-1][i] * 0.5; vol[SW] += g21[k][j-1][i] * 0.5;
2574 }
2575 }
2576 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j-1][i+1] > 0.1) {
2577 if (nvert[k][j][i-1] + nvert[k][j-1][i-1] < 0.1) {
2578 vol[CP] -= g21[k][j-1][i] * 0.5; vol[SP] -= g21[k][j-1][i] * 0.5;
2579 vol[WP] += g21[k][j-1][i] * 0.5; vol[SW] += g21[k][j-1][i] * 0.5;
2580 }
2581 }
2582 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j-1][i-1] > 0.1) {
2583 if (nvert[k][j][i+1] + nvert[k][j-1][i+1] < 0.1) {
2584 vol[EP] -= g21[k][j-1][i] * 0.5; vol[SE] -= g21[k][j-1][i] * 0.5;
2585 vol[CP] += g21[k][j-1][i] * 0.5; vol[SP] += g21[k][j-1][i] * 0.5;
2586 }
2587 }
2588 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j-1][i-1] > 0.1) {
2589 if (nvert[k][j][i+1] + nvert[k][j-1][i+1] < 0.1) {
2590 vol[EP] -= g21[k][j-1][i] * 0.5; vol[SE] -= g21[k][j-1][i] * 0.5;
2591 vol[CP] += g21[k][j-1][i] * 0.5; vol[SP] += g21[k][j-1][i] * 0.5;
2592 }
2593 }
2594 else {
2595 vol[EP] -= g21[k][j-1][i] * 0.25; vol[SE] -= g21[k][j-1][i] * 0.25;
2596 vol[WP] += g21[k][j-1][i] * 0.25; vol[SW] += g21[k][j-1][i] * 0.25;
2597 }
2598
2599 vol[CP] -= g22[k][j-1][i];
2600 vol[SP] += g22[k][j-1][i];
2601
2602 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j-1][i] > 0.1) {
2603 if (nvert[k-1][j][i] + nvert[k-1][j-1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
2604 vol[CP] -= g23[k][j-1][i] * 0.5; vol[SP] -= g23[k][j-1][i] * 0.5;
2605 vol[BP] += g23[k][j-1][i] * 0.5; vol[BS] += g23[k][j-1][i] * 0.5;
2606 }
2607 }
2608 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j-1][i] > 0.1) {
2609 if (nvert[k-1][j][i] + nvert[k-1][j-1][i] < 0.1 ) {
2610 vol[CP] -= g23[k][j-1][i] * 0.5; vol[SP] -= g23[k][j-1][i] * 0.5;
2611 vol[BP] += g23[k][j-1][i] * 0.5; vol[BS] += g23[k][j-1][i] * 0.5;
2612 }
2613 }
2614 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j-1][i] > 0.1) {
2615 if (nvert[k+1][j][i] + nvert[k+1][j-1][i] < 0.1) {
2616 vol[TP] -= g23[k][j-1][i] * 0.5; vol[TS] -= g23[k][j-1][i] * 0.5;
2617 vol[CP] += g23[k][j-1][i] * 0.5; vol[SP] += g23[k][j-1][i] * 0.5;
2618 }
2619 }
2620 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j-1][i] > 0.1) {
2621 if (nvert[k+1][j][i] + nvert[k+1][j-1][i] < 0.1) {
2622 vol[TP] -= g23[k][j-1][i] * 0.5; vol[TS] -= g23[k][j-1][i] * 0.5;
2623 vol[CP] += g23[k][j-1][i] * 0.5; vol[SP] += g23[k][j-1][i] * 0.5;
2624 }
2625 }
2626 else {
2627 vol[TP] -= g23[k][j-1][i] * 0.25; vol[TS] -= g23[k][j-1][i] * 0.25;
2628 vol[BP] += g23[k][j-1][i] * 0.25; vol[BS] += g23[k][j-1][i] * 0.25;
2629 }
2630 }
2631
2632 /************************************************************************
2633 * TOP FACE CONTRIBUTION (between k and k+1)
2634 ************************************************************************/
2635 if (nvert[k+1][j][i] < IBM_FLUID_THRESHOLD && k != z_end) {
2636 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
2637 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
2638 vol[CP] += g31[k][j][i] * 0.5; vol[TP] += g31[k][j][i] * 0.5;
2639 vol[WP] -= g31[k][j][i] * 0.5; vol[TW] -= g31[k][j][i] * 0.5;
2640 }
2641 }
2642 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
2643 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
2644 vol[CP] += g31[k][j][i] * 0.5; vol[TP] += g31[k][j][i] * 0.5;
2645 vol[WP] -= g31[k][j][i] * 0.5; vol[TW] -= g31[k][j][i] * 0.5;
2646 }
2647 }
2648 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
2649 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
2650 vol[EP] += g31[k][j][i] * 0.5; vol[TE] += g31[k][j][i] * 0.5;
2651 vol[CP] -= g31[k][j][i] * 0.5; vol[TP] -= g31[k][j][i] * 0.5;
2652 }
2653 }
2654 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
2655 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
2656 vol[EP] += g31[k][j][i] * 0.5; vol[TE] += g31[k][j][i] * 0.5;
2657 vol[CP] -= g31[k][j][i] * 0.5; vol[TP] -= g31[k][j][i] * 0.5;
2658 }
2659 }
2660 else {
2661 vol[EP] += g31[k][j][i] * 0.25; vol[TE] += g31[k][j][i] * 0.25;
2662 vol[WP] -= g31[k][j][i] * 0.25; vol[TW] -= g31[k][j][i] * 0.25;
2663 }
2664
2665 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
2666 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
2667 vol[CP] += g32[k][j][i] * 0.5; vol[TP] += g32[k][j][i] * 0.5;
2668 vol[SP] -= g32[k][j][i] * 0.5; vol[TS] -= g32[k][j][i] * 0.5;
2669 }
2670 }
2671 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
2672 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
2673 vol[CP] += g32[k][j][i] * 0.5; vol[TP] += g32[k][j][i] * 0.5;
2674 vol[SP] -= g32[k][j][i] * 0.5; vol[TS] -= g32[k][j][i] * 0.5;
2675 }
2676 }
2677 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
2678 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
2679 vol[NP] += g32[k][j][i] * 0.5; vol[TN] += g32[k][j][i] * 0.5;
2680 vol[CP] -= g32[k][j][i] * 0.5; vol[TP] -= g32[k][j][i] * 0.5;
2681 }
2682 }
2683 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
2684 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
2685 vol[NP] += g32[k][j][i] * 0.5; vol[TN] += g32[k][j][i] * 0.5;
2686 vol[CP] -= g32[k][j][i] * 0.5; vol[TP] -= g32[k][j][i] * 0.5;
2687 }
2688 }
2689 else {
2690 vol[NP] += g32[k][j][i] * 0.25; vol[TN] += g32[k][j][i] * 0.25;
2691 vol[SP] -= g32[k][j][i] * 0.25; vol[TS] -= g32[k][j][i] * 0.25;
2692 }
2693
2694 vol[CP] -= g33[k][j][i];
2695 vol[TP] += g33[k][j][i];
2696 }
2697
2698 /************************************************************************
2699 * BOTTOM FACE CONTRIBUTION (between k-1 and k)
2700 ************************************************************************/
2701 if (nvert[k-1][j][i] < IBM_FLUID_THRESHOLD && k != z_str) {
2702 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k-1][j][i+1] > 0.1) {
2703 if (nvert[k][j][i-1] + nvert[k-1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
2704 vol[CP] -= g31[k-1][j][i] * 0.5; vol[BP] -= g31[k-1][j][i] * 0.5;
2705 vol[WP] += g31[k-1][j][i] * 0.5; vol[BW] += g31[k-1][j][i] * 0.5;
2706 }
2707 }
2708 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k-1][j][i+1] > 0.1) {
2709 if (nvert[k][j][i-1] + nvert[k-1][j][i-1] < 0.1) {
2710 vol[CP] -= g31[k-1][j][i] * 0.5; vol[BP] -= g31[k-1][j][i] * 0.5;
2711 vol[WP] += g31[k-1][j][i] * 0.5; vol[BW] += g31[k-1][j][i] * 0.5;
2712 }
2713 }
2714 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k-1][j][i-1] > 0.1) {
2715 if (nvert[k][j][i+1] + nvert[k-1][j][i+1] < 0.1) {
2716 vol[EP] -= g31[k-1][j][i] * 0.5; vol[BE] -= g31[k-1][j][i] * 0.5;
2717 vol[CP] += g31[k-1][j][i] * 0.5; vol[BP] += g31[k-1][j][i] * 0.5;
2718 }
2719 }
2720 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k-1][j][i-1] > 0.1) {
2721 if (nvert[k][j][i+1] + nvert[k-1][j][i+1] < 0.1) {
2722 vol[EP] -= g31[k-1][j][i] * 0.5; vol[BE] -= g31[k-1][j][i] * 0.5;
2723 vol[CP] += g31[k-1][j][i] * 0.5; vol[BP] += g31[k-1][j][i] * 0.5;
2724 }
2725 }
2726 else {
2727 vol[EP] -= g31[k-1][j][i] * 0.25; vol[BE] -= g31[k-1][j][i] * 0.25;
2728 vol[WP] += g31[k-1][j][i] * 0.25; vol[BW] += g31[k-1][j][i] * 0.25;
2729 }
2730
2731 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k-1][j+1][i] > 0.1) {
2732 if (nvert[k][j-1][i] + nvert[k-1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
2733 vol[CP] -= g32[k-1][j][i] * 0.5; vol[BP] -= g32[k-1][j][i] * 0.5;
2734 vol[SP] += g32[k-1][j][i] * 0.5; vol[BS] += g32[k-1][j][i] * 0.5;
2735 }
2736 }
2737 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k-1][j+1][i] > 0.1) {
2738 if (nvert[k][j-1][i] + nvert[k-1][j-1][i] < 0.1) {
2739 vol[CP] -= g32[k-1][j][i] * 0.5; vol[BP] -= g32[k-1][j][i] * 0.5;
2740 vol[SP] += g32[k-1][j][i] * 0.5; vol[BS] += g32[k-1][j][i] * 0.5;
2741 }
2742 }
2743 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k-1][j-1][i] > 0.1) {
2744 if (nvert[k][j+1][i] + nvert[k-1][j+1][i] < 0.1) {
2745 vol[NP] -= g32[k-1][j][i] * 0.5; vol[BN] -= g32[k-1][j][i] * 0.5;
2746 vol[CP] += g32[k-1][j][i] * 0.5; vol[BP] += g32[k-1][j][i] * 0.5;
2747 }
2748 }
2749 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k-1][j-1][i] > 0.1) {
2750 if (nvert[k][j+1][i] + nvert[k-1][j+1][i] < 0.1) {
2751 vol[NP] -= g32[k-1][j][i] * 0.5; vol[BN] -= g32[k-1][j][i] * 0.5;
2752 vol[CP] += g32[k-1][j][i] * 0.5; vol[BP] += g32[k-1][j][i] * 0.5;
2753 }
2754 }
2755 else {
2756 vol[NP] -= g32[k-1][j][i] * 0.25; vol[BN] -= g32[k-1][j][i] * 0.25;
2757 vol[SP] += g32[k-1][j][i] * 0.25; vol[BS] += g32[k-1][j][i] * 0.25;
2758 }
2759
2760 vol[CP] -= g33[k-1][j][i];
2761 vol[BP] += g33[k-1][j][i];
2762 }
2763
2764 // --- Final scaling and insertion into the matrix ---
2765
2766 // Scale all stencil coefficients by the negative cell volume (-aj).
2767 for (PetscInt m = 0; m < 19; m++) {
2768 vol[m] *= -aj[k][j][i];
2769 }
2770
2771 // Set the global column indices for the 19 stencil points, handling periodic BCs.
2772 idx[CP] = Gidx(i, j, k, user);
2773 if (user->bctype[0]==7 && i==mx-2) idx[EP] = Gidx(1, j, k, user); else idx[EP] = Gidx(i+1, j, k, user);
2774 if (user->bctype[0]==7 && i==1) idx[WP] = Gidx(mx-2, j, k, user); else idx[WP] = Gidx(i-1, j, k, user);
2775 if (user->bctype[2]==7 && j==my-2) idx[NP] = Gidx(i, 1, k, user); else idx[NP] = Gidx(i, j+1, k, user);
2776 if (user->bctype[2]==7 && j==1) idx[SP] = Gidx(i, my-2, k, user); else idx[SP] = Gidx(i, j-1, k, user);
2777 if (user->bctype[4]==7 && k==mz-2) idx[TP] = Gidx(i, j, 1, user); else idx[TP] = Gidx(i, j, k+1, user);
2778 if (user->bctype[4]==7 && k==1) idx[BP] = Gidx(i, j, mz-2, user); else idx[BP] = Gidx(i, j, k-1, user);
2779 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);
2780 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);
2781 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);
2782 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);
2783 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);
2784 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);
2785 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);
2786 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);
2787 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);
2788 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);
2789 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);
2790 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);
2791
2792 // Insert the computed row into the matrix A.
2793 MatSetValues(user->A, 1, &row, 19, idx, vol, INSERT_VALUES);
2794 }
2795 }
2796 }
2797 }
2798
2799 //================================================================================
2800 // Section 4: Finalize Matrix and Cleanup
2801 //================================================================================
2802
2803 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Finalizing matrix assembly.\n");
2804 MatAssemblyBegin(user->A, MAT_FINAL_ASSEMBLY);
2805 MatAssemblyEnd(user->A, MAT_FINAL_ASSEMBLY);
2806
2807 PetscReal max_A;
2808
2809 ierr = MatNorm(user->A,NORM_INFINITY,&max_A);CHKERRQ(ierr);
2810
2811 LOG_ALLOW(GLOBAL,LOG_DEBUG," Max value in A matrix for level %d = %le.\n",user->thislevel,max_A);
2812
2813 // if (get_log_level() >= LOG_DEBUG) {
2814 // ierr = MatView(user->A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
2815 // }
2816
2817 // --- Restore access to all PETSc vectors and destroy temporary ones ---
2818 DMDAVecRestoreArray(da, G11, &g11); DMDAVecRestoreArray(da, G12, &g12); DMDAVecRestoreArray(da, G13, &g13);
2819 DMDAVecRestoreArray(da, G21, &g21); DMDAVecRestoreArray(da, G22, &g22); DMDAVecRestoreArray(da, G23, &g23);
2820 DMDAVecRestoreArray(da, G31, &g31); DMDAVecRestoreArray(da, G32, &g32); DMDAVecRestoreArray(da, G33, &g33);
2821
2822 VecDestroy(&G11); VecDestroy(&G12); VecDestroy(&G13);
2823 VecDestroy(&G21); VecDestroy(&G22); VecDestroy(&G23);
2824 VecDestroy(&G31); VecDestroy(&G32); VecDestroy(&G33);
2825
2826 DMDAVecRestoreArray(fda, user->lCsi, &csi); DMDAVecRestoreArray(fda, user->lEta, &eta); DMDAVecRestoreArray(fda, user->lZet, &zet);
2827 DMDAVecRestoreArray(fda, user->lICsi, &icsi); DMDAVecRestoreArray(fda, user->lIEta, &ieta); DMDAVecRestoreArray(fda, user->lIZet, &izet);
2828 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi); DMDAVecRestoreArray(fda, user->lJEta, &jeta); DMDAVecRestoreArray(fda, user->lJZet, &jzet);
2829 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi); DMDAVecRestoreArray(fda, user->lKEta, &keta); DMDAVecRestoreArray(fda, user->lKZet, &kzet);
2830 DMDAVecRestoreArray(da, user->lAj, &aj); DMDAVecRestoreArray(da, user->lIAj, &iaj); DMDAVecRestoreArray(da, user->lJAj, &jaj); DMDAVecRestoreArray(da, user->lKAj, &kaj);
2831 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2832
2833 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting PoissonLHSNew.\n");
2834 PetscFunctionReturn(0);
2835}
2836
2837PetscErrorCode PoissonRHS(UserCtx *user, Vec B)
2838{
2839 DMDALocalInfo info = user->info;
2840 PetscInt xs = info.xs, xe = info.xs + info.xm;
2841 PetscInt ys = info.ys, ye = info.ys + info.ym;
2842 PetscInt zs = info.zs, ze = info.zs + info.zm;
2843 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2844
2845 PetscInt i, j, k;
2846 PetscReal ***nvert, ***aj, ***rb, dt = user->simCtx->dt;
2847 struct Components{
2848 PetscReal x;
2849 PetscReal y;
2850 PetscReal z;
2851 } *** ucont;
2852
2853 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonRHS to compute pressure equation RHS.\n");
2854
2855 DMDAVecGetArray(user->da, B, &rb);
2856 DMDAVecGetArray(user->fda, user->lUcont, &ucont);
2857 DMDAVecGetArray(user->da, user->lNvert, &nvert);
2858 DMDAVecGetArray(user->da, user->lAj, &aj);
2859
2860
2861 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing RHS values for each cell.\n");
2862
2863 for (k=zs; k<ze; k++) {
2864 for (j=ys; j<ye; j++) {
2865 for (i=xs; i<xe; i++) {
2866
2867 if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1) {
2868 rb[k][j][i] = 0.;
2869 }
2870 else if (nvert[k][j][i] > 0.1) {
2871 rb[k][j][i] = 0;
2872 }
2873 else {
2874 rb[k][j][i] = -(ucont[k][j][i].x - ucont[k][j][i-1].x +
2875 ucont[k][j][i].y - ucont[k][j-1][i].y +
2876 ucont[k][j][i].z - ucont[k-1][j][i].z) / dt
2877 * aj[k][j][i] / user->simCtx->st * COEF_TIME_ACCURACY;
2878
2879 }
2880 }
2881 }
2882 }
2883
2884
2885 // --- Check the solvability condition for the Poisson equation ---
2886 // The global sum of the RHS (proportional to the total divergence) must be zero.
2887 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying solvability condition (sum of RHS terms).\n");
2888 PetscReal lsum=0., sum=0.;
2889
2890 for (k=zs; k<ze; k++) {
2891 for (j=ys; j<ye; j++) {
2892 for (i=xs; i<xe; i++) {
2893
2894 lsum += rb[k][j][i] / aj[k][j][i]* dt/COEF_TIME_ACCURACY;
2895
2896 }
2897 }
2898 }
2899
2900 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2901
2902 LOG_ALLOW(GLOBAL, LOG_INFO, "Global Sum of RHS (Divergence Check): %le\n", sum);
2903
2904 user->simCtx->summationRHS = sum;
2905
2906 DMDAVecRestoreArray(user->fda, user->lUcont, &ucont);
2907 DMDAVecRestoreArray(user->da, user->lNvert, &nvert);
2908 DMDAVecRestoreArray(user->da, user->lAj, &aj);
2909 DMDAVecRestoreArray(user->da, B, &rb);
2910
2911 return 0;
2912}
2913
2914PetscErrorCode VolumeFlux_rev(UserCtx *user, PetscReal *ibm_Flux,
2915 PetscReal *ibm_Area, PetscInt flg)
2916{
2917
2918 // --- CONTEXT ACQUISITION BLOCK ---
2919 // Get the master simulation context from the UserCtx.
2920 SimCtx *simCtx = user->simCtx;
2921
2922 // Create a local variable to mirror the legacy global for minimal code changes.
2923 const PetscInt NumberOfBodies = simCtx->NumberOfBodies;
2924 // --- END CONTEXT ACQUISITION BLOCK ---
2925
2926 DM da = user->da, fda = user->fda;
2927
2928 DMDALocalInfo info = user->info;
2929
2930 PetscInt xs = info.xs, xe = info.xs + info.xm;
2931 PetscInt ys = info.ys, ye = info.ys + info.ym;
2932 PetscInt zs = info.zs, ze = info.zs + info.zm;
2933 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2934
2935 PetscInt i, j, k;
2936 PetscInt lxs, lys, lzs, lxe, lye, lze;
2937
2938 lxs = xs; lxe = xe;
2939 lys = ys; lye = ye;
2940 lzs = zs; lze = ze;
2941
2942 if (xs==0) lxs = xs+1;
2943 if (ys==0) lys = ys+1;
2944 if (zs==0) lzs = zs+1;
2945
2946 if (xe==mx) lxe = xe-1;
2947 if (ye==my) lye = ye-1;
2948 if (ze==mz) lze = ze-1;
2949
2950 PetscReal ***nvert, ibmval=1.5;
2951 Cmpnts ***ucor, ***csi, ***eta, ***zet;
2952 DMDAVecGetArray(fda, user->Ucont, &ucor);
2953 DMDAVecGetArray(fda, user->lCsi, &csi);
2954 DMDAVecGetArray(fda, user->lEta, &eta);
2955 DMDAVecGetArray(fda, user->lZet, &zet);
2956 DMDAVecGetArray(da, user->lNvert, &nvert);
2957
2958 PetscReal libm_Flux, libm_area;
2959 libm_Flux = 0;
2960 libm_area = 0;
2961 for (k=lzs; k<lze; k++) {
2962 for (j=lys; j<lye; j++) {
2963 for (i=lxs; i<lxe; i++) {
2964 if (nvert[k][j][i] < 0.1) {
2965 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
2966 libm_Flux += ucor[k][j][i].x;
2967 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2968 csi[k][j][i].y * csi[k][j][i].y +
2969 csi[k][j][i].z * csi[k][j][i].z);
2970
2971 }
2972 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
2973 libm_Flux += ucor[k][j][i].y;
2974 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2975 eta[k][j][i].y * eta[k][j][i].y +
2976 eta[k][j][i].z * eta[k][j][i].z);
2977 }
2978 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
2979 libm_Flux += ucor[k][j][i].z;
2980 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
2981 zet[k][j][i].y * zet[k][j][i].y +
2982 zet[k][j][i].z * zet[k][j][i].z);
2983 }
2984 }
2985
2986 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
2987 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
2988 libm_Flux -= ucor[k][j][i].x;
2989 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
2990 csi[k][j][i].y * csi[k][j][i].y +
2991 csi[k][j][i].z * csi[k][j][i].z);
2992
2993 }
2994 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
2995 libm_Flux -= ucor[k][j][i].y;
2996 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
2997 eta[k][j][i].y * eta[k][j][i].y +
2998 eta[k][j][i].z * eta[k][j][i].z);
2999 }
3000 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
3001 libm_Flux -= ucor[k][j][i].z;
3002 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3003 zet[k][j][i].y * zet[k][j][i].y +
3004 zet[k][j][i].z * zet[k][j][i].z);
3005 }
3006 }
3007
3008 }
3009 }
3010 }
3011
3012 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3013 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3014
3015 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
3016/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
3017 PetscPrintf(PETSC_COMM_WORLD, "IBMFlux REV %le %le\n", *ibm_Flux, *ibm_Area);
3018
3019 PetscReal correction;
3020
3021 if (*ibm_Area > 1.e-15) {
3022 if (flg)
3023 correction = (*ibm_Flux + user->FluxIntpSum) / *ibm_Area;
3024 else
3025 correction = *ibm_Flux / *ibm_Area;
3026 }
3027 else {
3028 correction = 0;
3029 }
3030
3031 for (k=lzs; k<lze; k++) {
3032 for (j=lys; j<lye; j++) {
3033 for (i=lxs; i<lxe; i++) {
3034 if (nvert[k][j][i] < 0.1) {
3035 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
3036 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
3037 csi[k][j][i].y * csi[k][j][i].y +
3038 csi[k][j][i].z * csi[k][j][i].z) *
3039 correction;
3040
3041 }
3042 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
3043 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
3044 eta[k][j][i].y * eta[k][j][i].y +
3045 eta[k][j][i].z * eta[k][j][i].z) *
3046 correction;
3047 }
3048 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
3049 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
3050 zet[k][j][i].y * zet[k][j][i].y +
3051 zet[k][j][i].z * zet[k][j][i].z) *
3052 correction;
3053 }
3054 }
3055
3056 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
3057 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
3058 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3059 csi[k][j][i].y * csi[k][j][i].y +
3060 csi[k][j][i].z * csi[k][j][i].z) *
3061 correction;
3062
3063 }
3064 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
3065 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3066 eta[k][j][i].y * eta[k][j][i].y +
3067 eta[k][j][i].z * eta[k][j][i].z) *
3068 correction;
3069 }
3070 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
3071 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3072 zet[k][j][i].y * zet[k][j][i].y +
3073 zet[k][j][i].z * zet[k][j][i].z) *
3074 correction;
3075 }
3076 }
3077
3078 }
3079 }
3080 }
3081
3082
3083
3084 libm_Flux = 0;
3085 libm_area = 0;
3086 for (k=lzs; k<lze; k++) {
3087 for (j=lys; j<lye; j++) {
3088 for (i=lxs; i<lxe; i++) {
3089 if (nvert[k][j][i] < 0.1) {
3090 if (nvert[k][j][i+1] > ibmval-0.4 && nvert[k][j][i+1] < ibmval && i < mx-2) {
3091 libm_Flux += ucor[k][j][i].x;
3092 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3093 csi[k][j][i].y * csi[k][j][i].y +
3094 csi[k][j][i].z * csi[k][j][i].z);
3095
3096 }
3097 if (nvert[k][j+1][i] > ibmval-0.4 && nvert[k][j+1][i] < ibmval && j < my-2) {
3098 libm_Flux += ucor[k][j][i].y;
3099 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3100 eta[k][j][i].y * eta[k][j][i].y +
3101 eta[k][j][i].z * eta[k][j][i].z);
3102 }
3103 if (nvert[k+1][j][i] > ibmval-0.4 && nvert[k+1][j][i] < ibmval && k < mz-2) {
3104 libm_Flux += ucor[k][j][i].z;
3105 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3106 zet[k][j][i].y * zet[k][j][i].y +
3107 zet[k][j][i].z * zet[k][j][i].z);
3108 }
3109 }
3110
3111 if (nvert[k][j][i] > ibmval-0.4 && nvert[k][j][i] < ibmval) {
3112 if (nvert[k][j][i+1] < 0.1 && i < mx-2) {
3113 libm_Flux -= ucor[k][j][i].x;
3114 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3115 csi[k][j][i].y * csi[k][j][i].y +
3116 csi[k][j][i].z * csi[k][j][i].z);
3117
3118 }
3119 if (nvert[k][j+1][i] < 0.1 && j < my-2) {
3120 libm_Flux -= ucor[k][j][i].y;
3121 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3122 eta[k][j][i].y * eta[k][j][i].y +
3123 eta[k][j][i].z * eta[k][j][i].z);
3124 }
3125 if (nvert[k+1][j][i] < 0.1 && k < mz-2) {
3126 libm_Flux -= ucor[k][j][i].z;
3127 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3128 zet[k][j][i].y * zet[k][j][i].y +
3129 zet[k][j][i].z * zet[k][j][i].z);
3130 }
3131 }
3132
3133 }
3134 }
3135 }
3136
3137 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3138 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3139
3140 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
3141/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
3142 PetscPrintf(PETSC_COMM_WORLD, "IBMFlux22 REV %le %le\n", *ibm_Flux, *ibm_Area);
3143
3144 DMDAVecRestoreArray(da, user->lNvert, &nvert);
3145 DMDAVecRestoreArray(fda, user->lCsi, &csi);
3146 DMDAVecRestoreArray(fda, user->lEta, &eta);
3147 DMDAVecRestoreArray(fda, user->lZet, &zet);
3148 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
3149
3150 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3151 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3152 return 0;
3153}
3154
3155
3156PetscErrorCode VolumeFlux(UserCtx *user, PetscReal *ibm_Flux, PetscReal *ibm_Area, PetscInt flg)
3157{
3158 // --- CONTEXT ACQUISITION BLOCK ---
3159 // Get the master simulation context from the UserCtx.
3160 SimCtx *simCtx = user->simCtx;
3161
3162 // Create local variables to mirror the legacy globals for minimal code changes.
3163 const PetscInt NumberOfBodies = simCtx->NumberOfBodies;
3164 const PetscInt channelz = simCtx->channelz;
3165 const PetscInt fish = simCtx->fish;
3166 // --- END CONTEXT ACQUISITION BLOCK ---
3167
3168 DM da = user->da, fda = user->fda;
3169
3170 DMDALocalInfo info = user->info;
3171
3172 PetscInt xs = info.xs, xe = info.xs + info.xm;
3173 PetscInt ys = info.ys, ye = info.ys + info.ym;
3174 PetscInt zs = info.zs, ze = info.zs + info.zm;
3175 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3176
3177 PetscInt i, j, k,ibi;
3178 PetscInt lxs, lys, lzs, lxe, lye, lze;
3179
3180 PetscInt gxs, gxe, gys, gye, gzs, gze;
3181
3182 gxs = info.gxs; gxe = gxs + info.gxm;
3183 gys = info.gys; gye = gys + info.gym;
3184 gzs = info.gzs; gze = gzs + info.gzm;
3185
3186 lxs = xs; lxe = xe;
3187 lys = ys; lye = ye;
3188 lzs = zs; lze = ze;
3189
3190 if (xs==0) lxs = xs+1;
3191 if (ys==0) lys = ys+1;
3192 if (zs==0) lzs = zs+1;
3193
3194 if (xe==mx) lxe = xe-1;
3195 if (ye==my) lye = ye-1;
3196 if (ze==mz) lze = ze-1;
3197
3198 PetscReal epsilon=1.e-8;
3199 PetscReal ***nvert, ibmval=1.9999;
3200
3201 struct Components {
3202 PetscReal x;
3203 PetscReal y;
3204 PetscReal z;
3205 }***ucor, ***lucor,***csi, ***eta, ***zet;
3206
3207
3208 PetscInt xend=mx-2 ,yend=my-2,zend=mz-2;
3209
3210 if (user->bctype[0]==7) xend=mx-1;
3211 if (user->bctype[2]==7) yend=my-1;
3212 if (user->bctype[4]==7) zend=mz-1;
3213
3214 DMDAVecGetArray(fda, user->Ucont, &ucor);
3215 DMDAVecGetArray(fda, user->lCsi, &csi);
3216 DMDAVecGetArray(fda, user->lEta, &eta);
3217 DMDAVecGetArray(fda, user->lZet, &zet);
3218 DMDAVecGetArray(da, user->lNvert, &nvert);
3219
3220 PetscReal libm_Flux, libm_area, libm_Flux_abs=0., ibm_Flux_abs;
3221 libm_Flux = 0;
3222 libm_area = 0;
3223
3224 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering VolumeFlux to enforce no-penetration condition.\n");
3225
3226 //Mohsen March 2017
3227 PetscReal *lIB_Flux, *lIB_area,*IB_Flux,*IB_Area;
3228 if (NumberOfBodies > 1) {
3229
3230 lIB_Flux=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3231 lIB_area=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3232 IB_Flux=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3233 IB_Area=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3234
3235
3236 for (ibi=0; ibi<NumberOfBodies; ibi++) {
3237 lIB_Flux[ibi]=0.0;
3238 lIB_area[ibi]=0.0;
3239 IB_Flux[ibi]=0.0;
3240 IB_Area[ibi]=0.0;
3241 }
3242 }
3243
3244
3245 //================================================================================
3246 // PASS 1: Calculate Uncorrected Flux and Area
3247 // This pass measures the total fluid "leakage" across the immersed boundary.
3248 //================================================================================
3249 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 1: Measuring uncorrected flux and area.\n");
3250
3251 for (k=lzs; k<lze; k++) {
3252 for (j=lys; j<lye; j++) {
3253 for (i=lxs; i<lxe; i++) {
3254 if (nvert[k][j][i] < 0.1) {
3255 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] < ibmval && i < xend) {
3256
3257 if (fabs(ucor[k][j][i].x)>epsilon) {
3258 libm_Flux += ucor[k][j][i].x;
3259 if (flg==3)
3260 libm_Flux_abs += fabs(ucor[k][j][i].x)/sqrt(csi[k][j][i].x * csi[k][j][i].x +
3261 csi[k][j][i].y * csi[k][j][i].y +
3262 csi[k][j][i].z * csi[k][j][i].z);
3263 else
3264 libm_Flux_abs += fabs(ucor[k][j][i].x);
3265
3266 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3267 csi[k][j][i].y * csi[k][j][i].y +
3268 csi[k][j][i].z * csi[k][j][i].z);
3269
3270 if (NumberOfBodies > 1) {
3271
3272 ibi=(int)((nvert[k][j][i+1]-1.0)*1001);
3273 lIB_Flux[ibi] += ucor[k][j][i].x;
3274 lIB_area[ibi] += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3275 csi[k][j][i].y * csi[k][j][i].y +
3276 csi[k][j][i].z * csi[k][j][i].z);
3277 }
3278 } else
3279 ucor[k][j][i].x=0.;
3280
3281 }
3282 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
3283
3284 if (fabs(ucor[k][j][i].y)>epsilon) {
3285 libm_Flux += ucor[k][j][i].y;
3286 if (flg==3)
3287 libm_Flux_abs += fabs(ucor[k][j][i].y)/sqrt(eta[k][j][i].x * eta[k][j][i].x +
3288 eta[k][j][i].y * eta[k][j][i].y +
3289 eta[k][j][i].z * eta[k][j][i].z);
3290 else
3291 libm_Flux_abs += fabs(ucor[k][j][i].y);
3292 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3293 eta[k][j][i].y * eta[k][j][i].y +
3294 eta[k][j][i].z * eta[k][j][i].z);
3295 if (NumberOfBodies > 1) {
3296
3297 ibi=(int)((nvert[k][j+1][i]-1.0)*1001);
3298
3299 lIB_Flux[ibi] += ucor[k][j][i].y;
3300 lIB_area[ibi] += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3301 eta[k][j][i].y * eta[k][j][i].y +
3302 eta[k][j][i].z * eta[k][j][i].z);
3303 }
3304 } else
3305 ucor[k][j][i].y=0.;
3306 }
3307 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
3308
3309 if (fabs(ucor[k][j][i].z)>epsilon) {
3310 libm_Flux += ucor[k][j][i].z;
3311 if (flg==3)
3312 libm_Flux_abs += fabs(ucor[k][j][i].z)/sqrt(zet[k][j][i].x * zet[k][j][i].x +
3313 zet[k][j][i].y * zet[k][j][i].y +
3314 zet[k][j][i].z * zet[k][j][i].z);
3315 else
3316 libm_Flux_abs += fabs(ucor[k][j][i].z);
3317 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3318 zet[k][j][i].y * zet[k][j][i].y +
3319 zet[k][j][i].z * zet[k][j][i].z);
3320
3321 if (NumberOfBodies > 1) {
3322
3323 ibi=(int)((nvert[k+1][j][i]-1.0)*1001);
3324 lIB_Flux[ibi] += ucor[k][j][i].z;
3325 lIB_area[ibi] += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3326 zet[k][j][i].y * zet[k][j][i].y +
3327 zet[k][j][i].z * zet[k][j][i].z);
3328 }
3329 }else
3330 ucor[k][j][i].z=0.;
3331 }
3332 }
3333
3334 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
3335
3336 if (nvert[k][j][i+1] < 0.1 && i < xend) {
3337 if (fabs(ucor[k][j][i].x)>epsilon) {
3338 libm_Flux -= ucor[k][j][i].x;
3339 if (flg==3)
3340 libm_Flux_abs += fabs(ucor[k][j][i].x)/sqrt(csi[k][j][i].x * csi[k][j][i].x +
3341 csi[k][j][i].y * csi[k][j][i].y +
3342 csi[k][j][i].z * csi[k][j][i].z);
3343 else
3344 libm_Flux_abs += fabs(ucor[k][j][i].x);
3345 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3346 csi[k][j][i].y * csi[k][j][i].y +
3347 csi[k][j][i].z * csi[k][j][i].z);
3348 if (NumberOfBodies > 1) {
3349
3350 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3351 lIB_Flux[ibi] -= ucor[k][j][i].x;
3352 lIB_area[ibi] += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3353 csi[k][j][i].y * csi[k][j][i].y +
3354 csi[k][j][i].z * csi[k][j][i].z);
3355 }
3356
3357 }else
3358 ucor[k][j][i].x=0.;
3359 }
3360 if (nvert[k][j+1][i] < 0.1 && j < yend) {
3361 if (fabs(ucor[k][j][i].y)>epsilon) {
3362 libm_Flux -= ucor[k][j][i].y;
3363 if (flg==3)
3364 libm_Flux_abs += fabs(ucor[k][j][i].y)/ sqrt(eta[k][j][i].x * eta[k][j][i].x +
3365 eta[k][j][i].y * eta[k][j][i].y +
3366 eta[k][j][i].z * eta[k][j][i].z);
3367 else
3368 libm_Flux_abs += fabs(ucor[k][j][i].y);
3369 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3370 eta[k][j][i].y * eta[k][j][i].y +
3371 eta[k][j][i].z * eta[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].y;
3376 lIB_area[ibi] += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3377 eta[k][j][i].y * eta[k][j][i].y +
3378 eta[k][j][i].z * eta[k][j][i].z);
3379 }
3380 }else
3381 ucor[k][j][i].y=0.;
3382 }
3383 if (nvert[k+1][j][i] < 0.1 && k < zend) {
3384 if (fabs(ucor[k][j][i].z)>epsilon) {
3385 libm_Flux -= ucor[k][j][i].z;
3386 if (flg==3)
3387 libm_Flux_abs += fabs(ucor[k][j][i].z)/sqrt(zet[k][j][i].x * zet[k][j][i].x +
3388 zet[k][j][i].y * zet[k][j][i].y +
3389 zet[k][j][i].z * zet[k][j][i].z);
3390 else
3391 libm_Flux_abs += fabs(ucor[k][j][i].z);
3392 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3393 zet[k][j][i].y * zet[k][j][i].y +
3394 zet[k][j][i].z * zet[k][j][i].z);
3395 if (NumberOfBodies > 1) {
3396
3397 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3398 lIB_Flux[ibi] -= ucor[k][j][i].z;
3399 lIB_area[ibi] += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3400 zet[k][j][i].y * zet[k][j][i].y +
3401 zet[k][j][i].z * zet[k][j][i].z);
3402 }
3403 }else
3404 ucor[k][j][i].z=0.;
3405 }
3406 }
3407
3408 }
3409 }
3410 }
3411
3412 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3413 MPI_Allreduce(&libm_Flux_abs, &ibm_Flux_abs,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3414 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3415
3416 if (NumberOfBodies > 1) {
3417 MPI_Allreduce(lIB_Flux,IB_Flux,NumberOfBodies,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
3418 MPI_Allreduce(lIB_area,IB_Area,NumberOfBodies,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
3419 }
3420
3421 PetscReal correction;
3422
3423 PetscReal *Correction;
3424 if (NumberOfBodies > 1) {
3425 Correction=(PetscReal *)calloc(NumberOfBodies,sizeof(PetscReal));
3426 for (ibi=0; ibi<NumberOfBodies; ibi++) Correction[ibi]=0.0;
3427 }
3428
3429 if (*ibm_Area > 1.e-15) {
3430 if (flg>1)
3431 correction = (*ibm_Flux + user->FluxIntpSum)/ ibm_Flux_abs;
3432 else if (flg)
3433 correction = (*ibm_Flux + user->FluxIntpSum) / *ibm_Area;
3434 else
3435 correction = *ibm_Flux / *ibm_Area;
3436 if (NumberOfBodies > 1)
3437 for (ibi=0; ibi<NumberOfBodies; ibi++) if (IB_Area[ibi]>1.e-15) Correction[ibi] = IB_Flux[ibi] / IB_Area[ibi];
3438 }
3439 else {
3440 correction = 0;
3441 }
3442 // --- Log the uncorrected results and calculated correction ---
3443 LOG_ALLOW(GLOBAL, LOG_INFO, "IBM Uncorrected Flux: %g, Area: %g, Correction: %g\n", *ibm_Flux, *ibm_Area, correction);
3444 if (NumberOfBodies>1){
3445 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]);
3446 }
3447
3448 //================================================================================
3449 // PASS 2: Apply Correction to Velocity Field
3450 // This pass modifies the velocity at the boundary to enforce zero net flux.
3451 //================================================================================
3452 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 2: Applying velocity corrections at the boundary.\n");
3453
3454 for (k=lzs; k<lze; k++) {
3455 for (j=lys; j<lye; j++) {
3456 for (i=lxs; i<lxe; i++) {
3457 if (nvert[k][j][i] < 0.1) {
3458 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] <ibmval && i < xend) {
3459 if (fabs(ucor[k][j][i].x)>epsilon){
3460 if (flg==3)
3461 ucor[k][j][i].x -=correction*fabs(ucor[k][j][i].x)/
3462 sqrt(csi[k][j][i].x * csi[k][j][i].x +
3463 csi[k][j][i].y * csi[k][j][i].y +
3464 csi[k][j][i].z * csi[k][j][i].z);
3465 else if (flg==2)
3466 ucor[k][j][i].x -=correction*fabs(ucor[k][j][i].x);
3467 else if (NumberOfBodies > 1) {
3468 ibi=(int)((nvert[k][j][i+1]-1.0)*1001);
3469 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
3470 csi[k][j][i].y * csi[k][j][i].y +
3471 csi[k][j][i].z * csi[k][j][i].z) *
3472 Correction[ibi];
3473 }
3474 else
3475 ucor[k][j][i].x -= sqrt(csi[k][j][i].x * csi[k][j][i].x +
3476 csi[k][j][i].y * csi[k][j][i].y +
3477 csi[k][j][i].z * csi[k][j][i].z) *
3478 correction;
3479 }
3480 }
3481 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
3482 if (fabs(ucor[k][j][i].y)>epsilon) {
3483 if (flg==3)
3484 ucor[k][j][i].y -=correction*fabs(ucor[k][j][i].y)/
3485 sqrt(eta[k][j][i].x * eta[k][j][i].x +
3486 eta[k][j][i].y * eta[k][j][i].y +
3487 eta[k][j][i].z * eta[k][j][i].z);
3488 else if (flg==2)
3489 ucor[k][j][i].y -=correction*fabs(ucor[k][j][i].y);
3490 else if (NumberOfBodies > 1) {
3491 ibi=(int)((nvert[k][j+1][i]-1.0)*1001);
3492 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
3493 eta[k][j][i].y * eta[k][j][i].y +
3494 eta[k][j][i].z * eta[k][j][i].z) *
3495 Correction[ibi];
3496 }
3497 else
3498 ucor[k][j][i].y -= sqrt(eta[k][j][i].x * eta[k][j][i].x +
3499 eta[k][j][i].y * eta[k][j][i].y +
3500 eta[k][j][i].z * eta[k][j][i].z) *
3501 correction;
3502 }
3503 }
3504 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
3505 if (fabs(ucor[k][j][i].z)>epsilon) {
3506 if (flg==3)
3507 ucor[k][j][i].z -= correction*fabs(ucor[k][j][i].z)/
3508 sqrt(zet[k][j][i].x * zet[k][j][i].x +
3509 zet[k][j][i].y * zet[k][j][i].y +
3510 zet[k][j][i].z * zet[k][j][i].z);
3511 else if (flg==2)
3512 ucor[k][j][i].z -= correction*fabs(ucor[k][j][i].z);
3513 else if (NumberOfBodies > 1) {
3514 ibi=(int)((nvert[k+1][j][i]-1.0)*1001);
3515 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
3516 zet[k][j][i].y * zet[k][j][i].y +
3517 zet[k][j][i].z * zet[k][j][i].z) *
3518 Correction[ibi];
3519 }
3520 else
3521 ucor[k][j][i].z -= sqrt(zet[k][j][i].x * zet[k][j][i].x +
3522 zet[k][j][i].y * zet[k][j][i].y +
3523 zet[k][j][i].z * zet[k][j][i].z) *
3524 correction;
3525 }
3526 }
3527 }
3528
3529 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
3530 if (nvert[k][j][i+1] < 0.1 && i < xend) {
3531 if (fabs(ucor[k][j][i].x)>epsilon) {
3532 if (flg==3)
3533 ucor[k][j][i].x += correction*fabs(ucor[k][j][i].x)/
3534 sqrt(csi[k][j][i].x * csi[k][j][i].x +
3535 csi[k][j][i].y * csi[k][j][i].y +
3536 csi[k][j][i].z * csi[k][j][i].z);
3537 else if (flg==2)
3538 ucor[k][j][i].x += correction*fabs(ucor[k][j][i].x);
3539 else if (NumberOfBodies > 1) {
3540 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3541 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3542 csi[k][j][i].y * csi[k][j][i].y +
3543 csi[k][j][i].z * csi[k][j][i].z) *
3544 Correction[ibi];
3545 }
3546 else
3547 ucor[k][j][i].x += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3548 csi[k][j][i].y * csi[k][j][i].y +
3549 csi[k][j][i].z * csi[k][j][i].z) *
3550 correction;
3551 }
3552 }
3553 if (nvert[k][j+1][i] < 0.1 && j < yend) {
3554 if (fabs(ucor[k][j][i].y)>epsilon) {
3555 if (flg==3)
3556 ucor[k][j][i].y +=correction*fabs(ucor[k][j][i].y)/
3557 sqrt(eta[k][j][i].x * eta[k][j][i].x +
3558 eta[k][j][i].y * eta[k][j][i].y +
3559 eta[k][j][i].z * eta[k][j][i].z);
3560 else if (flg==2)
3561 ucor[k][j][i].y +=correction*fabs(ucor[k][j][i].y);
3562 else if (NumberOfBodies > 1) {
3563 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3564 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3565 eta[k][j][i].y * eta[k][j][i].y +
3566 eta[k][j][i].z * eta[k][j][i].z) *
3567 Correction[ibi];
3568 }
3569 else
3570 ucor[k][j][i].y += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3571 eta[k][j][i].y * eta[k][j][i].y +
3572 eta[k][j][i].z * eta[k][j][i].z) *
3573 correction;
3574 }
3575 }
3576 if (nvert[k+1][j][i] < 0.1 && k < zend) {
3577 if (fabs(ucor[k][j][i].z)>epsilon) {
3578 if (flg==3)
3579 ucor[k][j][i].z += correction*fabs(ucor[k][j][i].z)/
3580 sqrt(zet[k][j][i].x * zet[k][j][i].x +
3581 zet[k][j][i].y * zet[k][j][i].y +
3582 zet[k][j][i].z * zet[k][j][i].z);
3583 else if (flg==2)
3584 ucor[k][j][i].z += correction*fabs(ucor[k][j][i].z);
3585 else if (NumberOfBodies > 1) {
3586 ibi=(int)((nvert[k][j][i]-1.0)*1001);
3587 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3588 zet[k][j][i].y * zet[k][j][i].y +
3589 zet[k][j][i].z * zet[k][j][i].z) *
3590 Correction[ibi];
3591 }
3592 else
3593 ucor[k][j][i].z += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3594 zet[k][j][i].y * zet[k][j][i].y +
3595 zet[k][j][i].z * zet[k][j][i].z) *
3596 correction;
3597 }
3598 }
3599 }
3600
3601 }
3602 }
3603 }
3604
3605 //================================================================================
3606 // PASS 3: Verification
3607 // This optional pass recalculates the flux to confirm the correction was successful.
3608 //================================================================================
3609 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Pass 3: Verifying corrected flux.\n");
3610
3611 libm_Flux = 0;
3612 libm_area = 0;
3613 for (k=lzs; k<lze; k++) {
3614 for (j=lys; j<lye; j++) {
3615 for (i=lxs; i<lxe; i++) {
3616 if (nvert[k][j][i] < 0.1) {
3617 if (nvert[k][j][i+1] > 0.1 && nvert[k][j][i+1] < ibmval && i < xend) {
3618 libm_Flux += ucor[k][j][i].x;
3619 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3620 csi[k][j][i].y * csi[k][j][i].y +
3621 csi[k][j][i].z * csi[k][j][i].z);
3622
3623 }
3624 if (nvert[k][j+1][i] > 0.1 && nvert[k][j+1][i] < ibmval && j < yend) {
3625 libm_Flux += ucor[k][j][i].y;
3626 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3627 eta[k][j][i].y * eta[k][j][i].y +
3628 eta[k][j][i].z * eta[k][j][i].z);
3629 }
3630 if (nvert[k+1][j][i] > 0.1 && nvert[k+1][j][i] < ibmval && k < zend) {
3631 libm_Flux += ucor[k][j][i].z;
3632 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3633 zet[k][j][i].y * zet[k][j][i].y +
3634 zet[k][j][i].z * zet[k][j][i].z);
3635 }
3636 }
3637
3638 if (nvert[k][j][i] > 0.1 && nvert[k][j][i] < ibmval) {
3639 if (nvert[k][j][i+1] < 0.1 && i < xend) {
3640 libm_Flux -= ucor[k][j][i].x;
3641 libm_area += sqrt(csi[k][j][i].x * csi[k][j][i].x +
3642 csi[k][j][i].y * csi[k][j][i].y +
3643 csi[k][j][i].z * csi[k][j][i].z);
3644
3645 }
3646 if (nvert[k][j+1][i] < 0.1 && j < yend) {
3647 libm_Flux -= ucor[k][j][i].y;
3648 libm_area += sqrt(eta[k][j][i].x * eta[k][j][i].x +
3649 eta[k][j][i].y * eta[k][j][i].y +
3650 eta[k][j][i].z * eta[k][j][i].z);
3651 }
3652 if (nvert[k+1][j][i] < 0.1 && k < zend) {
3653 libm_Flux -= ucor[k][j][i].z;
3654 libm_area += sqrt(zet[k][j][i].x * zet[k][j][i].x +
3655 zet[k][j][i].y * zet[k][j][i].y +
3656 zet[k][j][i].z * zet[k][j][i].z);
3657 }
3658 }
3659
3660 }
3661 }
3662 }
3663
3664 MPI_Allreduce(&libm_Flux, ibm_Flux,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3665 MPI_Allreduce(&libm_area, ibm_Area,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
3666
3667 /* PetscGlobalSum(&libm_Flux, ibm_Flux, PETSC_COMM_WORLD); */
3668/* PetscGlobalSum(&libm_area, ibm_Area, PETSC_COMM_WORLD); */
3669 LOG_ALLOW(GLOBAL, LOG_INFO, "IBM Corrected (Verified) Flux: %g, Area: %g\n", *ibm_Flux, *ibm_Area);
3670
3671
3672 if (user->bctype[0]==7 || user->bctype[1]==7){
3673 if (xe==mx){
3674 i=mx-2;
3675 for (k=lzs; k<lze; k++) {
3676 for (j=lys; j<lye; j++) {
3677 // if(j>0 && k>0 && j<user->JM && k<user->KM){
3678 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;
3679
3680 // }
3681 }
3682 }
3683 }
3684 }
3685
3686 if (user->bctype[2]==7 || user->bctype[3]==7){
3687 if (ye==my){
3688 j=my-2;
3689 for (k=lzs; k<lze; k++) {
3690 for (i=lxs; i<lxe; i++) {
3691 // if(i>0 && k>0 && i<user->IM && k<user->KM){
3692 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;
3693 // }
3694 }
3695 }
3696 }
3697 }
3698
3699 if (user->bctype[4]==7 || user->bctype[5]==7){
3700 if (ze==mz){
3701 k=mz-2;
3702 for (j=lys; j<lye; j++) {
3703 for (i=lxs; i<lxe; i++) {
3704 // if(i>0 && j>0 && i<user->IM && j<user->JM){
3705 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;
3706 // }
3707 }
3708 }
3709 }
3710 }
3711
3712
3713 DMDAVecRestoreArray(da, user->lNvert, &nvert);
3714 DMDAVecRestoreArray(fda, user->lCsi, &csi);
3715 DMDAVecRestoreArray(fda, user->lEta, &eta);
3716 DMDAVecRestoreArray(fda, user->lZet, &zet);
3717 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
3718
3719 DMGlobalToLocalBegin(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3720 DMGlobalToLocalEnd(fda, user->Ucont, INSERT_VALUES, user->lUcont);
3721
3722 /* periodci boundary condition update corrected flux */
3723 //Mohsen Dec 2015
3724 DMDAVecGetArray(fda, user->lUcont, &lucor);
3725 DMDAVecGetArray(fda, user->Ucont, &ucor);
3726
3727 if (user->bctype[0]==7 || user->bctype[1]==7){
3728 if (xs==0){
3729 i=xs;
3730 for (k=zs; k<ze; k++) {
3731 for (j=ys; j<ye; j++) {
3732 if(j>0 && k>0 && j<user->JM && k<user->KM){
3733 ucor[k][j][i].x=lucor[k][j][i-2].x;
3734 }
3735 }
3736 }
3737 }
3738 }
3739 if (user->bctype[2]==7 || user->bctype[3]==7){
3740 if (ys==0){
3741 j=ys;
3742 for (k=zs; k<ze; k++) {
3743 for (i=xs; i<xe; i++) {
3744 if(i>0 && k>0 && i<user->IM && k<user->KM){
3745 ucor[k][j][i].y=lucor[k][j-2][i].y;
3746 }
3747 }
3748 }
3749 }
3750 }
3751 if (user->bctype[4]==7 || user->bctype[5]==7){
3752 if (zs==0){
3753 k=zs;
3754 for (j=ys; j<ye; j++) {
3755 for (i=xs; i<xe; i++) {
3756 if(i>0 && j>0 && i<user->IM && j<user->JM){
3757 ucor[k][j][i].z=lucor[k-2][j][i].z;
3758 }
3759 }
3760 }
3761 }
3762 }
3763 DMDAVecRestoreArray(fda, user->lUcont, &lucor);
3764 DMDAVecRestoreArray(fda, user->Ucont, &ucor);
3765
3766 /* DMLocalToGlobalBegin(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); */
3767/* DMLocalToGlobalEnd(user->fda, user->lUcont, INSERT_VALUES, user->Ucont); */
3768 DMGlobalToLocalBegin(user->fda, user->Ucont, INSERT_VALUES, user->lUcont);
3769 DMGlobalToLocalEnd(user->fda, user->Ucont, INSERT_VALUES, user->lUcont);
3770
3771 if (NumberOfBodies > 1) {
3772 free(lIB_Flux);
3773 free(lIB_area);
3774 free(IB_Flux);
3775 free(IB_Area);
3776 free(Correction);
3777 }
3778
3779 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Exiting VolumeFlux.\n");
3780
3781 return 0;
3782}
3783
3784PetscErrorCode FullyBlocked(UserCtx *user)
3785{
3786 DM da = user->da;
3787 Vec nNvert;
3788 DMDALocalInfo info = user->info;
3789
3790 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3791
3792 PetscInt i, j, k;
3793
3794 PetscInt *KSKE = user->KSKE;
3795 PetscReal ***nvert;
3796 PetscBool *Blocked;
3797
3798 DMDACreateNaturalVector(da, &nNvert);
3799 DMDAGlobalToNaturalBegin(da, user->Nvert, INSERT_VALUES, nNvert);
3800 DMDAGlobalToNaturalEnd(da, user->Nvert, INSERT_VALUES, nNvert);
3801
3802 VecScatter ctx;
3803 Vec Zvert;
3804 VecScatterCreateToZero(nNvert, &ctx, &Zvert);
3805
3806 VecScatterBegin(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
3807 VecScatterEnd(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
3808
3809 VecScatterDestroy(&ctx);
3810 VecDestroy(&nNvert);
3811
3812 PetscInt rank;
3813 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
3814
3815 if (!rank) {
3816
3817 VecGetArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
3818 PetscMalloc(mx*my*sizeof(PetscBool), &Blocked);
3819 for (j=1; j<my-1; j++) {
3820 for (i=1; i<mx-1; i++) {
3821 Blocked[j*mx+i] = PETSC_FALSE;
3822 for (k=0; k<mz; k++) {
3823 if (nvert[k][j][i] > 0.1) {
3824 if (!Blocked[j*mx+i]) {
3825 KSKE[2*(j*mx+i)] = k;
3826 Blocked[j*mx+i] = PETSC_TRUE;
3827 }
3828 else {
3829 KSKE[2*(j*mx+i)] = PetscMin(KSKE[2*(j*mx+i)], k);
3830 }
3831 }
3832 }
3833 }
3834 }
3835
3836
3837 user->multinullspace = PETSC_TRUE;
3838 for (j=1; j<my-1; j++) {
3839 for (i=1; i<mx-1; i++) {
3840 if (!Blocked[j*mx+i]) {
3841 user->multinullspace = PETSC_FALSE;
3842 break;
3843 }
3844 }
3845 }
3846 PetscFree(Blocked);
3847 VecRestoreArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
3848 MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD);
3849 if (user->multinullspace) {
3850 MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD);
3851
3852 }
3853 }
3854 else {
3855 MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD);
3856 if (user->multinullspace) {
3857 MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD);
3858 }
3859 }
3860
3861
3862
3863 VecDestroy(&Zvert);
3864 return 0;
3865}
3866
3867PetscErrorCode MyNvertRestriction(UserCtx *user_h, UserCtx *user_c)
3868{
3869 // DA da = user_c->da, fda = user_c->fda;
3870
3871
3872
3873 DMDALocalInfo info = user_c->info;
3874 PetscInt xs = info.xs, xe = info.xs + info.xm;
3875 PetscInt ys = info.ys, ye = info.ys + info.ym;
3876 PetscInt zs = info.zs, ze = info.zs + info.zm;
3877 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3878
3879 PetscInt i,j,k;
3880 PetscInt ih, jh, kh, ia, ja, ka;
3881 PetscInt lxs, lxe, lys, lye, lzs, lze;
3882
3883 PetscReal ***nvert, ***nvert_h;
3884
3885 DMDAVecGetArray(user_h->da, user_h->lNvert, &nvert_h);
3886 DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert);
3887
3888 lxs = xs; lxe = xe;
3889 lys = ys; lye = ye;
3890 lzs = zs; lze = ze;
3891
3892 if (xs==0) lxs = xs+1;
3893 if (ys==0) lys = ys+1;
3894 if (zs==0) lzs = zs+1;
3895
3896 if (xe==mx) lxe = xe-1;
3897 if (ye==my) lye = ye-1;
3898 if (ze==mz) lze = ze-1;
3899
3900 if ((user_c->isc)) ia = 0;
3901 else ia = 1;
3902
3903 if ((user_c->jsc)) ja = 0;
3904 else ja = 1;
3905
3906 if ((user_c->ksc)) ka = 0;
3907 else ka = 1;
3908
3909 VecSet(user_c->Nvert, 0.);
3910 if (user_c->thislevel > 0) {
3911 for (k=lzs; k<lze; k++) {
3912 for (j=lys; j<lye; j++) {
3913 for (i=lxs; i<lxe; i++) {
3914 GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
3915 if (nvert_h[kh ][jh ][ih ] *
3916 nvert_h[kh ][jh ][ih-ia] *
3917 nvert_h[kh ][jh-ja][ih ] *
3918 nvert_h[kh-ka][jh ][ih ] *
3919 nvert_h[kh ][jh-ja][ih-ia] *
3920 nvert_h[kh-ka][jh ][ih-ia] *
3921 nvert_h[kh-ka][jh-ja][ih ] *
3922 nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
3923 nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
3924 }
3925 }
3926 }
3927 }
3928 }
3929 else {
3930 for (k=lzs; k<lze; k++) {
3931 for (j=lys; j<lye; j++) {
3932 for (i=lxs; i<lxe; i++) {
3933 GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
3934 if (nvert_h[kh ][jh ][ih ] *
3935 nvert_h[kh ][jh ][ih-ia] *
3936 nvert_h[kh ][jh-ja][ih ] *
3937 nvert_h[kh-ka][jh ][ih ] *
3938 nvert_h[kh ][jh-ja][ih-ia] *
3939 nvert_h[kh-ka][jh ][ih-ia] *
3940 nvert_h[kh-ka][jh-ja][ih ] *
3941 nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
3942 nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
3943 }
3944 }
3945 }
3946 }
3947 }
3948 DMDAVecRestoreArray(user_h->da, user_h->lNvert, &nvert_h);
3949 DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert);
3950
3951 DMGlobalToLocalBegin(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3952 DMGlobalToLocalEnd(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3953 //Mohsen Dec 2015
3954 DMDAVecGetArray(user_c->da, user_c->lNvert, &nvert);
3955 DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert_h);
3956
3957 for (k=lzs; k<lze; k++) {
3958 for (j=lys; j<lye; j++) {
3959 for (i=lxs; i<lxe; i++) {
3960 if (nvert_h[k][j][i] < 0.1) {
3961 if (nvert[k][j][i+1] + nvert[k][j][i-1] > 1.1 &&
3962 nvert[k][j+1][i] + nvert[k][j-1][i] > 1.1 &&
3963 nvert[k+1][j][i] + nvert[k-1][j][i] > 1.1) {
3964 nvert_h[k][j][i] = 1.;
3965 }
3966 }
3967 }
3968 }
3969 }
3970
3971 DMDAVecRestoreArray(user_c->da, user_c->lNvert, &nvert);
3972 DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert_h);
3973 DMGlobalToLocalBegin(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3974 DMGlobalToLocalEnd(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3975 /* DMLocalToGlobalBegin(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
3976/* DMLocalToGlobalEnd(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
3977 return 0;
3978}
3979
3980#undef __FUNCT__
3981#define __FUNCT__ "PoissonSolver_MG"
3982PetscErrorCode PoissonSolver_MG(UserMG *usermg)
3983{
3984 // --- CONTEXT ACQUISITION BLOCK ---
3985 // Get the master simulation context from the first block's UserCtx on the finest level.
3986 // This provides access to all former global variables.
3987 SimCtx *simCtx = usermg->mgctx[0].user[0].simCtx;
3988
3989 // Create local variables to mirror the legacy globals for minimal code changes.
3990 const PetscInt block_number = simCtx->block_number;
3991 const PetscInt immersed = simCtx->immersed;
3992 const PetscInt MHV = simCtx->MHV;
3993 const PetscInt LV = simCtx->LV;
3994 PetscMPIInt rank = simCtx->rank;
3995 // --- END CONTEXT ACQUISITION BLOCK ---
3996
3997 PetscErrorCode ierr;
3998 PetscInt l, bi;
3999 MGCtx *mgctx = usermg->mgctx;
4000 KSP mgksp, subksp;
4001 PC mgpc, subpc;
4002 UserCtx *user;
4003
4004 PetscFunctionBeginUser; // Moved to after variable declarations
4005 PROFILE_FUNCTION_BEGIN;
4006 LOG_ALLOW(GLOBAL, LOG_INFO, "Starting Multigrid Poisson Solve...\n");
4007
4008 for (bi = 0; bi < block_number; bi++) {
4009
4010 // ====================================================================
4011 // SECTION: Immersed Boundary Specific Setup (Conditional)
4012 // ====================================================================
4013 if (immersed) {
4014 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Performing IBM pre-solve setup (Nvert restriction, etc.).\n", bi);
4015 for (l = usermg->mglevels - 1; l > 0; l--) {
4016 mgctx[l].user[bi].multinullspace = PETSC_FALSE;
4017 MyNvertRestriction(&mgctx[l].user[bi], &mgctx[l-1].user[bi]);
4018 }
4019 // Coarsest level check for disconnected domains due to IBM
4020 l = 0;
4021 user = mgctx[l].user;
4022 ierr = PetscMalloc1(user[bi].info.mx * user[bi].info.my * 2, &user[bi].KSKE); CHKERRQ(ierr);
4023 FullyBlocked(&user[bi]);
4024 }
4025
4026
4027 l = usermg->mglevels - 1;
4028 user = mgctx[l].user;
4029
4030 // --- 1. Compute RHS of the Poisson Equation ---
4031 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Computing Poisson RHS...\n", bi);
4032 ierr = VecDuplicate(user[bi].P, &user[bi].B); CHKERRQ(ierr);
4033
4034 PetscReal ibm_Flux, ibm_Area;
4035 PetscInt flg = immersed - 1;
4036
4037 // Calculate volume flux source terms (often from IBM)
4038 VolumeFlux(&user[bi], &ibm_Flux, &ibm_Area, flg);
4039 if (MHV || LV) {
4040 flg = ((MHV > 1 || LV) && bi == 0) ? 1 : 0;
4041 VolumeFlux_rev(&user[bi], &ibm_Flux, &ibm_Area, flg);
4042 }
4043 // Calculate the main divergence term
4044 PoissonRHS(&user[bi], user[bi].B);
4045
4046 // --- 2. Assemble LHS Matrix (Laplacian) on all MG levels ---
4047 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Assembling Poisson LHS on all levels...\n", bi);
4048 for (l = usermg->mglevels - 1; l >= 0; l--) {
4049 user = mgctx[l].user;
4050 LOG_ALLOW(GLOBAL,LOG_DEBUG," Calculating LHS for Level %d.\n",l);
4051 PoissonLHSNew(&user[bi]);
4052 }
4053
4054 // --- 3. Setup PETSc KSP and PCMG (Multigrid Preconditioner) ---
4055 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Configuring KSP and PCMG...\n", bi);
4056
4057 ierr = KSPCreate(PETSC_COMM_WORLD, &mgksp); CHKERRQ(ierr);
4058 ierr = KSPAppendOptionsPrefix(mgksp, "ps_"); CHKERRQ(ierr);
4059
4060 // =======================================================================
4061 DualMonitorCtx *monctx;
4062 char filen[128];
4063 PetscBool has_monitor_option;
4064
4065 // 1. Allocate the context and set it up.
4066 ierr = PetscNew(&monctx); CHKERRQ(ierr);
4067
4068 monctx->step = simCtx->step;
4069 monctx->block_id = bi;
4070 monctx->file_handle = NULL;
4071
4072 // Only rank 0 handles the file.
4073 if (!rank) {
4074 sprintf(filen, "logs/Poisson_Solver_Convergence_History_Block_%d.log", bi);
4075 // On the very first step of the entire simulation, TRUNCATE the file.
4076 if (simCtx->step == simCtx->StartStep) {
4077 monctx->file_handle = fopen(filen, "w");
4078 } else { // For all subsequent steps, APPEND to the file.
4079 monctx->file_handle = fopen(filen, "a");
4080 }
4081
4082 if (monctx->file_handle) {
4083 PetscFPrintf(PETSC_COMM_SELF, monctx->file_handle, "--- Convergence for Timestep %d, Block %d ---\n", (int)simCtx->step, bi);
4084 } else {
4085 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open KSP monitor log file: %s", filen);
4086 }
4087 }
4088
4089 ierr = PetscOptionsHasName(NULL, NULL, "-ps_ksp_pic_monitor_true_residual", &has_monitor_option); CHKERRQ(ierr);
4090 monctx->log_to_console = has_monitor_option;
4091
4092 ierr = KSPMonitorSet(mgksp, DualKSPMonitor, monctx, DualMonitorDestroy); CHKERRQ(ierr);
4093 // =======================================================================
4094
4095 ierr = KSPGetPC(mgksp, &mgpc); CHKERRQ(ierr);
4096 ierr = PCSetType(mgpc, PCMG); CHKERRQ(ierr);
4097
4098 ierr = PCMGSetLevels(mgpc, usermg->mglevels, PETSC_NULL); CHKERRQ(ierr);
4099 ierr = PCMGSetCycleType(mgpc, PC_MG_CYCLE_V); CHKERRQ(ierr);
4100 ierr = PCMGSetType(mgpc, PC_MG_MULTIPLICATIVE); CHKERRQ(ierr);
4101
4102 // --- 4. Define Restriction and Interpolation Operators for MG ---
4103 for (l = usermg->mglevels - 1; l > 0; l--) {
4104
4105 // Get stable pointers directly from the main mgctx array.
4106 // These pointers point to memory that will persist.
4107 UserCtx *fine_user_ctx = &mgctx[l].user[bi];
4108 UserCtx *coarse_user_ctx = &mgctx[l-1].user[bi];
4109
4110 // --- Configure the context pointers ---
4111 // The coarse UserCtx needs to know about the fine grid for restriction.
4112 coarse_user_ctx->da_f = &(fine_user_ctx->da);
4113 coarse_user_ctx->user_f = fine_user_ctx;
4114
4115 // The fine UserCtx needs to know about the coarse grid for interpolation.
4116 fine_user_ctx->da_c = &(coarse_user_ctx->da);
4117 fine_user_ctx->user_c = coarse_user_ctx;
4118 fine_user_ctx->lNvert_c = &(coarse_user_ctx->lNvert);
4119
4120 // --- Get matrix dimensions ---
4121 PetscInt m_c = (coarse_user_ctx->info.xm * coarse_user_ctx->info.ym * coarse_user_ctx->info.zm);
4122 PetscInt m_f = (fine_user_ctx->info.xm * fine_user_ctx->info.ym * fine_user_ctx->info.zm);
4123 PetscInt M_c = (coarse_user_ctx->info.mx * coarse_user_ctx->info.my * coarse_user_ctx->info.mz);
4124 PetscInt M_f = (fine_user_ctx->info.mx * fine_user_ctx->info.my * fine_user_ctx->info.mz);
4125
4126 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);
4127 // --- Create the MatShell objects ---
4128 // Pass the STABLE pointer coarse_user_ctx as the context for restriction.
4129 ierr = MatCreateShell(PETSC_COMM_WORLD, m_c, m_f, M_c, M_f, (void*)coarse_user_ctx, &fine_user_ctx->MR); CHKERRQ(ierr);
4130
4131 // Pass the STABLE pointer fine_user_ctx as the context for interpolation.
4132 ierr = MatCreateShell(PETSC_COMM_WORLD, m_f, m_c, M_f, M_c, (void*)fine_user_ctx, &fine_user_ctx->MP); CHKERRQ(ierr);
4133
4134 // --- Set the operations for the MatShells ---
4135 ierr = MatShellSetOperation(fine_user_ctx->MR, MATOP_MULT, (void(*)(void))RestrictResidual_SolidAware); CHKERRQ(ierr);
4136 ierr = MatShellSetOperation(fine_user_ctx->MP, MATOP_MULT, (void(*)(void))MyInterpolation); CHKERRQ(ierr);
4137
4138 // --- Register the operators with PCMG ---
4139 ierr = PCMGSetRestriction(mgpc, l, fine_user_ctx->MR); CHKERRQ(ierr);
4140 ierr = PCMGSetInterpolation(mgpc, l, fine_user_ctx->MP); CHKERRQ(ierr);
4141
4142 }
4143
4144 // --- 5. Configure Solvers on Each MG Level ---
4145 for (l = usermg->mglevels - 1; l >= 0; l--) {
4146 user = mgctx[l].user;
4147 if (l > 0) { // Smoother for fine levels
4148 ierr = PCMGGetSmoother(mgpc, l, &subksp); CHKERRQ(ierr);
4149 } else { // Direct or iterative solver for the coarsest level
4150 ierr = PCMGGetCoarseSolve(mgpc, &subksp); CHKERRQ(ierr);
4151 ierr = KSPSetTolerances(subksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, 40); CHKERRQ(ierr);
4152 }
4153
4154 ierr = KSPSetOperators(subksp, user[bi].A, user[bi].A); CHKERRQ(ierr);
4155 ierr = KSPSetFromOptions(subksp); CHKERRQ(ierr);
4156 ierr = KSPGetPC(subksp, &subpc); CHKERRQ(ierr);
4157 ierr = PCSetType(subpc, PCBJACOBI); CHKERRQ(ierr);
4158
4159 KSP *subsubksp;
4160 PC subsubpc;
4161 PetscInt nlocal;
4162
4163 // This logic is required for both the smoother and the coarse solve
4164 // since both use PCBJACOBI.
4165 ierr = KSPSetUp(subksp); CHKERRQ(ierr); // Set up KSP to allow access to sub-KSPs
4166 ierr = PCBJacobiGetSubKSP(subpc, &nlocal, NULL, &subsubksp); CHKERRQ(ierr);
4167
4168 for (PetscInt abi = 0; abi < nlocal; abi++) {
4169 ierr = KSPGetPC(subsubksp[abi], &subsubpc); CHKERRQ(ierr);
4170 // Add the critical shift amount
4171 ierr = PCFactorSetShiftAmount(subsubpc, 1.e-10); CHKERRQ(ierr);
4172 }
4173
4174 ierr = MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &user[bi].nullsp); CHKERRQ(ierr);
4175 ierr = MatNullSpaceSetFunction(user[bi].nullsp, PoissonNullSpaceFunction, &user[bi]); CHKERRQ(ierr);
4176 ierr = MatSetNullSpace(user[bi].A, user[bi].nullsp); CHKERRQ(ierr);
4177
4178 ierr = PCMGSetResidual(mgpc, l, PCMGResidualDefault, user[bi].A); CHKERRQ(ierr);
4179 ierr = KSPSetUp(subksp); CHKERRQ(ierr);
4180
4181 if (l < usermg->mglevels - 1) {
4182 ierr = MatCreateVecs(user[bi].A, &user[bi].R, PETSC_NULL); CHKERRQ(ierr);
4183 ierr = PCMGSetRhs(mgpc, l, user[bi].R); CHKERRQ(ierr);
4184 }
4185 }
4186
4187 // --- 6. Set Final KSP Operators and Solve ---
4188 l = usermg->mglevels - 1;
4189 user = mgctx[l].user;
4190
4191 LOG_ALLOW(LOCAL, LOG_DEBUG, "Block %d: Setting KSP operators and solving...\n", bi);
4192 ierr = KSPSetOperators(mgksp, user[bi].A, user[bi].A); CHKERRQ(ierr);
4193 ierr = MatSetNullSpace(user[bi].A, user[bi].nullsp); CHKERRQ(ierr);
4194 ierr = KSPSetFromOptions(mgksp); CHKERRQ(ierr);
4195 ierr = KSPSetUp(mgksp); CHKERRQ(ierr);
4196 ierr = KSPSolve(mgksp, user[bi].B, user[bi].Phi); CHKERRQ(ierr);
4197
4198 // --- 7. Cleanup for this block ---
4199 for (l = usermg->mglevels - 1; l >= 0; l--) {
4200 user = mgctx[l].user;
4201 MatNullSpaceDestroy(&user[bi].nullsp);
4202 MatDestroy(&user[bi].A);
4203 user[bi].assignedA = PETSC_FALSE;
4204 if (l > 0) {
4205 MatDestroy(&user[bi].MR);
4206 MatDestroy(&user[bi].MP);
4207 } else if (l==0 && immersed) {
4208 PetscFree(user[bi].KSKE);
4209 }
4210 if (l < usermg->mglevels - 1) {
4211 VecDestroy(&user[bi].R);
4212 }
4213 }
4214
4215 KSPDestroy(&mgksp);
4216 VecDestroy(&mgctx[usermg->mglevels-1].user[bi].B);
4217
4218 } // End of loop over blocks
4219
4220 LOG_ALLOW(GLOBAL, LOG_INFO, "Multigrid Poisson Solve complete.\n");
4221 PROFILE_FUNCTION_END;
4222 PetscFunctionReturn(0);
4223}
logging.h
Logging utilities and macros for PETSc-based applications.
DualMonitorDestroy
PetscErrorCode DualMonitorDestroy(void **ctx)
Destroys the DualMonitorCtx.
Definition logging.c:727
DualMonitorCtx::log_to_console
PetscBool log_to_console
Definition logging.h:65
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:207
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:767
DualMonitorCtx::step
PetscInt step
Definition logging.h:67
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:355
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:759
LOG_INFO
@ LOG_INFO
Informational messages about program execution.
Definition logging.h:32
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:758
DualMonitorCtx::file_handle
FILE * file_handle
Definition logging.h:64
DualMonitorCtx::block_id
PetscInt block_id
Definition logging.h:68
DualMonitorCtx
Context for a dual-purpose KSP monitor.
Definition logging.h:63
TW
#define TW
Definition poisson.c:946
MyNvertRestriction
PetscErrorCode MyNvertRestriction(UserCtx *user_h, UserCtx *user_c)
Definition poisson.c:3867
SE
#define SE
Definition poisson.c:937
BN
#define BN
Definition poisson.c:941
PoissonNullSpaceFunction
PetscErrorCode PoissonNullSpaceFunction(MatNullSpace nullsp, Vec X, void *ctx)
The callback function for PETSc's MatNullSpace object.
Definition poisson.c:1722
PoissonLHSNew
PetscErrorCode PoissonLHSNew(UserCtx *user)
Assembles the Left-Hand Side (LHS) matrix for the Pressure Poisson Equation.
Definition poisson.c:2240
WP
#define WP
Definition poisson.c:929
Projection
PetscErrorCode Projection(UserCtx *user)
Performs the projection step to enforce an incompressible velocity field.
Definition poisson.c:988
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:2914
Gidx
static PetscInt Gidx(PetscInt i, PetscInt j, PetscInt k, UserCtx *user)
Definition poisson.c:879
mymatmultadd
static PetscErrorCode mymatmultadd(Mat mat, Vec v1, Vec v2)
Definition poisson.c:1711
SW
#define SW
Definition poisson.c:939
BS
#define BS
Definition poisson.c:943
NE
#define NE
Definition poisson.c:936
MyRestriction
PetscErrorCode MyRestriction(Mat A, Vec X, Vec F)
The callback function for the multigrid restriction operator (MatShell).
Definition poisson.c:2097
GhostNodeVelocity
static PetscErrorCode GhostNodeVelocity(UserCtx *user)
Definition poisson.c:59
CP
#define CP
Definition poisson.c:926
BE
#define BE
Definition poisson.c:945
UpdatePressure
PetscErrorCode UpdatePressure(UserCtx *user)
Updates the pressure field and handles periodic boundary conditions.
Definition poisson.c:1542
BP
#define BP
Definition poisson.c:933
GridInterpolation
#define GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user)
Definition poisson.c:844
BW
#define BW
Definition poisson.c:947
RestrictResidual_SolidAware
static PetscErrorCode RestrictResidual_SolidAware(Mat A, Vec X, Vec F)
Definition poisson.c:2022
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:3156
TE
#define TE
Definition poisson.c:944
Calculate_Covariant_metrics
static void Calculate_Covariant_metrics(double g[3][3], double G[3][3])
Definition poisson.c:5
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:2837
MyInterpolation
PetscErrorCode MyInterpolation(Mat A, Vec X, Vec F)
The callback function for the multigrid interpolation operator (MatShell).
Definition poisson.c:1915
Calculate_normal
static void Calculate_normal(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3])
Definition poisson.c:24
TS
#define TS
Definition poisson.c:942
FullyBlocked
PetscErrorCode FullyBlocked(UserCtx *user)
Definition poisson.c:3784
NP
#define NP
Definition poisson.c:930
PoissonSolver_MG
PetscErrorCode PoissonSolver_MG(UserMG *usermg)
Solves the pressure-Poisson equation using a geometric multigrid method.
Definition poisson.c:3982
EP
#define EP
Definition poisson.c:928
TN
#define TN
Definition poisson.c:940
SP
#define SP
Definition poisson.c:931
TP
#define TP
Definition poisson.c:932
NW
#define NW
Definition poisson.c:938
GridRestriction
static PetscErrorCode GridRestriction(PetscInt i, PetscInt j, PetscInt k, PetscInt *ih, PetscInt *jh, PetscInt *kh, UserCtx *user)
Definition poisson.c:896
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:1785
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:779
SimCtx::MHV
PetscInt MHV
Definition variables.h:540
UserCtx::isc
PetscInt isc
Definition variables.h:643
MGCtx::user
UserCtx * user
Definition variables.h:418
SimCtx::rank
PetscMPIInt rank
Definition variables.h:516
SimCtx::fish
PetscInt fish
Definition variables.h:540
SimCtx::LV
PetscInt LV
Definition variables.h:540
SimCtx::block_number
PetscInt block_number
Definition variables.h:562
UserCtx::lIEta
Vec lIEta
Definition variables.h:675
UserCtx::KSKE
PetscInt * KSKE
Definition variables.h:665
UserCtx::lIZet
Vec lIZet
Definition variables.h:675
UserCtx::user_f
UserCtx * user_f
Definition variables.h:690
SimCtx::channelz
PetscInt channelz
Definition variables.h:562
UserCtx::lNvert
Vec lNvert
Definition variables.h:657
UserCtx::Phi
Vec Phi
Definition variables.h:657
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:633
SimCtx::rans
PetscInt rans
Definition variables.h:578
UserCtx::ksc
PetscInt ksc
Definition variables.h:643
UserCtx::KM
PetscInt KM
Definition variables.h:639
SimCtx::st
PetscReal st
Definition variables.h:549
UserCtx::lZet
Vec lZet
Definition variables.h:673
UserCtx::assignedA
PetscBool assignedA
Definition variables.h:669
UserCtx::lIAj
Vec lIAj
Definition variables.h:675
UserCtx::lKEta
Vec lKEta
Definition variables.h:677
SimCtx::dt
PetscReal dt
Definition variables.h:527
UserCtx::lUstar
Vec lUstar
Definition variables.h:652
UserCtx::jsc
PetscInt jsc
Definition variables.h:643
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:676
UserCtx::Ucont
Vec Ucont
Definition variables.h:657
SimCtx::StartStep
PetscInt StartStep
Definition variables.h:523
BCS::Ubcs
Vec Ubcs
Boundary condition velocity values. (Comment: "An ugly hack, waste of memory")
Definition variables.h:120
Cmpnts::x
PetscScalar x
Definition variables.h:100
UserCtx::Bcs
BCS Bcs
Definition variables.h:651
UserCtx::lPhi
Vec lPhi
Definition variables.h:657
UserCtx::da_c
DM * da_c
Definition variables.h:691
UserCtx::user_c
UserCtx * user_c
Definition variables.h:690
SimCtx::NumberOfBodies
PetscInt NumberOfBodies
Definition variables.h:555
UserCtx::lKZet
Vec lKZet
Definition variables.h:677
UserCtx::da_f
DM * da_f
Definition variables.h:691
UserCtx::lJEta
Vec lJEta
Definition variables.h:676
UserCtx::lCsi
Vec lCsi
Definition variables.h:673
UserCtx::thislevel
PetscInt thislevel
Definition variables.h:689
Cmpnts::z
PetscScalar z
Definition variables.h:100
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:677
UserCtx::Ucat
Vec Ucat
Definition variables.h:657
UserCtx::JM
PetscInt JM
Definition variables.h:639
SimCtx::wallfunction
PetscInt wallfunction
Definition variables.h:579
UserMG::mglevels
PetscInt mglevels
Definition variables.h:425
UserCtx::bctype
PetscInt bctype[6]
Definition variables.h:653
UserCtx::lJZet
Vec lJZet
Definition variables.h:676
SimCtx::U_bc
PetscReal U_bc
Definition variables.h:573
UserCtx::lUcont
Vec lUcont
Definition variables.h:657
SimCtx::step
PetscInt step
Definition variables.h:521
UserCtx::lAj
Vec lAj
Definition variables.h:673
UserCtx::lICsi
Vec lICsi
Definition variables.h:675
UserCtx::info
DMDALocalInfo info
Definition variables.h:637
UserCtx::lUcat
Vec lUcat
Definition variables.h:657
Cmpnts::y
PetscScalar y
Definition variables.h:100
UserCtx::IM
PetscInt IM
Definition variables.h:639
UserCtx::FluxIntpSum
PetscReal FluxIntpSum
Definition variables.h:654
UserCtx::lEta
Vec lEta
Definition variables.h:673
UserCtx::Cent
Vec Cent
Definition variables.h:673
SimCtx::les
PetscInt les
Definition variables.h:578
UserCtx::Nvert
Vec Nvert
Definition variables.h:657
UserMG::mgctx
MGCtx * mgctx
Definition variables.h:428
UserCtx::multinullspace
PetscBool multinullspace
Definition variables.h:666
UserCtx::lJAj
Vec lJAj
Definition variables.h:676
SimCtx::summationRHS
PetscReal summationRHS
Definition variables.h:605
UserCtx::lKAj
Vec lKAj
Definition variables.h:677
SimCtx::immersed
PetscInt immersed
Definition variables.h:535
COEF_TIME_ACCURACY
#define COEF_TIME_ACCURACY
Coefficient controlling the temporal accuracy scheme (e.g., 1.5 for BDF2).
Definition variables.h:57
UserCtx::lNvert_c
Vec * lNvert_c
Definition variables.h:692
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:99
MGCtx
Context for Multigrid operations.
Definition variables.h:417
SimCtx
The master context for the entire simulation.
Definition variables.h:513
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:630
UserMG
User-level context for managing the entire multigrid hierarchy.
Definition variables.h:424
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