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

Go to the source code of this file.

Macros

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

Functions

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

Macro Definition Documentation

◆ GridInterpolation

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

Definition at line 844 of file poisson.c.

845 { \
846 ic = i; \
847 ia = 0; \
848 } \
849 else { \
850 ic = (i+1) / 2; \
851 ia = (i - 2 * (ic)) == 0 ? 1 : -1; \
852 if (i==1 || i==mx-2) ia = 0; \
853 }\
854 if ((user->jsc)) { \
855 jc = j; \
856 ja = 0; \
857 } \
858 else { \
859 jc = (j+1) / 2; \
860 ja = (j - 2 * (jc)) == 0 ? 1 : -1; \
861 if (j==1 || j==my-2) ja = 0; \
862 } \
863 if ((user->ksc)) { \
864 kc = k; \
865 ka = 0; \
866 } \
867 else { \
868 kc = (k+1) / 2; \
869 ka = (k - 2 * (kc)) == 0 ? 1 : -1; \
870 if (k==1 || k==mz-2) ka = 0; \
871 } \
872 if (ka==-1 && nvert_c[kc-1][jc][ic] > 0.1) ka = 0; \
873 else if (ka==1 && nvert_c[kc+1][jc][ic] > 0.1) ka = 0; \
874 if (ja==-1 && nvert_c[kc][jc-1][ic] > 0.1) ja = 0; \
875 else if (ja==1 && nvert_c[kc][jc+1][ic] > 0.1) ja = 0; \
876 if (ia==-1 && nvert_c[kc][jc][ic-1] > 0.1) ia = 0; \
877 else if (ia==1 && nvert_c[kc][jc][ic+1] > 0.1) ia = 0;

◆ CP

#define CP   0

Definition at line 926 of file poisson.c.

◆ EP

#define EP   1

Definition at line 928 of file poisson.c.

◆ WP

#define WP   2

Definition at line 929 of file poisson.c.

◆ NP

#define NP   3

Definition at line 930 of file poisson.c.

◆ SP

#define SP   4

Definition at line 931 of file poisson.c.

◆ TP

#define TP   5

Definition at line 932 of file poisson.c.

◆ BP

#define BP   6

Definition at line 933 of file poisson.c.

◆ NE

#define NE   7

Definition at line 936 of file poisson.c.

◆ SE

#define SE   8

Definition at line 937 of file poisson.c.

◆ NW

#define NW   9

Definition at line 938 of file poisson.c.

◆ SW

#define SW   10

Definition at line 939 of file poisson.c.

◆ TN

#define TN   11

Definition at line 940 of file poisson.c.

◆ BN

#define BN   12

Definition at line 941 of file poisson.c.

◆ TS

#define TS   13

Definition at line 942 of file poisson.c.

◆ BS

#define BS   14

Definition at line 943 of file poisson.c.

◆ TE

#define TE   15

Definition at line 944 of file poisson.c.

◆ BE

#define BE   16

Definition at line 945 of file poisson.c.

◆ TW

#define TW   17

Definition at line 946 of file poisson.c.

◆ BW

#define BW   18

Definition at line 947 of file poisson.c.

◆ __FUNCT__ [1/3]

#define __FUNCT__   "Projection"

Definition at line 950 of file poisson.c.

◆ __FUNCT__ [2/3]

#define __FUNCT__   "UpdatePressure"

Definition at line 950 of file poisson.c.

◆ __FUNCT__ [3/3]

#define __FUNCT__   "PoissonSolver_MG"

Definition at line 950 of file poisson.c.

Function Documentation

◆ Calculate_Covariant_metrics()

static void Calculate_Covariant_metrics ( double  g[3][3],
double  G[3][3] 
)
static

Definition at line 5 of file poisson.c.

6{
7 /*
8 | csi.x csi.y csi.z |-1 | x.csi x.eta x.zet |
9 | eta.x eta.y eta.z | = | y.csi y.eta y.zet |
10 | zet.x zet.y zet.z | | z.csi z.eta z.zet |
11
12 */
13 const double a11=g[0][0], a12=g[0][1], a13=g[0][2];
14 const double a21=g[1][0], a22=g[1][1], a23=g[1][2];
15 const double a31=g[2][0], a32=g[2][1], a33=g[2][2];
16
17 double det= a11*(a33*a22-a32*a23)-a21*(a33*a12-a32*a13)+a31*(a23*a12-a22*a13);
18
19 G[0][0] = (a33*a22-a32*a23)/det, G[0][1] = - (a33*a12-a32*a13)/det, G[0][2] = (a23*a12-a22*a13)/det;
20 G[1][0] = -(a33*a21-a31*a23)/det, G[1][1] = (a33*a11-a31*a13)/det, G[1][2] = - (a23*a11-a21*a13)/det;
21 G[2][0] = (a32*a21-a31*a22)/det, G[2][1] = - (a32*a11-a31*a12)/det, G[2][2] = (a22*a11-a21*a12)/det;
22};
Here is the caller graph for this function:

◆ Calculate_normal()

static void Calculate_normal ( Cmpnts  csi,
Cmpnts  eta,
Cmpnts  zet,
double  ni[3],
double  nj[3],
double  nk[3] 
)
static

Definition at line 24 of file poisson.c.

25{
26 double g[3][3];
27 double G[3][3];
28
29 g[0][0]=csi.x, g[0][1]=csi.y, g[0][2]=csi.z;
30 g[1][0]=eta.x, g[1][1]=eta.y, g[1][2]=eta.z;
31 g[2][0]=zet.x, g[2][1]=zet.y, g[2][2]=zet.z;
32
33 Calculate_Covariant_metrics(g, G);
34 double xcsi=G[0][0], ycsi=G[1][0], zcsi=G[2][0];
35 double xeta=G[0][1], yeta=G[1][1], zeta=G[2][1];
36 double xzet=G[0][2], yzet=G[1][2], zzet=G[2][2];
37
38 double nx_i = xcsi, ny_i = ycsi, nz_i = zcsi;
39 double nx_j = xeta, ny_j = yeta, nz_j = zeta;
40 double nx_k = xzet, ny_k = yzet, nz_k = zzet;
41
42 double sum_i=sqrt(nx_i*nx_i+ny_i*ny_i+nz_i*nz_i);
43 double sum_j=sqrt(nx_j*nx_j+ny_j*ny_j+nz_j*nz_j);
44 double sum_k=sqrt(nx_k*nx_k+ny_k*ny_k+nz_k*nz_k);
45
46// *Ai = sqrt( g[0][0]*g[0][0] + g[0][1]*g[0][1] + g[0][2]*g[0][2] ); // area
47// *Aj = sqrt( g[1][0]*g[1][0] + g[1][1]*g[1][1] + g[1][2]*g[1][2] );
48// *Ak =sqrt( g[2][0]*g[2][0] + g[2][1]*g[2][1] + g[2][2]*g[2][2] );
49
50 nx_i /= sum_i, ny_i /= sum_i, nz_i /= sum_i;
51 nx_j /= sum_j, ny_j /= sum_j, nz_j /= sum_j;
52 nx_k /= sum_k, ny_k /= sum_k, nz_k /= sum_k;
53
54 ni[0] = nx_i, ni[1] = ny_i, ni[2] = nz_i;
55 nj[0] = nx_j, nj[1] = ny_j, nj[2] = nz_j;
56 nk[0] = nx_k, nk[1] = ny_k, nk[2] = nz_k;
57}
Calculate_Covariant_metrics
static void Calculate_Covariant_metrics(double g[3][3], double G[3][3])
Definition poisson.c:5
Cmpnts::x
PetscScalar x
Definition variables.h:100
Cmpnts::z
PetscScalar z
Definition variables.h:100
Cmpnts::y
PetscScalar y
Definition variables.h:100
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GhostNodeVelocity()

static PetscErrorCode GhostNodeVelocity ( UserCtx *  user)
static

Definition at line 59 of file poisson.c.

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

◆ Gidx()

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

Definition at line 879 of file poisson.c.

881{
882 PetscInt nidx;
883 DMDALocalInfo info = user->info;
884
885 PetscInt mx = info.mx, my = info.my;
886
887 AO ao;
888 DMDAGetAO(user->da, &ao);
889 nidx=i+j*mx+k*mx*my;
890
891 AOApplicationToPetsc(ao,1,&nidx);
892
893 return (nidx);
894}
Here is the caller graph for this function:

◆ GridRestriction()

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

Definition at line 896 of file poisson.c.

899{
900 if ((user->isc)) {
901 *ih = i;
902 }
903 else {
904 *ih = 2 * i;
905 }
906
907 if ((user->jsc)) {
908 *jh = j;
909 }
910 else {
911 *jh = 2 * j;
912 }
913
914 if ((user->ksc)) {
915 *kh = k;
916 }
917 else {
918 *kh = 2 * k;
919 }
920
921 return 0;
922}
UserCtx::isc
PetscInt isc
Definition variables.h:643
UserCtx::ksc
PetscInt ksc
Definition variables.h:643
UserCtx::jsc
PetscInt jsc
Definition variables.h:643
Here is the caller graph for this function:

◆ Projection()

PetscErrorCode Projection ( UserCtx *  user)

Performs the projection step to enforce an incompressible velocity field.

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

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

The main steps are:

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

Definition at line 988 of file poisson.c.

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

◆ UpdatePressure()

PetscErrorCode UpdatePressure ( UserCtx *  user)

Updates the pressure field and handles periodic boundary conditions.

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

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

The function performs two main operations:

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

Definition at line 1542 of file poisson.c.

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

◆ mymatmultadd()

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

Definition at line 1711 of file poisson.c.

1712{
1713
1714 Vec vt;
1715 VecDuplicate(v2, &vt);
1716 MatMult(mat, v1, vt);
1717 VecAYPX(v2, 1., vt);
1718 VecDestroy(&vt);
1719 return(0);
1720}

◆ PoissonNullSpaceFunction()

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

The callback function for PETSc's MatNullSpace object.

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

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

Definition at line 1722 of file poisson.c.

1723{
1724 UserCtx *user = (UserCtx*)ctx;
1725
1726 DM da = user->da;
1727
1728 DMDALocalInfo info = user->info;
1729 PetscInt xs = info.xs, xe = info.xs + info.xm;
1730 PetscInt ys = info.ys, ye = info.ys + info.ym;
1731 PetscInt zs = info.zs, ze = info.zs + info.zm;
1732 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1733 PetscInt lxs, lxe, lys, lye, lzs, lze;
1734
1735 PetscReal ***x, ***nvert;
1736 PetscInt i, j, k;
1737
1738/* /\* First remove a constant from the Vec field X *\/ */
1739
1740
1741 /* Then apply boundary conditions */
1742 DMDAVecGetArray(da, X, &x);
1743 DMDAVecGetArray(da, user->lNvert, &nvert);
1744
1745 lxs = xs; lxe = xe;
1746 lys = ys; lye = ye;
1747 lzs = zs; lze = ze;
1748
1749 if (xs==0) lxs = xs+1;
1750 if (ys==0) lys = ys+1;
1751 if (zs==0) lzs = zs+1;
1752
1753 if (xe==mx) lxe = xe-1;
1754 if (ye==my) lye = ye-1;
1755 if (ze==mz) lze = ze-1;
1756
1757 PetscReal lsum, sum;
1758 PetscReal lnum, num;
1759
1760 if (user->multinullspace) PetscPrintf(PETSC_COMM_WORLD, "MultiNullSpace!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");
1761 if (!user->multinullspace) {
1762 lsum = 0;
1763 lnum = 0;
1764 for (k=lzs; k<lze; k++) {
1765 for (j=lys; j<lye; j++) {
1766 for (i=lxs; i<lxe; i++) {
1767 if (nvert[k][j][i] < 0.1) {
1768 lsum += x[k][j][i];
1769 lnum ++;
1770 }
1771 }
1772 }
1773 }
1774
1775 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1776 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1777 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1778/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1779 sum = sum / (-1.0 * num);
1780
1781 for (k=lzs; k<lze; k++) {
1782 for (j=lys; j<lye; j++) {
1783 for (i=lxs; i<lxe; i++) {
1784 if (nvert[k][j][i] < 0.1) {
1785 x[k][j][i] +=sum;
1786 }
1787 }
1788 }
1789 }
1790 }
1791 else {
1792 lsum = 0;
1793 lnum = 0;
1794 for (j=lys; j<lye; j++) {
1795 for (i=lxs; i<lxe; i++) {
1796 for (k=lzs; k<lze; k++) {
1797 if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1798 lsum += x[k][j][i];
1799 lnum ++;
1800 }
1801 }
1802 }
1803 }
1804 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1805 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1806 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1807/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1808 sum /= -num;
1809 for (j=lys; j<lye; j++) {
1810 for (i=lxs; i<lxe; i++) {
1811 for (k=lzs; k<lze; k++) {
1812 if (k<user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1813 x[k][j][i] += sum;
1814 }
1815 }
1816 }
1817 }
1818
1819 lsum = 0;
1820 lnum = 0;
1821 for (j=lys; j<lye; j++) {
1822 for (i=lxs; i<lxe; i++) {
1823 for (k=lzs; k<lze; k++) {
1824 if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1825 lsum += x[k][j][i];
1826 lnum ++;
1827 }
1828 }
1829 }
1830 }
1831 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1832 MPI_Allreduce(&lnum,&num,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
1833 /* PetscGlobalSum(&lsum, &sum, PETSC_COMM_WORLD); */
1834/* PetscGlobalSum(&lnum, &num, PETSC_COMM_WORLD); */
1835 sum /= -num;
1836 for (j=lys; j<lye; j++) {
1837 for (i=lxs; i<lxe; i++) {
1838 for (k=lzs; k<lze; k++) {
1839 if (k>=user->KSKE[2*(j*mx+i)] && nvert[k][j][i]<0.1) {
1840 x[k][j][i] += sum;
1841 }
1842 }
1843 }
1844 }
1845
1846 } //if multinullspace
1847 if (zs == 0) {
1848 k = 0;
1849 for (j=ys; j<ye; j++) {
1850 for (i=xs; i<xe; i++) {
1851 x[k][j][i] = 0.;
1852 }
1853 }
1854 }
1855
1856 if (ze == mz) {
1857 k = mz-1;
1858 for (j=ys; j<ye; j++) {
1859 for (i=xs; i<xe; i++) {
1860 x[k][j][i] = 0.;
1861 }
1862 }
1863 }
1864
1865 if (ys == 0) {
1866 j = 0;
1867 for (k=zs; k<ze; k++) {
1868 for (i=xs; i<xe; i++) {
1869 x[k][j][i] = 0.;
1870 }
1871 }
1872 }
1873
1874 if (ye == my) {
1875 j = my-1;
1876 for (k=zs; k<ze; k++) {
1877 for (i=xs; i<xe; i++) {
1878 x[k][j][i] = 0.;
1879 }
1880 }
1881 }
1882
1883 if (xs == 0) {
1884 i = 0;
1885 for (k=zs; k<ze; k++) {
1886 for (j=ys; j<ye; j++) {
1887 x[k][j][i] = 0.;
1888 }
1889 }
1890 }
1891
1892 if (xe == mx) {
1893 i = mx-1;
1894 for (k=zs; k<ze; k++) {
1895 for (j=ys; j<ye; j++) {
1896 x[k][j][i] = 0.;
1897 }
1898 }
1899 }
1900
1901 for (k=zs; k<ze; k++) {
1902 for (j=ys; j<ye; j++) {
1903 for (i=xs; i<xe; i++) {
1904 if (nvert[k][j][i] > 0.1)
1905 x[k][j][i] = 0.;
1906 }
1907 }
1908 }
1909 DMDAVecRestoreArray(da, X, &x);
1910 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1911
1912 return 0;
1913}
UserCtx::KSKE
PetscInt * KSKE
Definition variables.h:665
UserCtx::multinullspace
PetscBool multinullspace
Definition variables.h:666
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:630
Here is the caller graph for this function:

◆ MyInterpolation()

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

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

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

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

Definition at line 1915 of file poisson.c.

1916{
1917 UserCtx *user;
1918
1919 MatShellGetContext(A, (void**)&user);
1920
1921
1922
1923 DM da = user->da;
1924
1925 DM da_c = *user->da_c;
1926
1927 DMDALocalInfo info = user->info;
1928 PetscInt xs = info.xs, xe = info.xs + info.xm;
1929 PetscInt ys = info.ys, ye = info.ys + info.ym;
1930 PetscInt zs = info.zs, ze = info.zs + info.zm;
1931 PetscInt mx = info.mx, my = info.my, mz = info.mz;
1932 PetscInt lxs, lxe, lys, lye, lzs, lze;
1933
1934 PetscReal ***f, ***x, ***nvert, ***nvert_c;
1935 PetscInt i, j, k, ic, jc, kc, ia, ja, ka;
1936
1937 lxs = xs; lxe = xe;
1938 lys = ys; lye = ye;
1939 lzs = zs; lze = ze;
1940
1941 if (xs==0) lxs = xs+1;
1942 if (ys==0) lys = ys+1;
1943 if (zs==0) lzs = zs+1;
1944
1945 if (xe==mx) lxe = xe-1;
1946 if (ye==my) lye = ye-1;
1947 if (ze==mz) lze = ze-1;
1948
1949
1950 DMDAVecGetArray(da, F, &f);
1951
1952
1953 Vec lX;
1954 DMCreateLocalVector(da_c, &lX);
1955
1956 DMGlobalToLocalBegin(da_c, X, INSERT_VALUES, lX);
1957 DMGlobalToLocalEnd(da_c, X, INSERT_VALUES, lX);
1958 DMDAVecGetArray(da_c, lX, &x);
1959
1960 DMDAVecGetArray(da, user->lNvert, &nvert);
1961 DMDAVecGetArray(da_c, *(user->lNvert_c), &nvert_c);
1962 for (k=lzs; k<lze; k++) {
1963 for (j=lys; j<lye; j++) {
1964 for (i=lxs; i<lxe; i++) {
1965
1966 GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user);
1967
1968 f[k][j][i] = (x[kc ][jc ][ic ] * 9 +
1969 x[kc ][jc+ja][ic ] * 3 +
1970 x[kc ][jc ][ic+ia] * 3 +
1971 x[kc ][jc+ja][ic+ia]) * 3./64. +
1972 (x[kc+ka][jc ][ic ] * 9 +
1973 x[kc+ka][jc+ja][ic ] * 3 +
1974 x[kc+ka][jc ][ic+ia] * 3 +
1975 x[kc+ka][jc+ja][ic+ia]) /64.;
1976 }
1977 }
1978 }
1979
1980 for (k=zs; k<ze; k++) {
1981 for (j=ys; j<ye; j++) {
1982 for (i=xs; i<xe; i++) {
1983
1984 if (i==0) {
1985 f[k][j][i] = 0.;//-f[k][j][i+1];
1986 }
1987 else if (i==mx-1) {
1988 f[k][j][i] = 0.;//-f[k][j][i-1];
1989 }
1990 else if (j==0) {
1991 f[k][j][i] = 0.;//-f[k][j+1][i];
1992 }
1993 else if (j==my-1) {
1994 f[k][j][i] = 0.;//-f[k][j-1][i];
1995 }
1996 else if (k==0) {
1997 f[k][j][i] = 0.;//-f[k+1][j][i];
1998 }
1999 else if (k==mz-1) {
2000 f[k][j][i] = 0.;//-f[k-1][j][i];
2001 }
2002 if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
2003
2004 }
2005 }
2006 }
2007
2008 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2009 DMDAVecRestoreArray(da_c, *(user->lNvert_c), &nvert_c);
2010
2011 DMDAVecRestoreArray(da_c, lX, &x);
2012
2013 VecDestroy(&lX);
2014 DMDAVecRestoreArray(da, F, &f);
2015
2016
2017
2018 return 0;
2019
2020}
GridInterpolation
#define GridInterpolation(i, j, k, ic, jc, kc, ia, ja, ka, user)
Definition poisson.c:844
UserCtx::da_c
DM * da_c
Definition variables.h:691
UserCtx::lNvert_c
Vec * lNvert_c
Definition variables.h:692
Here is the caller graph for this function:

