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