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

Go to the source code of this file.

Macros

#define __FUNCT__   "Convection"
 
#define __FUNCT__   "Viscous"
 
#define __FUNCT__   "ComputeRHS"
 

Functions

PetscErrorCode Convection (UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
 
PetscErrorCode Viscous (UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
 
PetscErrorCode ComputeRHS (UserCtx *user, Vec Rhs)
 Computes the Right-Hand-Side (RHS) of the momentum equations.
 

Macro Definition Documentation

◆ __FUNCT__ [1/3]

#define __FUNCT__   "Convection"

Definition at line 4 of file rhs.c.

◆ __FUNCT__ [2/3]

#define __FUNCT__   "Viscous"

Definition at line 4 of file rhs.c.

◆ __FUNCT__ [3/3]

#define __FUNCT__   "ComputeRHS"

Definition at line 4 of file rhs.c.

Function Documentation

◆ Convection()

PetscErrorCode Convection ( UserCtx *  user,
Vec  Ucont,
Vec  Ucat,
Vec  Conv 
)

Definition at line 5 of file rhs.c.

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 PetscInt 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->bctype[0]==7 || user->bctype[1]==7){
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->bctype[2]==7 || user->bctype[3]==7){
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->bctype[4]==7 || user->bctype[5]==7){
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->bctype[0]==7 && 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->bctype[0]==7 && 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->bctype[2]==7 && 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->bctype[2]==7 && 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->bctype[4]==7 && 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->bctype[4]==7 && 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}
GLOBAL
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:47
LOG_ALLOW
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:201
PROFILE_FUNCTION_END
#define PROFILE_FUNCTION_END
Marks the end of a profiled code block.
Definition logging.h:740
LOG_DEBUG
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:33
PROFILE_FUNCTION_BEGIN
#define PROFILE_FUNCTION_BEGIN
Marks the beginning of a profiled code block (typically a function).
Definition logging.h:731
UserCtx::lNvert
Vec lNvert
Definition variables.h:688
UserCtx::simCtx
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:664
Cmpnts::x
PetscScalar x
Definition variables.h:101
Cmpnts::z
PetscScalar z
Definition variables.h:101
SimCtx::central
PetscInt central
Definition variables.h:580
UserCtx::bctype
PetscInt bctype[6]
Definition variables.h:684
UserCtx::lAj
Vec lAj
Definition variables.h:705
Cmpnts::y
PetscScalar y
Definition variables.h:101
SimCtx::les
PetscInt les
Definition variables.h:609
Cmpnts
A 3D point or vector with PetscScalar components.
Definition variables.h:100
SimCtx
The master context for the entire simulation.
Definition variables.h:538
Here is the caller graph for this function:

◆ Viscous()

PetscErrorCode Viscous ( UserCtx *  user,
Vec  Ucont,
Vec  Ucat,
Vec  Visc 
)

Definition at line 504 of file rhs.c.

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 PetscInt 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->bctype[0]==1 && i==0) || (user->bctype[1]==1 && 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->bctype[2]==1 && j==0) || (user->bctype[3]==1 && 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->bctype[4]==1 && k==0) || (user->bctype[5]==1 && 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}
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
SimCtx::clark
PetscInt clark
Definition variables.h:610
UserCtx::lIEta
Vec lIEta
Definition variables.h:707
UserCtx::lIZet
Vec lIZet
Definition variables.h:707
SimCtx::rans
PetscInt rans
Definition variables.h:609
UserCtx::lZet
Vec lZet
Definition variables.h:705
SimCtx::ren
PetscReal ren
Definition variables.h:581
UserCtx::lIAj
Vec lIAj
Definition variables.h:707
UserCtx::lKEta
Vec lKEta
Definition variables.h:709
UserCtx::lJCsi
Vec lJCsi
Definition variables.h:708
UserCtx::lKZet
Vec lKZet
Definition variables.h:709
UserCtx::lNu_t
Vec lNu_t
Definition variables.h:712
UserCtx::lJEta
Vec lJEta
Definition variables.h:708
UserCtx::lCsi
Vec lCsi
Definition variables.h:705
UserCtx::lKCsi
Vec lKCsi
Definition variables.h:709
UserCtx::lJZet
Vec lJZet
Definition variables.h:708
SimCtx::step
PetscInt step
Definition variables.h:546
UserCtx::lICsi
Vec lICsi
Definition variables.h:707
UserCtx::lEta
Vec lEta
Definition variables.h:705
UserCtx::lJAj
Vec lJAj
Definition variables.h:708
UserCtx::lKAj
Vec lKAj
Definition variables.h:709
nu_t
double nu_t(double yplus)
Definition wallfunction.c:103
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ComputeRHS()

