PLUTO
Functions
viscosity.h File Reference

Go to the source code of this file.

Functions

void ViscousFlux (Data_Arr, double **, double **, double **, int, int, Grid *)
 
void Visc_nu (double *, double, double, double, double *, double *)
 Calculate first and second viscosity coefficients as functions of data and coordinates. More...
 

Function Documentation

void Visc_nu ( double *  v,
double  x1,
double  x2,
double  x3,
double *  nu1,
double *  nu2 
)

Calculate first and second viscosity coefficients as functions of data and coordinates.

Parameters
[in]vpointer to data array containing cell-centered quantities
[in]x1real, coordinate value
[in]x2real, coordinate value
[in]x3real, coordinate value
[in,out]nu1pointer to first viscous coefficient
[in,out]nu2pointer to second viscous coefficient
Returns
This function has no return value.

Definition at line 7 of file visc_nu.c.

20 {
21  *nu1 = 0.0;
22  *nu2 = 0.0;
23 }

Here is the caller graph for this function:

void ViscousFlux ( Data_Arr  V,
double **  ViF,
double **  ViS,
double **  dcoeff,
int  beg,
int  end,
Grid *  grid 
)
Parameters
[in]Vdata array containing cell-centered quantities
[in,out]ViFpointer to viscous fluxes
[in,out]ViSpointer to viscous source terms
[in,out]dcoeffpointer to diffusion coefficient for dt calculation
[in]beginteger, index for loop beg
[in]endinteger, index for loop end
[in]gridpointer to array of Grid structures
Returns
This function has no return value.

Calculates the stress tensor components (at the i+1/2 face of each cell) and adds explicit viscous terms to the energy and momentum equation. It is called in the during the sweep integrators. The stress tensor is given by

\[ T = \left( \begin{array}{ccc} T_{xx} & T_{xy} & T_{xz} \\ T_{yx} & T_{yy} & T_{yz} \\ T_{zx} & T_{zy} & T_{zz} \end{array}\right) \]

where $ T_{ij} = T_{ji}$ and the components are given by $ T_{ij} = 2\,m\,e(ij) + (l - 2/3 m) \nabla\cdot {\bf V}\, \delta_{ij} $ where $ e(ij)= e_{ij}/(h_ih_j)$ and $ e_{ij} = 0.5( V_{i;j} + V_{j;i} ) $ whereas m,l are the first and second parameter of viscosity respectively.

Definition at line 14 of file viscous_flux.c.