◆ RestrictResidual_SolidAware()

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

Definition at line 2022 of file poisson.c.

2023{
2024 UserCtx *user;
2025 MatShellGetContext(A, (void**)&user);
2026
2027 DM da = user->da;
2028 DM da_f = *user->da_f;
2029
2030 DMDALocalInfo info;
2031 DMDAGetLocalInfo(da, &info);
2032 PetscInt xs = info.xs, xe = info.xs + info.xm;
2033 PetscInt ys = info.ys, ye = info.ys + info.ym;
2034 PetscInt zs = info.zs, ze = info.zs + info.zm;
2035 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2036
2037 PetscReal ***f, ***x, ***nvert;
2038 PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
2039
2040 DMDAVecGetArray(da, F, &f);
2041
2042 Vec lX;
2043 DMCreateLocalVector(da_f, &lX);
2044 DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
2045 DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);
2046 DMDAVecGetArray(da_f, lX, &x);
2047
2048 DMDAVecGetArray(da, user->lNvert, &nvert);
2049
2050 PetscReal ***nvert_f;
2051 DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
2052
2053 if ((user->isc)) ia = 0;
2054 else ia = 1;
2055
2056 if ((user->jsc)) ja = 0;
2057 else ja = 1;
2058
2059 if ((user->ksc)) ka = 0;
2060 else ka = 1;
2061
2062 for (k=zs; k<ze; k++) {
2063 for (j=ys; j<ye; j++) {
2064 for (i=xs; i<xe; i++) {
2065 // --- CORRECTED LOGIC ---
2066 // First, check if the current point is a boundary point.
2067 // If it is, it does not contribute to the coarse grid residual.
2068 if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1 || nvert[k][j][i] > 0.1) {
2069 f[k][j][i] = 0.0;
2070 }
2071 // Only if it's a true interior fluid point, perform the restriction.
2072 else {
2073 GridRestriction(i, j, k, &ih, &jh, &kh, user);
2074 f[k][j][i] = 0.125 *
2075 (x[kh ][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih ]) +
2076 x[kh ][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih-ia]) +
2077 x[kh ][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih ]) +
2078 x[kh-ka][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih ]) +
2079 x[kh ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih-ia]) +
2080 x[kh-ka][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih ]) +
2081 x[kh-ka][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih-ia]) +
2082 x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
2083 }
2084 }
2085 }
2086 }
2087
2088 DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
2089 DMDAVecRestoreArray(da_f, lX, &x);
2090 VecDestroy(&lX);
2091 DMDAVecRestoreArray(da, F, &f);
2092 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2093
2094 return 0;
2095}
GridRestriction
static PetscErrorCode GridRestriction(PetscInt i, PetscInt j, PetscInt k, PetscInt *ih, PetscInt *jh, PetscInt *kh, UserCtx *user)
Definition poisson.c:896
UserCtx::user_f
UserCtx * user_f
Definition variables.h:690
UserCtx::da_f
DM * da_f
Definition variables.h:691
Here is the call graph for this function:
Here is the caller graph for this function:

◆ MyRestriction()

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

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

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

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

Definition at line 2097 of file poisson.c.