PetscErrorCode ComputeRHS ( UserCtx *  user,
Vec  Rhs 
)

Computes the Right-Hand-Side (RHS) of the momentum equations.

Computes the Right-Hand Side (RHS) of the momentum equations.

This function orchestrates the calculation of the spatial discretization of the convective and diffusive terms. It calls specialized helper functions (Convection, Viscous) to compute these terms and then combines them with the pressure gradient to form the final RHS vector.

This is a minimally-edited version of the legacy kernel. It operates on a single UserCtx (one grid block) and retrieves all global parameters (Re, rans, les, etc.) from the master SimCtx.

Parameters
userThe UserCtx for the specific block being computed.
RhsThe PETSc Vec where the calculated RHS will be stored.
Returns
PetscErrorCode 0 on success.

Definition at line 1159 of file rhs.c.

1160{
1161 PetscErrorCode ierr;
1162 SimCtx *simCtx = user->simCtx;
1163 DM da = user->da, fda = user->fda;
1164 DMDALocalInfo info = user->info;
1165 PetscInt i,j,k;
1166 PetscReal dpdc, dpde,dpdz;
1167 // --- Local Grid Indices and Parameters ---
1168 PetscInt xs = info.xs, xe = xs + info.xm, mx = info.mx;
1169 PetscInt ys = info.ys, ye = ys + info.ym, my = info.my;
1170 PetscInt zs = info.zs, ze = zs + info.zm, mz = info.mz;
1171 PetscInt lxs = (xs==0) ? xs+1 : xs;
1172 PetscInt lys = (ys==0) ? ys+1 : ys;
1173 PetscInt lzs = (zs==0) ? zs+1 : zs;
1174 PetscInt lxe = (xe==mx) ? xe-1 : xe;
1175 PetscInt lye = (ye==my) ? ye-1 : ye;
1176 PetscInt lze = (ze==mz) ? ze-1 : ze;
1177 PetscInt mz_end = (user->bctype[5] == 8) ? mz - 2 : mz - 3;
1178
1179 // --- Array Pointers ---
1180 Cmpnts ***csi, ***eta, ***zet, ***icsi, ***ieta, ***izet, ***jcsi, ***jeta, ***jzet, ***kcsi, ***keta, ***kzet;
1181 PetscReal ***p, ***iaj, ***jaj, ***kaj, ***aj, ***nvert, ***nvert_o;
1182 Cmpnts ***rhs, ***rc, ***rct,***ucont;
1183
1184 // --- Temporary Vectors ---
1185 Vec Conv, Visc, Rc, Rct;
1186
1187 PetscFunctionBeginUser;
1188 PROFILE_FUNCTION_BEGIN;
1189 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: Computing RHS (FormFunction1)...\n",
1190 simCtx->rank, user->_this);
1191
1192 // --- Get all necessary array pointers ---
1193 ierr = DMDAVecGetArrayRead(fda, user->lCsi, &csi); CHKERRQ(ierr);
1194 ierr = DMDAVecGetArrayRead(fda, user->lEta, &eta); CHKERRQ(ierr);
1195 ierr = DMDAVecGetArrayRead(fda, user->lZet, &zet); CHKERRQ(ierr);
1196 ierr = DMDAVecGetArrayRead(da, user->lAj, &aj); CHKERRQ(ierr);
1197 ierr = DMDAVecGetArrayRead(fda, user->lICsi, &icsi); CHKERRQ(ierr);
1198 ierr = DMDAVecGetArrayRead(fda, user->lIEta, &ieta); CHKERRQ(ierr);
1199 ierr = DMDAVecGetArrayRead(fda, user->lIZet, &izet); CHKERRQ(ierr);
1200 ierr = DMDAVecGetArrayRead(fda, user->lJCsi, &jcsi); CHKERRQ(ierr);
1201 ierr = DMDAVecGetArrayRead(fda, user->lJEta, &jeta); CHKERRQ(ierr);
1202 ierr = DMDAVecGetArrayRead(fda, user->lJZet, &jzet); CHKERRQ(ierr);
1203 ierr = DMDAVecGetArrayRead(fda, user->lKCsi, &kcsi); CHKERRQ(ierr);
1204 ierr = DMDAVecGetArrayRead(fda, user->lKEta, &keta); CHKERRQ(ierr);
1205 ierr = DMDAVecGetArrayRead(fda, user->lKZet, &kzet); CHKERRQ(ierr);
1206 ierr = DMDAVecGetArrayRead(da, user->lIAj, &iaj); CHKERRQ(ierr);
1207 ierr = DMDAVecGetArrayRead(da, user->lJAj, &jaj); CHKERRQ(ierr);
1208 ierr = DMDAVecGetArrayRead(da, user->lKAj, &kaj); CHKERRQ(ierr);
1209 ierr = DMDAVecGetArrayRead(da, user->lP, &p); CHKERRQ(ierr);
1210 ierr = DMDAVecGetArrayRead(da, user->lNvert, &nvert); CHKERRQ(ierr);
1211 ierr = DMDAVecGetArray(fda, Rhs, &rhs); CHKERRQ(ierr);
1212
1213 // --- Create temporary work vectors ---
1214 ierr = VecDuplicate(user->lUcont, &Rc); CHKERRQ(ierr);
1215 ierr = VecDuplicate(Rc, &Rct); CHKERRQ(ierr);
1216 ierr = VecDuplicate(Rct, &Conv); CHKERRQ(ierr);
1217 ierr = VecDuplicate(Rct, &Visc); CHKERRQ(ierr);
1218
1219 // ========================================================================
1220 // CORE LOGIC (UNCHANGED FROM LEGACY CODE)
1221 // ========================================================================
1222
1223 // 1. Obtain Cartesian velocity from Contravariant velocity
1224 ierr = Contra2Cart(user); CHKERRQ(ierr);
1225
1226 // 2. Compute Convective term
1227 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating convective terms...\n");
1228 if (simCtx->moveframe || simCtx->rotateframe) {
1229 // ierr = Convection_MV(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1230 } else {
1231 ierr = Convection(user, user->lUcont, user->lUcat, Conv); CHKERRQ(ierr);
1232 }
1233
1234 // 3. Compute Viscous term
1235 if (simCtx->invicid) {
1236 ierr = VecSet(Visc, 0.0); CHKERRQ(ierr);
1237 } else {
1238 LOG_ALLOW(LOCAL, LOG_DEBUG, " Calculating viscous terms...\n");
1239 ierr = Viscous(user, user->lUcont, user->lUcat, Visc); CHKERRQ(ierr);
1240 }
1241
1242 // 4. Combine terms to get Cartesian RHS: Rc = Visc - Conv
1243 ierr = VecWAXPY(Rc, -1.0, Conv, Visc); CHKERRQ(ierr);
1244
1245 // 5. Convert Cartesian RHS (Rc) to Contravariant RHS (Rct)
1246 LOG_ALLOW(LOCAL, LOG_DEBUG, " Converting Cartesian RHS to Contravariant RHS...\n");
1247 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1248 ierr = DMDAVecGetArray(fda, Rc, &rc); CHKERRQ(ierr);
1249
1250 for (k = lzs; k < lze; k++) {
1251 for (j = lys; j < lye; j++) {
1252 for (i = lxs; i < lxe; i++) {
1253 rct[k][j][i].x = aj[k][j][i] *
1254 (0.5 * (csi[k][j][i].x + csi[k][j][i-1].x) * rc[k][j][i].x +
1255 0.5 * (csi[k][j][i].y + csi[k][j][i-1].y) * rc[k][j][i].y +
1256 0.5 * (csi[k][j][i].z + csi[k][j][i-1].z) * rc[k][j][i].z);
1257 rct[k][j][i].y = aj[k][j][i] *
1258 (0.5 * (eta[k][j][i].x + eta[k][j-1][i].x) * rc[k][j][i].x +
1259 0.5 * (eta[k][j][i].y + eta[k][j-1][i].y) * rc[k][j][i].y +
1260 0.5 * (eta[k][j][i].z + eta[k][j-1][i].z) * rc[k][j][i].z);
1261 rct[k][j][i].z = aj[k][j][i] *
1262 (0.5 * (zet[k][j][i].x + zet[k-1][j][i].x) * rc[k][j][i].x +
1263 0.5 * (zet[k][j][i].y + zet[k-1][j][i].y) * rc[k][j][i].y +
1264 0.5 * (zet[k][j][i].z + zet[k-1][j][i].z) * rc[k][j][i].z);
1265 }
1266 }
1267 }
1268 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1269 ierr = DMDAVecRestoreArray(fda, Rc, &rc); CHKERRQ(ierr);
1270
1271 PetscBarrier(NULL);
1272
1273 DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct);
1274 DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct);
1275
1276 DMDAVecGetArray(fda, Rct, &rct);
1277
1278 if (user->bctype[0]==7 || user->bctype[1]==7){
1279 if (xs==0){
1280 i=xs;
1281 for (k=lzs; k<lze; k++) {
1282 for (j=lys; j<lye; j++) {
1283 rct[k][j][i].x=rct[k][j][i-2].x;
1284 }
1285 }
1286 }
1287
1288 if (xe==mx){
1289 i=mx-1;
1290 for (k=lzs; k<lze; k++) {
1291 for (j=lys; j<lye; j++) {
1292 rct[k][j][i].x=rct[k][j][i+2].x;
1293 }
1294 }
1295 }
1296 }
1297
1298 if (user->bctype[2]==7 || user->bctype[3]==7){
1299 if (ys==0){
1300 j=ys;
1301 for (k=lzs; k<lze; k++) {
1302 for (i=lxs; i<lxe; i++) {
1303 rct[k][j][i].y=rct[k][j-2][i].y;
1304 }
1305 }
1306 }
1307
1308 if (ye==my){
1309 j=my-1;
1310 for (k=lzs; k<lze; k++) {
1311 for (i=lxs; i<lxe; i++) {
1312 rct[k][j][i].y=rct[k][j+2][i].y;
1313 }
1314 }
1315 }
1316 }
1317 if (user->bctype[4]==7 || user->bctype[5]==7){
1318 if (zs==0){
1319 k=zs;
1320 for (j=lys; j<lye; j++) {
1321 for (i=lxs; i<lxe; i++) {
1322 rct[k][j][i].z=rct[k-2][j][i].z;
1323 }
1324 }
1325 }
1326 if (ze==mz){
1327 k=mz-1;
1328 for (j=lys; j<lye; j++) {
1329 for (i=lxs; i<lxe; i++) {
1330 rct[k][j][i].z=rct[k+2][j][i].z;
1331 }
1332 }
1333 }
1334 }
1335
1336 // 6. Add Pressure Gradient Term and Finalize RHS
1337 // This involves calculating pressure derivatives (dpdc, dpde, dpdz) and using
1338 // them to adjust the contravariant RHS. The full stencil logic is preserved.
1339 LOG_ALLOW(LOCAL, LOG_DEBUG, " Adding pressure gradient term to RHS...\n");
1340
1341 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1342
1343 ierr = DMLocalToLocalBegin(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1344 ierr = DMLocalToLocalEnd(fda, Rct, INSERT_VALUES, Rct); CHKERRQ(ierr);
1345
1346 ierr = DMDAVecGetArray(fda, Rct, &rct); CHKERRQ(ierr);
1347
1348 for (k = lzs; k < lze; k++) {
1349 for (j = lys; j < lye; j++) {
1350 for (i = lxs; i < lxe; i++) {
1351 PetscReal dpdc, dpde, dpdz;
1352 dpdc = p[k][j][i+1] - p[k][j][i];
1353
1354 if ((j==my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1355 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1356 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1357 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1358 }
1359 }
1360 else if ((j==my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1361 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1362 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1363 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1364 }
1365 }
1366 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1367 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1368 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1369 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1370 }
1371 }
1372 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1373 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1374 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1375 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1376 }
1377 }
1378 else {
1379 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1380 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1381 }
1382
1383 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1384 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1385 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1386 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1387 }
1388 }
1389 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1390 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1391 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1392 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1393 }
1394 }
1395 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1396 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1397 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1398 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1399 }
1400 }
1401 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1402 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1403 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1404 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1405 }
1406 }
1407 else {
1408 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1409 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1410 }
1411
1412 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1413
1414
1415 rhs[k][j][i].x -=
1416 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1417 icsi[k][j][i].y * icsi[k][j][i].y +
1418 icsi[k][j][i].z * icsi[k][j][i].z)+
1419 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1420 ieta[k][j][i].y * icsi[k][j][i].y +
1421 ieta[k][j][i].z * icsi[k][j][i].z)+
1422 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1423 izet[k][j][i].y * icsi[k][j][i].y +
1424 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1425
1426 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1427 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1428 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1429 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1430 }
1431 }
1432 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1433 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1434 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1435 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1436 }
1437 }
1438 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1439 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1440 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1441 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1442 }
1443 }
1444 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1445 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1446 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1447 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1448 }
1449 }
1450 else {
1451 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1452 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1453 }
1454
1455 dpde = p[k][j+1][i] - p[k][j][i];
1456
1457 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1458 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1459 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1460 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1461 }
1462 }
1463 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1464 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1465 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1466 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1467 }
1468 }
1469 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1470 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1471 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1472 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1473 }
1474 }
1475 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1476 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1477 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1478 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1479 }
1480 }
1481 else {
1482 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1483 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1484 }
1485
1486 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1487
1488
1489 rhs[k][j][i].y -=
1490 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1491 jcsi[k][j][i].y * jeta[k][j][i].y +
1492 jcsi[k][j][i].z * jeta[k][j][i].z) +
1493 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1494 jeta[k][j][i].y * jeta[k][j][i].y +
1495 jeta[k][j][i].z * jeta[k][j][i].z) +
1496 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1497 jzet[k][j][i].y * jeta[k][j][i].y +
1498 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1499
1500 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1501 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1502 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1503 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1504 }
1505 }
1506 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1507 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1508 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1509 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1510 }
1511 }
1512 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1513 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1514 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1515 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1516 }
1517 }
1518 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1519 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1520 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1521 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1522 }
1523 }
1524 else {
1525 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1526 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1527 }
1528
1529 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1530 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1531 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1532 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1533 }
1534 }
1535 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1536 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1537 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1538 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1539 }
1540 }
1541 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1542 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1543 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1544 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1545 }
1546 }
1547 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1548 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1549 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1550 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1551 }
1552 }
1553 else {
1554 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1555 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1556 }
1557
1558 dpdz = (p[k+1][j][i] - p[k][j][i]);
1559
1560 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1561
1562 rhs[k][j][i].z -=
1563 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1564 kcsi[k][j][i].y * kzet[k][j][i].y +
1565 kcsi[k][j][i].z * kzet[k][j][i].z) +
1566 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1567 keta[k][j][i].y * kzet[k][j][i].y +
1568 keta[k][j][i].z * kzet[k][j][i].z) +
1569 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1570 kzet[k][j][i].y * kzet[k][j][i].y +
1571 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1572
1573 }
1574 }
1575 }
1576
1577
1578 //Mohsen March 2012//
1579
1580 // rhs.x at boundaries for periodic bc at i direction//
1581 if (user->bctype[0]==7 && xs==0){
1582 for (k=lzs; k<lze; k++) {
1583 for (j=lys; j<lye; j++) {
1584 i=xs;
1585
1586 dpdc = p[k][j][i+1] - p[k][j][i];
1587
1588 if ((j==my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1589 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1590 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1591 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1592 }
1593 }
1594 else if ((j==my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i]+nvert[k][j+1][i+1] > 0.1) {
1595 if (nvert[k][j-1][i] + nvert[k][j-1][i+1] < 0.1) {
1596 dpde = (p[k][j ][i] + p[k][j ][i+1] -
1597 p[k][j-1][i] - p[k][j-1][i+1]) * 0.5;
1598 }
1599 }
1600 else if ((j == 1 && user->bctype[2]!=7) || nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1601 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1602 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1603 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1604 }
1605 }
1606 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k][j-1][i+1] > 0.1) {
1607 if (nvert[k][j+1][i] + nvert[k][j+1][i+1] < 0.1) {
1608 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1609 p[k][j ][i] - p[k][j ][i+1]) * 0.5;
1610 }
1611 }
1612 else {
1613 dpde = (p[k][j+1][i] + p[k][j+1][i+1] -
1614 p[k][j-1][i] - p[k][j-1][i+1]) * 0.25;
1615 }
1616
1617 if ((k == mz-2 && user->bctype[4]!=7) || nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1618 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1619 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1620 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1621 }
1622 }
1623 else if ((k == mz-2 || k==1) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j][i+1] > 0.1) {
1624 if (nvert[k-1][j][i] + nvert[k-1][j][i+1] < 0.1) {
1625 dpdz = (p[k ][j][i] + p[k ][j][i+1] -
1626 p[k-1][j][i] - p[k-1][j][i+1]) * 0.5;
1627 }
1628 }
1629 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1630 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1631 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1632 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1633 }
1634 }
1635 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j][i+1] > 0.1) {
1636 if (nvert[k+1][j][i] + nvert[k+1][j][i+1] < 0.1) {
1637 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1638 p[k ][j][i] - p[k ][j][i+1]) * 0.5;
1639 }
1640 }
1641 else {
1642 dpdz = (p[k+1][j][i] + p[k+1][j][i+1] -
1643 p[k-1][j][i] - p[k-1][j][i+1]) * 0.25;
1644 }
1645
1646 rhs[k][j][i].x =0.5 * (rct[k][j][i].x + rct[k][j][i+1].x);
1647 rhs[k][j][i].x -=
1648 (dpdc * (icsi[k][j][i].x * icsi[k][j][i].x +
1649 icsi[k][j][i].y * icsi[k][j][i].y +
1650 icsi[k][j][i].z * icsi[k][j][i].z)+
1651 dpde * (ieta[k][j][i].x * icsi[k][j][i].x +
1652 ieta[k][j][i].y * icsi[k][j][i].y +
1653 ieta[k][j][i].z * icsi[k][j][i].z)+
1654 dpdz * (izet[k][j][i].x * icsi[k][j][i].x +
1655 izet[k][j][i].y * icsi[k][j][i].y +
1656 izet[k][j][i].z * icsi[k][j][i].z)) * iaj[k][j][i];
1657 }
1658 }
1659 }
1660
1661// rhs.y at boundaries for periodic bc at j direction//
1662 if (user->bctype[2]==7 && ys==0){
1663 for (k=lzs; k<lze; k++) {
1664 for (i=lxs; i<lxe; i++) {
1665
1666 j=ys;
1667
1668 if ((i == mx-2 && user->bctype[0]!=7) || nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1669 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1670 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1671 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1672 }
1673 }
1674 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k][j+1][i+1] > 0.1) {
1675 if (nvert[k][j][i-1] + nvert[k][j+1][i-1] < 0.1) {
1676 dpdc = (p[k][j][i ] + p[k][j+1][i ] -
1677 p[k][j][i-1] - p[k][j+1][i-1]) * 0.5;
1678 }
1679 }
1680 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1681 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1682 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1683 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1684 }
1685 }
1686 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k][j+1][i-1] > 0.1) {
1687 if (nvert[k][j][i+1] + nvert[k][j+1][i+1] < 0.1) {
1688 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1689 p[k][j][i ] - p[k][j+1][i ]) * 0.5;
1690 }
1691 }
1692 else {
1693 dpdc = (p[k][j][i+1] + p[k][j+1][i+1] -
1694 p[k][j][i-1] - p[k][j+1][i-1]) * 0.25;
1695 }
1696
1697 dpde = p[k][j+1][i] - p[k][j][i];
1698
1699 if ((k == mz-2 && user->bctype[4]!=7)|| nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1700 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1 && (k!=1 || (k==1 && user->bctype[4]==7))) {
1701 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1702 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1703 }
1704 }
1705 else if ((k == mz-2 || k==1 ) && user->bctype[4]==7 && nvert[k+1][j][i] + nvert[k+1][j+1][i] > 0.1) {
1706 if (nvert[k-1][j][i] + nvert[k-1][j+1][i] < 0.1) {
1707 dpdz = (p[k ][j][i] + p[k ][j+1][i] -
1708 p[k-1][j][i] - p[k-1][j+1][i]) * 0.5;
1709 }
1710 }
1711 else if ((k == 1 && user->bctype[4]!=7)|| nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1712 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1713 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1714 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1715 }
1716 }
1717 else if ((k == 1 || k==mz-2) && user->bctype[4]==7 && nvert[k-1][j][i] + nvert[k-1][j+1][i] > 0.1) {
1718 if (nvert[k+1][j][i] + nvert[k+1][j+1][i] < 0.1) {
1719 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1720 p[k ][j][i] - p[k ][j+1][i]) * 0.5;
1721 }
1722 }
1723 else {
1724 dpdz = (p[k+1][j][i] + p[k+1][j+1][i] -
1725 p[k-1][j][i] - p[k-1][j+1][i]) * 0.25;
1726 }
1727
1728 rhs[k][j][i].y =0.5 * (rct[k][j][i].y + rct[k][j+1][i].y);
1729
1730 rhs[k][j][i].y -=
1731 (dpdc * (jcsi[k][j][i].x * jeta[k][j][i].x +
1732 jcsi[k][j][i].y * jeta[k][j][i].y +
1733 jcsi[k][j][i].z * jeta[k][j][i].z)+
1734 dpde * (jeta[k][j][i].x * jeta[k][j][i].x +
1735 jeta[k][j][i].y * jeta[k][j][i].y +
1736 jeta[k][j][i].z * jeta[k][j][i].z)+
1737 dpdz * (jzet[k][j][i].x * jeta[k][j][i].x +
1738 jzet[k][j][i].y * jeta[k][j][i].y +
1739 jzet[k][j][i].z * jeta[k][j][i].z)) * jaj[k][j][i];
1740
1741 }
1742 }
1743 }
1744
1745 // rhs.z at boundaries for periodic bc at k direction//
1746 if (user->bctype[4]==7&& zs==0){
1747 for (j=lys; j<lye; j++) {
1748 for (i=lxs; i<lxe; i++) {
1749
1750 k=zs;
1751
1752 if ((i == mx-2 && user->bctype[0]!=7)|| nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1753 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1 && (i!=1 || (i==1 && user->bctype[0]==7))) {
1754 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1755 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1756 }
1757 }
1758 else if ((i == mx-2 || i==1) && user->bctype[0]==7 && nvert[k][j][i+1] + nvert[k+1][j][i+1] > 0.1) {
1759 if (nvert[k][j][i-1] + nvert[k+1][j][i-1] < 0.1) {
1760 dpdc = (p[k][j][i ] + p[k+1][j][i ] -
1761 p[k][j][i-1] - p[k+1][j][i-1]) * 0.5;
1762 }
1763 }
1764 else if ((i == 1 && user->bctype[0]!=7)|| nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1765 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1766 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1767 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1768 }
1769 }
1770 else if ((i == 1 || i==mx-2) && user->bctype[0]==7 && nvert[k][j][i-1] + nvert[k+1][j][i-1] > 0.1) {
1771 if (nvert[k][j][i+1] + nvert[k+1][j][i+1] < 0.1) {
1772 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1773 p[k][j][i ] - p[k+1][j][i ]) * 0.5;
1774 }
1775 }
1776 else {
1777 dpdc = (p[k][j][i+1] + p[k+1][j][i+1] -
1778 p[k][j][i-1] - p[k+1][j][i-1]) * 0.25;
1779 }
1780
1781 if ((j == my-2 && user->bctype[2]!=7)|| nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1782 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1 && (j!=1 || (j==1 && user->bctype[2]==7))) {
1783 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1784 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1785 }
1786 }
1787 else if ((j == my-2 || j==1) && user->bctype[2]==7 && nvert[k][j+1][i] + nvert[k+1][j+1][i] > 0.1) {
1788 if (nvert[k][j-1][i] + nvert[k+1][j-1][i] < 0.1) {
1789 dpde = (p[k][j ][i] + p[k+1][j ][i] -
1790 p[k][j-1][i] - p[k+1][j-1][i]) * 0.5;
1791 }
1792 }
1793 else if ((j == 1 && user->bctype[2]!=7)|| nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1794 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1795 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1796 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1797 }
1798 }
1799 else if ((j == 1 || j==my-2) && user->bctype[2]==7 && nvert[k][j-1][i] + nvert[k+1][j-1][i] > 0.1) {
1800 if (nvert[k][j+1][i] + nvert[k+1][j+1][i] < 0.1) {
1801 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1802 p[k][j ][i] - p[k+1][j ][i]) * 0.5;
1803 }
1804 }
1805 else {
1806 dpde = (p[k][j+1][i] + p[k+1][j+1][i] -
1807 p[k][j-1][i] - p[k+1][j-1][i]) * 0.25;
1808 }
1809
1810 dpdz = (p[k+1][j][i] - p[k][j][i]);
1811
1812 rhs[k][j][i].z =0.5 * (rct[k][j][i].z + rct[k+1][j][i].z);
1813
1814 rhs[k][j][i].z -=
1815 (dpdc * (kcsi[k][j][i].x * kzet[k][j][i].x +
1816 kcsi[k][j][i].y * kzet[k][j][i].y +
1817 kcsi[k][j][i].z * kzet[k][j][i].z)+
1818 dpde * (keta[k][j][i].x * kzet[k][j][i].x +
1819 keta[k][j][i].y * kzet[k][j][i].y +
1820 keta[k][j][i].z * kzet[k][j][i].z)+
1821 dpdz * (kzet[k][j][i].x * kzet[k][j][i].x +
1822 kzet[k][j][i].y * kzet[k][j][i].y +
1823 kzet[k][j][i].z * kzet[k][j][i].z)) * kaj[k][j][i];
1824
1825 }
1826 }
1827 }
1828
1829 ierr = DMDAVecRestoreArray(fda, Rct, &rct); CHKERRQ(ierr);
1830
1831 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure Gradient added to RHS .\n");
1832 PetscInt TwoD = simCtx->TwoD;
1833
1834 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases initiated .\n");
1835
1836 // 7. Final clean-up for immersed boundaries and 2D cases
1837 for (k=lzs; k<lze; k++) {
1838 for (j=lys; j<lye; j++) {
1839 for (i=lxs; i<lxe; i++) {
1840 if (TwoD==1)
1841 rhs[k][j][i].x =0.;
1842 else if (TwoD==2)
1843 rhs[k][j][i].y =0.;
1844 else if (TwoD==3)
1845 rhs[k][j][i].z =0.;
1846
1847 if (nvert[k][j][i]>0.1) {
1848 rhs[k][j][i].x = 0;
1849 rhs[k][j][i].y = 0;
1850 rhs[k][j][i].z = 0;
1851 }
1852 if (nvert[k][j][i+1]>0.1) {
1853 rhs[k][j][i].x=0;
1854 }
1855 if (nvert[k][j+1][i]>0.1) {
1856 rhs[k][j][i].y=0;
1857 }
1858 if (nvert[k+1][j][i]>0.1) {
1859 rhs[k][j][i].z=0;
1860 }
1861 }
1862 }
1863 }
1864 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Final cleanup and edge-cases complete .\n");
1865
1866 // ========================================================================
1867
1868 // --- Restore all PETSc array pointers ---
1869 // DMDAVecRestoreArray(fda, user->lUcont, &ucont);
1870 DMDAVecRestoreArray(fda, Rhs, &rhs);
1871 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Rhs restored successfully! .\n");
1872
1873 DMDAVecRestoreArray(fda, user->lCsi, &csi);
1874 DMDAVecRestoreArray(fda, user->lEta, &eta);
1875 DMDAVecRestoreArray(fda, user->lZet, &zet);
1876 DMDAVecRestoreArray(da, user->lAj, &aj);
1877 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Face metrics restored successfully! .\n");
1878
1879 DMDAVecRestoreArray(fda, user->lICsi, &icsi);
1880 DMDAVecRestoreArray(fda, user->lIEta, &ieta);
1881 DMDAVecRestoreArray(fda, user->lIZet, &izet);
1882 DMDAVecRestoreArray(da, user->lIAj, &iaj);
1883 LOG_ALLOW(GLOBAL,LOG_DEBUG,"I Face metrics restored successfully! .\n");
1884
1885 DMDAVecRestoreArray(fda, user->lJCsi, &jcsi);
1886 DMDAVecRestoreArray(fda, user->lJEta, &jeta);
1887 DMDAVecRestoreArray(fda, user->lJZet, &jzet);
1888 DMDAVecRestoreArray(da, user->lJAj, &jaj);
1889 LOG_ALLOW(GLOBAL,LOG_DEBUG,"J Face metrics restored successfully! .\n");
1890
1891 DMDAVecRestoreArray(fda, user->lKCsi, &kcsi);
1892 DMDAVecRestoreArray(fda, user->lKEta, &keta);
1893 DMDAVecRestoreArray(fda, user->lKZet, &kzet);
1894 DMDAVecRestoreArray(da, user->lKAj, &kaj);
1895 LOG_ALLOW(GLOBAL,LOG_DEBUG,"K Face metrics restored successfully! .\n");
1896
1897 DMDAVecRestoreArray(da, user->lP, &p);
1898 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Pressure restored successfully! .\n");
1899
1900 DMDAVecRestoreArray(da, user->lNvert, &nvert);
1901 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Nvert restored successfully! .\n");
1902
1903 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Cell Centered scalars restored successfully! .\n");
1904
1905 // --- Destroy temporary work vectors ---
1906 ierr = VecDestroy(&Conv); CHKERRQ(ierr);
1907 ierr = VecDestroy(&Visc); CHKERRQ(ierr);
1908 ierr = VecDestroy(&Rc); CHKERRQ(ierr);
1909 ierr = VecDestroy(&Rct); CHKERRQ(ierr);
1910 LOG_ALLOW(GLOBAL,LOG_DEBUG,"Temporary work vectors destroyed successfully! .\n");
1911
1912 LOG_ALLOW(LOCAL, LOG_DEBUG, "Rank %d, Block %d: RHS computation complete.\n",
1913 simCtx->rank, user->_this);
1914
1915
1916 PROFILE_FUNCTION_END;
1917 PetscFunctionReturn(0);
1918}
LOCAL
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:46
Viscous
PetscErrorCode Viscous(UserCtx *user, Vec Ucont, Vec Ucat, Vec Visc)
Definition rhs.c:504
Convection
PetscErrorCode Convection(UserCtx *user, Vec Ucont, Vec Ucat, Vec Conv)
Definition rhs.c:5
Contra2Cart
PetscErrorCode Contra2Cart(UserCtx *user)
Reconstructs Cartesian velocity (Ucat) at cell centers from contravariant velocity (Ucont) defined on...
Definition setup.c:2084
SimCtx::moveframe
PetscInt moveframe
Definition variables.h:568
SimCtx::TwoD
PetscInt TwoD
Definition variables.h:568
SimCtx::rank
PetscMPIInt rank
Definition variables.h:541
UserCtx::_this
PetscInt _this
Definition variables.h:674
SimCtx::invicid
PetscInt invicid
Definition variables.h:568
UserCtx::lUcont
Vec lUcont
Definition variables.h:688
UserCtx::info
DMDALocalInfo info
Definition variables.h:668
UserCtx::lUcat
Vec lUcat
Definition variables.h:688
SimCtx::rotateframe
PetscInt rotateframe
Definition variables.h:568
Here is the call graph for this function:
Here is the caller graph for this function: