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#undef __FUNCT__
4#define __FUNCT__ "Convection"
5PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
6{
7
8 // --- CONTEXT ACQUISITION BLOCK ---
9 // Get the master simulation context from the UserCtx.
10 SimCtx *simCtx = user->simCtx;
11
12 // Create local variables to mirror the legacy globals for minimal code changes.
13 const LESModelType les = simCtx->les;
14 const PetscInt central = simCtx->central; // Get this from SimCtx now
15 // --- END CONTEXT ACQUISITION BLOCK ---
16
17 Cmpnts ***ucont, ***ucat;
18 DM da = user->da, fda = user->fda;
19 DMDALocalInfo info;
20 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
21 PetscInt mx, my, mz; // Dimensions in three directions
22 PetscInt i, j, k;
23 Vec Fp1, Fp2, Fp3;
24 Cmpnts ***fp1, ***fp2, ***fp3;
25 Cmpnts ***conv;
26
27 PetscReal ucon, up, um;
28 PetscReal coef = 0.125, innerblank=7.;
29
30 PetscInt lxs, lxe, lys, lye, lzs, lze, gxs, gxe, gys, gye, gzs,gze;
31
32 PetscReal ***nvert,***aj;
33
34 PROFILE_FUNCTION_BEGIN;
35
36 DMDAGetLocalInfo(da, &info);
37 mx = info.mx; my = info.my; mz = info.mz;
38 xs = info.xs; xe = xs + info.xm;
39 ys = info.ys; ye = ys + info.ym;
40 zs = info.zs; ze = zs + info.zm;
41 gxs = info.gxs; gxe = gxs + info.gxm;
42 gys = info.gys; gye = gys + info.gym;
43 gzs = info.gzs; gze = gzs + info.gzm;
44
45 DMDAVecGetArray(fda, Ucont, &ucont);
46 DMDAVecGetArray(fda, Ucat, &ucat);
47 DMDAVecGetArray(fda, Conv, &conv);
48 DMDAVecGetArray(da, user->lAj, &aj);
49
50 VecDuplicate(Ucont, &Fp1);
51 VecDuplicate(Ucont, &Fp2);
52 VecDuplicate(Ucont, &Fp3);
53
54 DMDAVecGetArray(fda, Fp1, &fp1);
55 DMDAVecGetArray(fda, Fp2, &fp2);
56 DMDAVecGetArray(fda, Fp3, &fp3);
57
58 DMDAVecGetArray(da, user->lNvert, &nvert);
59
60
61 /* We have two different sets of node: 1. grid node, the physical points
62 where grid lines intercross; 2. storage node, where we store variables.
63 All node without explicitly specified as "grid node" refers to
64 storage node.
65
66 The integer node is defined at cell center while half node refers to
67 the actual grid node. (The reason to choose this arrangement is we need
68 ghost node, which is half node away from boundaries, to specify boundary
69 conditions. By using this storage arrangement, the actual storage need
70 is (IM+1) * (JM + 1) * (KM+1) where IM, JM, & KM refer to the number of
71 grid nodes along i, j, k directions.)
72
73 DA, the data structure used to define the storage of 3D arrays, is defined
74 as mx * my * mz. mx = IM+1, my = JM+1, mz = KM+1.
75
76 Staggered grid arrangement is used in this solver.
77 Pressure is stored at interger node (hence the cell center) and volume
78 fluxes defined on the center of each surface of a given control volume
79 is stored on the cloest upper integer node. */
80
81 /* First we calculate the flux on cell surfaces. Stored on the upper integer
82 node. For example, along i direction, the flux are stored at node 0:mx-2*/
83
84 lxs = xs; lxe = xe;
85 lys = ys; lye = ye;
86 lzs = zs; lze = ze;
87
88 if (xs==0) lxs = xs+1;
89 if (ys==0) lys = ys+1;
90 if (zs==0) lzs = zs+1;
91
92 if (xe==mx) lxe=xe-1;
93 if (ye==my) lye=ye-1;
94 if (ze==mz) lze=ze-1;
95
96 //Mohsen Sep 2012//
97/* First update the computational ghost points velocity for periodic boundary conditions
98 just for this subroutine because of Quick scheme for velocity deravatives */
99 /* for (k=zs; k<ze; k++) { */
100/* for (j=ys; j<ye; j++) { */
101/* for (i=xs; i<xe; i++) { */
102/* if (i==1 && (j==1) && (k==0 || k==1)) */
103/* 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 ); */
104/* } */
105/* } */
106/* } */
107 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
108 if (xs==0){
109 i=xs-1;
110 for (k=gzs; k<gze; k++) {
111 for (j=gys; j<gye; j++) {
112 ucat[k][j][i].x=ucat[k][j][i-2].x;
113 ucat[k][j][i].y=ucat[k][j][i-2].y;
114 ucat[k][j][i].z=ucat[k][j][i-2].z;
115 nvert[k][j][i]=nvert[k][j][i-2];
116 }
117 }
118 }
119 if (xe==mx){
120 i=mx;
121 for (k=gzs; k<gze; k++) {
122 for (j=gys; j<gye; j++) {
123 ucat[k][j][i].x=ucat[k][j][i+2].x;
124 ucat[k][j][i].y=ucat[k][j][i+2].y;
125 ucat[k][j][i].z=ucat[k][j][i+2].z;
126 nvert[k][j][i]=nvert[k][j][i+2];
127 }
128 }
129 }
130 }
131 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
132 if (ys==0){
133 j=ys-1;
134 for (k=gzs; k<gze; k++) {
135 for (i=gxs; i<gxe; i++) {
136 ucat[k][j][i].x=ucat[k][j-2][i].x;
137 ucat[k][j][i].y=ucat[k][j-2][i].y;
138 ucat[k][j][i].z=ucat[k][j-2][i].z;
139 nvert[k][j][i]=nvert[k][j-2][i];
140 }
141 }
142 }
143 if (ye==my){
144 j=my;
145 for (k=gzs; k<gze; k++) {
146 for (i=gxs; i<gxe; i++) {
147 ucat[k][j][i].x=ucat[k][j+2][i].x;
148 ucat[k][j][i].y=ucat[k][j+2][i].y;
149 ucat[k][j][i].z=ucat[k][j+2][i].z;
150 nvert[k][j][i]=nvert[k][j+2][i];
151 }
152 }
153 }
154 }
155 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
156 if (zs==0){
157 k=zs-1;
158 for (j=gys; j<gye; j++) {
159 for (i=gxs; i<gxe; i++) {
160 ucat[k][j][i].x=ucat[k-2][j][i].x;
161 ucat[k][j][i].y=ucat[k-2][j][i].y;
162 ucat[k][j][i].z=ucat[k-2][j][i].z;
163 nvert[k][j][i]=nvert[k-2][j][i];
164 }
165 }
166 }
167 if (ze==mz){
168 k=mz;
169 for (j=gys; j<gye; j++) {
170 for (i=gxs; i<gxe; i++) {
171 ucat[k][j][i].x=ucat[k+2][j][i].x;
172 ucat[k][j][i].y=ucat[k+2][j][i].y;
173 ucat[k][j][i].z=ucat[k+2][j][i].z;
174 nvert[k][j][i]=nvert[k+2][j][i];
175 }
176 }
177 }
178 }
179
180 VecSet(Conv, 0.0);
181
182 /* Calculating the convective terms on cell centers.
183 First calcualte the contribution from i direction
184 The flux is evaluated by QUICK scheme */
185
186 for (k=lzs; k<lze; k++){
187 for (j=lys; j<lye; j++){
188 for (i=lxs-1; i<lxe; i++){
189
190
191 ucon = ucont[k][j][i].x * 0.5;
192
193 up = ucon + fabs(ucon);
194 um = ucon - fabs(ucon);
195
196 if (i>0 && i<mx-2 &&
197 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
198 (nvert[k][j][i-1] < 0.1 || nvert[k][j][i-1]>innerblank)) { // interial nodes
199 if ((les || central)) {
200 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
201 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
202 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
203
204 } else {
205 fp1[k][j][i].x =
206 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) +
207 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);
208 fp1[k][j][i].y =
209 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) +
210 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);
211 fp1[k][j][i].z =
212 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) +
213 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);
214 }
215 }
216 else if ((les || central) && (i==0 || i==mx-2) &&
217 (nvert[k][j][i+1] < 0.1 || nvert[k][j][i+1]>innerblank) &&
218 (nvert[k][j][i ] < 0.1 || nvert[k][j][i ]>innerblank))
219 {
220 fp1[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j][i+1].x );
221 fp1[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j][i+1].y );
222 fp1[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j][i+1].z );
223 }
224 else if (i==0 ||(nvert[k][j][i-1] > 0.1) ) {
225 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==0 && (nvert[k][j][i-1]<0.1 && nvert[k][j][i+1]<0.1)){//Mohsen Feb 12
226 fp1[k][j][i].x =
227 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) +
228 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);
229 fp1[k][j][i].y =
230 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) +
231 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);
232 fp1[k][j][i].z =
233 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) +
234 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);
235 }else{
236 fp1[k][j][i].x =
237 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) +
238 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);
239 fp1[k][j][i].y =
240 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) +
241 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);
242 fp1[k][j][i].z =
243 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) +
244 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);
245 }
246 }
247 else if (i==mx-2 ||(nvert[k][j][i+1]) > 0.1) {
248 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && i==mx-2 &&(nvert[k][j][i-1]<0.1 && nvert[k][j][i+1]<0.1)){//Mohsen Feb 12
249 fp1[k][j][i].x =
250 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) +
251 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);
252 fp1[k][j][i].y =
253 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) +
254 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);
255 fp1[k][j][i].z =
256 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) +
257 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);
258 }else{
259 fp1[k][j][i].x =
260 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) +
261 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);
262 fp1[k][j][i].y =
263 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) +
264 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);
265 fp1[k][j][i].z =
266 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) +
267 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);
268 }
269 }
270 }
271 }
272 }
273
274 /* j direction */
275 for (k=lzs; k<lze; k++) {
276 for(j=lys-1; j<lye; j++) {
277 for(i=lxs; i<lxe; i++) {
278 ucon = ucont[k][j][i].y * 0.5;
279
280 up = ucon + fabs(ucon);
281 um = ucon - fabs(ucon);
282
283 if (j>0 && j<my-2 &&
284 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i] > innerblank) &&
285 (nvert[k][j-1][i] < 0.1 || nvert[k][j-1][i] > innerblank)) {
286 if ((les || central)) {
287 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
288 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
289 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
290
291 } else {
292 fp2[k][j][i].x =
293 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) +
294 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);
295 fp2[k][j][i].y =
296 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) +
297 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);
298 fp2[k][j][i].z =
299 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) +
300 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);
301 }
302 }
303 else if ((les || central) && (j==0 || i==my-2) &&
304 (nvert[k][j+1][i] < 0.1 || nvert[k][j+1][i]>innerblank) &&
305 (nvert[k][j ][i] < 0.1 || nvert[k][j ][i]>innerblank))
306 {
307 fp2[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k][j+1][i].x );
308 fp2[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k][j+1][i].y );
309 fp2[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k][j+1][i].z );
310 }
311 else if (j==0 || (nvert[k][j-1][i]) > 0.1) {
312 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==0 && (nvert[k][j-1][i]<0.1 && nvert[k][j+1][i]<0.1 )){//Mohsen Feb 12 //
313 fp2[k][j][i].x =
314 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) +
315 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);
316 fp2[k][j][i].y =
317 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) +
318 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);
319 fp2[k][j][i].z =
320 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) +
321 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);
322 }else{
323 fp2[k][j][i].x =
324 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) +
325 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);
326 fp2[k][j][i].y =
327 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) +
328 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);
329 fp2[k][j][i].z =
330 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) +
331 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);
332 }
333 }
334 else if (j==my-2 ||(nvert[k][j+1][i]) > 0.1) {
335 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && j==my-2 && (nvert[k][j-1][i]<0.1 && nvert[k][j+1][i]<0.1 )){//Mohsen Feb 12//
336 fp2[k][j][i].x =
337 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) +
338 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);
339 fp2[k][j][i].y =
340 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) +
341 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);
342 fp2[k][j][i].z =
343 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) +
344 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);
345 }else{
346 fp2[k][j][i].x =
347 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) +
348 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);
349 fp2[k][j][i].y =
350 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) +
351 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);
352 fp2[k][j][i].z =
353 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) +
354 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);
355 }
356 }
357 }
358 }
359 }
360
361
362 /* k direction */
363 for (k=lzs-1; k<lze; k++) {
364 for(j=lys; j<lye; j++) {
365 for(i=lxs; i<lxe; i++) {
366 ucon = ucont[k][j][i].z * 0.5;
367
368 up = ucon + fabs(ucon);
369 um = ucon - fabs(ucon);
370
371 if (k>0 && k<mz-2 &&
372 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i] > innerblank) &&
373 (nvert[k-1][j][i] < 0.1 || nvert[k-1][j][i] > innerblank)) {
374 if ((les || central)) {
375 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
376 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
377 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
378
379 } else {
380 fp3[k][j][i].x =
381 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) +
382 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);
383 fp3[k][j][i].y =
384 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) +
385 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);
386 fp3[k][j][i].z =
387 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) +
388 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);
389 }
390 }
391 else if ((les || central) && (k==0 || k==mz-2) &&
392 (nvert[k+1][j][i] < 0.1 || nvert[k+1][j][i]>innerblank) &&
393 (nvert[k ][j][i] < 0.1 || nvert[k ][j][i]>innerblank))
394 {
395 fp3[k][j][i].x = ucon * ( ucat[k][j][i].x + ucat[k+1][j][i].x );
396 fp3[k][j][i].y = ucon * ( ucat[k][j][i].y + ucat[k+1][j][i].y );
397 fp3[k][j][i].z = ucon * ( ucat[k][j][i].z + ucat[k+1][j][i].z );
398 }
399 else if (k<mz-2 && (k==0 ||(nvert[k-1][j][i]) > 0.1)) {
400 if(user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==0 && (nvert[k-1][j][i]<0.1 && nvert[k+1][j][i]<0.1)){//Mohsen Feb 12//
401 fp3[k][j][i].x =
402 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) +
403 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);
404 fp3[k][j][i].y =
405 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) +
406 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);
407 fp3[k][j][i].z =
408 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) +
409 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);
410 }else{
411 fp3[k][j][i].x =
412 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) +
413 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);
414 fp3[k][j][i].y =
415 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) +
416 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);
417 fp3[k][j][i].z =
418 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) +
419 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);
420 }
421 }
422 else if (k>0 && (k==mz-2 ||(nvert[k+1][j][i]) > 0.1)) {
423 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && k==mz-2 && (nvert[k-1][j][i]<0.1 && nvert[k+1][j][i]<0.1)){//Mohsen Feb 12//
424 fp3[k][j][i].x =
425 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) +
426 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);
427 fp3[k][j][i].y =
428 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) +
429 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);
430 fp3[k][j][i].z =
431 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) +
432 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);
433 }else{
434 fp3[k][j][i].x =
435 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) +
436 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);
437 fp3[k][j][i].y =
438 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) +
439 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);
440 fp3[k][j][i].z =
441 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) +
442 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);
443 }
444 }
445 }
446 }
447 }
448
449 /* Calculate the convective terms under cartesian coordinates */
450
451 for (k=lzs; k<lze; k++) {
452 for (j=lys; j<lye; j++) {
453 for (i=lxs; i<lxe; i++) {
454 conv[k][j][i].x =
455 fp1[k][j][i].x - fp1[k][j][i-1].x +
456 fp2[k][j][i].x - fp2[k][j-1][i].x +
457 fp3[k][j][i].x - fp3[k-1][j][i].x;
458
459 conv[k][j][i].y =
460 fp1[k][j][i].y - fp1[k][j][i-1].y +
461 fp2[k][j][i].y - fp2[k][j-1][i].y +
462 fp3[k][j][i].y - fp3[k-1][j][i].y;
463
464 conv[k][j][i].z =
465 fp1[k][j][i].z - fp1[k][j][i-1].z +
466 fp2[k][j][i].z - fp2[k][j-1][i].z +
467 fp3[k][j][i].z - fp3[k-1][j][i].z;
468 }
469 }
470 }
471 /* for (k=zs; k<ze; k++) { */
472/* for (j=ys; j<ye; j++) { */
473/* for (i=xs; i<xe; i++) { */
474/* if (i==1 && (j==1) && (k==1 || k==21 || k==22|| k==200)) */
475/* 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); */
476/* } */
477/* } */
478/* } */
479
480 DMDAVecRestoreArray(fda, Ucont, &ucont);
481 DMDAVecRestoreArray(fda, Ucat, &ucat);
482 DMDAVecRestoreArray(fda, Conv, &conv);
483 DMDAVecRestoreArray(da, user->lAj, &aj);
484
485 DMDAVecRestoreArray(fda, Fp1, &fp1);
486 DMDAVecRestoreArray(fda, Fp2, &fp2);
487 DMDAVecRestoreArray(fda, Fp3, &fp3);
488 DMDAVecRestoreArray(da, user->lNvert, &nvert);
489
490 VecDestroy(&Fp1);
491 VecDestroy(&Fp2);
492 VecDestroy(&Fp3);
493
494
495 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Convective term calculated .\n");
496
497 PROFILE_FUNCTION_END;
498 return (0);
499}
500
501
502#undef __FUNCT__
503#define __FUNCT__ "Viscous"
504PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
505{
506
507 Vec Csi = user->lCsi, Eta = user->lEta, Zet = user->lZet;
508
509 Cmpnts ***ucont, ***ucat;
510
511 Cmpnts ***csi, ***eta, ***zet;
512 Cmpnts ***icsi, ***ieta, ***izet;
513 Cmpnts ***jcsi, ***jeta, ***jzet;
514 Cmpnts ***kcsi, ***keta, ***kzet;
515
516 PetscReal ***nvert;
517
518 DM da = user->da, fda = user->fda;
519 DMDALocalInfo info;
520 PetscInt xs, xe, ys, ye, zs, ze; // Local grid information
521 PetscInt mx, my, mz; // Dimensions in three directions
522 PetscInt i, j, k;
523 Vec Fp1, Fp2, Fp3;
524 Cmpnts ***fp1, ***fp2, ***fp3;
525 Cmpnts ***visc;
526 PetscReal ***aj, ***iaj, ***jaj, ***kaj;
527
528 PetscInt lxs, lxe, lys, lye, lzs, lze;
529
530 PetscReal ajc;
531
532 PetscReal dudc, dude, dudz, dvdc, dvde, dvdz, dwdc, dwde, dwdz;
533 PetscReal csi0, csi1, csi2, eta0, eta1, eta2, zet0, zet1, zet2;
534 PetscReal g11, g21, g31;
535 PetscReal r11, r21, r31, r12, r22, r32, r13, r23, r33;
536
537 PetscScalar solid,innerblank;
538
539 // --- CONTEXT ACQUISITION BLOCK ---
540 // Get the master simulation context from the UserCtx.
541 SimCtx *simCtx = user->simCtx;
542
543 // Create local variables to mirror the legacy globals for minimal code changes.
544 const LESModelType les = simCtx->les;
545 const PetscInt rans = simCtx->rans;
546 const PetscInt ti = simCtx->step; // Assuming simCtx->step is the new integer time counter
547 const PetscReal ren = simCtx->ren;
548 const PetscInt clark = simCtx->clark;
549 solid = 0.5;
550 innerblank = 7.;
551
552 PROFILE_FUNCTION_BEGIN;
553
554 DMDAVecGetArray(fda, Ucont, &ucont);
555 DMDAVecGetArray(fda, Ucat, &ucat);
556 DMDAVecGetArray(fda, Visc, &visc);
557
558 DMDAVecGetArray(fda, Csi, &csi);
559 DMDAVecGetArray(fda, Eta, &eta);
560 DMDAVecGetArray(fda, Zet, &zet);
561
562 DMDAVecGetArray(fda, user->lICsi, &icsi);
563 DMDAVecGetArray(fda, user->lIEta, &ieta);
564 DMDAVecGetArray(fda, user->lIZet, &izet);
565
566 DMDAVecGetArray(fda, user->lJCsi, &jcsi);
567 DMDAVecGetArray(fda, user->lJEta, &jeta);
568 DMDAVecGetArray(fda, user->lJZet, &jzet);
569
570 DMDAVecGetArray(fda, user->lKCsi, &kcsi);
571 DMDAVecGetArray(fda, user->lKEta, &keta);
572 DMDAVecGetArray(fda, user->lKZet, &kzet);
573
574 DMDAVecGetArray(da, user->lNvert, &nvert);
575
576 VecDuplicate(Ucont, &Fp1);
577 VecDuplicate(Ucont, &Fp2);
578 VecDuplicate(Ucont, &Fp3);
579
580 DMDAVecGetArray(fda, Fp1, &fp1);
581 DMDAVecGetArray(fda, Fp2, &fp2);
582 DMDAVecGetArray(fda, Fp3, &fp3);
583
584 DMDAVecGetArray(da, user->lAj, &aj);
585
586 DMDAGetLocalInfo(da, &info);
587
588 mx = info.mx; my = info.my; mz = info.mz;
589 xs = info.xs; xe = xs + info.xm;
590 ys = info.ys; ye = ys + info.ym;
591 zs = info.zs; ze = zs + info.zm;
592
593 /* First we calculate the flux on cell surfaces. Stored on the upper integer
594 node. For example, along i direction, the flux are stored at node 0:mx-2*/
595 lxs = xs; lxe = xe;
596 lys = ys; lye = ye;
597 lzs = zs; lze = ze;
598
599 if (xs==0) lxs = xs+1;
600 if (ys==0) lys = ys+1;
601 if (zs==0) lzs = zs+1;
602
603
604 if (xe==mx) lxe=xe-1;
605 if (ye==my) lye=ye-1;
606 if (ze==mz) lze=ze-1;
607
608 VecSet(Visc,0.0);
609
610 PetscReal ***lnu_t;
611
612 if(les) {
613 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
614 } else if (rans) {
615
616 DMDAVecGetArray(da, user->lNu_t, &lnu_t);
617 }
618
619 /* The visc flux on each surface center is stored at previous integer node */
620
621 DMDAVecGetArray(da, user->lIAj, &iaj);
622 /* for (k=zs; k<ze; k++) { */
623/* for (j=ys; j<ye; j++) { */
624/* for (i=xs; i<xe; i++) { */
625/* if (i==1 && (j==0 ||j==1 || j==2) && (k==21 || k==22|| k==20)) */
626/* 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 ); */
627/* } */
628/* } */
629/* } */
630 // i direction
631 for (k=lzs; k<lze; k++) {
632 for (j=lys; j<lye; j++) {
633 for (i=lxs-1; i<lxe; i++) {
634
635 dudc = ucat[k][j][i+1].x - ucat[k][j][i].x;
636 dvdc = ucat[k][j][i+1].y - ucat[k][j][i].y;
637 dwdc = ucat[k][j][i+1].z - ucat[k][j][i].z;
638
639 if ((nvert[k][j+1][i ]> solid && nvert[k][j+1][i ]<innerblank) ||
640 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
641 dude = (ucat[k][j ][i+1].x + ucat[k][j ][i].x -
642 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.5;
643 dvde = (ucat[k][j ][i+1].y + ucat[k][j ][i].y -
644 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.5;
645 dwde = (ucat[k][j ][i+1].z + ucat[k][j ][i].z -
646 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.5;
647 }
648 else if ((nvert[k][j-1][i ]> solid && nvert[k][j-1][i ]<innerblank) ||
649 (nvert[k][j-1][i+1]> solid && nvert[k][j-1][i+1]<innerblank)) {
650 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
651 ucat[k][j ][i+1].x - ucat[k][j ][i].x) * 0.5;
652 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
653 ucat[k][j ][i+1].y - ucat[k][j ][i].y) * 0.5;
654 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
655 ucat[k][j ][i+1].z - ucat[k][j ][i].z) * 0.5;
656 }
657 else {
658 dude = (ucat[k][j+1][i+1].x + ucat[k][j+1][i].x -
659 ucat[k][j-1][i+1].x - ucat[k][j-1][i].x) * 0.25;
660 dvde = (ucat[k][j+1][i+1].y + ucat[k][j+1][i].y -
661 ucat[k][j-1][i+1].y - ucat[k][j-1][i].y) * 0.25;
662 dwde = (ucat[k][j+1][i+1].z + ucat[k][j+1][i].z -
663 ucat[k][j-1][i+1].z - ucat[k][j-1][i].z) * 0.25;
664 }
665
666 if ((nvert[k+1][j][i ]> solid && nvert[k+1][j][i ]<innerblank)||
667 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
668 dudz = (ucat[k ][j][i+1].x + ucat[k ][j][i].x -
669 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.5;
670 dvdz = (ucat[k ][j][i+1].y + ucat[k ][j][i].y -
671 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.5;
672 dwdz = (ucat[k ][j][i+1].z + ucat[k ][j][i].z -
673 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.5;
674 }
675 else if ((nvert[k-1][j][i ]> solid && nvert[k-1][j][i ]<innerblank) ||
676 (nvert[k-1][j][i+1]> solid && nvert[k-1][j][i+1]<innerblank)) {
677
678 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
679 ucat[k ][j][i+1].x - ucat[k ][j][i].x) * 0.5;
680 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
681 ucat[k ][j][i+1].y - ucat[k ][j][i].y) * 0.5;
682 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
683 ucat[k ][j][i+1].z - ucat[k ][j][i].z) * 0.5;
684 }
685 else {
686 dudz = (ucat[k+1][j][i+1].x + ucat[k+1][j][i].x -
687 ucat[k-1][j][i+1].x - ucat[k-1][j][i].x) * 0.25;
688 dvdz = (ucat[k+1][j][i+1].y + ucat[k+1][j][i].y -
689 ucat[k-1][j][i+1].y - ucat[k-1][j][i].y) * 0.25;
690 dwdz = (ucat[k+1][j][i+1].z + ucat[k+1][j][i].z -
691 ucat[k-1][j][i+1].z - ucat[k-1][j][i].z) * 0.25;
692 }
693
694 csi0 = icsi[k][j][i].x;
695 csi1 = icsi[k][j][i].y;
696 csi2 = icsi[k][j][i].z;
697
698 eta0 = ieta[k][j][i].x;
699 eta1 = ieta[k][j][i].y;
700 eta2 = ieta[k][j][i].z;
701
702 zet0 = izet[k][j][i].x;
703 zet1 = izet[k][j][i].y;
704 zet2 = izet[k][j][i].z;
705
706 g11 = csi0 * csi0 + csi1 * csi1 + csi2 * csi2;
707 g21 = eta0 * csi0 + eta1 * csi1 + eta2 * csi2;
708 g31 = zet0 * csi0 + zet1 * csi1 + zet2 * csi2;
709
710 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
711 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
712 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
713
714 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
715 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
716 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
717
718 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
719 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
720 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
721
722 ajc = iaj[k][j][i];
723
724 double nu = 1./ren, nu_t=0;
725
726 if( les || (rans && ti>0) ) {
727 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j][i+1]) ), 2.0) * Sabs;
728 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j][i+1]);
729 if ( (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == WALL && i==0) || (user->boundary_faces[BC_FACE_POS_X].mathematical_type == WALL && i==mx-2) ) nu_t=0;
730 fp1[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * csi0 + r21 * csi1 + r31 * csi2) * ajc * (nu_t);
731 fp1[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * csi0 + r22 * csi1 + r32 * csi2) * ajc * (nu_t);
732 fp1[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * csi0 + r23 * csi1 + r33 * csi2) * ajc * (nu_t);
733 }
734 else {
735 fp1[k][j][i].x = 0;
736 fp1[k][j][i].y = 0;
737 fp1[k][j][i].z = 0;
738 }
739
740 fp1[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * csi0 + r21 * csi1 + r31 * csi2 ) * ajc * (nu);
741 fp1[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * csi0 + r22 * csi1 + r32 * csi2 ) * ajc * (nu);
742 fp1[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * csi0 + r23 * csi1 + r33 * csi2 ) * ajc * (nu);
743
744
745 if(clark) {
746 double dc, de, dz;
747 ComputeCellCharacteristicLengthScale (ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
748 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
749
750 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
751 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
752 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
753 double t21 = t12;
754 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
755 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
756 double t31 = t13;
757 double t32 = t23;
758 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
759
760 fp1[k][j][i].x -= ( t11 * csi0 + t12 * csi1 + t13 * csi2 ) / 12.;
761 fp1[k][j][i].y -= ( t21 * csi0 + t22 * csi1 + t23 * csi2 ) / 12.;
762 fp1[k][j][i].z -= ( t31 * csi0 + t32 * csi1 + t33 * csi2 ) / 12.;
763 }
764
765 }
766 }
767 }
768 DMDAVecRestoreArray(da, user->lIAj, &iaj);
769
770
771 // j direction
772 DMDAVecGetArray(da, user->lJAj, &jaj);
773 for (k=lzs; k<lze; k++) {
774 for (j=lys-1; j<lye; j++) {
775 for (i=lxs; i<lxe; i++) {
776
777 if ((nvert[k][j ][i+1]> solid && nvert[k][j ][i+1]<innerblank)||
778 (nvert[k][j+1][i+1]> solid && nvert[k][j+1][i+1]<innerblank)) {
779 dudc = (ucat[k][j+1][i ].x + ucat[k][j][i ].x -
780 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.5;
781 dvdc = (ucat[k][j+1][i ].y + ucat[k][j][i ].y -
782 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.5;
783 dwdc = (ucat[k][j+1][i ].z + ucat[k][j][i ].z -
784 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.5;
785 }
786 else if ((nvert[k][j ][i-1]> solid && nvert[k][j ][i-1]<innerblank) ||
787 (nvert[k][j+1][i-1]> solid && nvert[k][j+1][i-1]<innerblank)) {
788 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
789 ucat[k][j+1][i ].x - ucat[k][j][i ].x) * 0.5;
790 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
791 ucat[k][j+1][i ].y - ucat[k][j][i ].y) * 0.5;
792 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
793 ucat[k][j+1][i ].z - ucat[k][j][i ].z) * 0.5;
794 }
795 else {
796 dudc = (ucat[k][j+1][i+1].x + ucat[k][j][i+1].x -
797 ucat[k][j+1][i-1].x - ucat[k][j][i-1].x) * 0.25;
798 dvdc = (ucat[k][j+1][i+1].y + ucat[k][j][i+1].y -
799 ucat[k][j+1][i-1].y - ucat[k][j][i-1].y) * 0.25;
800 dwdc = (ucat[k][j+1][i+1].z + ucat[k][j][i+1].z -
801 ucat[k][j+1][i-1].z - ucat[k][j][i-1].z) * 0.25;
802 }
803
804 dude = ucat[k][j+1][i].x - ucat[k][j][i].x;
805 dvde = ucat[k][j+1][i].y - ucat[k][j][i].y;
806 dwde = ucat[k][j+1][i].z - ucat[k][j][i].z;
807
808 if ((nvert[k+1][j ][i]> solid && nvert[k+1][j ][i]<innerblank)||
809 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
810 dudz = (ucat[k ][j+1][i].x + ucat[k ][j][i].x -
811 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.5;
812 dvdz = (ucat[k ][j+1][i].y + ucat[k ][j][i].y -
813 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.5;
814 dwdz = (ucat[k ][j+1][i].z + ucat[k ][j][i].z -
815 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.5;
816 }
817 else if ((nvert[k-1][j ][i]> solid && nvert[k-1][j ][i]<innerblank)||
818 (nvert[k-1][j+1][i]> solid && nvert[k-1][j+1][i]<innerblank)) {
819 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
820 ucat[k ][j+1][i].x - ucat[k ][j][i].x) * 0.5;
821 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
822 ucat[k ][j+1][i].y - ucat[k ][j][i].y) * 0.5;
823 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
824 ucat[k ][j+1][i].z - ucat[k ][j][i].z) * 0.5;
825 }
826 else {
827 dudz = (ucat[k+1][j+1][i].x + ucat[k+1][j][i].x -
828 ucat[k-1][j+1][i].x - ucat[k-1][j][i].x) * 0.25;
829 dvdz = (ucat[k+1][j+1][i].y + ucat[k+1][j][i].y -
830 ucat[k-1][j+1][i].y - ucat[k-1][j][i].y) * 0.25;
831 dwdz = (ucat[k+1][j+1][i].z + ucat[k+1][j][i].z -
832 ucat[k-1][j+1][i].z - ucat[k-1][j][i].z) * 0.25;
833 }
834
835 csi0 = jcsi[k][j][i].x;
836 csi1 = jcsi[k][j][i].y;
837 csi2 = jcsi[k][j][i].z;
838
839 eta0 = jeta[k][j][i].x;
840 eta1 = jeta[k][j][i].y;
841 eta2 = jeta[k][j][i].z;
842
843 zet0 = jzet[k][j][i].x;
844 zet1 = jzet[k][j][i].y;
845 zet2 = jzet[k][j][i].z;
846
847
848 g11 = csi0 * eta0 + csi1 * eta1 + csi2 * eta2;
849 g21 = eta0 * eta0 + eta1 * eta1 + eta2 * eta2;
850 g31 = zet0 * eta0 + zet1 * eta1 + zet2 * eta2;
851
852 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
853 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
854 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
855
856 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
857 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
858 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
859
860 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
861 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
862 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
863
864 // 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);
865 // 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);
866 // 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);
867 // 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);
868
869
870
871 ajc = jaj[k][j][i];
872
873 double nu = 1./ren, nu_t = 0;
874
875 if( les || (rans && ti>0) ) {
876 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k][j+1][i]) ), 2.0) * Sabs;
877 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k][j+1][i]);
878 if ( (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == WALL && j==0) || (user->boundary_faces[BC_FACE_POS_Y].mathematical_type == WALL && j==my-2) ) nu_t=0;
879
880 fp2[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * eta0 + r21 * eta1 + r31 * eta2) * ajc * (nu_t);
881 fp2[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * eta0 + r22 * eta1 + r32 * eta2) * ajc * (nu_t);
882 fp2[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * eta0 + r23 * eta1 + r33 * eta2) * ajc * (nu_t);
883 }
884 else {
885 fp2[k][j][i].x = 0;
886 fp2[k][j][i].y = 0;
887 fp2[k][j][i].z = 0;
888 }
889
890 fp2[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz+ r11 * eta0 + r21 * eta1 + r31 * eta2 ) * ajc * (nu);
891 fp2[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz+ r12 * eta0 + r22 * eta1 + r32 * eta2 ) * ajc * (nu);
892 fp2[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz+ r13 * eta0 + r23 * eta1 + r33 * eta2 ) * ajc * (nu);
893
894 if(clark) {
895 double dc, de, dz;
896 ComputeCellCharacteristicLengthScale(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
897 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
898
899 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
900 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
901 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
902 double t21 = t12;
903 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
904 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
905 double t31 = t13;
906 double t32 = t23;
907 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
908
909 fp2[k][j][i].x -= ( t11 * eta0 + t12 * eta1 + t13 * eta2 ) / 12.;
910 fp2[k][j][i].y -= ( t21 * eta0 + t22 * eta1 + t23 * eta2 ) / 12.;
911 fp2[k][j][i].z -= ( t31 * eta0 + t32 * eta1 + t33 * eta2 ) / 12.;
912 }
913 }
914 }
915 }
916
917 DMDAVecRestoreArray(da, user->lJAj, &jaj);
918 // k direction
919
920 DMDAVecGetArray(da, user->lKAj, &kaj);
921 for (k=lzs-1; k<lze; k++) {
922 for (j=lys; j<lye; j++) {
923 for (i=lxs; i<lxe; i++) {
924 if ((nvert[k ][j][i+1]> solid && nvert[k ][j][i+1]<innerblank)||
925 (nvert[k+1][j][i+1]> solid && nvert[k+1][j][i+1]<innerblank)) {
926 dudc = (ucat[k+1][j][i ].x + ucat[k][j][i ].x -
927 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.5;
928 dvdc = (ucat[k+1][j][i ].y + ucat[k][j][i ].y -
929 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.5;
930 dwdc = (ucat[k+1][j][i ].z + ucat[k][j][i ].z -
931 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.5;
932 }
933 else if ((nvert[k ][j][i-1]> solid && nvert[k ][j][i-1]<innerblank) ||
934 (nvert[k+1][j][i-1]> solid && nvert[k+1][j][i-1]<innerblank)) {
935 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
936 ucat[k+1][j][i ].x - ucat[k][j][i ].x) * 0.5;
937 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
938 ucat[k+1][j][i ].y - ucat[k][j][i ].y) * 0.5;
939 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
940 ucat[k+1][j][i ].z - ucat[k][j][i ].z) * 0.5;
941 }
942 else {
943 dudc = (ucat[k+1][j][i+1].x + ucat[k][j][i+1].x -
944 ucat[k+1][j][i-1].x - ucat[k][j][i-1].x) * 0.25;
945 dvdc = (ucat[k+1][j][i+1].y + ucat[k][j][i+1].y -
946 ucat[k+1][j][i-1].y - ucat[k][j][i-1].y) * 0.25;
947 dwdc = (ucat[k+1][j][i+1].z + ucat[k][j][i+1].z -
948 ucat[k+1][j][i-1].z - ucat[k][j][i-1].z) * 0.25;
949 }
950
951 if ((nvert[k ][j+1][i]> solid && nvert[k ][j+1][i]<innerblank)||
952 (nvert[k+1][j+1][i]> solid && nvert[k+1][j+1][i]<innerblank)) {
953 dude = (ucat[k+1][j ][i].x + ucat[k][j ][i].x -
954 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.5;
955 dvde = (ucat[k+1][j ][i].y + ucat[k][j ][i].y -
956 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.5;
957 dwde = (ucat[k+1][j ][i].z + ucat[k][j ][i].z -
958 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.5;
959 }
960 else if ((nvert[k ][j-1][i]> solid && nvert[k ][j-1][i]<innerblank) ||
961 (nvert[k+1][j-1][i]> solid && nvert[k+1][j-1][i]<innerblank)){
962 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
963 ucat[k+1][j ][i].x - ucat[k][j ][i].x) * 0.5;
964 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
965 ucat[k+1][j ][i].y - ucat[k][j ][i].y) * 0.5;
966 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
967 ucat[k+1][j ][i].z - ucat[k][j ][i].z) * 0.5;
968 }
969 else {
970 dude = (ucat[k+1][j+1][i].x + ucat[k][j+1][i].x -
971 ucat[k+1][j-1][i].x - ucat[k][j-1][i].x) * 0.25;
972 dvde = (ucat[k+1][j+1][i].y + ucat[k][j+1][i].y -
973 ucat[k+1][j-1][i].y - ucat[k][j-1][i].y) * 0.25;
974 dwde = (ucat[k+1][j+1][i].z + ucat[k][j+1][i].z -
975 ucat[k+1][j-1][i].z - ucat[k][j-1][i].z) * 0.25;
976 }
977
978 dudz = ucat[k+1][j][i].x - ucat[k][j][i].x;
979 dvdz = ucat[k+1][j][i].y - ucat[k][j][i].y;
980 dwdz = ucat[k+1][j][i].z - ucat[k][j][i].z;
981
982
983 csi0 = kcsi[k][j][i].x;
984 csi1 = kcsi[k][j][i].y;
985 csi2 = kcsi[k][j][i].z;
986
987 eta0 = keta[k][j][i].x;
988 eta1 = keta[k][j][i].y;
989 eta2 = keta[k][j][i].z;
990
991 zet0 = kzet[k][j][i].x;
992 zet1 = kzet[k][j][i].y;
993 zet2 = kzet[k][j][i].z;
994
995
996 g11 = csi0 * zet0 + csi1 * zet1 + csi2 * zet2;
997 g21 = eta0 * zet0 + eta1 * zet1 + eta2 * zet2;
998 g31 = zet0 * zet0 + zet1 * zet1 + zet2 * zet2;
999
1000 r11 = dudc * csi0 + dude * eta0 + dudz * zet0;
1001 r21 = dvdc * csi0 + dvde * eta0 + dvdz * zet0;
1002 r31 = dwdc * csi0 + dwde * eta0 + dwdz * zet0;
1003
1004 r12 = dudc * csi1 + dude * eta1 + dudz * zet1;
1005 r22 = dvdc * csi1 + dvde * eta1 + dvdz * zet1;
1006 r32 = dwdc * csi1 + dwde * eta1 + dwdz * zet1;
1007
1008 r13 = dudc * csi2 + dude * eta2 + dudz * zet2;
1009 r23 = dvdc * csi2 + dvde * eta2 + dvdz * zet2;
1010 r33 = dwdc * csi2 + dwde * eta2 + dwdz * zet2;
1011
1012 ajc = kaj[k][j][i];
1013
1014 double nu = 1./ren, nu_t =0;
1015
1016 if( les || (rans && ti>0) ) {
1017 //nu_t = pow( 0.5 * ( sqrt(lnu_t[k][j][i]) + sqrt(lnu_t[k+1][j][i]) ), 2.0) * Sabs;
1018 nu_t = 0.5 * (lnu_t[k][j][i] + lnu_t[k+1][j][i]);
1019 if ( (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == WALL && k==0) || (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == WALL && k==mz-2) ) nu_t=0;
1020
1021 fp3[k][j][i].x = (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu_t);
1022 fp3[k][j][i].y = (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu_t);
1023 fp3[k][j][i].z = (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu_t);
1024 }
1025 else {
1026 fp3[k][j][i].x = 0;
1027 fp3[k][j][i].y = 0;
1028 fp3[k][j][i].z = 0;
1029 }
1030 fp3[k][j][i].x += (g11 * dudc + g21 * dude + g31 * dudz + r11 * zet0 + r21 * zet1 + r31 * zet2) * ajc * (nu);//
1031 fp3[k][j][i].y += (g11 * dvdc + g21 * dvde + g31 * dvdz + r12 * zet0 + r22 * zet1 + r32 * zet2) * ajc * (nu);//
1032 fp3[k][j][i].z += (g11 * dwdc + g21 * dwde + g31 * dwdz + r13 * zet0 + r23 * zet1 + r33 * zet2) * ajc * (nu);//
1033
1034 if(clark) {
1035 double dc, de, dz;
1036 ComputeCellCharacteristicLengthScale(ajc, csi[k][j][i], eta[k][j][i], zet[k][j][i], &dc, &de, &dz);
1037 double dc2=dc*dc, de2=de*de, dz2=dz*dz;
1038
1039 double t11 = ( dudc * dudc * dc2 + dude * dude * de2 + dudz * dudz * dz2 );
1040 double t12 = ( dudc * dvdc * dc2 + dude * dvde * de2 + dudz * dvdz * dz2 );
1041 double t13 = ( dudc * dwdc * dc2 + dude * dwde * de2 + dudz * dwdz * dz2 );
1042 double t21 = t12;
1043 double t22 = ( dvdc * dvdc * dc2 + dvde * dvde * de2 + dvdz * dvdz * dz2 );
1044 double t23 = ( dvdc * dwdc * dc2 + dvde * dwde * de2 + dvdz * dwdz * dz2 );
1045 double t31 = t13;
1046 double t32 = t23;
1047 double t33 = ( dwdc * dwdc * dc2 + dwde * dwde * de2 + dwdz * dwdz * dz2 );
1048
1049 fp3[k][j][i].x -= ( t11 * zet0 + t12 * zet1 + t13 * zet2 ) / 12.;
1050 fp3[k][j][i].y -= ( t21 * zet0 + t22 * zet1 + t23 * zet2 ) / 12.;
1051 fp3[k][j][i].z -= ( t31 * zet0 + t32 * zet1 + t33 * zet2 ) / 12.;
1052 }
1053 }
1054 }
1055 }
1056
1057 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1058
1059 for (k=lzs; k<lze; k++) {
1060 for (j=lys; j<lye; j++) {
1061 for (i=lxs; i<lxe; i++) {
1062 visc[k][j][i].x =
1063 (fp1[k][j][i].x - fp1[k][j][i-1].x +
1064 fp2[k][j][i].x - fp2[k][j-1][i].x +
1065 fp3[k][j][i].x - fp3[k-1][j][i].x);
1066
1067 visc[k][j][i].y =
1068 (fp1[k][j][i].y - fp1[k][j][i-1].y +
1069 fp2[k][j][i].y - fp2[k][j-1][i].y +
1070 fp3[k][j][i].y - fp3[k-1][j][i].y);
1071
1072 visc[k][j][i].z =
1073 (fp1[k][j][i].z - fp1[k][j][i-1].z +
1074 fp2[k][j][i].z - fp2[k][j-1][i].z +
1075 fp3[k][j][i].z - fp3[k-1][j][i].z);
1076
1077 }
1078 }
1079 }
1080/* for (k=zs; k<ze; k++) { */
1081/* for (j=ys; j<ye; j++) { */
1082/* for (i=xs; i<xe; i++) { */
1083/* 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); */
1084/* 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); */
1085/* 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); */
1086/* 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); */
1087/* 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); */
1088/* 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); */
1089
1090/* } */
1091/* } */
1092/* } */
1093 DMDAVecRestoreArray(fda, Ucont, &ucont);
1094 DMDAVecRestoreArray(fda, Ucat, &ucat);
1095 DMDAVecRestoreArray(fda, Visc, &visc);
1096
1097 DMDAVecRestoreArray(fda, Csi, &csi);
1098 DMDAVecRestoreArray(fda, Eta, &eta);
1099 DMDAVecRestoreArray(fda, Zet, &zet);
1100
1101 DMDAVecRestoreArray(fda, Fp1, &fp1);
1102 DMDAVecRestoreArray(fda, Fp2, &fp2);
1103 DMDAVecRestoreArray(fda, Fp3, &fp3);
1104
1105 DMDAVecRestoreArray(da, user->lAj, &aj);
1106
1107 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1108 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1109 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1110
1111 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1112 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1113 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1114
1115 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1116 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1117 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1118
1119 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1120
1121 if(les) {
1122 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1123 } else if (rans) {
1124
1125 DMDAVecRestoreArray(da, user->lNu_t, &lnu_t);
1126 }
1127
1128
1129 VecDestroy(&Fp1);
1130 VecDestroy(&Fp2);
1131 VecDestroy(&Fp3);
1132
1133
1134 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Viscous terms calculated .\n");
1135
1136 PROFILE_FUNCTION_END;
1137
1138 return(0);
1139}
1140
1141#undef __FUNCT__
1142#define __FUNCT__ "ComputeBodyForces"
1143/**
1144 * @brief General dispatcher for applying all active body forces (momentum sources).
1145 *
1146 * This function serves as a central hub for adding momentum source terms to the
1147 * contravariant right-hand-side (Rct) of the momentum equations. It is called once per RHS
1148 * evaluation (e.g., once per Runge-Kutta stage).
1149 *
1150 * It introspects the simulation configuration to determine which, if any, body
1151 * forces are active and calls their specific implementation functions.
1152 *
1153 * @param user The UserCtx containing the simulation state for a single block.
1154 * @param Rhs The PETSc Vec for the contravariant RHS, which will be modified in-place.
1155 * @return PetscErrorCode 0 on success.
1156 */
1157PetscErrorCode ComputeBodyForces(UserCtx *user, Vec Rct)
1158{
1159 PetscErrorCode ierr;
1160 PetscFunctionBeginUser;
1161
1162 // --- 1. Apply momentum source for driven channel/pipe flows ---
1163 // This function will internally check if a driven flow BC is active.
1164 ierr = ComputeDrivenChannelFlowSource(user, Rct); CHKERRQ(ierr);
1165
1166 // --- 2. (Future Extension) Apply gravitational force ---
1167 // if (user->simCtx->gravityEnabled) {
1168 // ierr = ApplyGravitationalForce(user, Rhs); CHKERRQ(ierr);
1169 // }
1170
1171 PetscFunctionReturn(0);
1172}
1173
1174#undef __FUNCT__
1175#define __FUNCT__ "ComputeRHS"
1176/**
1177 * @brief Computes the Right-Hand-Side (RHS) of the momentum equations.
1178 *
1179 * This function orchestrates the calculation of the spatial discretization of the
1180 * convective and diffusive terms. It calls specialized helper functions
1181 * (`Convection`, `Viscous`) to compute these terms and then combines them with
1182 * the pressure gradient to form the final RHS vector.
1183 *
1184 * This is a minimally-edited version of the legacy kernel. It operates on a
1185 * single UserCtx (one grid block) and retrieves all global parameters
1186 * (Re, rans, les, etc.) from the master SimCtx.
1187 *
1188 * @param user The UserCtx for the specific block being computed.
1189 * @param Rhs The PETSc Vec where the calculated RHS will be stored.
1190 * @return PetscErrorCode 0 on success.
1191 */
1192PetscErrorCode ComputeRHS(UserCtx *user, Vec Rhs)
1193{
1194 PetscErrorCode ierr;
1195 SimCtx *simCtx = user->simCtx;
1196 DM da = user->da, fda = user->fda;
1197 DMDALocalInfo info = user->info;
1198 PetscInt i,j,k;
1199 PetscReal dpdc, dpde,dpdz;
1200 // --- Local Grid Indices and Parameters ---
1201 PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
1202 PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
1203 PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
1204 PetscInt lxs = (xs==0) ? xs+1 : xs;
1205 PetscInt lys = (ys==0) ? ys+1 : ys;
1206 PetscInt lzs = (zs==0) ? zs+1 : zs;
1207 PetscInt lxe = (xe==mx) ? xe-1 : xe;
1208 PetscInt lye = (ye==my) ? ye-1 : ye;
1209 PetscInt lze = (ze==mz) ? ze-1 : ze;
1210 PetscInt mz_end = (user->boundary_faces[BC_FACE_POS_Z].mathematical_type == CHARACTERISTIC_BC) ? mz - 2 : mz - 3;
1211
1212 // --- Array Pointers ---
1213 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1214 PetscReal ***p, ***iaj, ***jaj, ***kaj, ***aj, ***nvert, ***nvert_o;
1215 Cmpnts ***rhs, ***rc, ***rct,***ucont;
1216
1217 // --- Temporary Vectors ---
1218 Vec Conv, Visc, Rc, Rct;
1219
1220 PetscFunctionBeginUser;
1221 PROFILE_FUNCTION_BEGIN;
1222 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Computing RHS (FormFunction1)...\n",
1223 simCtx->rank, user->_this);
1224
1225 // --- Get all necessary array pointers ---
1226 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1227 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1228 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1229 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1230 ierr = DMDAVecGetArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1231 ierr = DMDAVecGetArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1232 ierr = DMDAVecGetArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1233 ierr = DMDAVecGetArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1234 ierr = DMDAVecGetArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1235 ierr = DMDAVecGetArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1236 ierr = DMDAVecGetArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1237 ierr = DMDAVecGetArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1238 ierr = DMDAVecGetArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1239 ierr = DMDAVecGetArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1240 ierr = DMDAVecGetArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1241 ierr = DMDAVecGetArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1242 ierr = DMDAVecGetArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1243 ierr = DMDAVecGetArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1244 ierr = DMDAVecGetArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1245
1246 // --- Create temporary work vectors ---
1247 ierr = VecDuplicate(user->lUcont, &Rc); CHKERRQ(ierr);
1248 ierr = VecDuplicate(Rc, &Rct); CHKERRQ(ierr);
1249 ierr = VecDuplicate(Rct, &Conv); CHKERRQ(ierr);
1250 ierr = VecDuplicate(Rct, &Visc); CHKERRQ(ierr);
1251
1252 // ========================================================================
1253 // CORE LOGIC (UNCHANGED FROM LEGACY CODE)
1254 // ========================================================================
1255
1256 // 1. Obtain Cartesian velocity from Contravariant velocity
1257 ierr = Contra2Cart(user); CHKERRQ(ierr);
1258 ierr = UpdateLocalGhosts(user,"Ucat"); CHKERRQ(ierr);
1259
1260 // 2. Compute Convective term
1261 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating convective terms...\n");
1262 if (simCtx->moveframe || simCtx->rotateframe) {
1263 // ierr = Convection_MV(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1264 } else {
1265 ierr = Convection(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1266 }
1267
1268 // 3. Compute Viscous term
1269 if (simCtx->invicid) {
1270 ierr = VecSet(Visc, 0.0); CHKERRQ(ierr);
1271 } else {
1272 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating viscous terms...\n");
1273 ierr = Viscous(user, user->lUcont, user->lUcat, Visc); CHKERRQ(ierr);
1274 }
1275
1276 // 4. Combine terms to get Cartesian RHS: Rc = Visc - Conv
1277 ierr = VecWAXPY(Rc, -1.0, Conv, Visc); CHKERRQ(ierr);
1278
1279 // 5. Convert Cartesian RHS (Rc) to Contravariant RHS (Rct)
1280 LOG_ALLOW(LOCAL, LOG_DEBUG, " Converting Cartesian RHS to Contravariant RHS...\n");
1281 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1282 ierr = DMDAVecGetArray(fda, Rc, &rc); CHKERRQ(ierr);
1283
1284 for (k = lzs; k < lze; k++) {
1285 for (j = lys; j < lye; j++) {
1286 for (i = lxs; i < lxe; i++) {
1287 rct[k][j][i].x = aj[k][j][i] *
1288 (0.5 * (csi[k][j][i].x + csi[k][j][i-1].x) * rc[k][j][i].x +
1289 0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
1290 0.5 * (csi[k][j][i].z + csi[k][j][i-1].z) * rc[k][j][i].z);
1291 rct[k][j][i].y = aj[k][j][i] *
1292 (0.5 * (eta[k][j][i].x + eta[k][j-1][i].x) * rc[k][j][i].x +
1293 0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
1294 0.5 * (eta[k][j][i].z + eta[k][j-1][i].z) * rc[k][j][i].z);
1295 rct[k][j][i].z = aj[k][j][i] *
1296 (0.5 * (zet[k][j][i].x + zet[k-1][j][i].x) * rc[k][j][i].x +
1297 0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
1298 0.5 * (zet[k][j][i].z + zet[k-1][j][i].z) * rc[k][j][i].z);
1299 }
1300 }
1301 }
1302 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1303 ierr = DMDAVecRestoreArray(fda, Rc, &rc); CHKERRQ(ierr);
1304
1305 PetscBarrier(NULL);
1306
1307 DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct);
1308 DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct);
1309
1310 // Compute and Add Body Force term if applicable.
1311 ierr = ComputeBodyForces(user,Rct);
1312
1313 DMDAVecGetArray(fda, Rct, &rct);
1314
1315
1316 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_X].mathematical_type == PERIODIC){
1317 if (xs==0){
1318 i=xs;
1319 for (k=lzs; k<lze; k++) {
1320 for (j=lys; j<lye; j++) {
1321 rct[k][j][i].x=rct[k][j][i-2].x;
1322 }
1323 }
1324 }
1325
1326 if (xe==mx){
1327 i=mx-1;
1328 for (k=lzs; k<lze; k++) {
1329 for (j=lys; j<lye; j++) {
1330 rct[k][j][i].x=rct[k][j][i+2].x;
1331 }
1332 }
1333 }
1334 }
1335
1336 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Y].mathematical_type == PERIODIC){
1337 if (ys==0){
1338 j=ys;
1339 for (k=lzs; k<lze; k++) {
1340 for (i=lxs; i<lxe; i++) {
1341 rct[k][j][i].y=rct[k][j-2][i].y;
1342 }
1343 }
1344 }
1345
1346 if (ye==my){
1347 j=my-1;
1348 for (k=lzs; k<lze; k++) {
1349 for (i=lxs; i<lxe; i++) {
1350 rct[k][j][i].y=rct[k][j+2][i].y;
1351 }
1352 }
1353 }
1354 }
1355 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC || user->boundary_faces[BC_FACE_POS_Z].mathematical_type == PERIODIC){
1356 if (zs==0){
1357 k=zs;
1358 for (j=lys; j<lye; j++) {
1359 for (i=lxs; i<lxe; i++) {
1360 rct[k][j][i].z=rct[k-2][j][i].z;
1361 }
1362 }
1363 }
1364 if (ze==mz){
1365 k=mz-1;
1366 for (j=lys; j<lye; j++) {
1367 for (i=lxs; i<lxe; i++) {
1368 rct[k][j][i].z=rct[k+2][j][i].z;
1369 }
1370 }
1371 }
1372 }
1373
1374 // 6. Add Pressure Gradient Term and Finalize RHS
1375 // This involves calculating pressure derivatives (dpdc, dpde, dpdz) and using
1376 // them to adjust the contravariant RHS. The full stencil logic is preserved.
1377 LOG_ALLOW(LOCAL, LOG_DEBUG, " Adding pressure gradient term to RHS...\n");
1378
1379 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1380
1381 ierr = DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1382 ierr = DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1383
1384 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1385
1386 for (k = lzs; k < lze; k++) {
1387 for (j = lys; j < lye; j++) {
1388 for (i = lxs; i < lxe; i++) {
1389 PetscReal dpdc, dpde, dpdz;
1390 dpdc = p[k][j][i+1] - p[k][j][i];
1391
1392 if ((j==my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1393 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1394 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1395 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1396 }
1397 }
1398 else if ((j==my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1399 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1400 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1401 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1402 }
1403 }
1404 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1405 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1406 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1407 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1408 }
1409 }
1410 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1411 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1412 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1413 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1414 }
1415 }
1416 else {
1417 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1418 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1419 }
1420
1421 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1422 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1423 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1424 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1425 }
1426 }
1427 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1428 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1429 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1430 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1431 }
1432 }
1433 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1434 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1435 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1436 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1437 }
1438 }
1439 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1440 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1441 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1442 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1443 }
1444 }
1445 else {
1446 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1447 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1448 }
1449
1450 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1451
1452
1453 rhs[k][j][i].x -=
1454 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1455 icsi[k][j][i].y * icsi[k][j][i].y +
1456 icsi[k][j][i].z * icsi[k][j][i].z)+
1457 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1458 ieta[k][j][i].y * icsi[k][j][i].y +
1459 ieta[k][j][i].z * icsi[k][j][i].z)+
1460 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1461 izet[k][j][i].y * icsi[k][j][i].y +
1462 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1463
1464 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1465 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1466 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1467 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1468 }
1469 }
1470 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1471 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1472 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1473 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1474 }
1475 }
1476 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1477 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1478 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1479 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1480 }
1481 }
1482 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1483 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1484 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1485 p[k][j][i] - p[k][j+1][i]) * 0.5;
1486 }
1487 }
1488 else {
1489 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1490 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1491 }
1492
1493 dpde = p[k][j+1][i] - p[k][j][i];
1494
1495 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1496 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1497 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1498 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1499 }
1500 }
1501 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1502 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1503 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1504 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1505 }
1506 }
1507 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1508 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1509 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1510 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1511 }
1512 }
1513 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1514 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1515 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1516 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1517 }
1518 }
1519 else {
1520 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1521 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1522 }
1523
1524 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1525
1526
1527 rhs[k][j][i].y -=
1528 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1529 jcsi[k][j][i].y * jeta[k][j][i].y +
1530 jcsi[k][j][i].z * jeta[k][j][i].z) +
1531 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1532 jeta[k][j][i].y * jeta[k][j][i].y +
1533 jeta[k][j][i].z * jeta[k][j][i].z) +
1534 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1535 jzet[k][j][i].y * jeta[k][j][i].y +
1536 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1537
1538 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1539 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1540 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1541 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1542 }
1543 }
1544 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1545 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1546 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1547 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1548 }
1549 }
1550 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1551 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1552 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1553 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1554 }
1555 }
1556 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1557 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1558 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1559 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1560 }
1561 }
1562 else {
1563 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1564 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1565 }
1566
1567 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1568 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1569 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1570 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1571 }
1572 }
1573 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1574 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1575 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1576 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1577 }
1578 }
1579 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1580 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1581 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1582 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1583 }
1584 }
1585 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1586 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1587 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1588 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1589 }
1590 }
1591 else {
1592 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1593 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1594 }
1595
1596 dpdz = (p[k+1][j][i] - p[k][j][i]);
1597
1598 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1599
1600 rhs[k][j][i].z -=
1601 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1602 kcsi[k][j][i].y * kzet[k][j][i].y +
1603 kcsi[k][j][i].z * kzet[k][j][i].z) +
1604 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1605 keta[k][j][i].y * kzet[k][j][i].y +
1606 keta[k][j][i].z * kzet[k][j][i].z) +
1607 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1608 kzet[k][j][i].y * kzet[k][j][i].y +
1609 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1610
1611 }
1612 }
1613 }
1614
1615
1616 //Mohsen March 2012//
1617
1618 // rhs.x at boundaries for periodic bc at i direction//
1619 if (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && xs==0){
1620 for (k=lzs; k<lze; k++) {
1621 for (j=lys; j<lye; j++) {
1622 i=xs;
1623
1624 dpdc = p[k][j][i+1] - p[k][j][i];
1625
1626 if ((j==my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1627 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1628 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1629 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1630 }
1631 }
1632 else if ((j==my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1633 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1634 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1635 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1636 }
1637 }
1638 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1639 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1640 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1641 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1642 }
1643 }
1644 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1645 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1646 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1647 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1648 }
1649 }
1650 else {
1651 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1652 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1653 }
1654
1655 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1656 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1657 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1658 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1659 }
1660 }
1661 else if ((k == mz-2 || k==1) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1662 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1663 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1664 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1665 }
1666 }
1667 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1668 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1669 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1670 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1671 }
1672 }
1673 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1674 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1675 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1676 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1677 }
1678 }
1679 else {
1680 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1681 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1682 }
1683
1684 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1685 rhs[k][j][i].x -=
1686 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1687 icsi[k][j][i].y * icsi[k][j][i].y +
1688 icsi[k][j][i].z * icsi[k][j][i].z)+
1689 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1690 ieta[k][j][i].y * icsi[k][j][i].y +
1691 ieta[k][j][i].z * icsi[k][j][i].z)+
1692 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1693 izet[k][j][i].y * icsi[k][j][i].y +
1694 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1695 }
1696 }
1697 }
1698
1699// rhs.y at boundaries for periodic bc at j direction//
1700 if (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && ys==0){
1701 for (k=lzs; k<lze; k++) {
1702 for (i=lxs; i<lxe; i++) {
1703
1704 j=ys;
1705
1706 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1707 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1708 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1709 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1710 }
1711 }
1712 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1713 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1714 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1715 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1716 }
1717 }
1718 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1719 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1720 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1721 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1722 }
1723 }
1724 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1725 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1726 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1727 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1728 }
1729 }
1730 else {
1731 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1732 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1733 }
1734
1735 dpde = p[k][j+1][i] - p[k][j][i];
1736
1737 if ((k == mz-2 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1738 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC))) {
1739 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1740 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1741 }
1742 }
1743 else if ((k == mz-2 || k==1 ) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1744 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1745 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1746 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1747 }
1748 }
1749 else if ((k == 1 && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type != PERIODIC)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1750 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1751 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1752 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1753 }
1754 }
1755 else if ((k == 1 || k==mz-2) && user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1756 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1757 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1758 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1759 }
1760 }
1761 else {
1762 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1763 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1764 }
1765
1766 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1767
1768 rhs[k][j][i].y -=
1769 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1770 jcsi[k][j][i].y * jeta[k][j][i].y +
1771 jcsi[k][j][i].z * jeta[k][j][i].z)+
1772 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1773 jeta[k][j][i].y * jeta[k][j][i].y +
1774 jeta[k][j][i].z * jeta[k][j][i].z)+
1775 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1776 jzet[k][j][i].y * jeta[k][j][i].y +
1777 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1778
1779 }
1780 }
1781 }
1782
1783 // rhs.z at boundaries for periodic bc at k direction//
1784 if (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC&& zs==0){
1785 for (j=lys; j<lye; j++) {
1786 for (i=lxs; i<lxe; i++) {
1787
1788 k=zs;
1789
1790 if ((i == mx-2 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1791 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC))) {
1792 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1793 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1794 }
1795 }
1796 else if ((i == mx-2 || i==1) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1797 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1798 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1799 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1800 }
1801 }
1802 else if ((i == 1 && user->boundary_faces[BC_FACE_NEG_X].mathematical_type != PERIODIC)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1803 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1804 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1805 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1806 }
1807 }
1808 else if ((i == 1 || i==mx-2) && user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1809 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1810 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1811 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1812 }
1813 }
1814 else {
1815 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1816 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1817 }
1818
1819 if ((j == my-2 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1820 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC))) {
1821 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1822 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1823 }
1824 }
1825 else if ((j == my-2 || j==1) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1826 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1827 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1828 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1829 }
1830 }
1831 else if ((j == 1 && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type != PERIODIC)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1832 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1833 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1834 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1835 }
1836 }
1837 else if ((j == 1 || j==my-2) && user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1838 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1839 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1840 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1841 }
1842 }
1843 else {
1844 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1845 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1846 }
1847
1848 dpdz = (p[k+1][j][i] - p[k][j][i]);
1849
1850 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1851
1852 rhs[k][j][i].z -=
1853 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1854 kcsi[k][j][i].y * kzet[k][j][i].y +
1855 kcsi[k][j][i].z * kzet[k][j][i].z)+
1856 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1857 keta[k][j][i].y * kzet[k][j][i].y +
1858 keta[k][j][i].z * kzet[k][j][i].z)+
1859 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1860 kzet[k][j][i].y * kzet[k][j][i].y +
1861 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1862
1863 }
1864 }
1865 }
1866
1867 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1868
1869 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure Gradient added to RHS .\n");
1870 PetscInt TwoD = simCtx->TwoD;
1871
1872 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases initiated .\n");
1873
1874 // 7. Final clean-up for immersed boundaries and 2D cases
1875 for (k=lzs; k<lze; k++) {
1876 for (j=lys; j<lye; j++) {
1877 for (i=lxs; i<lxe; i++) {
1878 if (TwoD==1)
1879 rhs[k][j][i].x =0.;
1880 else if (TwoD==2)
1881 rhs[k][j][i].y =0.;
1882 else if (TwoD==3)
1883 rhs[k][j][i].z =0.;
1884
1885 if (nvert[k][j][i]>0.1) {
1886 rhs[k][j][i].x = 0;
1887 rhs[k][j][i].y = 0;
1888 rhs[k][j][i].z = 0;
1889 }
1890 if (nvert[k][j][i+1]>0.1) {
1891 rhs[k][j][i].x=0;
1892 }
1893 if (nvert[k][j+1][i]>0.1) {
1894 rhs[k][j][i].y=0;
1895 }
1896 if (nvert[k+1][j][i]>0.1) {
1897 rhs[k][j][i].z=0;
1898 }
1899 }
1900 }
1901 }
1902 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases complete .\n");
1903
1904 // ========================================================================
1905
1906 // --- Restore all PETSc array pointers ---
1907 // DMDAVecRestoreArray(fda, user->lUcont, &ucont);
1908 DMDAVecRestoreArray(fda, Rhs, &rhs);
1909 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Rhs restored successfully! .\n");
1910
1911 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1912 DMDAVecRestoreArray(fda, user->lEta, &eta);
1913 DMDAVecRestoreArray(fda, user->lZet, &zet);
1914 DMDAVecRestoreArray(da, user->lAj, &aj);
1915 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Face metrics restored successfully! .\n");
1916
1917 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1918 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1919 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1920 DMDAVecRestoreArray(da, user->lIAj, &iaj);
1921 LOG_ALLOW(GLOBAL,LOG_DEBUG,"I Face metrics restored successfully! .\n");
1922
1923 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1924 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1925 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1926 DMDAVecRestoreArray(da, user->lJAj, &jaj);
1927 LOG_ALLOW(GLOBAL,LOG_DEBUG,"J Face metrics restored successfully! .\n");
1928
1929 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1930 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1931 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1932 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1933 LOG_ALLOW(GLOBAL,LOG_DEBUG,"K Face metrics restored successfully! .\n");
1934
1935 DMDAVecRestoreArray(da, user->lP, &p);
1936 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure restored successfully! .\n");
1937
1938 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1939 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Nvert restored successfully! .\n");
1940
1941 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Cell Centered scalars restored successfully! .\n");
1942
1943 // --- Destroy temporary work vectors ---
1944 ierr = VecDestroy(&Conv); CHKERRQ(ierr);
1945 ierr = VecDestroy(&Visc); CHKERRQ(ierr);
1946 ierr = VecDestroy(&Rc); CHKERRQ(ierr);
1947 ierr = VecDestroy(&Rct); CHKERRQ(ierr);
1948 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Temporary work vectors destroyed successfully! .\n");
1949
1950 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: RHS computation complete.\n",
1951 simCtx->rank, user->_this);
1952
1953 PROFILE_FUNCTION_END;
1954 PetscFunctionReturn(0);
1955}
1956
1957#undef __FUNCT__
1958#define __FUNCT__ "ComputeEulerianDiffusivity"
1959/**
1960 * @brief Computes the effective diffusivity scalar field ($\Gamma_{eff}$) on the Eulerian grid.
1961 *
1962 * This function calculates the total diffusivity used to drive the stochastic
1963 * motion of particles (Scalar FDF). It combines molecular diffusion and
1964 * turbulent diffusion.
1965 *
1966 * **Formula:**
1967 * \f[
1968 * \Gamma_{eff} = \underbrace{\frac{\nu}{Sc}}_{\text{Molecular}} + \underbrace{\frac{\nu_t}{Sc_t}}_{\text{Turbulent}}
1969 * \f]
1970 *
1971 * Where:
1972 * - \f$ \nu = 1/Re \f$ (Kinematic Viscosity)
1973 * - \f$ \nu_t \f$ (Eddy Viscosity from LES/RANS model)
1974 * - \f$ Sc \f$ (Molecular Schmidt Number, user-defined)
1975 * - \f$ Sc_t \f$ (Turbulent Schmidt Number, user-defined)
1976 *
1977 * @note If turbulence models are disabled, \f$ \nu_t \f$ is assumed to be 0.
1978 * @note This function updates the local ghost values of lDiffusivity at the end
1979 * to ensure gradients can be computed correctly at subdomain boundaries.
1980 *
1981 * @param[in,out] user Pointer to the user context containing grid data and simulation parameters.
1982 *
1983 * @return PetscErrorCode 0 on success.
1984 */
1985PetscErrorCode ComputeEulerianDiffusivity(UserCtx *user)
1986{
1987 PetscErrorCode ierr;
1988 DM da = user->da;
1989 PetscInt i, j, k, xs, ys, zs, xm, ym, zm, xe, ye, ze;
1990 PetscInt lxs, lys, lzs, lxe, lye, lze;
1991
1992 // Pointers for 3D grid access
1993 PetscReal ***diff_arr; // Output: Diffusivity field
1994 PetscReal ***nut_arr; // Input: Eddy Viscosity field (optional)
1995
1996 // Physics parameters
1997 PetscReal nu_molecular, gamma_molecular;
1998 PetscReal nu_turbulent, gamma_turbulent;
1999 PetscReal Sc, Sct;
2000 PetscBool use_turbulence_model;
2001
2002 PetscFunctionBeginUser;
2003 PROFILE_FUNCTION_BEGIN;
2004
2005 // ------------------------------------------------------------------------
2006 // 1. Parameter Setup & Safety Checks
2007 // ------------------------------------------------------------------------
2008
2009 // Determine Molecular Viscosity (nu = 1/Re)
2010 // Guard against division by zero if Re is not set or infinite (inviscid)
2011 if (user->simCtx->ren > 1.0e-12) {
2012 nu_molecular = 1.0 / user->simCtx->ren;
2013 } else {
2014 nu_molecular = 0.0;
2015 }
2016
2017 // Set Schmidt Numbers (Default to 1.0 if not provided to prevent NaN)
2018 Sc = (user->simCtx->schmidt_number > 1.0e-6) ? user->simCtx->schmidt_number : 1.0;
2019 Sct = (user->simCtx->Turbulent_schmidt_number > 1.0e-6) ? user->simCtx->Turbulent_schmidt_number : 0.7;
2020
2021 // Pre-calculate molecular component
2022 gamma_molecular = nu_molecular / Sc;
2023
2024 // Check if a turbulence model is active (LES or RANS)
2025 use_turbulence_model = (user->simCtx->les || user->simCtx->rans) ? PETSC_TRUE : PETSC_FALSE;
2026
2027 // ------------------------------------------------------------------------
2028 // 2. Data Access
2029 // ------------------------------------------------------------------------
2030
2031 // Get local grid boundaries
2032 DMDALocalInfo info;
2033 ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr);
2034
2035 xs = info.xs; ys = info.ys; zs = info.zs;
2036 xm = info.xm; ym = info.ym; zm = info.zm;
2037 xe = xs + xm; ye = ys + ym; ze = zs + zm;
2038
2039 lxs = (xs == 0)? xs + 1 : xs;
2040 lys = (ys == 0)? ys + 1 : ys;
2041 lzs = (zs == 0)? zs + 1 : zs;
2042 lxe = (xe == info.mx)? xe - 1 : xe;
2043 lye = (ye == info.my)? ye - 1 : ye;
2044 lze = (ze == info.mz)? ze - 1 : ze;
2045
2046 // Get write access to the output Diffusivity array
2047 ierr = DMDAVecGetArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
2048
2049 // Get read access to Eddy Viscosity only if turbulence is active
2050 if (use_turbulence_model) {
2051 ierr = DMDAVecGetArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
2052 }
2053
2054 // ------------------------------------------------------------------------
2055 // 3. Calculation Loop
2056 // ------------------------------------------------------------------------
2057
2058 for (k = lzs; k < lze; k++) {
2059 for (j = lys; j < lye; j++) {
2060 for (i = lxs; i < lxe; i++) {
2061
2062 gamma_turbulent = 0.0;
2063
2064 if (use_turbulence_model) {
2065 // Fetch local eddy viscosity
2066 nu_turbulent = nut_arr[k][j][i];
2067
2068 // NUMERICAL SAFETY:
2069 // Some turbulence models (dynamic SGS) can locally produce
2070 // slightly negative viscosity. We clamp this to 0 to prevent
2071 // negative diffusivity, which crashes the Langevin sqrt().
2072 if (nu_turbulent < 0.0) {
2073 nu_turbulent = 0.0;
2074 }
2075
2076 gamma_turbulent = nu_turbulent / Sct;
2077 }
2078
2079 // Sum components
2080 diff_arr[k][j][i] = gamma_molecular + gamma_turbulent;
2081 }
2082 }
2083 }
2084
2085 // ------------------------------------------------------------------------
2086 // 4. Cleanup & Synchronization
2087 // ------------------------------------------------------------------------
2088
2089 // Restore arrays
2090 if (use_turbulence_model) {
2091 ierr = DMDAVecRestoreArrayRead(da, user->Nu_t, &nut_arr); CHKERRQ(ierr);
2092 }
2093 ierr = DMDAVecRestoreArray(da, user->Diffusivity, &diff_arr); CHKERRQ(ierr);
2094
2095 // Update Ghost Points
2096 // This is required because downstream operations (Drift Gradient Calculation
2097 // and Particle Interpolation) will need access to the halo regions of this field.
2098 ierr = UpdateLocalGhosts(user,"Diffusivity"); CHKERRQ(ierr);
2099 PROFILE_FUNCTION_END;
2100 PetscFunctionReturn(0);
2101}
2102
2103#undef __FUNCT__
2104#define __FUNCT__ "ComputeEulerianDiffusivityGradient"
2105/**
2106 * @brief Computes the gradient of the scalar Diffusivity field (Drift Vector).
2107 *
2108 * Logic:
2109 * 1. Iterates over the PHYSICAL cell centers [1 ... M-2].
2110 * 2. Uses 2nd-Order Central Difference for interior cells.
2111 * 3. Uses 2nd-Order One-Sided Differences at physical boundaries
2112 * (Forward at start, Backward at end) unless the boundary is Periodic.
2113 * 4. Transforms results to Cartesian coordinates using the Chain Rule.
2114 */
2115PetscErrorCode ComputeEulerianDiffusivityGradient(UserCtx *user)
2116{
2117 PetscErrorCode ierr;
2118 DM da = user->da, fda = user->fda;
2119 DMDALocalInfo info = user->info;
2120
2121 // 1. Determine Global Dimensions
2122 PetscInt mx = info.mx, my = info.my, mz = info.mz;
2123
2124 // 2. Determine Local Loop Bounds (Skip Ghosts/Unused Indices)
2125 // Grid uses indices 1 to M-2 for physical cells.
2126 // Index 0 and M-1 are ghost/boundary holders.
2127
2128 // Start: If we own the global start (0), skip it and start at 1.
2129 PetscInt lxs = (info.xs == 0) ? 1 : info.xs;
2130 // End: If we own the global end (mx), stop before it (mx-1), so loop covers mx-2.
2131 PetscInt lxe = (info.xs + info.xm == mx) ? mx - 1 : info.xs + info.xm;
2132
2133 PetscInt lys = (info.ys == 0) ? 1 : info.ys;
2134 PetscInt lye = (info.ys + info.ym == my) ? my - 1 : info.ys + info.ym;
2135
2136 PetscInt lzs = (info.zs == 0) ? 1 : info.zs;
2137 PetscInt lze = (info.zs + info.zm == mz) ? mz - 1 : info.zs + info.zm;
2138
2139 PetscInt i, j, k;
2140
2141 // Pointers
2142 PetscReal ***diff; // Input (Scalar)
2143 Cmpnts ***grad_diff; // Output (Vector)
2144 Cmpnts ***csi, ***eta, ***zet;
2145 PetscReal ***aj;
2146
2147 // Boundary Flags (Check if Periodic)
2148 PetscBool p_x = (user->boundary_faces[BC_FACE_NEG_X].mathematical_type == PERIODIC);
2149 PetscBool p_y = (user->boundary_faces[BC_FACE_NEG_Y].mathematical_type == PERIODIC);
2150 PetscBool p_z = (user->boundary_faces[BC_FACE_NEG_Z].mathematical_type == PERIODIC);
2151
2152 PetscFunctionBeginUser;
2153 PROFILE_FUNCTION_BEGIN;
2154
2155 // 3. Update Ghosts for Diffusivity
2156 // Required so that Central Differences at i=2 can read i=1,
2157 // and Central Differences at Periodic Boundaries work correctly.
2158 ierr = UpdateLocalGhosts(user, "Diffusivity"); CHKERRQ(ierr);
2159
2160 // 4. Get Arrays (Read Only)
2161 ierr = DMDAVecGetArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2162 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2163 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2164 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2165 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2166
2167 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity and Metrics arrays accessed successfully! .\n");
2168
2169 // 5. Get Output Array (Read/Write)
2170 ierr = DMDAVecGetArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2171
2172 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Diffusivity Gradient array accessed successfully! .\n");
2173
2174 // 6. Loop over Physical Domain
2175 for (k = lzs; k < lze; k++) {
2176 for (j = lys; j < lye; j++) {
2177 for (i = lxs; i < lxe; i++) {
2178
2179 PetscReal dGdCsi, dGdEta, dGdZet;
2180
2181 // ---------------------------------------------------------
2182 // I-Direction (Csi)
2183 // ---------------------------------------------------------
2184 if (!p_x && i == 1) {
2185 // Physical Start: 2nd Order Forward Difference
2186 // Stencil: [-3, 4, -1] / 2
2187 dGdCsi = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j][i+1] - diff[k][j][i+2]) * 0.5;
2188 }
2189 else if (!p_x && i == mx - 2) {
2190 // Physical End: 2nd Order Backward Difference
2191 // Stencil: [1, -4, 3] / 2
2192 dGdCsi = (3.0 * diff[k][j][i] - 4.0 * diff[k][j][i-1] + diff[k][j][i-2]) * 0.5;
2193 }
2194 else {
2195 // Interior / Periodic: Central Difference
2196 // Stencil: [-1, 0, 1] / 2
2197 dGdCsi = (diff[k][j][i+1] - diff[k][j][i-1]) * 0.5;
2198 }
2199
2200 // ---------------------------------------------------------
2201 // J-Direction (Eta)
2202 // ---------------------------------------------------------
2203 if (!p_y && j == 1) {
2204 // Physical Start: Forward
2205 dGdEta = (-3.0 * diff[k][j][i] + 4.0 * diff[k][j+1][i] - diff[k][j+2][i]) * 0.5;
2206 }
2207 else if (!p_y && j == my - 2) {
2208 // Physical End: Backward
2209 dGdEta = (3.0 * diff[k][j][i] - 4.0 * diff[k][j-1][i] + diff[k][j-2][i]) * 0.5;
2210 }
2211 else {
2212 // Interior: Central
2213 dGdEta = (diff[k][j+1][i] - diff[k][j-1][i]) * 0.5;
2214 }
2215
2216 // ---------------------------------------------------------
2217 // K-Direction (Zet)
2218 // ---------------------------------------------------------
2219 if (!p_z && k == 1) {
2220 // Physical Start: Forward
2221 dGdZet = (-3.0 * diff[k][j][i] + 4.0 * diff[k+1][j][i] - diff[k+2][j][i]) * 0.5;
2222 }
2223 else if (!p_z && k == mz - 2) {
2224 // Physical End: Backward
2225 dGdZet = (3.0 * diff[k][j][i] - 4.0 * diff[k-1][j][i] + diff[k-2][j][i]) * 0.5;
2226 }
2227 else {
2228 // Interior: Central
2229 dGdZet = (diff[k+1][j][i] - diff[k-1][j][i]) * 0.5;
2230 }
2231
2232 // ---------------------------------------------------------
2233 // Transform to Physical Space (Cartesian Gradient)
2234 // ---------------------------------------------------------
2235 TransformScalarDerivativesToPhysical(aj[k][j][i],
2236 csi[k][j][i], eta[k][j][i], zet[k][j][i],
2237 dGdCsi, dGdEta, dGdZet,
2238 &grad_diff[k][j][i]);
2239 }
2240 }
2241 }
2242
2243 // 7. Restore Arrays
2244 ierr = DMDAVecRestoreArrayRead(da, user->lDiffusivity, &diff); CHKERRQ(ierr);
2245 ierr = DMDAVecRestoreArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
2246 ierr = DMDAVecRestoreArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
2247 ierr = DMDAVecRestoreArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
2248 ierr = DMDAVecRestoreArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
2249 ierr = DMDAVecRestoreArray(fda, user->DiffusivityGradient, &grad_diff); CHKERRQ(ierr);
2250
2251 // 8. Update Ghosts for the Result
2252 // Important: Particles near subdomain boundaries will need to interpolate
2253 // this gradient vector, so the ghosts of the vector field must be filled.
2254 ierr = UpdateLocalGhosts(user, "DiffusivityGradient"); CHKERRQ(ierr);
2255
2256 PROFILE_FUNCTION_END;
2257 PetscFunctionReturn(0);
2258}
ComputeDrivenChannelFlowSource
PetscErrorCode ComputeDrivenChannelFlowSource(UserCtx *user, Vec Rct)
Applies a momentum source term to drive flow in a periodic channel or pipe.
Definition BodyForces.c:29
ComputeCellCharacteristicLengthScale
PetscErrorCode ComputeCellCharacteristicLengthScale(PetscReal ajc, Cmpnts csi, Cmpnts eta, Cmpnts zet, double *dx, double *dy, double *dz)
Computes characteristic length scales (dx, dy, dz) for a curvilinear cell.
Definition Metric.c:313
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
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:200
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:746
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:737
ComputeBodyForces
PetscErrorCode ComputeBodyForces(UserCtx *user, Vec Rct)
General dispatcher for applying all active body forces (momentum sources).
Definition rhs.c:1157
Viscous
PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
Definition rhs.c:504
ComputeEulerianDiffusivity
PetscErrorCode ComputeEulerianDiffusivity(UserCtx *user)
Computes the effective diffusivity scalar field ($\Gamma_{eff}$) on the Eulerian grid.
Definition rhs.c:1985
ComputeEulerianDiffusivityGradient
PetscErrorCode ComputeEulerianDiffusivityGradient(UserCtx *user)
Computes the gradient of the scalar Diffusivity field (Drift Vector).
Definition rhs.c:2115
Convection
PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
Definition rhs.c:5
ComputeRHS
PetscErrorCode ComputeRHS(UserCtx *user, Vec Rhs)
Computes the Right-Hand-Side (RHS) of the momentum equations.
Definition rhs.c:1192
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2177
TransformScalarDerivativesToPhysical
void TransformScalarDerivativesToPhysical(PetscReal jacobian, Cmpnts csi_metrics, Cmpnts eta_metrics, Cmpnts zet_metrics, PetscReal dPhi_dcsi, PetscReal dPhi_deta, PetscReal dPhi_dzet, Cmpnts *gradPhi)
Transforms scalar derivatives from computational space to physical space using the chain rule.
Definition setup.c:2779
UpdateLocalGhosts
PetscErrorCode UpdateLocalGhosts(UserCtx *user, const char *fieldName)
Updates the local vector (including ghost points) from its corresponding global vector.
Definition setup.c:1157
LESModelType
LESModelType
Identifies the six logical faces of a structured computational block.
Definition variables.h:447
SimCtx::clark
PetscInt clark
Definition variables.h:670
CHARACTERISTIC_BC
@ CHARACTERISTIC_BC
Definition variables.h:220
PERIODIC
@ PERIODIC
Definition variables.h:219
WALL
@ WALL
Definition variables.h:213
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:615
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:615
SimCtx::schmidt_number
PetscReal schmidt_number
Definition variables.h:646
SimCtx::rank
PetscMPIInt rank
Definition variables.h:588
SimCtx::Turbulent_schmidt_number
PetscReal Turbulent_schmidt_number
Definition variables.h:646
UserCtx::boundary_faces
BoundaryFaceConfig boundary_faces[6]
Definition variables.h:746
UserCtx::lIEta
Vec lIEta
Definition variables.h:778
UserCtx::lIZet
Vec lIZet
Definition variables.h:778
UserCtx::lNvert
Vec lNvert
Definition variables.h:755
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:731
SimCtx::rans
PetscInt rans
Definition variables.h:669
UserCtx::lZet
Vec lZet
Definition variables.h:776
SimCtx::ren
PetscReal ren
Definition variables.h:632
UserCtx::lIAj
Vec lIAj
Definition variables.h:778
UserCtx::_this
PetscInt _this
Definition variables.h:741
UserCtx::lKEta
Vec lKEta
Definition variables.h:780
UserCtx::DiffusivityGradient
Vec DiffusivityGradient
Definition variables.h:759
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:779
Cmpnts::x
PetscScalar x
Definition variables.h:101
SimCtx::invicid
PetscInt invicid
Definition variables.h:615
UserCtx::lKZet
Vec lKZet
Definition variables.h:780
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:783
UserCtx::Nu_t
Vec Nu_t
Definition variables.h:783
UserCtx::lJEta
Vec lJEta
Definition variables.h:779
UserCtx::lCsi
Vec lCsi
Definition variables.h:776
Cmpnts::z
PetscScalar z
Definition variables.h:101
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:780
SimCtx::central
PetscInt central
Definition variables.h:630
UserCtx::lJZet
Vec lJZet
Definition variables.h:779
UserCtx::lUcont
Vec lUcont
Definition variables.h:755
SimCtx::step
PetscInt step
Definition variables.h:593
UserCtx::Diffusivity
Vec Diffusivity
Definition variables.h:758
UserCtx::lAj
Vec lAj
Definition variables.h:776
UserCtx::lICsi
Vec lICsi
Definition variables.h:778
UserCtx::info
DMDALocalInfo info
Definition variables.h:735
UserCtx::lUcat
Vec lUcat
Definition variables.h:755
Cmpnts::y
PetscScalar y
Definition variables.h:101
UserCtx::lEta
Vec lEta
Definition variables.h:776
SimCtx::les
PetscInt les
Definition variables.h:669
UserCtx::lDiffusivity
Vec lDiffusivity
Definition variables.h:758
BoundaryFaceConfig::mathematical_type
BCType mathematical_type
Definition variables.h:295
UserCtx::lJAj
Vec lJAj
Definition variables.h:779
SimCtx::rotateframe
PetscInt rotateframe
Definition variables.h:615
UserCtx::lKAj
Vec lKAj
Definition variables.h:780
BC_FACE_NEG_X
@ BC_FACE_NEG_X
Definition variables.h:204
BC_FACE_POS_Z
@ BC_FACE_POS_Z
Definition variables.h:206
BC_FACE_POS_Y
@ BC_FACE_POS_Y
Definition variables.h:205
BC_FACE_NEG_Z
@ BC_FACE_NEG_Z
Definition variables.h:206
BC_FACE_POS_X
@ BC_FACE_POS_X
Definition variables.h:204
BC_FACE_NEG_Y
@ BC_FACE_NEG_Y
Definition variables.h:205
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:585
UserCtx
User-defined context containing data specific to a single computational grid level.
Definition variables.h:728
nu_t
double nu_t(double yplus)
Computes turbulent eddy viscosity ratio (ν_t / ν)
Definition wallfunction.c:423