PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
rhs.c
Go to the documentation of this file.
1#include "rhs.h"
2
3#include "solvers.h" // Or a new "rhs.h"
4#include "logging.h"
5
6#include "rhs.h"
7#include "logging.h"
8
9// Forward declarations for the legacy helper functions if they are not in a header.
10// It is best practice to move these to rhs.h if they are in rhs.c.
11
12
13static void Calculate_Covariant_metrics(double g[3][3], double G[3][3])
14{
15 /*
16 | csi.x csi.y csi.z |-1 | x.csi x.eta x.zet |
17 | eta.x eta.y eta.z | = | y.csi y.eta y.zet |
18 | zet.x zet.y zet.z | | z.csi z.eta z.zet |
19
20 */
21 const double a11=g[0][0], a12=g[0][1], a13=g[0][2];
22 const double a21=g[1][0], a22=g[1][1], a23=g[1][2];
23 const double a31=g[2][0], a32=g[2][1], a33=g[2][2];
24
25 double det= a11*(a33*a22-a32*a23)-a21*(a33*a12-a32*a13)+a31*(a23*a12-a22*a13);
26
27 G[0][0] = (a33*a22-a32*a23)/det, G[0][1] = - (a33*a12-a32*a13)/det, G[0][2] = (a23*a12-a22*a13)/det;
28 G[1][0] = -(a33*a21-a31*a23)/det, G[1][1] = (a33*a11-a31*a13)/det, G[1][2] = - (a23*a11-a21*a13)/det;
29 G[2][0] = (a32*a21-a31*a22)/det, G[2][1] = - (a32*a11-a31*a12)/det, G[2][2] = (a22*a11-a21*a12)/det;
30};
31
32static void Calculate_normal_and_area(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
33{
34 double g[3][3];
35 double G[3][3];
36
37 g[0][0]=csi.x, g[0][1]=csi.y, g[0][2]=csi.z;
38 g[1][0]=eta.x, g[1][1]=eta.y, g[1][2]=eta.z;
39 g[2][0]=zet.x, g[2][1]=zet.y, g[2][2]=zet.z;
40
41 Calculate_Covariant_metrics(g, G);
42 double xcsi=G[0][0], ycsi=G[1][0], zcsi=G[2][0];
43 double xeta=G[0][1], yeta=G[1][1], zeta=G[2][1];
44 double xzet=G[0][2], yzet=G[1][2], zzet=G[2][2];
45
46 double nx_i = xcsi, ny_i = ycsi, nz_i = zcsi;
47 double nx_j = xeta, ny_j = yeta, nz_j = zeta;
48 double nx_k = xzet, ny_k = yzet, nz_k = zzet;
49
50 double sum_i=sqrt(nx_i*nx_i+ny_i*ny_i+nz_i*nz_i);
51 double sum_j=sqrt(nx_j*nx_j+ny_j*ny_j+nz_j*nz_j);
52 double sum_k=sqrt(nx_k*nx_k+ny_k*ny_k+nz_k*nz_k);
53
54 *Ai = sqrt( g[0][0]*g[0][0] + g[0][1]*g[0][1] + g[0][2]*g[0][2] ); // area
55 *Aj = sqrt( g[1][0]*g[1][0] + g[1][1]*g[1][1] + g[1][2]*g[1][2] );
56 *Ak =sqrt( g[2][0]*g[2][0] + g[2][1]*g[2][1] + g[2][2]*g[2][2] );
57
58 nx_i /= sum_i, ny_i /= sum_i, nz_i /= sum_i;
59 nx_j /= sum_j, ny_j /= sum_j, nz_j /= sum_j;
60 nx_k /= sum_k, ny_k /= sum_k, nz_k /= sum_k;
61
62 ni[0] = nx_i, ni[1] = ny_i, ni[2] = nz_i;
63 nj[0] = nx_j, nj[1] = ny_j, nj[2] = nz_j;
64 nk[0] = nx_k, nk[1] = ny_k, nk[2] = nz_k;
65}
66
67
68static void Calculate_dxdydz(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
69{
70 double ni[3], nj[3], nk[3];
71 double Li, Lj, Lk;
72 double Ai, Aj, Ak;
73 double vol = 1./ajc;
74
75 Calculate_normal_and_area(csi, eta, zet, ni, nj, nk, &Ai, &Aj, &Ak);
76 Li = vol / Ai;
77 Lj = vol / Aj;
78 Lk = vol / Ak;
79
80 // Length scale vector = di * ni_vector + dj * nj_vector + dk * nk_vector
81 *dx = fabs( Li * ni[0] + Lj * nj[0] + Lk * nk[0] );
82 *dy = fabs( Li * ni[1] + Lj * nj[1] + Lk * nk[1] );
83 *dz = fabs( Li * ni[2] + Lj * nj[2] + Lk * nk[2] );
84}
85
86
87PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
88{
89
90 // --- CONTEXT ACQUISITION BLOCK ---
91 // Get the master simulation context from the UserCtx.
92 SimCtx *simCtx = user->simCtx;
93
94 // Create local variables to mirror the legacy globals for minimal code changes.
95 const PetscInt les = simCtx->les;
96 const PetscInt central = simCtx->central; // Get this from SimCtx now
97 // --- END CONTEXT ACQUISITION BLOCK ---
98
99 Cmpnts ***ucont, ***ucat;
100 DM da = user->da, fda = user->fda;
101 DMDALocalInfo info;
102 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
103 PetscInt mx, my, mz; // Dimensions in three directions
104 PetscInt i, j, k;
105 Vec Fp1, Fp2, Fp3;
106 Cmpnts ***fp1, ***fp2, ***fp3;
107 Cmpnts ***conv;
108
109 PetscReal ucon, up, um;
110 PetscReal coef = 0.125, innerblank=7.;
111
112 PetscInt lxs, lxe, lys, lye, lzs, lze, gxs, gxe, gys, gye, gzs,gze;
113
114 PetscReal ***nvert,***aj;
115
116 DMDAGetLocalInfo(da, &info);
117 mx = info.mx; my = info.my; mz = info.mz;
118 xs = info.xs; xe = xs + info.xm;
119 ys = info.ys; ye = ys + info.ym;
120 zs = info.zs; ze = zs + info.zm;
121 gxs = info.gxs; gxe = gxs + info.gxm;
122 gys = info.gys; gye = gys + info.gym;
123 gzs = info.gzs; gze = gzs + info.gzm;
124
125 DMDAVecGetArray(fda, Ucont, &ucont);
126 DMDAVecGetArray(fda, Ucat, &ucat);
127 DMDAVecGetArray(fda, Conv, &conv);
128 DMDAVecGetArray(da, user->lAj, &aj);
129
130 VecDuplicate(Ucont, &Fp1);
131 VecDuplicate(Ucont, &Fp2);
132 VecDuplicate(Ucont, &Fp3);
133
134 DMDAVecGetArray(fda, Fp1, &fp1);
135 DMDAVecGetArray(fda, Fp2, &fp2);
136 DMDAVecGetArray(fda, Fp3, &fp3);
137
138 DMDAVecGetArray(da, user->lNvert, &nvert);
139
140
141 /* We have two different sets of node: 1. grid node, the physical points
142 where grid lines intercross; 2. storage node, where we store variables.
143 All node without explicitly specified as "grid node" refers to
144 storage node.
145
146 The integer node is defined at cell center while half node refers to
147 the actual grid node. (The reason to choose this arrangement is we need
148 ghost node, which is half node away from boundaries, to specify boundary
149 conditions. By using this storage arrangement, the actual storage need
150 is (IM+1) * (JM + 1) * (KM+1) where IM, JM, & KM refer to the number of
151 grid nodes along i, j, k directions.)
152
153 DA, the data structure used to define the storage of 3D arrays, is defined
154 as mx * my * mz. mx = IM+1, my = JM+1, mz = KM+1.
155
156 Staggered grid arrangement is used in this solver.
157 Pressure is stored at interger node (hence the cell center) and volume
158 fluxes defined on the center of each surface of a given control volume
159 is stored on the cloest upper integer node. */
160
161 /* First we calculate the flux on cell surfaces. Stored on the upper integer
162 node. For example, along i direction, the flux are stored at node 0:mx-2*/
163
164 lxs = xs; lxe = xe;
165 lys = ys; lye = ye;
166 lzs = zs; lze = ze;
167
168 if (xs==0) lxs = xs+1;
169 if (ys==0) lys = ys+1;
170 if (zs==0) lzs = zs+1;
171
172 if (xe==mx) lxe=xe-1;
173 if (ye==my) lye=ye-1;
174 if (ze==mz) lze=ze-1;
175
176 //Mohsen Sep 2012//
177/* First update the computational ghost points velocity for periodic boundary conditions
178 just for this subroutine because of Quick scheme for velocity deravatives */
179 /* for (k=zs; k<ze; k++) { */
180/* for (j=ys; j<ye; j++) { */
181/* for (i=xs; i<xe; i++) { */
182/* if (i==1 && (j==1) && (k==0 || k==1)) */
183/* PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d u is %.15le v is %.15le w is %.15le \n",i,j,k,ucat[k][j][i].x,ucat[k][j][i].y,ucat[k][j][i].z ); */
184/* } */
185/* } */
186/* } */
187 if (user->bctype[0]==7 || user->bctype[1]==7){
188 if (xs==0){
189 i=xs-1;
190 for (k=gzs; k<gze; k++) {
191 for (j=gys; j<gye; j++) {
192 ucat[k][j][i].x=ucat[k][j][i-2].x;
193 ucat[k][j][i].y=ucat[k][j][i-2].y;
194 ucat[k][j][i].z=ucat[k][j][i-2].z;
195 nvert[k][j][i]=nvert[k][j][i-2];
196 }
197 }
198 }
199 if (xe==mx){
200 i=mx;
201 for (k=gzs; k<gze; k++) {
202 for (j=gys; j<gye; j++) {
203 ucat[k][j][i].x=ucat[k][j][i+2].x;
204 ucat[k][j][i].y=ucat[k][j][i+2].y;
205 ucat[k][j][i].z=ucat[k][j][i+2].z;
206 nvert[k][j][i]=nvert[k][j][i+2];
207 }
208 }
209 }
210 }
211 if (user->bctype[2]==7 || user->bctype[3]==7){
212 if (ys==0){
213 j=ys-1;
214 for (k=gzs; k<gze; k++) {
215 for (i=gxs; i<gxe; i++) {
216 ucat[k][j][i].x=ucat[k][j-2][i].x;
217 ucat[k][j][i].y=ucat[k][j-2][i].y;
218 ucat[k][j][i].z=ucat[k][j-2][i].z;
219 nvert[k][j][i]=nvert[k][j-2][i];
220 }
221 }
222 }
223 if (ye==my){
224 j=my;
225 for (k=gzs; k<gze; k++) {
226 for (i=gxs; i<gxe; i++) {
227 ucat[k][j][i].x=ucat[k][j+2][i].x;
228 ucat[k][j][i].y=ucat[k][j+2][i].y;
229 ucat[k][j][i].z=ucat[k][j+2][i].z;
230 nvert[k][j][i]=nvert[k][j+2][i];
231 }
232 }
233 }
234 }
235 if (user->bctype[4]==7 || user->bctype[5]==7){
236 if (zs==0){
237 k=zs-1;
238 for (j=gys; j<gye; j++) {
239 for (i=gxs; i<gxe; i++) {
240 ucat[k][j][i].x=ucat[k-2][j][i].x;
241 ucat[k][j][i].y=ucat[k-2][j][i].y;
242 ucat[k][j][i].z=ucat[k-2][j][i].z;
243 nvert[k][j][i]=nvert[k-2][j][i];
244 }
245 }
246 }
247 if (ze==mz){
248 k=mz;
249 for (j=gys; j<gye; j++) {
250 for (i=gxs; i<gxe; i++) {
251 ucat[k][j][i].x=ucat[k+2][j][i].x;
252 ucat[k][j][i].y=ucat[k+2][j][i].y;
253 ucat[k][j][i].z=ucat[k+2][j][i].z;
254 nvert[k][j][i]=nvert[k+2][j][i];
255 }
256 }
257 }
258 }
259
260 VecSet(Conv, 0.0);
261
262 /* Calculating the convective terms on cell centers.
263 First calcualte the contribution from i direction
264 The flux is evaluated by QUICK scheme */
265
266 for (k=lzs; k<lze; k++){
267 for (j=lys; j<lye; j++){
268 for (i=lxs-1; i<lxe; i++){
269
270
271 ucon = ucont[k][j][i].x * 0.5;
272
273 up = ucon + fabs(ucon);
274 um = ucon - fabs(ucon);
275
276 if (i>0 && i<mx-2 &&
277 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
278 (nvert[k][j][i-1] < 0.1 || nvert[k][j][i-1]>innerblank)) { // interial nodes
279 if ((les || central)) {
280 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
281 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
282 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
283
284 } else {
285 fp1[k][j][i].x =
286 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
287 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
288 fp1[k][j][i].y =
289 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
290 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
291 fp1[k][j][i].z =
292 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
293 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
294 }
295 }
296 else if ((les || central) && (i==0 || i==mx-2) &&
297 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
298 (nvert[k][j][i ] < 0.1 || nvert[k][j][i ]>innerblank))
299 {
300 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
301 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
302 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
303 }
304 else if (i==0 ||(nvert[k][j][i-1] > 0.1) ) {
305 if (user->bctype[0]==7 && i==0 && (nvert[k][j][i-1]<0.1 && nvert[k][j][i+1]<0.1)){//Mohsen Feb 12
306 fp1[k][j][i].x =
307 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
308 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
309 fp1[k][j][i].y =
310 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
311 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
312 fp1[k][j][i].z =
313 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
314 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
315 }else{
316 fp1[k][j][i].x =
317 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
318 up * (coef * (-ucat[k][j][i ].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
319 fp1[k][j][i].y =
320 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
321 up * (coef * (-ucat[k][j][i ].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
322 fp1[k][j][i].z =
323 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
324 up * (coef * (-ucat[k][j][i ].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
325 }
326 }
327 else if (i==mx-2 ||(nvert[k][j][i+1]) > 0.1) {
328 if (user->bctype[0]==7 && i==mx-2 &&(nvert[k][j][i-1]<0.1 && nvert[k][j][i+1]<0.1)){//Mohsen Feb 12
329 fp1[k][j][i].x =
330 um * (coef * (-ucat[k][j][i+2].x -2.* ucat[k][j][i+1].x +3.* ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
331 up * (coef * (-ucat[k][j][i-1].x -2.* ucat[k][j][i ].x +3.* ucat[k][j][i+1].x) +ucat[k][j][i ].x);
332 fp1[k][j][i].y =
333 um * (coef * (-ucat[k][j][i+2].y -2.* ucat[k][j][i+1].y +3.* ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
334 up * (coef * (-ucat[k][j][i-1].y -2.* ucat[k][j][i ].y +3.* ucat[k][j][i+1].y) +ucat[k][j][i ].y);
335 fp1[k][j][i].z =
336 um * (coef * (-ucat[k][j][i+2].z -2.* ucat[k][j][i+1].z +3.* ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
337 up * (coef * (-ucat[k][j][i-1].z -2.* ucat[k][j][i ].z +3.* ucat[k][j][i+1].z) +ucat[k][j][i ].z);
338 }else{
339 fp1[k][j][i].x =
340 um * (coef * (-ucat[k][j][i+1].x -2. * ucat[k][j][i+1].x +3. * ucat[k][j][i ].x) +ucat[k][j][i+1].x) +
341 up * (coef * (-ucat[k][j][i-1].x -2. * ucat[k][j][i ].x +3. * ucat[k][j][i+1].x) +ucat[k][j][i ].x);
342 fp1[k][j][i].y =
343 um * (coef * (-ucat[k][j][i+1].y -2. * ucat[k][j][i+1].y +3. * ucat[k][j][i ].y) +ucat[k][j][i+1].y) +
344 up * (coef * (-ucat[k][j][i-1].y -2. * ucat[k][j][i ].y +3. * ucat[k][j][i+1].y) +ucat[k][j][i ].y);
345 fp1[k][j][i].z =
346 um * (coef * (-ucat[k][j][i+1].z -2. * ucat[k][j][i+1].z +3. * ucat[k][j][i ].z) +ucat[k][j][i+1].z) +
347 up * (coef * (-ucat[k][j][i-1].z -2. * ucat[k][j][i ].z +3. * ucat[k][j][i+1].z) +ucat[k][j][i ].z);
348 }
349 }
350 }
351 }
352 }
353
354 /* j direction */
355 for (k=lzs; k<lze; k++) {
356 for(j=lys-1; j<lye; j++) {
357 for(i=lxs; i<lxe; i++) {
358 ucon = ucont[k][j][i].y * 0.5;
359
360 up = ucon + fabs(ucon);
361 um = ucon - fabs(ucon);
362
363 if (j>0 && j<my-2 &&
364 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i] > innerblank) &&
365 (nvert[k][j-1][i] < 0.1 || nvert[k][j-1][i] > innerblank)) {
366 if ((les || central)) {
367 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
368 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
369 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
370
371 } else {
372 fp2[k][j][i].x =
373 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
374 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
375 fp2[k][j][i].y =
376 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
377 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
378 fp2[k][j][i].z =
379 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
380 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
381 }
382 }
383 else if ((les || central) && (j==0 || i==my-2) &&
384 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i]>innerblank) &&
385 (nvert[k][j ][i] < 0.1 || nvert[k][j ][i]>innerblank))
386 {
387 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
388 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
389 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
390 }
391 else if (j==0 || (nvert[k][j-1][i]) > 0.1) {
392 if (user->bctype[2]==7 && j==0 && (nvert[k][j-1][i]<0.1 && nvert[k][j+1][i]<0.1 )){//Mohsen Feb 12 //
393 fp2[k][j][i].x =
394 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
395 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
396 fp2[k][j][i].y =
397 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
398 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
399 fp2[k][j][i].z =
400 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
401 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
402 }else{
403 fp2[k][j][i].x =
404 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
405 up * (coef * (-ucat[k][j ][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
406 fp2[k][j][i].y =
407 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
408 up * (coef * (-ucat[k][j ][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
409 fp2[k][j][i].z =
410 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
411 up * (coef * (-ucat[k][j ][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
412 }
413 }
414 else if (j==my-2 ||(nvert[k][j+1][i]) > 0.1) {
415 if (user->bctype[2]==7 && j==my-2 && (nvert[k][j-1][i]<0.1 && nvert[k][j+1][i]<0.1 )){//Mohsen Feb 12//
416 fp2[k][j][i].x =
417 um * (coef * (-ucat[k][j+2][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
418 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
419 fp2[k][j][i].y =
420 um * (coef * (-ucat[k][j+2][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
421 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j ][i].y);
422 fp2[k][j][i].z =
423 um * (coef * (-ucat[k][j+2][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
424 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j ][i].z);
425 }else{
426 fp2[k][j][i].x =
427 um * (coef * (-ucat[k][j+1][i].x -2. * ucat[k][j+1][i].x +3. * ucat[k][j ][i].x) +ucat[k][j+1][i].x) +
428 up * (coef * (-ucat[k][j-1][i].x -2. * ucat[k][j ][i].x +3. * ucat[k][j+1][i].x) +ucat[k][j ][i].x);
429 fp2[k][j][i].y =
430 um * (coef * (-ucat[k][j+1][i].y -2. * ucat[k][j+1][i].y +3. * ucat[k][j ][i].y) +ucat[k][j+1][i].y) +
431 up * (coef * (-ucat[k][j-1][i].y -2. * ucat[k][j ][i].y +3. * ucat[k][j+1][i].y) +ucat[k][j][i ].y);
432 fp2[k][j][i].z =
433 um * (coef * (-ucat[k][j+1][i].z -2. * ucat[k][j+1][i].z +3. * ucat[k][j ][i].z) +ucat[k][j+1][i].z) +
434 up * (coef * (-ucat[k][j-1][i].z -2. * ucat[k][j ][i].z +3. * ucat[k][j+1][i].z) +ucat[k][j][i ].z);
435 }
436 }
437 }
438 }
439 }
440
441
442 /* k direction */
443 for (k=lzs-1; k<lze; k++) {
444 for(j=lys; j<lye; j++) {
445 for(i=lxs; i<lxe; i++) {
446 ucon = ucont[k][j][i].z * 0.5;
447
448 up = ucon + fabs(ucon);
449 um = ucon - fabs(ucon);
450
451 if (k>0 && k<mz-2 &&
452 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i] > innerblank) &&
453 (nvert[k-1][j][i] < 0.1 || nvert[k-1][j][i] > innerblank)) {
454 if ((les || central)) {
455 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
456 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
457 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
458
459 } else {
460 fp3[k][j][i].x =
461 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
462 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
463 fp3[k][j][i].y =
464 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
465 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
466 fp3[k][j][i].z =
467 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
468 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
469 }
470 }
471 else if ((les || central) && (k==0 || k==mz-2) &&
472 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i]>innerblank) &&
473 (nvert[k ][j][i] < 0.1 || nvert[k ][j][i]>innerblank))
474 {
475 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
476 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
477 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
478 }
479 else if (k<mz-2 && (k==0 ||(nvert[k-1][j][i]) > 0.1)) {
480 if(user->bctype[4]==7 && k==0 && (nvert[k-1][j][i]<0.1 && nvert[k+1][j][i]<0.1)){//Mohsen Feb 12//
481 fp3[k][j][i].x =
482 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
483 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
484 fp3[k][j][i].y =
485 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
486 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
487 fp3[k][j][i].z =
488 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
489 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
490 }else{
491 fp3[k][j][i].x =
492 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
493 up * (coef * (-ucat[k ][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k][j][i ].x);
494 fp3[k][j][i].y =
495 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
496 up * (coef * (-ucat[k ][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k][j][i ].y);
497 fp3[k][j][i].z =
498 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
499 up * (coef * (-ucat[k ][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k][j][i ].z);
500 }
501 }
502 else if (k>0 && (k==mz-2 ||(nvert[k+1][j][i]) > 0.1)) {
503 if (user->bctype[4]==7 && k==mz-2 && (nvert[k-1][j][i]<0.1 && nvert[k+1][j][i]<0.1)){//Mohsen Feb 12//
504 fp3[k][j][i].x =
505 um * (coef * (-ucat[k+2][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
506 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k ][j][i].x);
507 fp3[k][j][i].y =
508 um * (coef * (-ucat[k+2][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
509 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k ][j][i].y);
510 fp3[k][j][i].z =
511 um * (coef * (-ucat[k+2][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
512 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k ][j][i].z);
513 }else{
514 fp3[k][j][i].x =
515 um * (coef * (-ucat[k+1][j][i].x -2. * ucat[k+1][j][i].x +3. * ucat[k ][j][i].x) +ucat[k+1][j][i].x) +
516 up * (coef * (-ucat[k-1][j][i].x -2. * ucat[k ][j][i].x +3. * ucat[k+1][j][i].x) +ucat[k][j][i ].x);
517 fp3[k][j][i].y =
518 um * (coef * (-ucat[k+1][j][i].y -2. * ucat[k+1][j][i].y +3. * ucat[k ][j][i].y) +ucat[k+1][j][i].y) +
519 up * (coef * (-ucat[k-1][j][i].y -2. * ucat[k ][j][i].y +3. * ucat[k+1][j][i].y) +ucat[k][j][i ].y);
520 fp3[k][j][i].z =
521 um * (coef * (-ucat[k+1][j][i].z -2. * ucat[k+1][j][i].z +3. * ucat[k ][j][i].z) +ucat[k+1][j][i].z) +
522 up * (coef * (-ucat[k-1][j][i].z -2. * ucat[k ][j][i].z +3. * ucat[k+1][j][i].z) +ucat[k][j][i ].z);
523 }
524 }
525 }
526 }
527 }
528
529 /* Calculate the convective terms under cartesian coordinates */
530
531 for (k=lzs; k<lze; k++) {
532 for (j=lys; j<lye; j++) {
533 for (i=lxs; i<lxe; i++) {
534 conv[k][j][i].x =
535 fp1[k][j][i].x - fp1[k][j][i-1].x +
536 fp2[k][j][i].x - fp2[k][j-1][i].x +
537 fp3[k][j][i].x - fp3[k-1][j][i].x;
538
539 conv[k][j][i].y =
540 fp1[k][j][i].y - fp1[k][j][i-1].y +
541 fp2[k][j][i].y - fp2[k][j-1][i].y +
542 fp3[k][j][i].y - fp3[k-1][j][i].y;
543
544 conv[k][j][i].z =
545 fp1[k][j][i].z - fp1[k][j][i-1].z +
546 fp2[k][j][i].z - fp2[k][j-1][i].z +
547 fp3[k][j][i].z - fp3[k-1][j][i].z;
548 }
549 }
550 }
551 /* for (k=zs; k<ze; k++) { */
552/* for (j=ys; j<ye; j++) { */
553/* for (i=xs; i<xe; i++) { */
554/* if (i==1 && (j==1) && (k==1 || k==21 || k==22|| k==200)) */
555/* PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d conv.y is %.15le conv.z is %.15le \n",i,j,k,conv[k][j][i].y,conv[k][j][i].z); */
556/* } */
557/* } */
558/* } */
559
560 DMDAVecRestoreArray(fda, Ucont, &ucont);
561 DMDAVecRestoreArray(fda, Ucat, &ucat);
562 DMDAVecRestoreArray(fda, Conv, &conv);
563 DMDAVecRestoreArray(da, user->lAj, &aj);
564
565 DMDAVecRestoreArray(fda, Fp1, &fp1);
566 DMDAVecRestoreArray(fda, Fp2, &fp2);
567 DMDAVecRestoreArray(fda, Fp3, &fp3);
568 DMDAVecRestoreArray(da, user->lNvert, &nvert);
569
570 VecDestroy(&Fp1);
571 VecDestroy(&Fp2);
572 VecDestroy(&Fp3);
573
574
575 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Convective term calculated .\n");
576
577 return (0);
578}
579
580PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
581{
582
583 Vec Csi = user->lCsi, Eta = user->lEta, Zet = user->lZet;
584
585 Cmpnts ***ucont, ***ucat;
586
587 Cmpnts ***csi, ***eta, ***zet;
588 Cmpnts ***icsi, ***ieta, ***izet;
589 Cmpnts ***jcsi, ***jeta, ***jzet;
590 Cmpnts ***kcsi, ***keta, ***kzet;
591
592 PetscReal ***nvert;
593
594 DM da = user->da, fda = user->fda;
595 DMDALocalInfo info;
596 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
597 PetscInt mx, my, mz; // Dimensions in three directions
598 PetscInt i, j, k;
599 Vec Fp1, Fp2, Fp3;
600 Cmpnts ***fp1, ***fp2, ***fp3;
601 Cmpnts ***visc;
602 PetscReal ***aj, ***iaj, ***jaj, ***kaj;
603
604 PetscInt lxs, lxe, lys, lye, lzs, lze;
605
606 PetscReal ajc;
607
608 PetscReal dudc, dude, dudz, dvdc, dvde, dvdz, dwdc, dwde, dwdz;
609 PetscReal csi0, csi1, csi2, eta0, eta1, eta2, zet0, zet1, zet2;
610 PetscReal g11, g21, g31;
611 PetscReal r11, r21, r31, r12, r22, r32, r13, r23, r33;
612
613 PetscScalar solid,innerblank;
614
615 // --- CONTEXT ACQUISITION BLOCK ---
616 // Get the master simulation context from the UserCtx.
617 SimCtx *simCtx = user->simCtx;
618
619 // Create local variables to mirror the legacy globals for minimal code changes.
620 const PetscInt les = simCtx->les;
621 const PetscInt rans = simCtx->rans;
622 const PetscInt ti = simCtx->step; // Assuming simCtx->step is the new integer time counter
623 const PetscReal ren = simCtx->ren;
624 const PetscInt clark = simCtx->clark;
625 solid = 0.5;
626 innerblank = 7.;
627
628 DMDAVecGetArray(fda, Ucont, &ucont);
629 DMDAVecGetArray(fda, Ucat, &ucat);
630 DMDAVecGetArray(fda, Visc, &visc);
631
632 DMDAVecGetArray(fda, Csi, &csi);
633 DMDAVecGetArray(fda, Eta, &eta);
634 DMDAVecGetArray(fda, Zet, &zet);
635
636 DMDAVecGetArray(fda, user->lICsi, &icsi);
637 DMDAVecGetArray(fda, user->lIEta, &ieta);
638 DMDAVecGetArray(fda, user->lIZet, &izet);
639
640 DMDAVecGetArray(fda, user->lJCsi, &jcsi);
641 DMDAVecGetArray(fda, user->lJEta, &jeta);
642 DMDAVecGetArray(fda, user->lJZet, &jzet);
643
644 DMDAVecGetArray(fda, user->lKCsi, &kcsi);
645 DMDAVecGetArray(fda, user->lKEta, &keta);
646 DMDAVecGetArray(fda, user->lKZet, &kzet);
647
648 DMDAVecGetArray(da, user->lNvert, &nvert);
649
650 VecDuplicate(Ucont, &Fp1);
651 VecDuplicate(Ucont, &Fp2);
652 VecDuplicate(Ucont, &Fp3);
653
654 DMDAVecGetArray(fda, Fp1, &fp1);
655 DMDAVecGetArray(fda, Fp2, &fp2);
656 DMDAVecGetArray(fda, Fp3, &fp3);
657
658 DMDAVecGetArray(da, user->lAj, &aj);
659
660 DMDAGetLocalInfo(da, &info);
661
662 mx = info.mx; my = info.my; mz = info.mz;
663 xs = info.xs; xe = xs + info.xm;
664 ys = info.ys; ye = ys + info.ym;
665 zs = info.zs; ze = zs + info.zm;
666
667 /* First we calculate the flux on cell surfaces. Stored on the upper integer
668 node. For example, along i direction, the flux are stored at node 0:mx-2*/
669 lxs = xs; lxe = xe;
670 lys = ys; lye = ye;
671 lzs = zs; lze = ze;
672
673 if (xs==0) lxs = xs+1;
674 if (ys==0) lys = ys+1;
675 if (zs==0) lzs = zs+1;
676
677
678 if (xe==mx) lxe=xe-1;
679 if (ye==my) lye=ye-1;
680 if (ze==mz) lze=ze-1;
681
682 VecSet(Visc,0.0);
683
684 PetscReal ***lnu_t;
685
686 if(les) {
687 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
688 } else if (rans) {
689
690 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
691 }
692
693 /* The visc flux on each surface center is stored at previous integer node */
694
695 DMDAVecGetArray(da, user->lIAj, &iaj);
696 /* for (k=zs; k<ze; k++) { */
697/* for (j=ys; j<ye; j++) { */
698/* for (i=xs; i<xe; i++) { */
699/* if (i==1 && (j==0 ||j==1 || j==2) && (k==21 || k==22|| k==20)) */
700/* PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d u is %.15le v is %.15le w is %.15le \n",i,j,k,ucat[k][j][i].x,ucat[k][j][i].y,ucat[k][j][i].z ); */
701/* } */
702/* } */
703/* } */
704 // i direction
705 for (k=lzs; k<lze; k++) {
706 for (j=lys; j<lye; j++) {
707 for (i=lxs-1; i<lxe; i++) {
708
709 dudc = ucat[k][j][i+1].x - ucat[k][j][i].x;
710 dvdc = ucat[k][j][i+1].y - ucat[k][j][i].y;
711 dwdc = ucat[k][j][i+1].z - ucat[k][j][i].z;
712
713 if ((nvert[k][j+1][i ]> solid && nvert[k][j+1][i ]<innerblank) ||
714 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
715 dude = (ucat[k][j ][i+1].x + ucat[k][j ][i].x -
716 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.5;
717 dvde = (ucat[k][j ][i+1].y + ucat[k][j ][i].y -
718 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.5;
719 dwde = (ucat[k][j ][i+1].z + ucat[k][j ][i].z -
720 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.5;
721 }
722 else if ((nvert[k][j-1][i ]> solid && nvert[k][j-1][i ]<innerblank) ||
723 (nvert[k][j-1][i+1]> solid && nvert[k][j-1][i+1]<innerblank)) {
724 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
725 ucat[k][j ][i+1].x - ucat[k][j ][i].x) * 0.5;
726 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
727 ucat[k][j ][i+1].y - ucat[k][j ][i].y) * 0.5;
728 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
729 ucat[k][j ][i+1].z - ucat[k][j ][i].z) * 0.5;
730 }
731 else {
732 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
733 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.25;
734 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
735 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.25;
736 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
737 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.25;
738 }
739
740 if ((nvert[k+1][j][i ]> solid && nvert[k+1][j][i ]<innerblank)||
741 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
742 dudz = (ucat[k ][j][i+1].x + ucat[k ][j][i].x -
743 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.5;
744 dvdz = (ucat[k ][j][i+1].y + ucat[k ][j][i].y -
745 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.5;
746 dwdz = (ucat[k ][j][i+1].z + ucat[k ][j][i].z -
747 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.5;
748 }
749 else if ((nvert[k-1][j][i ]> solid && nvert[k-1][j][i ]<innerblank) ||
750 (nvert[k-1][j][i+1]> solid && nvert[k-1][j][i+1]<innerblank)) {
751
752 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
753 ucat[k ][j][i+1].x - ucat[k ][j][i].x) * 0.5;
754 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
755 ucat[k ][j][i+1].y - ucat[k ][j][i].y) * 0.5;
756 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
757 ucat[k ][j][i+1].z - ucat[k ][j][i].z) * 0.5;
758 }
759 else {
760 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
761 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.25;
762 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
763 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.25;
764 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
765 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.25;
766 }
767
768 csi0 = icsi[k][j][i].x;
769 csi1 = icsi[k][j][i].y;
770 csi2 = icsi[k][j][i].z;
771
772 eta0 = ieta[k][j][i].x;
773 eta1 = ieta[k][j][i].y;
774 eta2 = ieta[k][j][i].z;
775
776 zet0 = izet[k][j][i].x;
777 zet1 = izet[k][j][i].y;
778 zet2 = izet[k][j][i].z;
779
780 g11 = csi0 * csi0 + csi1 * csi1 + csi2 * csi2;
781 g21 = eta0 * csi0 + eta1 * csi1 + eta2 * csi2;
782 g31 = zet0 * csi0 + zet1 * csi1 + zet2 * csi2;
783
784 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
785 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
786 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
787
788 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
789 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
790 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
791
792 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
793 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
794 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
795
796 ajc = iaj[k][j][i];
797
798 double nu = 1./ren, nu_t=0;
799
800 if( les || (rans && ti>0) ) {
801 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j][i+1]) ), 2.0) * Sabs;
802 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j][i+1]);
803 if ( (user->bctype[0]==1 && i==0) || (user->bctype[1]==1 && i==mx-2) ) nu_t=0;
804 fp1[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * csi0 + r21 * csi1 + r31 * csi2) * ajc * (nu_t);
805 fp1[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * csi0 + r22 * csi1 + r32 * csi2) * ajc * (nu_t);
806 fp1[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * csi0 + r23 * csi1 + r33 * csi2) * ajc * (nu_t);
807 }
808 else {
809 fp1[k][j][i].x = 0;
810 fp1[k][j][i].y = 0;
811 fp1[k][j][i].z = 0;
812 }
813
814 fp1[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * csi0 + r21 * csi1 + r31 * csi2 ) * ajc * (nu);
815 fp1[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * csi0 + r22 * csi1 + r32 * csi2 ) * ajc * (nu);
816 fp1[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * csi0 + r23 * csi1 + r33 * csi2 ) * ajc * (nu);
817
818
819 if(clark) {
820 double dc, de, dz;
821 Calculate_dxdydz (ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
822 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
823
824 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
825 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
826 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
827 double t21 = t12;
828 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
829 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
830 double t31 = t13;
831 double t32 = t23;
832 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
833
834 fp1[k][j][i].x -= ( t11 * csi0 + t12 * csi1 + t13 * csi2 ) / 12.;
835 fp1[k][j][i].y -= ( t21 * csi0 + t22 * csi1 + t23 * csi2 ) / 12.;
836 fp1[k][j][i].z -= ( t31 * csi0 + t32 * csi1 + t33 * csi2 ) / 12.;
837 }
838
839 }
840 }
841 }
842 DMDAVecRestoreArray(da, user->lIAj, &iaj);
843
844
845 // j direction
846 DMDAVecGetArray(da, user->lJAj, &jaj);
847 for (k=lzs; k<lze; k++) {
848 for (j=lys-1; j<lye; j++) {
849 for (i=lxs; i<lxe; i++) {
850
851 if ((nvert[k][j ][i+1]> solid && nvert[k][j ][i+1]<innerblank)||
852 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
853 dudc = (ucat[k][j+1][i ].x + ucat[k][j][i ].x -
854 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.5;
855 dvdc = (ucat[k][j+1][i ].y + ucat[k][j][i ].y -
856 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.5;
857 dwdc = (ucat[k][j+1][i ].z + ucat[k][j][i ].z -
858 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.5;
859 }
860 else if ((nvert[k][j ][i-1]> solid && nvert[k][j ][i-1]<innerblank) ||
861 (nvert[k][j+1][i-1]> solid && nvert[k][j+1][i-1]<innerblank)) {
862 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
863 ucat[k][j+1][i ].x - ucat[k][j][i ].x) * 0.5;
864 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
865 ucat[k][j+1][i ].y - ucat[k][j][i ].y) * 0.5;
866 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
867 ucat[k][j+1][i ].z - ucat[k][j][i ].z) * 0.5;
868 }
869 else {
870 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
871 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.25;
872 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
873 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.25;
874 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
875 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.25;
876 }
877
878 dude = ucat[k][j+1][i].x - ucat[k][j][i].x;
879 dvde = ucat[k][j+1][i].y - ucat[k][j][i].y;
880 dwde = ucat[k][j+1][i].z - ucat[k][j][i].z;
881
882 if ((nvert[k+1][j ][i]> solid && nvert[k+1][j ][i]<innerblank)||
883 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
884 dudz = (ucat[k ][j+1][i].x + ucat[k ][j][i].x -
885 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.5;
886 dvdz = (ucat[k ][j+1][i].y + ucat[k ][j][i].y -
887 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.5;
888 dwdz = (ucat[k ][j+1][i].z + ucat[k ][j][i].z -
889 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.5;
890 }
891 else if ((nvert[k-1][j ][i]> solid && nvert[k-1][j ][i]<innerblank)||
892 (nvert[k-1][j+1][i]> solid && nvert[k-1][j+1][i]<innerblank)) {
893 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
894 ucat[k ][j+1][i].x - ucat[k ][j][i].x) * 0.5;
895 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
896 ucat[k ][j+1][i].y - ucat[k ][j][i].y) * 0.5;
897 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
898 ucat[k ][j+1][i].z - ucat[k ][j][i].z) * 0.5;
899 }
900 else {
901 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
902 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.25;
903 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
904 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.25;
905 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
906 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.25;
907 }
908
909 csi0 = jcsi[k][j][i].x;
910 csi1 = jcsi[k][j][i].y;
911 csi2 = jcsi[k][j][i].z;
912
913 eta0 = jeta[k][j][i].x;
914 eta1 = jeta[k][j][i].y;
915 eta2 = jeta[k][j][i].z;
916
917 zet0 = jzet[k][j][i].x;
918 zet1 = jzet[k][j][i].y;
919 zet2 = jzet[k][j][i].z;
920
921
922 g11 = csi0 * eta0 + csi1 * eta1 + csi2 * eta2;
923 g21 = eta0 * eta0 + eta1 * eta1 + eta2 * eta2;
924 g31 = zet0 * eta0 + zet1 * eta1 + zet2 * eta2;
925
926 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
927 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
928 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
929
930 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
931 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
932 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
933
934 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
935 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
936 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
937
938 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d dvdc is %.15le dvde is %.15le dvdz is %.15le \n",i,j,k,dvdc,dvde,dvdz);
939 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d dwdc is %.15le dwde is %.15le dwdz is %.15le \n",i,j,k,dwdc,dwde,dwdz);
940 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d jcsi is %.15le jeta is %.15le jzet is %.15le \n",i,j,k,jcsi[k][j][i].z,jeta[k][j][i].z,jzet[k][j][i].z);
941 // if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i=%d j=%d k=%d r13 is %.15le r23 is %.15le r33 is %.15le \n",i,j,k,r13,r23,r33);
942
943
944
945 ajc = jaj[k][j][i];
946
947 double nu = 1./ren, nu_t = 0;
948
949 if( les || (rans && ti>0) ) {
950 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j+1][i]) ), 2.0) * Sabs;
951 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j+1][i]);
952 if ( (user->bctype[2]==1 && j==0) || (user->bctype[3]==1 && j==my-2) ) nu_t=0;
953
954 fp2[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * eta0 + r21 * eta1 + r31 * eta2) * ajc * (nu_t);
955 fp2[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * eta0 + r22 * eta1 + r32 * eta2) * ajc * (nu_t);
956 fp2[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * eta0 + r23 * eta1 + r33 * eta2) * ajc * (nu_t);
957 }
958 else {
959 fp2[k][j][i].x = 0;
960 fp2[k][j][i].y = 0;
961 fp2[k][j][i].z = 0;
962 }
963
964 fp2[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * eta0 + r21 * eta1 + r31 * eta2 ) * ajc * (nu);
965 fp2[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * eta0 + r22 * eta1 + r32 * eta2 ) * ajc * (nu);
966 fp2[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * eta0 + r23 * eta1 + r33 * eta2 ) * ajc * (nu);
967
968 if(clark) {
969 double dc, de, dz;
970 Calculate_dxdydz(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
971 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
972
973 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
974 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
975 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
976 double t21 = t12;
977 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
978 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
979 double t31 = t13;
980 double t32 = t23;
981 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
982
983 fp2[k][j][i].x -= ( t11 * eta0 + t12 * eta1 + t13 * eta2 ) / 12.;
984 fp2[k][j][i].y -= ( t21 * eta0 + t22 * eta1 + t23 * eta2 ) / 12.;
985 fp2[k][j][i].z -= ( t31 * eta0 + t32 * eta1 + t33 * eta2 ) / 12.;
986 }
987 }
988 }
989 }
990
991 DMDAVecRestoreArray(da, user->lJAj, &jaj);
992 // k direction
993
994 DMDAVecGetArray(da, user->lKAj, &kaj);
995 for (k=lzs-1; k<lze; k++) {
996 for (j=lys; j<lye; j++) {
997 for (i=lxs; i<lxe; i++) {
998 if ((nvert[k ][j][i+1]> solid && nvert[k ][j][i+1]<innerblank)||
999 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
1000 dudc = (ucat[k+1][j][i ].x + ucat[k][j][i ].x -
1001 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.5;
1002 dvdc = (ucat[k+1][j][i ].y + ucat[k][j][i ].y -
1003 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.5;
1004 dwdc = (ucat[k+1][j][i ].z + ucat[k][j][i ].z -
1005 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.5;
1006 }
1007 else if ((nvert[k ][j][i-1]> solid && nvert[k ][j][i-1]<innerblank) ||
1008 (nvert[k+1][j][i-1]> solid && nvert[k+1][j][i-1]<innerblank)) {
1009 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
1010 ucat[k+1][j][i ].x - ucat[k][j][i ].x) * 0.5;
1011 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
1012 ucat[k+1][j][i ].y - ucat[k][j][i ].y) * 0.5;
1013 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
1014 ucat[k+1][j][i ].z - ucat[k][j][i ].z) * 0.5;
1015 }
1016 else {
1017 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
1018 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.25;
1019 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
1020 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.25;
1021 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
1022 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.25;
1023 }
1024
1025 if ((nvert[k ][j+1][i]> solid && nvert[k ][j+1][i]<innerblank)||
1026 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
1027 dude = (ucat[k+1][j ][i].x + ucat[k][j ][i].x -
1028 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.5;
1029 dvde = (ucat[k+1][j ][i].y + ucat[k][j ][i].y -
1030 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.5;
1031 dwde = (ucat[k+1][j ][i].z + ucat[k][j ][i].z -
1032 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.5;
1033 }
1034 else if ((nvert[k ][j-1][i]> solid && nvert[k ][j-1][i]<innerblank) ||
1035 (nvert[k+1][j-1][i]> solid && nvert[k+1][j-1][i]<innerblank)){
1036 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
1037 ucat[k+1][j ][i].x - ucat[k][j ][i].x) * 0.5;
1038 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
1039 ucat[k+1][j ][i].y - ucat[k][j ][i].y) * 0.5;
1040 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
1041 ucat[k+1][j ][i].z - ucat[k][j ][i].z) * 0.5;
1042 }
1043 else {
1044 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
1045 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.25;
1046 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
1047 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.25;
1048 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
1049 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.25;
1050 }
1051
1052 dudz = ucat[k+1][j][i].x - ucat[k][j][i].x;
1053 dvdz = ucat[k+1][j][i].y - ucat[k][j][i].y;
1054 dwdz = ucat[k+1][j][i].z - ucat[k][j][i].z;
1055
1056
1057 csi0 = kcsi[k][j][i].x;
1058 csi1 = kcsi[k][j][i].y;
1059 csi2 = kcsi[k][j][i].z;
1060
1061 eta0 = keta[k][j][i].x;
1062 eta1 = keta[k][j][i].y;
1063 eta2 = keta[k][j][i].z;
1064
1065 zet0 = kzet[k][j][i].x;
1066 zet1 = kzet[k][j][i].y;
1067 zet2 = kzet[k][j][i].z;
1068
1069
1070 g11 = csi0 * zet0 + csi1 * zet1 + csi2 * zet2;
1071 g21 = eta0 * zet0 + eta1 * zet1 + eta2 * zet2;
1072 g31 = zet0 * zet0 + zet1 * zet1 + zet2 * zet2;
1073
1074 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
1075 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
1076 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
1077
1078 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
1079 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
1080 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
1081
1082 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
1083 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
1084 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
1085
1086 ajc = kaj[k][j][i];
1087
1088 double nu = 1./ren, nu_t =0;
1089
1090 if( les || (rans && ti>0) ) {
1091 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k+1][j][i]) ), 2.0) * Sabs;
1092 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k+1][j][i]);
1093 if ( (user->bctype[4]==1 && k==0) || (user->bctype[5]==1 && k==mz-2) ) nu_t=0;
1094
1095 fp3[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu_t);
1096 fp3[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu_t);
1097 fp3[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu_t);
1098 }
1099 else {
1100 fp3[k][j][i].x = 0;
1101 fp3[k][j][i].y = 0;
1102 fp3[k][j][i].z = 0;
1103 }
1104 fp3[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu);//
1105 fp3[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu);//
1106 fp3[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu);//
1107
1108 if(clark) {
1109 double dc, de, dz;
1110 Calculate_dxdydz(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
1111 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
1112
1113 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
1114 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
1115 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
1116 double t21 = t12;
1117 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
1118 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
1119 double t31 = t13;
1120 double t32 = t23;
1121 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
1122
1123 fp3[k][j][i].x -= ( t11 * zet0 + t12 * zet1 + t13 * zet2 ) / 12.;
1124 fp3[k][j][i].y -= ( t21 * zet0 + t22 * zet1 + t23 * zet2 ) / 12.;
1125 fp3[k][j][i].z -= ( t31 * zet0 + t32 * zet1 + t33 * zet2 ) / 12.;
1126 }
1127 }
1128 }
1129 }
1130
1131 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1132
1133 for (k=lzs; k<lze; k++) {
1134 for (j=lys; j<lye; j++) {
1135 for (i=lxs; i<lxe; i++) {
1136 visc[k][j][i].x =
1137 (fp1[k][j][i].x - fp1[k][j][i-1].x +
1138 fp2[k][j][i].x - fp2[k][j-1][i].x +
1139 fp3[k][j][i].x - fp3[k-1][j][i].x);
1140
1141 visc[k][j][i].y =
1142 (fp1[k][j][i].y - fp1[k][j][i-1].y +
1143 fp2[k][j][i].y - fp2[k][j-1][i].y +
1144 fp3[k][j][i].y - fp3[k-1][j][i].y);
1145
1146 visc[k][j][i].z =
1147 (fp1[k][j][i].z - fp1[k][j][i-1].z +
1148 fp2[k][j][i].z - fp2[k][j-1][i].z +
1149 fp3[k][j][i].z - fp3[k-1][j][i].z);
1150
1151 }
1152 }
1153 }
1154/* for (k=zs; k<ze; k++) { */
1155/* for (j=ys; j<ye; j++) { */
1156/* for (i=xs; i<xe; i++) { */
1157/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp1.z is %.15le \n",i,j,k,fp1[k][j][i].z); */
1158/* if (i==0 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp1.z is %.15le \n",i,j,k,fp1[k][j][i].z); */
1159/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp2.z is %.15le \n",i,j,k,fp2[k][j][i].z); */
1160/* if (i==1 && j==0 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp2.z is %.15le \n",i,j,k,fp2[k][j][i].z); */
1161/* if (i==1 && j==1 && k==21) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp3.z is %.15le \n",i,j,k,fp3[k][j][i].z); */
1162/* if (i==1 && j==1 && k==20) PetscPrintf(PETSC_COMM_SELF, "@ i= %d j=%d k=%d fp3.z is %.15le \n",i,j,k,fp3[k][j][i].z); */
1163
1164/* } */
1165/* } */
1166/* } */
1167 DMDAVecRestoreArray(fda, Ucont, &ucont);
1168 DMDAVecRestoreArray(fda, Ucat, &ucat);
1169 DMDAVecRestoreArray(fda, Visc, &visc);
1170
1171 DMDAVecRestoreArray(fda, Csi, &csi);
1172 DMDAVecRestoreArray(fda, Eta, &eta);
1173 DMDAVecRestoreArray(fda, Zet, &zet);
1174
1175 DMDAVecRestoreArray(fda, Fp1, &fp1);
1176 DMDAVecRestoreArray(fda, Fp2, &fp2);
1177 DMDAVecRestoreArray(fda, Fp3, &fp3);
1178
1179 DMDAVecRestoreArray(da, user->lAj, &aj);
1180
1181 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1182 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1183 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1184
1185 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1186 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1187 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1188
1189 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1190 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1191 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1192
1193 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1194
1195 if(les) {
1196 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1197 } else if (rans) {
1198
1199 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1200 }
1201
1202
1203 VecDestroy(&Fp1);
1204 VecDestroy(&Fp2);
1205 VecDestroy(&Fp3);
1206
1207
1208 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Viscous terms calculated .\n");
1209
1210 return(0);
1211}
1212
1213#undef __FUNCT__
1214#define __FUNCT__ "FormFunction1"
1215/**
1216 * @brief Computes the Right-Hand-Side (RHS) of the momentum equations.
1217 *
1218 * This function orchestrates the calculation of the spatial discretization of the
1219 * convective and diffusive terms. It calls specialized helper functions
1220 * (`Convection`, `Viscous`) to compute these terms and then combines them with
1221 * the pressure gradient to form the final RHS vector.
1222 *
1223 * This is a minimally-edited version of the legacy kernel. It operates on a
1224 * single UserCtx (one grid block) and retrieves all global parameters
1225 * (Re, rans, les, etc.) from the master SimCtx.
1226 *
1227 * @param user The UserCtx for the specific block being computed.
1228 * @param Rhs The PETSc Vec where the calculated RHS will be stored.
1229 * @return PetscErrorCode 0 on success.
1230 */
1231PetscErrorCode FormFunction1(UserCtx *user, Vec Rhs)
1232{
1233 PetscErrorCode ierr;
1234 SimCtx *simCtx = user->simCtx;
1235 DM da = user->da, fda = user->fda;
1236 DMDALocalInfo info = user->info;
1237 PetscInt i,j,k;
1238 PetscReal dpdc, dpde,dpdz;
1239 // --- Local Grid Indices and Parameters ---
1240 PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
1241 PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
1242 PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
1243 PetscInt lxs = (xs==0) ? xs+1 : xs;
1244 PetscInt lys = (ys==0) ? ys+1 : ys;
1245 PetscInt lzs = (zs==0) ? zs+1 : zs;
1246 PetscInt lxe = (xe==mx) ? xe-1 : xe;
1247 PetscInt lye = (ye==my) ? ye-1 : ye;
1248 PetscInt lze = (ze==mz) ? ze-1 : ze;
1249 PetscInt mz_end = (user->bctype[5] == 8) ? mz - 2 : mz - 3;
1250
1251 // --- Array Pointers ---
1252 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1253 PetscReal ***p, ***iaj, ***jaj, ***kaj, ***aj, ***nvert, ***nvert_o;
1254 Cmpnts ***rhs, ***rc, ***rct,***ucont;
1255
1256 // --- Temporary Vectors ---
1257 Vec Conv, Visc, Rc, Rct;
1258
1259 PetscFunctionBeginUser;
1260 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Computing RHS (FormFunction1)...\n",
1261 simCtx->rank, user->_this);
1262
1263 // --- Get all necessary array pointers ---
1264 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1265 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1266 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1267 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1268 ierr = DMDAVecGetArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1269 ierr = DMDAVecGetArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1270 ierr = DMDAVecGetArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1271 ierr = DMDAVecGetArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1272 ierr = DMDAVecGetArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1273 ierr = DMDAVecGetArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1274 ierr = DMDAVecGetArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1275 ierr = DMDAVecGetArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1276 ierr = DMDAVecGetArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1277 ierr = DMDAVecGetArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1278 ierr = DMDAVecGetArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1279 ierr = DMDAVecGetArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1280 ierr = DMDAVecGetArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1281 ierr = DMDAVecGetArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1282 ierr = DMDAVecGetArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1283
1284 // --- Create temporary work vectors ---
1285 ierr = VecDuplicate(user->lUcont, &Rc); CHKERRQ(ierr);
1286 ierr = VecDuplicate(Rc, &Rct); CHKERRQ(ierr);
1287 ierr = VecDuplicate(Rct, &Conv); CHKERRQ(ierr);
1288 ierr = VecDuplicate(Rct, &Visc); CHKERRQ(ierr);
1289
1290 // ========================================================================
1291 // CORE LOGIC (UNCHANGED FROM LEGACY CODE)
1292 // ========================================================================
1293
1294 // 1. Obtain Cartesian velocity from Contravariant velocity
1295 ierr = Contra2Cart(user); CHKERRQ(ierr);
1296
1297 // 2. Compute Convective term
1298 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating convective terms...\n");
1299 if (simCtx->moveframe || simCtx->rotateframe) {
1300 // ierr = Convection_MV(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1301 } else {
1302 ierr = Convection(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1303 }
1304
1305 // 3. Compute Viscous term
1306 if (simCtx->invicid) {
1307 ierr = VecSet(Visc, 0.0); CHKERRQ(ierr);
1308 } else {
1309 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating viscous terms...\n");
1310 ierr = Viscous(user, user->lUcont, user->lUcat, Visc); CHKERRQ(ierr);
1311 }
1312
1313 // 4. Combine terms to get Cartesian RHS: Rc = Visc - Conv
1314 ierr = VecWAXPY(Rc, -1.0, Conv, Visc); CHKERRQ(ierr);
1315
1316 // 5. Convert Cartesian RHS (Rc) to Contravariant RHS (Rct)
1317 LOG_ALLOW(LOCAL, LOG_DEBUG, " Converting Cartesian RHS to Contravariant RHS...\n");
1318 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1319 ierr = DMDAVecGetArray(fda, Rc, &rc); CHKERRQ(ierr);
1320
1321 for (k = lzs; k < lze; k++) {
1322 for (j = lys; j < lye; j++) {
1323 for (i = lxs; i < lxe; i++) {
1324 rct[k][j][i].x = aj[k][j][i] *
1325 (0.5 * (csi[k][j][i].x + csi[k][j][i-1].x) * rc[k][j][i].x +
1326 0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
1327 0.5 * (csi[k][j][i].z + csi[k][j][i-1].z) * rc[k][j][i].z);
1328 rct[k][j][i].y = aj[k][j][i] *
1329 (0.5 * (eta[k][j][i].x + eta[k][j-1][i].x) * rc[k][j][i].x +
1330 0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
1331 0.5 * (eta[k][j][i].z + eta[k][j-1][i].z) * rc[k][j][i].z);
1332 rct[k][j][i].z = aj[k][j][i] *
1333 (0.5 * (zet[k][j][i].x + zet[k-1][j][i].x) * rc[k][j][i].x +
1334 0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
1335 0.5 * (zet[k][j][i].z + zet[k-1][j][i].z) * rc[k][j][i].z);
1336 }
1337 }
1338 }
1339 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1340 ierr = DMDAVecRestoreArray(fda, Rc, &rc); CHKERRQ(ierr);
1341
1342 PetscBarrier(NULL);
1343
1344 DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct);
1345 DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct);
1346
1347 DMDAVecGetArray(fda, Rct, &rct);
1348
1349 if (user->bctype[0]==7 || user->bctype[1]==7){
1350 if (xs==0){
1351 i=xs;
1352 for (k=lzs; k<lze; k++) {
1353 for (j=lys; j<lye; j++) {
1354 rct[k][j][i].x=rct[k][j][i-2].x;
1355 }
1356 }
1357 }
1358
1359 if (xe==mx){
1360 i=mx-1;
1361 for (k=lzs; k<lze; k++) {
1362 for (j=lys; j<lye; j++) {
1363 rct[k][j][i].x=rct[k][j][i+2].x;
1364 }
1365 }
1366 }
1367 }
1368
1369 if (user->bctype[2]==7 || user->bctype[3]==7){
1370 if (ys==0){
1371 j=ys;
1372 for (k=lzs; k<lze; k++) {
1373 for (i=lxs; i<lxe; i++) {
1374 rct[k][j][i].y=rct[k][j-2][i].y;
1375 }
1376 }
1377 }
1378
1379 if (ye==my){
1380 j=my-1;
1381 for (k=lzs; k<lze; k++) {
1382 for (i=lxs; i<lxe; i++) {
1383 rct[k][j][i].y=rct[k][j+2][i].y;
1384 }
1385 }
1386 }
1387 }
1388 if (user->bctype[4]==7 || user->bctype[5]==7){
1389 if (zs==0){
1390 k=zs;
1391 for (j=lys; j<lye; j++) {
1392 for (i=lxs; i<lxe; i++) {
1393 rct[k][j][i].z=rct[k-2][j][i].z;
1394 }
1395 }
1396 }
1397 if (ze==mz){
1398 k=mz-1;
1399 for (j=lys; j<lye; j++) {
1400 for (i=lxs; i<lxe; i++) {
1401 rct[k][j][i].z=rct[k+2][j][i].z;
1402 }
1403 }
1404 }
1405 }
1406
1407 // 6. Add Pressure Gradient Term and Finalize RHS
1408 // This involves calculating pressure derivatives (dpdc, dpde, dpdz) and using
1409 // them to adjust the contravariant RHS. The full stencil logic is preserved.
1410 LOG_ALLOW(LOCAL, LOG_DEBUG, " Adding pressure gradient term to RHS...\n");
1411
1412 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1413
1414 ierr = DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1415 ierr = DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1416
1417 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1418
1419 for (k = lzs; k < lze; k++) {
1420 for (j = lys; j < lye; j++) {
1421 for (i = lxs; i < lxe; i++) {
1422 PetscReal dpdc, dpde, dpdz;
1423 dpdc = p[k][j][i+1] - p[k][j][i];
1424
1425 if ((j==my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1426 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1427 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1428 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1429 }
1430 }
1431 else if ((j==my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1432 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1433 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1434 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1435 }
1436 }
1437 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1438 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1439 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1440 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1441 }
1442 }
1443 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1444 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1445 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1446 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1447 }
1448 }
1449 else {
1450 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1451 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1452 }
1453
1454 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1455 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1456 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1457 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1458 }
1459 }
1460 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1461 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1462 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1463 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1464 }
1465 }
1466 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1467 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1468 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1469 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1470 }
1471 }
1472 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1473 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1474 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1475 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1476 }
1477 }
1478 else {
1479 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1480 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1481 }
1482
1483 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1484
1485
1486 rhs[k][j][i].x -=
1487 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1488 icsi[k][j][i].y * icsi[k][j][i].y +
1489 icsi[k][j][i].z * icsi[k][j][i].z)+
1490 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1491 ieta[k][j][i].y * icsi[k][j][i].y +
1492 ieta[k][j][i].z * icsi[k][j][i].z)+
1493 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1494 izet[k][j][i].y * icsi[k][j][i].y +
1495 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1496
1497 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1498 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1499 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1500 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1501 }
1502 }
1503 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1504 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1505 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1506 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1507 }
1508 }
1509 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1510 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1511 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1512 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1513 }
1514 }
1515 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1516 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1517 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1518 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1519 }
1520 }
1521 else {
1522 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1523 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1524 }
1525
1526 dpde = p[k][j+1][i] - p[k][j][i];
1527
1528 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1529 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1530 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1531 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1532 }
1533 }
1534 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1535 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1536 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1537 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1538 }
1539 }
1540 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1541 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1542 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1543 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1544 }
1545 }
1546 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1547 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1548 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1549 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1550 }
1551 }
1552 else {
1553 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1554 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1555 }
1556
1557 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1558
1559
1560 rhs[k][j][i].y -=
1561 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1562 jcsi[k][j][i].y * jeta[k][j][i].y +
1563 jcsi[k][j][i].z * jeta[k][j][i].z) +
1564 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1565 jeta[k][j][i].y * jeta[k][j][i].y +
1566 jeta[k][j][i].z * jeta[k][j][i].z) +
1567 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1568 jzet[k][j][i].y * jeta[k][j][i].y +
1569 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1570
1571 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1572 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1573 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1574 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1575 }
1576 }
1577 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1578 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1579 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1580 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1581 }
1582 }
1583 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1584 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1585 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1586 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1587 }
1588 }
1589 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1590 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1591 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1592 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1593 }
1594 }
1595 else {
1596 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1597 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1598 }
1599
1600 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1601 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1602 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1603 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1604 }
1605 }
1606 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1607 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1608 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1609 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1610 }
1611 }
1612 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1613 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1614 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1615 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1616 }
1617 }
1618 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1619 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1620 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1621 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1622 }
1623 }
1624 else {
1625 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1626 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1627 }
1628
1629 dpdz = (p[k+1][j][i] - p[k][j][i]);
1630
1631 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1632
1633 rhs[k][j][i].z -=
1634 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1635 kcsi[k][j][i].y * kzet[k][j][i].y +
1636 kcsi[k][j][i].z * kzet[k][j][i].z) +
1637 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1638 keta[k][j][i].y * kzet[k][j][i].y +
1639 keta[k][j][i].z * kzet[k][j][i].z) +
1640 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1641 kzet[k][j][i].y * kzet[k][j][i].y +
1642 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1643
1644 }
1645 }
1646 }
1647
1648
1649 //Mohsen March 2012//
1650
1651 // rhs.x at boundaries for periodic bc at i direction//
1652 if (user->bctype[0]==7 && xs==0){
1653 for (k=lzs; k<lze; k++) {
1654 for (j=lys; j<lye; j++) {
1655 i=xs;
1656
1657 dpdc = p[k][j][i+1] - p[k][j][i];
1658
1659 if ((j==my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1660 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1661 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1662 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1663 }
1664 }
1665 else if ((j==my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1666 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1667 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1668 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1669 }
1670 }
1671 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1672 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1673 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1674 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1675 }
1676 }
1677 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1678 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1679 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1680 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1681 }
1682 }
1683 else {
1684 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1685 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1686 }
1687
1688 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1689 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1690 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1691 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1692 }
1693 }
1694 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1695 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1696 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1697 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1698 }
1699 }
1700 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1701 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1702 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1703 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1704 }
1705 }
1706 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1707 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1708 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1709 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1710 }
1711 }
1712 else {
1713 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1714 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1715 }
1716
1717 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1718 rhs[k][j][i].x -=
1719 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1720 icsi[k][j][i].y * icsi[k][j][i].y +
1721 icsi[k][j][i].z * icsi[k][j][i].z)+
1722 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1723 ieta[k][j][i].y * icsi[k][j][i].y +
1724 ieta[k][j][i].z * icsi[k][j][i].z)+
1725 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1726 izet[k][j][i].y * icsi[k][j][i].y +
1727 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1728 }
1729 }
1730 }
1731
1732// rhs.y at boundaries for periodic bc at j direction//
1733 if (user->bctype[2]==7 && ys==0){
1734 for (k=lzs; k<lze; k++) {
1735 for (i=lxs; i<lxe; i++) {
1736
1737 j=ys;
1738
1739 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1740 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1741 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1742 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1743 }
1744 }
1745 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1746 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1747 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1748 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1749 }
1750 }
1751 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1752 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1753 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1754 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1755 }
1756 }
1757 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1758 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1759 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1760 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1761 }
1762 }
1763 else {
1764 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1765 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1766 }
1767
1768 dpde = p[k][j+1][i] - p[k][j][i];
1769
1770 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1771 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1772 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1773 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1774 }
1775 }
1776 else if ((k == mz-2 || k==1 ) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1777 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1778 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1779 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1780 }
1781 }
1782 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1783 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1784 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1785 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1786 }
1787 }
1788 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1789 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1790 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1791 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1792 }
1793 }
1794 else {
1795 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1796 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1797 }
1798
1799 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1800
1801 rhs[k][j][i].y -=
1802 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1803 jcsi[k][j][i].y * jeta[k][j][i].y +
1804 jcsi[k][j][i].z * jeta[k][j][i].z)+
1805 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1806 jeta[k][j][i].y * jeta[k][j][i].y +
1807 jeta[k][j][i].z * jeta[k][j][i].z)+
1808 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1809 jzet[k][j][i].y * jeta[k][j][i].y +
1810 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1811
1812 }
1813 }
1814 }
1815
1816 // rhs.z at boundaries for periodic bc at k direction//
1817 if (user->bctype[4]==7&& zs==0){
1818 for (j=lys; j<lye; j++) {
1819 for (i=lxs; i<lxe; i++) {
1820
1821 k=zs;
1822
1823 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1824 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1825 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1826 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1827 }
1828 }
1829 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1830 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1831 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1832 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1833 }
1834 }
1835 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1836 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1837 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1838 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1839 }
1840 }
1841 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1842 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1843 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1844 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1845 }
1846 }
1847 else {
1848 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1849 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1850 }
1851
1852 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1853 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1854 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1855 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1856 }
1857 }
1858 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1859 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1860 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1861 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1862 }
1863 }
1864 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1865 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1866 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1867 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1868 }
1869 }
1870 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1871 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1872 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1873 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1874 }
1875 }
1876 else {
1877 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1878 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1879 }
1880
1881 dpdz = (p[k+1][j][i] - p[k][j][i]);
1882
1883 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1884
1885 rhs[k][j][i].z -=
1886 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1887 kcsi[k][j][i].y * kzet[k][j][i].y +
1888 kcsi[k][j][i].z * kzet[k][j][i].z)+
1889 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1890 keta[k][j][i].y * kzet[k][j][i].y +
1891 keta[k][j][i].z * kzet[k][j][i].z)+
1892 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1893 kzet[k][j][i].y * kzet[k][j][i].y +
1894 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1895
1896 }
1897 }
1898 }
1899
1900 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1901
1902 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure Gradient added to RHS .\n");
1903 PetscInt TwoD = simCtx->TwoD;
1904
1905 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases initiated .\n");
1906
1907 // 7. Final clean-up for immersed boundaries and 2D cases
1908 for (k=lzs; k<lze; k++) {
1909 for (j=lys; j<lye; j++) {
1910 for (i=lxs; i<lxe; i++) {
1911 if (TwoD==1)
1912 rhs[k][j][i].x =0.;
1913 else if (TwoD==2)
1914 rhs[k][j][i].y =0.;
1915 else if (TwoD==3)
1916 rhs[k][j][i].z =0.;
1917
1918 if (nvert[k][j][i]>0.1) {
1919 rhs[k][j][i].x = 0;
1920 rhs[k][j][i].y = 0;
1921 rhs[k][j][i].z = 0;
1922 }
1923 if (nvert[k][j][i+1]>0.1) {
1924 rhs[k][j][i].x=0;
1925 }
1926 if (nvert[k][j+1][i]>0.1) {
1927 rhs[k][j][i].y=0;
1928 }
1929 if (nvert[k+1][j][i]>0.1) {
1930 rhs[k][j][i].z=0;
1931 }
1932 }
1933 }
1934 }
1935 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases complete .\n");
1936
1937 // ========================================================================
1938
1939 // --- Restore all PETSc array pointers ---
1940 // DMDAVecRestoreArray(fda, user->lUcont, &ucont);
1941 DMDAVecRestoreArray(fda, Rhs, &rhs);
1942 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Rhs restored successfully! .\n");
1943
1944 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1945 DMDAVecRestoreArray(fda, user->lEta, &eta);
1946 DMDAVecRestoreArray(fda, user->lZet, &zet);
1947 DMDAVecRestoreArray(da, user->lAj, &aj);
1948 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Face metrics restored successfully! .\n");
1949
1950 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1951 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1952 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1953 DMDAVecRestoreArray(da, user->lIAj, &iaj);
1954 LOG_ALLOW(GLOBAL,LOG_DEBUG,"I Face metrics restored successfully! .\n");
1955
1956 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1957 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1958 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1959 DMDAVecRestoreArray(da, user->lJAj, &jaj);
1960 LOG_ALLOW(GLOBAL,LOG_DEBUG,"J Face metrics restored successfully! .\n");
1961
1962 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1963 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1964 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1965 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1966 LOG_ALLOW(GLOBAL,LOG_DEBUG,"K Face metrics restored successfully! .\n");
1967
1968 DMDAVecRestoreArray(da, user->lP, &p);
1969 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure restored successfully! .\n");
1970
1971 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1972 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Nvert restored successfully! .\n");
1973
1974 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Cell Centered scalars restored successfully! .\n");
1975
1976 // --- Destroy temporary work vectors ---
1977 ierr = VecDestroy(&Conv); CHKERRQ(ierr);
1978 ierr = VecDestroy(&Visc); CHKERRQ(ierr);
1979 ierr = VecDestroy(&Rc); CHKERRQ(ierr);
1980 ierr = VecDestroy(&Rct); CHKERRQ(ierr);
1981 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Temporary work vectors destroyed successfully! .\n");
1982
1983 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: RHS computation complete.\n",
1984 simCtx->rank, user->_this);
1985 PetscFunctionReturn(0);
1986}
logging.h
Logging utilities and macros for PETSc-based applications.
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:44
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:45
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:207
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
Calculate_normal_and_area
static void Calculate_normal_and_area(Cmpnts csi, Cmpnts eta, Cmpnts zet, double ni[3], double nj[3], double nk[3], double *Ai, double *Aj, double *Ak)
Definition rhs.c:32
Viscous
PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
Definition rhs.c:580
Convection
PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
Definition rhs.c:87
Calculate_Covariant_metrics
static void Calculate_Covariant_metrics(double g[3][3], double G[3][3])
Definition rhs.c:13
Calculate_dxdydz
static void Calculate_dxdydz(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
Definition rhs.c:68
FormFunction1
PetscErrorCode FormFunction1(UserCtx *user, Vec Rhs)
Computes the Right-Hand-Side (RHS) of the momentum equations.
Definition rhs.c:1231
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:1785
SimCtx::clark
PetscInt clark
Definition variables.h:579
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:536
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:536
SimCtx::rank
PetscMPIInt rank
Definition variables.h:516
UserCtx::lIEta
Vec lIEta
Definition variables.h:675
UserCtx::lIZet
Vec lIZet
Definition variables.h:675
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
SimCtx::ren
PetscReal ren
Definition variables.h:549
UserCtx::lIAj
Vec lIAj
Definition variables.h:675
UserCtx::_this
PetscInt _this
Definition variables.h:643
UserCtx::lKEta
Vec lKEta
Definition variables.h:677
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:676
Cmpnts::x
PetscScalar x
Definition variables.h:100
SimCtx::invicid
PetscInt invicid
Definition variables.h:536
UserCtx::lKZet
Vec lKZet
Definition variables.h:677
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:680
UserCtx::lJEta
Vec lJEta
Definition variables.h:676
UserCtx::lCsi
Vec lCsi
Definition variables.h:673
Cmpnts::z
PetscScalar z
Definition variables.h:100
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:677
SimCtx::central
PetscInt central
Definition variables.h:548
UserCtx::bctype
PetscInt bctype[6]
Definition variables.h:653
UserCtx::lJZet
Vec lJZet
Definition variables.h:676
UserCtx::lUcont
Vec lUcont
Definition variables.h:657
SimCtx::step
PetscInt step
Definition variables.h:521
UserCtx::lAj
Vec lAj
Definition variables.h:673
UserCtx::lICsi
Vec lICsi
Definition variables.h:675
UserCtx::info
DMDALocalInfo info
Definition variables.h:637
UserCtx::lUcat
Vec lUcat
Definition variables.h:657
Cmpnts::y
PetscScalar y
Definition variables.h:100
UserCtx::lEta
Vec lEta
Definition variables.h:673
SimCtx::les
PetscInt les
Definition variables.h:578
UserCtx::lJAj
Vec lJAj
Definition variables.h:676
SimCtx::rotateframe
PetscInt rotateframe
Definition variables.h:536
UserCtx::lKAj
Vec lKAj
Definition variables.h:677
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
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:630
nu_t
double nu_t(double yplus)
Definition wallfunction.c:103