1# -*- coding: utf-8 -*-
2<%namespace module='pyfr.backends.base.makoutil' name='pyfr'/>
3
4% if ndims == 2:
5<%pyfr:macro name='viscous_flux_add' params='uin, grad_uin, fout'>
6    fpdtype_t rho = uin[0], rhou = uin[1], rhov = uin[2], E = uin[3];
7
8    fpdtype_t rcprho = 1.0/rho;
9    fpdtype_t u = rcprho*rhou, v = rcprho*rhov;
10
11    fpdtype_t rho_x = grad_uin[0][0];
12    fpdtype_t rho_y = grad_uin[1][0];
13
14    // Velocity derivatives (rho*grad[u,v])
15    fpdtype_t u_x = grad_uin[0][1] - u*rho_x;
16    fpdtype_t u_y = grad_uin[1][1] - u*rho_y;
17    fpdtype_t v_x = grad_uin[0][2] - v*rho_x;
18    fpdtype_t v_y = grad_uin[1][2] - v*rho_y;
19
20    fpdtype_t E_x = grad_uin[0][3];
21    fpdtype_t E_y = grad_uin[1][3];
22
23% if visc_corr == 'sutherland':
24    // Compute the temperature and viscosity
25    fpdtype_t cpT = ${c['gamma']}*(rcprho*E - 0.5*(u*u + v*v));
26    fpdtype_t Trat = ${1/c['cpTref']}*cpT;
27    fpdtype_t mu_c = ${c['mu']*(c['cpTref'] + c['cpTs'])}*Trat*sqrt(Trat)
28                   / (cpT + ${c['cpTs']});
29% else:
30    fpdtype_t mu_c = ${c['mu']};
31% endif
32
33    // Compute temperature derivatives (c_v*dT/d[x,y])
34    fpdtype_t T_x = rcprho*(E_x - (rcprho*rho_x*E + u*u_x + v*v_x));
35    fpdtype_t T_y = rcprho*(E_y - (rcprho*rho_y*E + u*u_y + v*v_y));
36
37    // Negated stress tensor elements
38    fpdtype_t t_xx = -2*mu_c*rcprho*(u_x - ${1.0/3.0}*(u_x + v_y));
39    fpdtype_t t_yy = -2*mu_c*rcprho*(v_y - ${1.0/3.0}*(u_x + v_y));
40    fpdtype_t t_xy = -mu_c*rcprho*(v_x + u_y);
41
42    fout[0][1] += t_xx;     fout[1][1] += t_xy;
43    fout[0][2] += t_xy;     fout[1][2] += t_yy;
44
45    fout[0][3] += u*t_xx + v*t_xy + -mu_c*${c['gamma']/c['Pr']}*T_x;
46    fout[1][3] += u*t_xy + v*t_yy + -mu_c*${c['gamma']/c['Pr']}*T_y;
47</%pyfr:macro>
48% elif ndims == 3:
49<%pyfr:macro name='viscous_flux_add' params='uin, grad_uin, fout'>
50    fpdtype_t rho  = uin[0];
51    fpdtype_t rhou = uin[1], rhov = uin[2], rhow = uin[3];
52    fpdtype_t E    = uin[4];
53
54    fpdtype_t rcprho = 1.0/rho;
55    fpdtype_t u = rcprho*rhou, v = rcprho*rhov, w = rcprho*rhow;
56
57    fpdtype_t rho_x = grad_uin[0][0];
58    fpdtype_t rho_y = grad_uin[1][0];
59    fpdtype_t rho_z = grad_uin[2][0];
60
61    // Velocity derivatives (rho*grad[u,v,w])
62    fpdtype_t u_x = grad_uin[0][1] - u*rho_x;
63    fpdtype_t u_y = grad_uin[1][1] - u*rho_y;
64    fpdtype_t u_z = grad_uin[2][1] - u*rho_z;
65    fpdtype_t v_x = grad_uin[0][2] - v*rho_x;
66    fpdtype_t v_y = grad_uin[1][2] - v*rho_y;
67    fpdtype_t v_z = grad_uin[2][2] - v*rho_z;
68    fpdtype_t w_x = grad_uin[0][3] - w*rho_x;
69    fpdtype_t w_y = grad_uin[1][3] - w*rho_y;
70    fpdtype_t w_z = grad_uin[2][3] - w*rho_z;
71
72    fpdtype_t E_x = grad_uin[0][4];
73    fpdtype_t E_y = grad_uin[1][4];
74    fpdtype_t E_z = grad_uin[2][4];
75
76% if visc_corr == 'sutherland':
77    // Compute the temperature and viscosity
78    fpdtype_t cpT = ${c['gamma']}*(rcprho*E - 0.5*(u*u + v*v + w*w));
79    fpdtype_t Trat = ${1/c['cpTref']}*cpT;
80    fpdtype_t mu_c = ${c['mu']*(c['cpTref'] + c['cpTs'])}*Trat*sqrt(Trat)
81                   / (cpT + ${c['cpTs']});
82% else:
83    fpdtype_t mu_c = ${c['mu']};
84% endif
85
86    // Compute temperature derivatives (c_v*dT/d[x,y,z])
87    fpdtype_t T_x = rcprho*(E_x - (rcprho*rho_x*E + u*u_x + v*v_x + w*w_x));
88    fpdtype_t T_y = rcprho*(E_y - (rcprho*rho_y*E + u*u_y + v*v_y + w*w_y));
89    fpdtype_t T_z = rcprho*(E_z - (rcprho*rho_z*E + u*u_z + v*v_z + w*w_z));
90
91    // Negated stress tensor elements
92    fpdtype_t t_xx = -2*mu_c*rcprho*(u_x - ${1.0/3.0}*(u_x + v_y + w_z));
93    fpdtype_t t_yy = -2*mu_c*rcprho*(v_y - ${1.0/3.0}*(u_x + v_y + w_z));
94    fpdtype_t t_zz = -2*mu_c*rcprho*(w_z - ${1.0/3.0}*(u_x + v_y + w_z));
95    fpdtype_t t_xy = -mu_c*rcprho*(v_x + u_y);
96    fpdtype_t t_xz = -mu_c*rcprho*(u_z + w_x);
97    fpdtype_t t_yz = -mu_c*rcprho*(w_y + v_z);
98
99    fout[0][1] += t_xx;     fout[1][1] += t_xy;     fout[2][1] += t_xz;
100    fout[0][2] += t_xy;     fout[1][2] += t_yy;     fout[2][2] += t_yz;
101    fout[0][3] += t_xz;     fout[1][3] += t_yz;     fout[2][3] += t_zz;
102
103    fout[0][4] += u*t_xx + v*t_xy + w*t_xz + -mu_c*${c['gamma']/c['Pr']}*T_x;
104    fout[1][4] += u*t_xy + v*t_yy + w*t_yz + -mu_c*${c['gamma']/c['Pr']}*T_y;
105    fout[2][4] += u*t_xz + v*t_yz + w*t_zz + -mu_c*${c['gamma']/c['Pr']}*T_z;
106</%pyfr:macro>
107% endif
108