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