69 {
70  int i,j,k,n,nv;
71  double nu1,nu2;
72  double div_v;
73  double dx_1, dy_1, dz_1;
74  double dx1, dx2, dx3;
75  double *x1, *x2, *x3;
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;
82  double scrh;
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;
87  double vc[NVAR], vi[NVAR]; /* Center and interface values */
88  static double *one_dVr, *one_dmu; /*auxillary volume components for r_1 singularity @ cylindrical and spherical*/
89 
90  EXPAND(Vx = V[VX1]; ,
91  Vy = V[VX2]; ,
92  Vz = V[VX3];)
93 
94  if (tau_xx == NULL){
95  tau_xx = ARRAY_1D(NMAX_POINT, double);
96  tau_xy = ARRAY_1D(NMAX_POINT, double);
97  tau_xz = ARRAY_1D(NMAX_POINT, double);
98  tau_yx = ARRAY_1D(NMAX_POINT, double);
99  tau_yy = ARRAY_1D(NMAX_POINT, double);
100  tau_yz = ARRAY_1D(NMAX_POINT, double);
101  tau_zx = ARRAY_1D(NMAX_POINT, double);
102  tau_zy = ARRAY_1D(NMAX_POINT, double);
103  tau_zz = ARRAY_1D(NMAX_POINT, double);
104  }
105 
106  #if GEOMETRY != CARTESIAN
107  r = grid[IDIR].x; th = grid[JDIR].x;
108  if (one_dVr == NULL){
109  one_dVr = ARRAY_1D(NX1_TOT, double); /* -- intercell (i) and (i+1) volume -- */
110  one_dmu = ARRAY_1D(NX2_TOT, double); /* -- intercell (j) and (j+1) volume -- */
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];
114  }
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];
118  }
119  }
120  #endif
121 
122 /* -- set pointers to coordinates and grid indices --*/
123 
124  x1 = grid[IDIR].x; x1r = grid[IDIR].xr; i = g_i;
125  x2 = grid[JDIR].x; x2r = grid[JDIR].xr; j = g_j;
126  x3 = grid[KDIR].x; x3r = grid[KDIR].xr; k = g_k;
127 
128  if (g_dir == IDIR){
129 
130  /* ---------------------------------------------------------
131  Loop on X1 direction
132  --------------------------------------------------------- */
133 
134  dy_1 = 1.0/grid[JDIR].dx[j];
135  dz_1 = 1.0/grid[KDIR].dx[k];
136  for (i = beg; i <= end; i++){
137  dx_1 = 1.0/grid[IDIR].dx[i];
138 
139  /* -- compute face- and cell-centered values */
140 
141  VAR_LOOP(nv) {
142  vi[nv] = 0.5*(V[nv][k][j][i] + V[nv][k][j][i+1]);
143  vc[nv] = V[nv][k][j][i];
144  }
145 
146  /* -- compute viscosity (face center) -- */
147 
148  Visc_nu(vi, x1r[i], x2[j], x3[k], &nu1, &nu2);
149 
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]=
152  tau_zz[i] = 0.0;
153 
154  dVxi = dVxj = dVxk = dVyi = dVyj = dVyk = dVzi = dVzj = dVzk = 0.0;
155  dxVx = dxVy = dxVz = dyVx = dyVy = dyVz = dzVx = dzVy = dzVz = 0.0;
156 
157  EXPAND (dVxi = D_DX_I(Vx);, dVyi = D_DX_I(Vy);, dVzi = D_DX_I(Vz);)
158  #if DIMENSIONS >= 2
159  EXPAND (dVxj = D_DY_I(Vx);, dVyj = D_DY_I(Vy);, dVzj = D_DY_I(Vz);)
160  #if DIMENSIONS == 3
161  dVxk = D_DZ_I(Vx); dVyk = D_DZ_I(Vy); dVzk = D_DZ_I(Vz);
162  #endif
163  #endif
164 
165  EXPAND(dxVx = dVxi*dx_1; , dxVy = dVyi*dx_1; , dxVz = dVzi*dx_1; )
166  #if DIMENSIONS >= 2
167  EXPAND(dyVx = dVxj*dy_1; , dyVy = dVyj*dy_1; , dyVz = dVzj*dy_1; )
168  #if DIMENSIONS == 3
169  dzVx = dVxk*dz_1; dzVy = dVyk*dz_1; dzVz = dVzk*dz_1;
170  #endif
171  #endif
172 
173  #if GEOMETRY == CARTESIAN
174  div_v = D_EXPAND(dxVx, + dyVy , + dzVz);
175 
176  /* -- stress tensor components (only some needed for flux computation) -- */
177 
178  tau_xx[i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v; /*2eta1 dxVx + (eta2 - 2/3 eta1) divV*/
179  tau_xy[i] = nu1*(dyVx + dxVy); /*eta1 (dyVx + dxVy)*/
180  tau_xz[i] = nu1*(dzVx + dxVz); /*eta1 (dzVx + dxVz)*/
181  tau_yx[i] = tau_xy[i];
182  tau_zx[i] = tau_xz[i];
183 
184  /* -- compute source terms -- */
185 
186  EXPAND(ViS[i][MX1] = 0.0; ,
187  ViS[i][MX2] = 0.0; ,
188  ViS[i][MX3] = 0.0; )
189 
190  #elif GEOMETRY == CYLINDRICAL
191 
192  r = grid[IDIR].x; dr = grid[IDIR].dx[i];
193  dr_1 = 1.0/dr; dz_1 = 0.0;
194  r_1 = 1.0/grid[IDIR].x[i];
195 
196  /* -- calculate the div(v) (trick @ axis for dxVx) -- */
197 
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);
200  dxVx = dVxi*dx_1;
201 
202  /* -- stress tensor components (only some needed for flux computation) -- */
203 
204  tau_xx[i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v; /* tau_rr = 2 eta1 drVr + (eta2 - 2/3 eta1) divV */
205  tau_xy[i] = nu1*(dyVx + dxVy); /*tau_rz = eta1 (dzVr + drVz)*/
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]);)
210  /*tau_rphi = eta1 (1/r dphiVr + drVphi - 1/r Vphi)= eta1 (r dr(1/r Vphi) )*/
211  tau_yx[i] = tau_xy[i]; /*tau_zr*/
212  tau_zx[i] = tau_xz[i]; /*tau_phir*/
213 
214  /* -- we calculate at the center cause we don't need it for flux
215  but for src, avoiding 1/r->inf at r=r_f=0 -- */
216 
217  Visc_nu(vc, x1[i], x2[j], x3[k], &nu1, &nu2);
218 
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,
221  + 0.0 );
222 
223  tau_zz[i] = 2.0*nu1*r_1*Vx[k][j][i]
224  + (nu2 - (2.0/3.0)*nu1)*div_v;
225 
226  /* -- compute source terms -- */
227 
228  EXPAND(ViS[i][MX1] = -tau_zz[i]*r_1; ,
229  ViS[i][MX2] = 0.0; ,
230  ViS[i][MX3] = 0.0; )
231 
232  #elif GEOMETRY == POLAR
233 
234  r = grid[IDIR].xr; th = grid[JDIR].x;
235  dr = grid[IDIR].dx[i]; dr_1 = 1.0/dr;
236  r_1 = 1.0/grid[IDIR].xr[i];
237 
238  /* -- calculate div(v) -- */
239 
240  div_v = D_EXPAND(r_1*vi[VX1] + dVxi*dr_1, + r_1*dVyj*dy_1, + dVzk*dz_1);
241 
242  /* -- stress tensor components (only some needed for flux computation) -- */
243 
244  tau_xx[i] = 2.0*nu1*dxVx + (nu2 - (2.0/3.0)*nu1)*div_v; /*tau_rr as in cylindrical */
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]);)
248  /*tau_rphi = eta1 (1/r dphiVr + drVphi - 1/r Vphi) */
249  tau_xz[i] = nu1*(dzVx + dxVz); /*tau_rz as in cylindrical*/
250  tau_yx[i] = tau_xy[i]; /*tau_phir = tau_rphi*/
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;)
256  /*tau_phiphi = 2 eta1 (1/r dphiVphi + 1/r Vr) + (eta2 - 2/3 eta1) divV*/
257  tau_zx[i] = tau_xz[i]; /*tau_zr = tau_rz*/
258 
259  /* -- compute source terms -- */
260 
261  r_1 = 1.0/grid[IDIR].x[i];
262 
263  EXPAND(ViS[i][MX1] = -0.5*(tau_yy[i - 1] + tau_yy[i])*r_1; ,
264  ViS[i][MX2] = 0.0; ,
265  ViS[i][MX3] = 0.0;)
266 
267  #elif GEOMETRY == SPHERICAL
268 
269  r = grid[IDIR].xr; th = grid[JDIR].x;
270  dr_1 = 1.0/grid[IDIR].dx[i]; r_1 = 1.0/grid[IDIR].xr[i];
271  tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
272 
273  /* -- calculate div(v) -- */
274 
275  div_v = D_EXPAND( 2.0*r_1*vi[VX1] + dVxi*dr_1 ,
276  + r_1*dVyj*dy_1 + r_1*tan_1*vi[VX2],
277  + r_1*s_1*dVzk*dz_1);
278 
279  /* -- stress tensor components (only some needed for flux computation) -- */
280 
281  /* tau_rr = 2eta1 drVr + + (eta2 - 2/3 eta1) divV */
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]);)
286  /*tau_rthita = eta1 (1/r dthitaVr + drVthita -1/r Vthita)*/
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]);)
290  /*tau_rphi = eta1 (1/r dthitaVr + drVthita -1/r Vthita)*/
291  tau_yx[i] = tau_xy[i]; /*tau_thitar*/
292  tau_yy[i] = 2.0*nu1*(r_1*dyVy + r_1*vi[VX1])
293  + (nu2 - (2.0/3.0)*nu1)*div_v;
294  /*tau_thitathita= 2 eta1 (1/r dthitaVthita + 1/r Vr) + (eta2 - 2/3 eta1)divV*/
295  tau_zx[i] = tau_xz[i]; /*tau_phir*/
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]
299  + tan_1*r_1*vi[VX2])
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]
302  + tan_1*r_1*vi[VX2])
303  + (nu2 - (2.0/3.0)*nu1)*div_v;)
304  /*tau_phiphi= 2 eta1 (1/rs dphiVphi + 1/r Vr + 1/r cot Vthita) + (eta2 - 2/3 eta1)divV*/
305 
306  /* -- compute source terms -- */
307 
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;,
311  ViS[i][MX2] = 0.0;,
312  ViS[i][MX3] = 0.0;)
313 
314  #endif /* -- end #if GEOMETRY -- */
315 
316  /* -- compute fluxes -- */
317 
318  EXPAND(ViF[i][MX1] = tau_xx[i]; ,
319  ViF[i][MX2] = tau_xy[i]; ,
320  ViF[i][MX3] = tau_xz[i]; )
321 
322  #if EOS != ISOTHERMAL
323  ViF[i][ENG] = EXPAND( vi[VX1]*tau_xx[i],
324  + vi[VX2]*tau_yx[i],
325  + vi[VX3]*tau_zx[i]);
326  #endif
327 
328  dcoeff[i][MX1] = MAX(nu1, nu2);
329  dcoeff[i][MX1] /= vi[RHO];
330  }
331 
332  }else if (g_dir == JDIR){
333 
334  /* ---------------------------------------------------------
335  Loop on X2 direction
336  --------------------------------------------------------- */
337 
338  dx_1 = 1.0/grid[IDIR].dx[i];
339  dz_1 = 1.0/grid[KDIR].dx[k];
340  for (j = beg ; j <= end; j++){
341  dy_1 = 1.0/grid[JDIR].dx[j];
342 
343  /* -- compute face- and cell-centered values */
344 
345  VAR_LOOP(nv) {
346  vi[nv] = 0.5*(V[nv][k][j+1][i] + V[nv][k][j][i]);
347  vc[nv] = V[nv][k][j][i];
348  }
349 
350  /* -- compute viscosity (face center) -- */
351 
352  Visc_nu(vi, x1[i], x2r[j], x3[k], &nu1, &nu2);
353 
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;
356 
357  dVxi = dVxj = dVxk = dVyi = dVyj = dVyk = dVzi = dVzj = dVzk = 0.0;
358  dxVx = dxVy = dxVz = dyVx = dyVy = dyVz = dzVx = dzVy = dzVz = 0.0;
359 
360  /* ------ CALCULATE DERIVATIVES ---------- */
361 
362  EXPAND (dVxi = D_DX_J(Vx);, dVyi = D_DX_J(Vy);, dVzi = D_DX_J(Vz);)
363  EXPAND (dVxj = D_DY_J(Vx);, dVyj = D_DY_J(Vy);, dVzj = D_DY_J(Vz);)
364  #if DIMENSIONS == 3
365  dVxk = D_DZ_J(Vx); dVyk = D_DZ_J(Vy); dVzk = D_DZ_J(Vz);
366  #endif
367 
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; )
370  #if DIMENSIONS == 3
371  dzVx = dVxk*dz_1; dzVy = dVyk*dz_1; dzVz = dVzk*dz_1;
372  #endif
373 
374  #if GEOMETRY == CARTESIAN
375 
376  /* -- calculate div(v) -- */
377 
378  div_v = D_EXPAND(dxVx, + dyVy, + dzVz);
379 
380  /* -- stress tensor components (only some needed for flux computation) -- */
381 
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];
387 
388  /* -- compute source terms -- */
389 
390  EXPAND(ViS[j][MX1] = 0.0; ,
391  ViS[j][MX2] = 0.0; ,
392  ViS[j][MX3] = 0.0; )
393 
394  #elif GEOMETRY == CYLINDRICAL
395 
396  r = grid[IDIR].x; th = grid[KDIR].x;
397  dr = grid[IDIR].dx[i]; dr_1 = 1.0/dr;
398  r_1 = 1.0/grid[IDIR].x[i]; dz_1 = 0.0;
399 
400  /* -- calculate div(v) -- */
401 
402  div_v = D_EXPAND(r_1*vi[VX1] + dVxi*dr_1, + dVyj*dy_1, + 0.0);
403 
404  /* -- stress tensor components (only some needed for flux computation) -- */
405 
406  tau_xy[j] = nu1*(dyVx + dxVy); /*tau_rz = eta1 (dzVr + drVz)*/
407  tau_yx[j] = tau_xy[j]; /*tau_zr */
408  tau_yy[j] = 2.0*nu1*dyVy + (nu2 - (2.0/3.0)*nu1)*div_v; /*tau_zz = 2eta1 dzVz + (eta2 - 2/3 eta1)divV*/
409  EXPAND(tau_yz[j] = 0.0; ,
410  tau_yz[j] = 0.0; ,
411  tau_yz[j] = nu1*(dyVz);) /*tau_zphi = eta1(1/r dphiVz + dzVphi)*/
412  tau_zy[j] = tau_yz[j]; /*tau_phiz*/
413 
414  /* -- compute source terms -- */
415 
416  EXPAND(ViS[j][MX1] = 0.0; ,
417  ViS[j][MX2] = 0.0; ,
418  ViS[j][MX3] = 0.0; )
419 
420  #elif GEOMETRY == POLAR
421 
422  r = grid[IDIR].x; th = grid[JDIR].xr; dr = grid[IDIR].dx[i];
423  dr_1 = 1.0/dr; r_1 = 1.0/grid[IDIR].x[i];
424 
425  /* -- calculate div(v) -- */
426 
427  div_v = D_EXPAND(r_1*vi[VX1] + dVxi*dr_1, + r_1*dVyj*dy_1, + dVzk*dz_1);
428 
429  /* -- stress tensor components (only some needed for flux computation) -- */
430 
431  tau_xy[j] = nu1*(r_1*dyVx + dxVy - r_1*vi[VX2]);
432  /*tau_rphi = eta1(1/r dphiVr + drVphi -1/r Vphi)*/
433  tau_yx[j] = tau_xy[j]; /*tau_phir*/
434  tau_yy[j] = 2.0*nu1*(r_1*dyVy + r_1*vi[VX1])
435  + (nu2 - (2.0/3.0)*nu1)*div_v;
436  /*tau_phiphi = 2 eta1 (1/r dphiVphi + 1/rVr) + (eta2 - 2/3 eta1)divV*/
437  tau_yz[j] = nu1*(r_1*dyVz + dzVy); /* tau_phiz = eta1(dzVphi + 1/r dphiVz)*/
438  tau_zy[j] = tau_yz[j]; /*tau_zphi*/
439 
440  r_1 = 1.0/grid[IDIR].x[i];
441 
442  /* -- compute source terms -- */
443 
444  EXPAND(ViS[j][MX1] = 0.0; ,
445  ViS[j][MX2] = 0.0; ,
446  ViS[j][MX3] = 0.0; )
447 
448  #elif GEOMETRY == SPHERICAL
449 
450  r = grid[IDIR].x; th = grid[JDIR].xr;
451  dr_1 = 1.0/grid[IDIR].dx[i]; r_1 = 1.0/grid[IDIR].x[i];
452  tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
453 
454  /*------------------------------------
455  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456  Hotfix for tan_1 at the axis:
457  since we have terms tan_1*smth that
458  cannot be treated with the volume
459  trick, we set an if condition
460  for this to go to zero, as it is
461  the correct behaviour of the
462  term in question there.
463  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
464  --------------------------------------*/
465  if (fabs(tan(th[j]))< 1.e-12) tan_1 = 0.0;
466 
467  /* -- calculate div(v) -- */
468 
469  th = grid[JDIR].x;
470 
471  div_v = D_EXPAND( 2.0*r_1*vi[VX1] + dVxi*dr_1 ,
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 );
474 
475  th = grid[JDIR].xr;
476 
477  /* -- stress tensor components (only some needed for flux computation) -- */
478 
479  tau_xy[j] = nu1*(r_1*dyVx + dxVy - r_1*vi[VX2]);
480  /*tau_rthita = eta1(1/r dthitaVr + drVthita -1/r Vthita)*/
481  tau_yx[j] = tau_xy[j]; /*tau_thitar*/
482  tau_yy[j] = 2.0*nu1*(r_1*dyVy + r_1*vi[VX1])
483  + (nu2 - (2.0/3.0)*nu1)*div_v;
484  /*tau_thitathita= 2 eta1 (1/r dthitaVthita + 1/r Vr) + (eta2 - 2/3 eta1)divV*/
485  #if DIMENSIONS <= 2
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]));
489  #endif
490  #if DIMENSIONS == 3
491  tau_yz[j] = nu1*(s_1*r_1*dzVy + r_1*dyVz - tan_1*r_1*vi[VX3]);
492  #endif /*tau_thitaphi= eta1 (1/rs dphiVthita + 1/r dthitaVphi -1/r cot Vphi)*/
493  tau_zy[j] = tau_yz[j]; /*tau_phithita*/
494  #if DIMENSIONS <= 2
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;)
501  #endif
502  #if DIMENSIONS == 3
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;
505  #endif /*tau_phiphi = 2eta1(1/rs dphiVphi + 1/r Vr + 1/r cot Vthita) + (eta2 -2/3 eta1)divV */
506 
507  /* -- compute source terms -- */
508 
509  th = grid[JDIR].x; tan_1= 1.0/tan(th[j]);
510 
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; ,
514  ViS[j][MX3] = 0.0; )
515 
516  #endif /* -- end #if GEOMETRY -- */
517 
518  /* -- compute fluxes -- */
519 
520  EXPAND(ViF[j][MX1] = tau_yx[j]; ,
521  ViF[j][MX2] = tau_yy[j]; ,
522  ViF[j][MX3] = tau_yz[j]; )
523 
524  #if EOS != ISOTHERMAL
525  ViF[j][ENG] = EXPAND( vi[VX1]*tau_xy[j],
526  + vi[VX2]*tau_yy[j],
527  + vi[VX3]*tau_zy[j]);
528  #endif
529 
530  dcoeff[j][MX1] = MAX(nu1,nu2);
531  dcoeff[j][MX1] /= vi[RHO];
532  }
533 
534  }else if (g_dir == KDIR){
535 
536  /* ---------------------------------------------------------
537  Loop on X3 direction
538  --------------------------------------------------------- */
539 
540  dx_1 = 1.0/grid[IDIR].dx[i];
541  dy_1 = 1.0/grid[JDIR].dx[j];
542  for (k = beg; k <= end; k++){
543  dz_1 = 1.0/grid[KDIR].dx[k];
544 
545  /* -- compute face- and cell-centered values */
546 
547  VAR_LOOP(nv) {
548  vi[nv] = 0.5*(V[nv][k][j][i] + V[nv][k+1][j][i]);
549  vc[nv] = V[nv][k][j][i];
550  }
551 
552  /* -- compute viscosity (face center) -- */
553 
554  Visc_nu(vi, x1[i], x2[j], x3r[k], &nu1, &nu2);
555 
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;
559 
560  dVxi = D_DX_K(Vx); dVyi = D_DX_K(Vy); dVzi = D_DX_K(Vz);
561  dVxj = D_DY_K(Vx); dVyj = D_DY_K(Vy); dVzj = D_DY_K(Vz);
562  dVxk = D_DZ_K(Vx); dVyk = D_DZ_K(Vy); dVzk = D_DZ_K(Vz);
563 
564  /* ------ CALCULATE DERIVATIVES ---------- */
565 
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;
569 
570  #if GEOMETRY == CARTESIAN
571 
572  /* -- calculate div(v) -- */
573 
574  div_v = dxVx + dyVy + dzVz;
575 
576  /* -- stress tensor components (only some needed for flux computation) -- */
577 
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;
583 
584  /* -- compute source terms -- */
585 
586  EXPAND(ViS[k][MX1] = 0.0; ,
587  ViS[k][MX2] = 0.0; ,
588  ViS[k][MX3] = 0.0; )
589 
590  #elif GEOMETRY == POLAR
591 
592  r = grid[IDIR].x; th = grid[JDIR].x;
593  dr = grid[IDIR].dx[i]; dr_1 = 1.0/dr;
594  r_1 = 1.0/grid[IDIR].x[i];
595 
596  /* -- calculate the div U -- */
597 
598  div_v = r_1*vi[VX1] + dVxi*dr_1 + r_1*dVyj*dy_1 + dVzk*dz_1;
599 
600  /* -- stress tensor components (only some needed for flux computation) -- */
601 
602  tau_xz[k] = nu1*(dzVx + dxVz); /*tau_rz = eta1(drVz + dzVr)*/
603  tau_yy[k] = 2.0*nu1*(r_1*dyVy + r_1*vi[VX1])
604  + (nu2 - (2.0/3.0)*nu1)*div_v;
605  /*tau_phiphi = 2eta1(1/r dphiVphi + 1/r V) + (eta2- 2/3 eta1)divV*/
606  tau_yz[k] = nu1*(r_1*dyVz + dzVy); /*tau_phiz = eta1(dzVphi + 1/r dphiVz)*/
607  tau_zx[k] = tau_xz[k]; /* tau_zr */
608  tau_zy[k] = tau_yz[k]; /* tau_zphi */
609  tau_zz[k] = 2.0*nu1*dzVz + (nu2 - (2.0/3.0)*nu1)*div_v;
610  /*tau_zz = 2eta1(dzVz) + (eta2- 2/3 eta1)divV*/
611 
612  EXPAND(ViS[k][MX1] = 0.0; ,
613  ViS[k][MX2] = 0.0; ,
614  ViS[k][MX3] = 0.0; )
615 
616  #elif GEOMETRY == SPHERICAL
617  r = grid[IDIR].x; th = grid[JDIR].x;
618  dr_1 = 1.0/grid[IDIR].dx[i];
619  r_1 = 1.0/grid[IDIR].x[i]; tan_1= 1.0/tan(th[j]); s_1 = 1.0/sin(th[j]);
620 
621  /* -- calculate div(v) -- */
622 
623  div_v = 2.0*r_1*vi[VX1] + dVxi*dr_1 + r_1*dVyj*dy_1 + r_1*tan_1*vi[VX2]
624  + r_1*s_1*dVzk*dz_1;
625 
626  /* -- stress tensor components (only some needed for flux computation) -- */
627 
628  tau_xz[k] = nu1*(r_1*s_1*dzVx + dxVz - r_1*vi[VX3]);
629  /*tau_rphi = eta1(drVphi + 1/rs dphiVr - 1/r Vphi)*/
630  tau_yz[k] = nu1*(s_1*r_1*dzVy + r_1*dyVz - tan_1*r_1*vi[VX3]);
631  /*tau_thitaphi = eta1(1/rs dphiVthita + 1/r dthitaVphi - 1/r cot Vphi)*/
632  tau_zx[k] = tau_xz[k]; /* tau_phir */
633  tau_zy[k] = tau_yz[k]; /* tau_phithita */
634  tau_zz[k] = 2.0*nu1*(r_1*s_1*dzVz + r_1*vi[VX1]
635  + tan_1*r_1*vi[VX2])
636  + (nu2 - (2.0/3.0)*nu1)*div_v;
637  /*tau_phiphi = 2eta1(1/rs dphiVphi +1/r Vr +1/r cot Vthita) + (eta2 - 2/3 eta1)divV*/
638 
639  /* -- compute source terms -- */
640 
641  EXPAND(ViS[k][MX1] = 0.0; ,
642  ViS[k][MX2] = 0.0; ,
643  ViS[k][MX3] = 0.0; )
644 
645  #endif /* -- end #if GEOMETRY -- */
646 
647  /* -- compute fluxes -- */
648 
649  ViF[k][MX1] = tau_zx[k];
650  ViF[k][MX2] = tau_zy[k];
651  ViF[k][MX3] = tau_zz[k];
652 
653  #if EOS != ISOTHERMAL
654  ViF[k][ENG] = vi[VX1]*tau_xz[k]
655  + vi[VX2]*tau_yz[k]
656  + vi[VX3]*tau_zz[k];
657  #endif
658  dcoeff[k][MX1] = MAX(nu1, nu2);
659  dcoeff[k][MX1] /= vi[RHO];
660  }/*loop*/
661  }/*sweep*/
662 
663 }/*function*/
MX3
#define MX3
Definition: mod_defs.h:22
EOS
#define EOS
Definition: pluto.h:341
MAX
#define MAX(a, b)
Definition: macros.h:101
GRID::end
int end
Global end index for the local array.
Definition: structs.h:116
MX1
#define MX1
Definition: mod_defs.h:20
VX2
#define VX2
Definition: mod_defs.h:29
n
static int n
Definition: analysis.c:3
GRID::xr
double * xr
Definition: structs.h:81
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
SPHERICAL
#define SPHERICAL
Definition: pluto.h:36
D_DX_K
#define D_DX_K(q)
D_DY_K
#define D_DY_K(q)
D_DZ_K
#define D_DZ_K(q)
GRID::dx
double * dx
Definition: structs.h:83
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
POLAR
#define POLAR
Definition: pluto.h:35
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
GRID::beg
int beg
Global start index for the local array.
Definition: structs.h:115
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
D_DY_J
#define D_DY_J(q)
j
int j
Definition: analysis.c:2
MX2
#define MX2
Definition: mod_defs.h:21
D_DZ_J
#define D_DZ_J(q)
D_DZ_I
#define D_DZ_I(q)
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
g_k
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
GEOMETRY
#define GEOMETRY
Definition: definitions_01.h:4
ISOTHERMAL
#define ISOTHERMAL
Definition: pluto.h:49
GRID::r_1
double * r_1
Geometrical factor 1/r.
Definition: structs.h:88
CYLINDRICAL
#define CYLINDRICAL
Definition: pluto.h:34
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
D_EXPAND
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;)
Definition: analysis.c:27
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
D_DX_I
#define D_DX_I(q)
i
int i
Definition: analysis.c:2
Visc_nu
void Visc_nu(double *v, double x1, double x2, double x3, double *nu1, double *nu2)
Definition: visc_nu.c:7
VX3
#define VX3
Definition: mod_defs.h:30
D_DX_J
#define D_DX_J(q)
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
D_DY_I
#define D_DY_I(q)
DIMENSIONS
#define DIMENSIONS
Definition: definitions_01.h:2

Here is the call graph for this function:

Here is the caller graph for this function: