15 double **dcoeff,
int beg,
int end,
Grid *grid)
48 #define D_DX_I(q) (q[k][j][i + 1] - q[k][j][i])
49 #define D_DY_J(q) (q[k][j + 1][i] - q[k][j][i])
50 #define D_DZ_K(q) (q[k + 1][j][i] - q[k][j][i])
52 #define D_DY_I(q) ( 0.25*(q[k][j + 1][i] + q[k][j + 1][i + 1]) \
53 - 0.25*(q[k][j - 1][i] + q[k][j - 1][i + 1]))
55 #define D_DZ_I(q) ( 0.25*(q[k + 1][j][i] + q[k + 1][j][i + 1]) \
56 - 0.25*(q[k - 1][j][i] + q[k - 1][j][i + 1]))
58 #define D_DX_J(q) ( 0.25*(q[k][j][i + 1] + q[k][j + 1][i + 1]) \
59 - 0.25*(q[k][j][i - 1] + q[k][j + 1][i - 1]))
61 #define D_DZ_J(q) ( 0.25*(q[k + 1][j][i] + q[k + 1][j + 1][i]) \
62 - 0.25*(q[k - 1][j][i] + q[k - 1][j + 1][i]))
64 #define D_DX_K(q) ( 0.25*(q[k][j][i + 1] + q[k + 1][j][i + 1]) \
65 - 0.25*(q[k][j][i - 1] + q[k + 1][j][i - 1]))
67 #define D_DY_K(q) ( 0.25*(q[k][j + 1][i] + q[k + 1][j + 1][i]) \
68 - 0.25*(q[k][j - 1][i] + q[k + 1][j - 1][i]))
73 double dx_1, dy_1, dz_1;
76 double *x1r, *x2r, *x3r;
77 double dxVx, dxVy, dxVz, dyVx, dyVy, dyVz, dzVx, dzVy, dzVz;
78 static double *tau_xx, *tau_xy, *tau_xz,
79 *tau_yx, *tau_yy, *tau_yz,
80 *tau_zx, *tau_zy, *tau_zz;
81 double ***Vx, ***Vy, ***Vz;
83 double *r, *th, r_1, dr, dr_1, s_1, tan_1, dy;
84 double dVxi,dVyi,dVzi;
85 double dVxj,dVyj,dVzj;
86 double dVxk,dVyk,dVzk;
88 static double *one_dVr, *one_dmu;
106 #if GEOMETRY != CARTESIAN
108 if (one_dVr == NULL){
111 for (i = 0; i <
NX1_TOT - 1; i++){
112 one_dVr[
i] = r[i+1]*fabs(r[i + 1]) - r[
i]*fabs(r[i]);
113 one_dVr[
i] = 2.0/one_dVr[
i];
115 for (j = 0; j <
NX2_TOT - 1; j++){
116 one_dmu[
j] = 1.0 - cos(th[j + 1]) - (1.0 - cos(th[j]))*(th[j] > 0.0 ? 1.0:-1.0);
117 one_dmu[
j] = 1.0/one_dmu[
j];
136 for (i = beg; i <= end; i++){
142 vi[nv] = 0.5*(V[nv][
k][
j][
i] + V[nv][
k][
j][i+1]);
143 vc[nv] = V[nv][
k][
j][
i];
148 Visc_nu(vi, x1r[i], x2[j], x3[k], &nu1, &nu2);
150 tau_xx[
i] = tau_xy[
i]= tau_xz[
i]= tau_yx[
i]=
151 tau_yy[
i] = tau_yz[
i]= tau_zx[
i]= tau_zy[
i]=
154 dVxi = dVxj = dVxk = dVyi = dVyj = dVyk = dVzi = dVzj = dVzk = 0.0;
155 dxVx = dxVy = dxVz = dyVx = dyVy = dyVz = dzVx = dzVy = dzVz = 0.0;
165 EXPAND(dxVx = dVxi*dx_1; , dxVy = dVyi*dx_1; , dxVz = dVzi*dx_1; )
167 EXPAND(dyVx = dVxj*dy_1; , dyVy = dVyj*dy_1; , dyVz = dVzj*dy_1; )
169 dzVx = dVxk*dz_1; dzVy = dVyk*dz_1; dzVz = dVzk*dz_1;
173 #if GEOMETRY == CARTESIAN
174 div_v =
D_EXPAND(dxVx, + dyVy , + dzVz);
178 tau_xx[
i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v;
179 tau_xy[
i] = nu1*(dyVx + dxVy);
180 tau_xz[
i] = nu1*(dzVx + dxVz);
181 tau_yx[
i] = tau_xy[
i];
182 tau_zx[
i] = tau_xz[
i];
186 EXPAND(ViS[i][
MX1] = 0.0; ,
193 dr_1 = 1.0/dr; dz_1 = 0.0;
194 r_1 = 1.0/grid[
IDIR].
x[
i];
198 dxVx = (Vx[
k][
j][i+1]*r[i+1] - Vx[
k][
j][
i]*fabs(r[i]))*one_dVr[i];
199 div_v =
D_EXPAND( dxVx, + dVyj*dy_1, + 0.0);
204 tau_xx[
i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v;
205 tau_xy[
i] = nu1*(dyVx + dxVy);
206 EXPAND(tau_xz[i] = nu1*(0.0);,
207 tau_xz[
i] = nu1*(0.0);,
208 tau_xz[
i] = nu1*0.5*(r[
i]+r[i+1])*dr_1*((1./r[i+1])*Vz[
k][
j][i+1]
209 - (1./r[
i])*Vz[k][j][i]);)
211 tau_yx[i] = tau_xy[i];
212 tau_zx[
i] = tau_xz[
i];
217 Visc_nu(vc, x1[i], x2[j], x3[k], &nu1, &nu2);
219 div_v =
D_EXPAND( 0.5*(Vx[k][j][i+1]-Vx[k][j][i-1])*dx_1 + Vx[k][j][i]*r_1,
220 + 0.5*(Vy[k][j + 1][i]-Vy[k][j - 1][i])*dy_1,
223 tau_zz[
i] = 2.0*nu1*r_1*Vx[
k][
j][
i]
224 + (nu2 - (2.0/3.0)*nu1)*div_v;
228 EXPAND(ViS[i][
MX1] = -tau_zz[i]*r_1; ,
235 dr = grid[
IDIR].
dx[
i]; dr_1 = 1.0/dr;
240 div_v =
D_EXPAND(r_1*vi[
VX1] + dVxi*dr_1, + r_1*dVyj*dy_1, + dVzk*dz_1);
244 tau_xx[
i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v;
245 EXPAND(tau_xy[i] = nu1*dxVy; ,
246 tau_xy[
i] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]); ,
247 tau_xy[
i] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]);)
249 tau_xz[i] = nu1*(dzVx + dxVz);
250 tau_yx[
i] = tau_xy[
i];
251 EXPAND(tau_yy[i] = 2.0*nu1*r_1*dyVy + (nu2 - (2.0/3.0)*nu1)*div_v;,
252 tau_yy[
i] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
253 + (nu2 - (2.0/3.0)*nu1)*div_v;,
254 tau_yy[
i] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
255 + (nu2 - (2.0/3.0)*nu1)*div_v;)
257 tau_zx[i] = tau_xz[i];
261 r_1 = 1.0/grid[
IDIR].
x[
i];
263 EXPAND(ViS[i][
MX1] = -0.5*(tau_yy[i - 1] + tau_yy[i])*r_1; ,
271 tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
276 + r_1*dVyj*dy_1 + r_1*tan_1*vi[
VX2],
277 + r_1*s_1*dVzk*dz_1);
282 tau_xx[
i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v;
283 EXPAND(tau_xy[i] = nu1*(r_1*dyVx + dxVy ); ,
284 tau_xy[
i] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]); ,
285 tau_xy[
i] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]);)
287 EXPAND(tau_xz[i] = nu1*(r_1*s_1*dzVx + dxVz) ; ,
288 tau_xz[
i] = nu1*(r_1*s_1*dzVx + dxVz ); ,
289 tau_xz[
i] = nu1*(r_1*s_1*dzVx + dxVz - r_1*vi[
VX3]);)
291 tau_yx[i] = tau_xy[i];
292 tau_yy[
i] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
293 + (nu2 - (2.0/3.0)*nu1)*div_v;
295 tau_zx[
i] = tau_xz[
i];
296 EXPAND(tau_zz[i] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1] )
297 + (nu2 - (2.0/3.0)*nu1)*div_v;,
298 tau_zz[
i] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1]
300 + (nu2 - (2.0/3.0)*nu1)*div_v;,
301 tau_zz[
i] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1]
303 + (nu2 - (2.0/3.0)*nu1)*div_v;)
308 r_1 = 1.0/grid[
IDIR].x[i];
309 EXPAND(ViS[i][
MX1] = -0.5*( tau_yy[i - 1] + tau_yy[i]
310 + tau_zz[i - 1] + tau_zz[i])*r_1;,
318 EXPAND(ViF[i][
MX1] = tau_xx[i]; ,
319 ViF[
i][
MX2] = tau_xy[
i]; ,
320 ViF[
i][
MX3] = tau_xz[
i]; )
323 ViF[i][ENG] = EXPAND( vi[VX1]*tau_xx[i],
325 + vi[
VX3]*tau_zx[i]);
328 dcoeff[
i][
MX1] =
MAX(nu1, nu2);
340 for (j = beg ; j <= end; j++){
346 vi[nv] = 0.5*(V[nv][
k][j+1][
i] + V[nv][
k][
j][
i]);
347 vc[nv] = V[nv][
k][
j][
i];
352 Visc_nu(vi, x1[i], x2r[j], x3[k], &nu1, &nu2);
354 tau_xx[
j]= tau_xy[
j]= tau_xz[
j]= tau_yx[
j] = tau_yy[
j]=
355 tau_yz[
j]= tau_zx[
j]= tau_zy[
j]= tau_zz[
j] =0.0;
357 dVxi = dVxj = dVxk = dVyi = dVyj = dVyk = dVzi = dVzj = dVzk = 0.0;
358 dxVx = dxVy = dxVz = dyVx = dyVy = dyVz = dzVx = dzVy = dzVz = 0.0;
368 EXPAND(dxVx = dVxi*dx_1; , dxVy = dVyi*dx_1; , dxVz = dVzi*dx_1; )
369 EXPAND(dyVx = dVxj*dy_1; , dyVy = dVyj*dy_1; , dyVz = dVzj*dy_1; )
371 dzVx = dVxk*dz_1; dzVy = dVyk*dz_1; dzVz = dVzk*dz_1;
374 #if GEOMETRY == CARTESIAN
378 div_v =
D_EXPAND(dxVx, + dyVy, + dzVz);
382 tau_xy[
j] = nu1*(dyVx + dxVy);
383 tau_yx[
j] = tau_xy[
j];
384 tau_yy[
j] = 2.0*nu1*dyVy + (nu2 - (2.0/3.0)*nu1)*div_v;
385 tau_yz[
j] = nu1*(dyVz + dzVy);
386 tau_zy[
j] = tau_yz[
j];
390 EXPAND(ViS[j][
MX1] = 0.0; ,
397 dr = grid[
IDIR].
dx[
i]; dr_1 = 1.0/dr;
398 r_1 = 1.0/grid[
IDIR].
x[
i]; dz_1 = 0.0;
402 div_v =
D_EXPAND(r_1*vi[
VX1] + dVxi*dr_1, + dVyj*dy_1, + 0.0);
406 tau_xy[
j] = nu1*(dyVx + dxVy);
407 tau_yx[
j] = tau_xy[
j];
408 tau_yy[
j] = 2.0*nu1*dyVy + (nu2 - (2.0/3.0)*nu1)*div_v;
409 EXPAND(tau_yz[j] = 0.0; ,
411 tau_yz[
j] = nu1*(dyVz);)
412 tau_zy[j] = tau_yz[j];
416 EXPAND(ViS[j][
MX1] = 0.0; ,
423 dr_1 = 1.0/dr; r_1 = 1.0/grid[
IDIR].
x[
i];
427 div_v =
D_EXPAND(r_1*vi[
VX1] + dVxi*dr_1, + r_1*dVyj*dy_1, + dVzk*dz_1);
431 tau_xy[
j] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]);
433 tau_yx[
j] = tau_xy[
j];
434 tau_yy[
j] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
435 + (nu2 - (2.0/3.0)*nu1)*div_v;
437 tau_yz[
j] = nu1*(r_1*dyVz + dzVy);
438 tau_zy[
j] = tau_yz[
j];
440 r_1 = 1.0/grid[
IDIR].
x[
i];
444 EXPAND(ViS[j][
MX1] = 0.0; ,
452 tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
465 if (fabs(tan(th[j]))< 1.e-12) tan_1 = 0.0;
472 + ( sin(th[j + 1])*Vy[k][j + 1][i] - fabs(sin(th[j]))*Vy[k][j][i])*r_1*one_dmu[j],
473 + r_1*s_1*dVzk*dz_1 );
479 tau_xy[
j] = nu1*(r_1*dyVx + dxVy - r_1*vi[
VX2]);
481 tau_yx[
j] = tau_xy[
j];
482 tau_yy[
j] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
483 + (nu2 - (2.0/3.0)*nu1)*div_v;
486 EXPAND(tau_yz[j] = nu1*r_1*dyVz; ,
487 tau_yz[
j] = nu1*r_1*dyVz; ,
488 tau_yz[
j] = nu1*(r_1*dyVz - tan_1*r_1*vi[
VX3]));
491 tau_yz[
j] = nu1*(s_1*r_1*dzVy + r_1*dyVz - tan_1*r_1*vi[
VX3]);
493 tau_zy[
j] = tau_yz[
j];
495 EXPAND(tau_zz[j] = 2.0*nu1*r_1*vi[
VX1]
496 + (nu2 - (2.0/3.0)*nu1)*div_v; ,
497 tau_zz[
j] = 2.0*nu1*(r_1*vi[
VX1] + tan_1*r_1*vi[
VX2])
498 + (nu2 - (2.0/3.0)*nu1)*div_v;,
499 tau_zz[
j] = 2.0*nu1*(r_1*vi[
VX1] + tan_1*r_1*vi[
VX2])
500 + (nu2 - (2.0/3.0)*nu1)*div_v;)
503 tau_zz[j] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1] + tan_1*r_1*vi[
VX2])
504 + (nu2 - (2.0/3.0)*nu1)*div_v;
509 th = grid[
JDIR].
x; tan_1= 1.0/tan(th[j]);
511 EXPAND(ViS[j][
MX1] = 0.0; ,
512 ViS[
j][
MX2] = 0.5*( tau_yx[j - 1] + tau_yx[
j])*r_1
513 - tan_1*0.5*(tau_zz[j - 1] + tau_zz[j])*r_1; ,
520 EXPAND(ViF[j][
MX1] = tau_yx[j]; ,
521 ViF[
j][
MX2] = tau_yy[
j]; ,
522 ViF[
j][
MX3] = tau_yz[
j]; )
525 ViF[j][ENG] = EXPAND( vi[
VX1]*tau_xy[j],
527 + vi[
VX3]*tau_zy[j]);
542 for (k = beg; k <= end; k++){
548 vi[nv] = 0.5*(V[nv][
k][
j][
i] + V[nv][k+1][
j][
i]);
549 vc[nv] = V[nv][
k][
j][
i];
554 Visc_nu(vi, x1[i], x2[j], x3r[k], &nu1, &nu2);
556 tau_xx[
k] = tau_xy[
k] = tau_xz[
k]= tau_yx[
k]= tau_yy[
k]= tau_yz[
k]=
557 tau_zx[
k] = tau_zy[
k] = tau_zz[
k]=0.0;
558 dVxi = dVxj = dVxk = dVyi = dVyj = dVyk = dVzi = dVzj = dVzk = 0.0;
566 dxVx = dVxi*dx_1; dxVy = dVyi*dx_1; dxVz = dVzi*dx_1;
567 dyVx = dVxj*dy_1; dyVy = dVyj*dy_1; dyVz = dVzj*dy_1;
568 dzVx = dVxk*dz_1; dzVy = dVyk*dz_1; dzVz = dVzk*dz_1;
570 #if GEOMETRY == CARTESIAN
574 div_v = dxVx + dyVy + dzVz;
578 tau_xz[
k] = nu1*(dzVx + dxVz);
579 tau_yz[
k] = nu1*(dyVz + dzVy);
580 tau_zx[
k] = tau_xz[
k];
581 tau_zy[
k] = tau_yz[
k];
582 tau_zz[
k] = 2.0*nu1*dzVz + (nu2 - (2.0/3.0)*nu1)*div_v;
586 EXPAND(ViS[k][
MX1] = 0.0; ,
593 dr = grid[
IDIR].
dx[
i]; dr_1 = 1.0/dr;
594 r_1 = 1.0/grid[
IDIR].
x[
i];
598 div_v = r_1*vi[
VX1] + dVxi*dr_1 + r_1*dVyj*dy_1 + dVzk*dz_1;
602 tau_xz[
k] = nu1*(dzVx + dxVz);
603 tau_yy[
k] = 2.0*nu1*(r_1*dyVy + r_1*vi[
VX1])
604 + (nu2 - (2.0/3.0)*nu1)*div_v;
606 tau_yz[
k] = nu1*(r_1*dyVz + dzVy);
607 tau_zx[
k] = tau_xz[
k];
608 tau_zy[
k] = tau_yz[
k];
609 tau_zz[
k] = 2.0*nu1*dzVz + (nu2 - (2.0/3.0)*nu1)*div_v;
612 EXPAND(ViS[k][
MX1] = 0.0; ,
619 r_1 = 1.0/grid[
IDIR].
x[
i]; tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
623 div_v = 2.0*r_1*vi[
VX1] + dVxi*dr_1 + r_1*dVyj*dy_1 + r_1*tan_1*vi[
VX2]
628 tau_xz[
k] = nu1*(r_1*s_1*dzVx + dxVz - r_1*vi[
VX3]);
630 tau_yz[
k] = nu1*(s_1*r_1*dzVy + r_1*dyVz - tan_1*r_1*vi[
VX3]);
632 tau_zx[
k] = tau_xz[
k];
633 tau_zy[
k] = tau_yz[
k];
634 tau_zz[
k] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[
VX1]
636 + (nu2 - (2.0/3.0)*nu1)*div_v;
641 EXPAND(ViS[k][
MX1] = 0.0; ,
649 ViF[k][
MX1] = tau_zx[k];
650 ViF[
k][
MX2] = tau_zy[
k];
651 ViF[
k][
MX3] = tau_zz[
k];
653 #if EOS != ISOTHERMAL
654 ViF[
k][ENG] = vi[
VX1]*tau_xz[
k]
658 dcoeff[
k][
MX1] =
MAX(nu1, nu2);
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
void ViscousFlux(Data_Arr V, double **ViF, double **ViS, double **dcoeff, int beg, int end, Grid *grid)
int g_dir
Specifies the current sweep or direction of integration.
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
#define ARRAY_1D(nx, type)
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
void Visc_nu(double *v, double x1, double x2, double x3, double *nu1, double *nu2)
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.