2098{
2099 UserCtx *user;
2100
2101 MatShellGetContext(A, (void**)&user);
2102
2103
2104 DM da = user->da;
2105
2106 DM da_f = *user->da_f;
2107
2108 DMDALocalInfo info;
2109 DMDAGetLocalInfo(da, &info);
2110 PetscInt xs = info.xs, xe = info.xs + info.xm;
2111 PetscInt ys = info.ys, ye = info.ys + info.ym;
2112 PetscInt zs = info.zs, ze = info.zs + info.zm;
2113 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2114 // PetscInt lxs, lxe, lys, lye, lzs, lze;
2115
2116 PetscReal ***f, ***x, ***nvert;
2117 PetscInt i, j, k, ih, jh, kh, ia, ja, ka;
2118
2119 DMDAVecGetArray(da, F, &f);
2120
2121 Vec lX;
2122
2123 DMCreateLocalVector(da_f, &lX);
2124 DMGlobalToLocalBegin(da_f, X, INSERT_VALUES, lX);
2125 DMGlobalToLocalEnd(da_f, X, INSERT_VALUES, lX);
2126 DMDAVecGetArray(da_f, lX, &x);
2127
2128 DMDAVecGetArray(da, user->lNvert, &nvert);
2129
2130 PetscReal ***nvert_f;
2131 DMDAVecGetArray(da_f, user->user_f->lNvert, &nvert_f);
2132
2133 if ((user->isc)) ia = 0;
2134 else ia = 1;
2135
2136 if ((user->jsc)) ja = 0;
2137 else ja = 1;
2138
2139 if ((user->ksc)) ka = 0;
2140 else ka = 1;
2141
2142 for (k=zs; k<ze; k++) {
2143 for (j=ys; j<ye; j++) {
2144 for (i=xs; i<xe; i++) {
2145 if (k==0) {
2146 f[k][j][i] = 0.;
2147 }
2148 else if (k==mz-1) {
2149 f[k][j][i] = 0.;
2150 }
2151 else if (j==0) {
2152 f[k][j][i] = 0.;
2153 }
2154 else if (j==my-1) {
2155 f[k][j][i] = 0.;
2156 }
2157 else if (i==0) {
2158 f[k][j][i] = 0.;
2159 }
2160 else if (i==mx-1) {
2161 f[k][j][i] = 0.;
2162 }
2163 else {
2164 GridRestriction(i, j, k, &ih, &jh, &kh, user);
2165 f[k][j][i] = 0.125 *
2166 (x[kh ][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih ]) +
2167 x[kh ][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh ][ih-ia]) +
2168 x[kh ][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih ]) +
2169 x[kh-ka][jh ][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih ]) +
2170 x[kh ][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh ][jh-ja][ih-ia]) +
2171 x[kh-ka][jh-ja][ih ] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih ]) +
2172 x[kh-ka][jh ][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh ][ih-ia]) +
2173 x[kh-ka][jh-ja][ih-ia] * PetscMax(0., 1 - nvert_f[kh-ka][jh-ja][ih-ia]));
2174
2175
2176
2177 if (nvert[k][j][i] > 0.1) f[k][j][i] = 0.;
2178 }
2179 }
2180 }
2181 }
2182
2183
2184 DMDAVecRestoreArray(da_f, user->user_f->lNvert, &nvert_f);
2185
2186 DMDAVecRestoreArray(da_f, lX, &x);
2187 VecDestroy(&lX);
2188
2189 DMDAVecRestoreArray(da, F, &f);
2190 DMDAVecRestoreArray(da, user->lNvert, &nvert);
2191
2192
2193 return 0;
2194}
Here is the call graph for this function:

◆ PoissonLHSNew()

PetscErrorCode PoissonLHSNew ( UserCtx *  user)

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

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

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

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

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

Definition at line 2240 of file poisson.c.

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

◆ PoissonRHS()

PetscErrorCode PoissonRHS ( UserCtx *  user,
Vec  B 
)

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

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

Definition at line 2837 of file poisson.c.

2838{
2839 DMDALocalInfo info = user->info;
2840 PetscInt xs = info.xs, xe = info.xs + info.xm;
2841 PetscInt ys = info.ys, ye = info.ys + info.ym;
2842 PetscInt zs = info.zs, ze = info.zs + info.zm;
2843 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2844
2845 PetscInt i, j, k;
2846 PetscReal ***nvert, ***aj, ***rb, dt = user->simCtx->dt;
2847 struct Components{
2848 PetscReal x;
2849 PetscReal y;
2850 PetscReal z;
2851 } *** ucont;
2852
2853 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Entering PoissonRHS to compute pressure equation RHS.\n");
2854
2855 DMDAVecGetArray(user->da, B, &rb);
2856 DMDAVecGetArray(user->fda, user->lUcont, &ucont);
2857 DMDAVecGetArray(user->da, user->lNvert, &nvert);
2858 DMDAVecGetArray(user->da, user->lAj, &aj);
2859
2860
2861 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Computing RHS values for each cell.\n");
2862
2863 for (k=zs; k<ze; k++) {
2864 for (j=ys; j<ye; j++) {
2865 for (i=xs; i<xe; i++) {
2866
2867 if (i==0 || i==mx-1 || j==0 || j==my-1 || k==0 || k==mz-1) {
2868 rb[k][j][i] = 0.;
2869 }
2870 else if (nvert[k][j][i] > 0.1) {
2871 rb[k][j][i] = 0;
2872 }
2873 else {
2874 rb[k][j][i] = -(ucont[k][j][i].x - ucont[k][j][i-1].x +
2875 ucont[k][j][i].y - ucont[k][j-1][i].y +
2876 ucont[k][j][i].z - ucont[k-1][j][i].z) / dt
2877 * aj[k][j][i] / user->simCtx->st * COEF_TIME_ACCURACY;
2878
2879 }
2880 }
2881 }
2882 }
2883
2884
2885 // --- Check the solvability condition for the Poisson equation ---
2886 // The global sum of the RHS (proportional to the total divergence) must be zero.
2887 LOG_ALLOW(GLOBAL, LOG_DEBUG, "Verifying solvability condition (sum of RHS terms).\n");
2888 PetscReal lsum=0., sum=0.;
2889
2890 for (k=zs; k<ze; k++) {
2891 for (j=ys; j<ye; j++) {
2892 for (i=xs; i<xe; i++) {
2893
2894 lsum += rb[k][j][i] / aj[k][j][i]* dt/COEF_TIME_ACCURACY;
2895
2896 }
2897 }
2898 }
2899
2900 MPI_Allreduce(&lsum,&sum,1,MPI_DOUBLE,MPI_SUM, PETSC_COMM_WORLD);
2901
2902 LOG_ALLOW(GLOBAL, LOG_INFO, "Global Sum of RHS (Divergence Check): %le\n", sum);
2903
2904 user->simCtx->summationRHS = sum;
2905
2906 DMDAVecRestoreArray(user->fda, user->lUcont, &ucont);
2907 DMDAVecRestoreArray(user->da, user->lNvert, &nvert);
2908 DMDAVecRestoreArray(user->da, user->lAj, &aj);
2909 DMDAVecRestoreArray(user->da, B, &rb);
2910
2911 return 0;
2912}
SimCtx::summationRHS
PetscReal summationRHS
Definition variables.h:605
Here is the caller graph for this function:

◆ VolumeFlux_rev()

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

A specialized version of VolumeFlux, likely for reversed normals.

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

Definition at line 2914 of file poisson.c.

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

◆ VolumeFlux()

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

Calculates the net flux across the immersed boundary surface.

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

Definition at line 3156 of file poisson.c.

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

◆ FullyBlocked()

PetscErrorCode FullyBlocked ( UserCtx *  user)

Definition at line 3784 of file poisson.c.

3785{
3786 DM da = user->da;
3787 Vec nNvert;
3788 DMDALocalInfo info = user->info;
3789
3790 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3791
3792 PetscInt i, j, k;
3793
3794 PetscInt *KSKE = user->KSKE;
3795 PetscReal ***nvert;
3796 PetscBool *Blocked;
3797
3798 DMDACreateNaturalVector(da, &nNvert);
3799 DMDAGlobalToNaturalBegin(da, user->Nvert, INSERT_VALUES, nNvert);
3800 DMDAGlobalToNaturalEnd(da, user->Nvert, INSERT_VALUES, nNvert);
3801
3802 VecScatter ctx;
3803 Vec Zvert;
3804 VecScatterCreateToZero(nNvert, &ctx, &Zvert);
3805
3806 VecScatterBegin(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
3807 VecScatterEnd(ctx, nNvert, Zvert, INSERT_VALUES, SCATTER_FORWARD);
3808
3809 VecScatterDestroy(&ctx);
3810 VecDestroy(&nNvert);
3811
3812 PetscInt rank;
3813 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
3814
3815 if (!rank) {
3816
3817 VecGetArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
3818 PetscMalloc(mx*my*sizeof(PetscBool), &Blocked);
3819 for (j=1; j<my-1; j++) {
3820 for (i=1; i<mx-1; i++) {
3821 Blocked[j*mx+i] = PETSC_FALSE;
3822 for (k=0; k<mz; k++) {
3823 if (nvert[k][j][i] > 0.1) {
3824 if (!Blocked[j*mx+i]) {
3825 KSKE[2*(j*mx+i)] = k;
3826 Blocked[j*mx+i] = PETSC_TRUE;
3827 }
3828 else {
3829 KSKE[2*(j*mx+i)] = PetscMin(KSKE[2*(j*mx+i)], k);
3830 }
3831 }
3832 }
3833 }
3834 }
3835
3836
3837 user->multinullspace = PETSC_TRUE;
3838 for (j=1; j<my-1; j++) {
3839 for (i=1; i<mx-1; i++) {
3840 if (!Blocked[j*mx+i]) {
3841 user->multinullspace = PETSC_FALSE;
3842 break;
3843 }
3844 }
3845 }
3846 PetscFree(Blocked);
3847 VecRestoreArray3d(Zvert, mz, my, mx, 0, 0, 0, &nvert);
3848 MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD);
3849 if (user->multinullspace) {
3850 MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD);
3851
3852 }
3853 }
3854 else {
3855 MPI_Bcast(&user->multinullspace, 1, MPI_INT, 0, PETSC_COMM_WORLD);
3856 if (user->multinullspace) {
3857 MPI_Bcast(user->KSKE, 2*mx*my, MPI_INT, 0, PETSC_COMM_WORLD);
3858 }
3859 }
3860
3861
3862
3863 VecDestroy(&Zvert);
3864 return 0;
3865}
UserCtx::Nvert
Vec Nvert
Definition variables.h:657
Here is the caller graph for this function:

◆ MyNvertRestriction()

PetscErrorCode MyNvertRestriction ( UserCtx *  user_h,
UserCtx *  user_c 
)

Definition at line 3867 of file poisson.c.

3868{
3869 // DA da = user_c->da, fda = user_c->fda;
3870
3871
3872
3873 DMDALocalInfo info = user_c->info;
3874 PetscInt xs = info.xs, xe = info.xs + info.xm;
3875 PetscInt ys = info.ys, ye = info.ys + info.ym;
3876 PetscInt zs = info.zs, ze = info.zs + info.zm;
3877 PetscInt mx = info.mx, my = info.my, mz = info.mz;
3878
3879 PetscInt i,j,k;
3880 PetscInt ih, jh, kh, ia, ja, ka;
3881 PetscInt lxs, lxe, lys, lye, lzs, lze;
3882
3883 PetscReal ***nvert, ***nvert_h;
3884
3885 DMDAVecGetArray(user_h->da, user_h->lNvert, &nvert_h);
3886 DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert);
3887
3888 lxs = xs; lxe = xe;
3889 lys = ys; lye = ye;
3890 lzs = zs; lze = ze;
3891
3892 if (xs==0) lxs = xs+1;
3893 if (ys==0) lys = ys+1;
3894 if (zs==0) lzs = zs+1;
3895
3896 if (xe==mx) lxe = xe-1;
3897 if (ye==my) lye = ye-1;
3898 if (ze==mz) lze = ze-1;
3899
3900 if ((user_c->isc)) ia = 0;
3901 else ia = 1;
3902
3903 if ((user_c->jsc)) ja = 0;
3904 else ja = 1;
3905
3906 if ((user_c->ksc)) ka = 0;
3907 else ka = 1;
3908
3909 VecSet(user_c->Nvert, 0.);
3910 if (user_c->thislevel > 0) {
3911 for (k=lzs; k<lze; k++) {
3912 for (j=lys; j<lye; j++) {
3913 for (i=lxs; i<lxe; i++) {
3914 GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
3915 if (nvert_h[kh ][jh ][ih ] *
3916 nvert_h[kh ][jh ][ih-ia] *
3917 nvert_h[kh ][jh-ja][ih ] *
3918 nvert_h[kh-ka][jh ][ih ] *
3919 nvert_h[kh ][jh-ja][ih-ia] *
3920 nvert_h[kh-ka][jh ][ih-ia] *
3921 nvert_h[kh-ka][jh-ja][ih ] *
3922 nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
3923 nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
3924 }
3925 }
3926 }
3927 }
3928 }
3929 else {
3930 for (k=lzs; k<lze; k++) {
3931 for (j=lys; j<lye; j++) {
3932 for (i=lxs; i<lxe; i++) {
3933 GridRestriction(i, j, k, &ih, &jh, &kh, user_c);
3934 if (nvert_h[kh ][jh ][ih ] *
3935 nvert_h[kh ][jh ][ih-ia] *
3936 nvert_h[kh ][jh-ja][ih ] *
3937 nvert_h[kh-ka][jh ][ih ] *
3938 nvert_h[kh ][jh-ja][ih-ia] *
3939 nvert_h[kh-ka][jh ][ih-ia] *
3940 nvert_h[kh-ka][jh-ja][ih ] *
3941 nvert_h[kh-ka][jh-ja][ih-ia] > 0.1) {
3942 nvert[k][j][i] = PetscMax(1., nvert[k][j][i]);
3943 }
3944 }
3945 }
3946 }
3947 }
3948 DMDAVecRestoreArray(user_h->da, user_h->lNvert, &nvert_h);
3949 DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert);
3950
3951 DMGlobalToLocalBegin(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3952 DMGlobalToLocalEnd(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3953 //Mohsen Dec 2015
3954 DMDAVecGetArray(user_c->da, user_c->lNvert, &nvert);
3955 DMDAVecGetArray(user_c->da, user_c->Nvert, &nvert_h);
3956
3957 for (k=lzs; k<lze; k++) {
3958 for (j=lys; j<lye; j++) {
3959 for (i=lxs; i<lxe; i++) {
3960 if (nvert_h[k][j][i] < 0.1) {
3961 if (nvert[k][j][i+1] + nvert[k][j][i-1] > 1.1 &&
3962 nvert[k][j+1][i] + nvert[k][j-1][i] > 1.1 &&
3963 nvert[k+1][j][i] + nvert[k-1][j][i] > 1.1) {
3964 nvert_h[k][j][i] = 1.;
3965 }
3966 }
3967 }
3968 }
3969 }
3970
3971 DMDAVecRestoreArray(user_c->da, user_c->lNvert, &nvert);
3972 DMDAVecRestoreArray(user_c->da, user_c->Nvert, &nvert_h);
3973 DMGlobalToLocalBegin(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3974 DMGlobalToLocalEnd(user_c->da, user_c->Nvert, INSERT_VALUES, user_c->lNvert);
3975 /* DMLocalToGlobalBegin(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
3976/* DMLocalToGlobalEnd(user_c->da, user_c->lNvert, INSERT_VALUES, user_c->Nvert); */
3977 return 0;
3978}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PoissonSolver_MG()

PetscErrorCode PoissonSolver_MG ( UserMG *  usermg)

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

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

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

Definition at line 3982 of file poisson.c.

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