!******************************************************************************** !* Flux Jacobians by automatic differentiation and hand-coded subroutine (2019), !* !* written by Dr. Katate Masatsuka (info[at]cfdbooks.com), !* !* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com). !* !* !* This program runs 2 test problems for verifying hand-coded flux Jacobians: !* !* 1. Compute Jacobians for the Roe numerical flux for the 3D Euler equations. !* -> d(Roe)/dUL and d(Roe)/dUR !* !* 2. Compute Jacobians for the alpha viscous numerical flux for the 3D NS equations. !* -> d(alpha_viscous)/dUL and d(alpha_viscous)/dUR !* !* !* In each test problem, the AD results are compared with the analytical formulas !* (i.e., hand-coded Jacobians). The difference will be printed, which must be zero. !* !* !* Compile and run it. See screen output for further details. !* !* !* NOTE: Hand-coded Jacobians are complicated as you can expect. !* But it can be done and more computationally efficient than AD, and !* so widely used in practical CFD codes. This program shows one way to !* verify the hand-coded Jacobians: by comparing with those computed by !* AD. That's how I typically verify my hand-coded Jacobian subroutines. !* !* !* NOTE: This is written and made available for an educational purpose. !* It is not complete and may not be the best implementation. !* It only attempts to let beginners understand how derivatives can be !* computed by the automatic differentiation method and by hand-coded !* subroutines. If it did help beginners understand the idea of automatic !* differentiation and hand-coded Jacobians, the author would feel happy. !* !* Note: Please let me (Hiro) know at sunmasen@hotmail.com if you find bugs. !* I'll greatly appreciate it and fix the bugs. !* !* !* This file may be updated in future. !* !* Katate Masatsuka, April 2019. http://www.cfdbooks.com !******************************************************************************** ! Must include the following module to enable automatic differentiation. include 'module_ddt5.f90' !required by the inviscid and viscous Jacobians. ! Main program begins here. program test_flux_jacobians implicit none integer , parameter :: p2 = selected_real_kind(15) !for double rpecision. real(p2), dimension(5) :: wL, wR !Primitive variables (rho,u,v,w,p) for testing. real(p2), dimension(5) :: uL, uR !Conservative variables from wL and wR. real(p2), dimension(5,3) :: gradwL, gradwR !Left and right gradients for viscous flux. real(p2), dimension(3) :: njk !Unit face normal vector for testing. real(p2) :: mag_njk !Magnitude of face normal vector for testing. real(p2), dimension(3) :: ejk !Unit edge vector for testing. real(p2) :: mag_ejk !Magnitude of edge vector for testing. write(*,*) write(*,*) !-------------------------------------------------------------------------------- ! Now test the inviscid and viscous flux Jacobians. ! Some random values for the left and right states (primitive variables) wL = (/ 0.993_p2, 0.495_p2, 0.121_p2, -5.62_p2, 0.714_p2 /) wR = (/ 1.011_p2, 0.396_p2, -0.639_p2, 4.31_p2, 0.731_p2 /) ! Compute the corresponding conservative variables for input to flux functions. ! u(1)=rho, u(2)=rho*u, u(3)=rho*v, u(4)=rho*w, u(5)=rho*E=p/(gamma-1)+rho/2*(u^2+v^2+w^2). uL(1) = wL(1) uL(2) = wL(1)*wL(2) uL(3) = wL(1)*wL(3) uL(4) = wL(1)*wL(4) uL(5) = wL(5)/(1.4_p2-1.0_p2) + 0.5_p2*wL(1)*(wL(2)**2+wL(3)**2+wL(4)**2) uR(1) = wR(1) uR(2) = wR(1)*wR(2) uR(3) = wR(1)*wR(3) uR(4) = wR(1)*wR(4) uR(5) = wR(5)/(1.4_p2-1.0_p2) + 0.5_p2*wR(1)*(wR(2)**2+wR(3)**2+wR(4)**2) write(*,*) write(*,*) write(*,*) "*******************************************************************" write(*,*) " Test problem 01: Jacobians of the Roe inviscid flux" write(*,*) "*******************************************************************" write(*,*) ! Define a randomely chosen normal vector. njk(1) = 5.2_p2 njk(2) = -2.7_p2 njk(3) = 0.9_p2 mag_njk = sqrt( njk(1)**2 + njk(2)**2 + njk(3)**2 ) njk = njk / mag_njk !Call the inviscid flux jacobian testing subroutine. call test_inviscid_roe_flux_jacobian(uL,uR,njk,mag_njk) write(*,*) write(*,*) write(*,*) "*******************************************************************" write(*,*) " Test problem 02: Jacobians of the alpha-damping viscous flux" write(*,*) "*******************************************************************" write(*,*) ! Viscous flux requires an edge vector also. Here is an randomly chosen vector: ejk(1) = 0.2_p2 ejk(2) = 2.9_p2 ejk(3) = -0.9_p2 mag_ejk = sqrt( ejk(1)**2 + ejk(2)**2 + ejk(3)**2 ) ejk = ejk / mag_ejk ! Viscous flux requries gradients also, but here they are set to be zero ! since the viscous flux Jacobian is computed with zero gradients. ! Zero gradients mean that the viscous flux has only the damping terms, ! which leads to a compact Jacobian, depending only on the negihbors. gradwL = 0.0_p2 gradwR = 0.0_p2 !Call the viscous flux jacobian testing subroutine. call test_viscous_alpha_flux_jacobian(uL,uR,gradwL,gradwR,njk,ejk,mag_njk,mag_ejk) write(*,*) write(*,*) "*******************************************************************" write(*,*) " Congratulations. All tests have been performed successfully. " write(*,*) "*******************************************************************" write(*,*) contains !******************************************************************************** !******************************************************************************** !******************************************************************************** !******************************************************************************** !* Example of automatic differentiation (AD) and verifying hand-coded Jacobian. !* !* Computation of Jacobians for the Roe flux for the 3D Euler. !* !******************************************************************************** !******************************************************************************** !******************************************************************************** !******************************************************************************** subroutine test_inviscid_roe_flux_jacobian(uL_real, uR_real, njk,mag_njk) use derivative_data_df5 implicit none integer , parameter :: p2 = selected_real_kind(15) real(p2) , dimension(5), intent(in) :: uL_real, uR_real real(p2) , dimension(3), intent(in) :: njk real(p2) , intent(in) :: mag_njk type(derivative_data_type_df5), dimension(5) :: flux, uL, uR real(p2) , dimension(5,5) :: jacobian real(p2) , dimension(5,5) :: jacobianL, jacobianR integer :: i, j write(*,*) write(*,'(a,5es25.15)') " uL = ", uL_real(:) write(*,'(a,5es25.15)') " uR = ", uR_real(:) write(*,'(a,3es25.15)') " njk = ", njk write(*,*) ! Some random values for the left and right states (conservative variables) write(*,*) write(*,*) "--- Check the left Jacobian ---" write(*,*) write(*,*) " Left Jacobian: Jacobian for the Roe flux with respect to the left state (uL)" ! Initialization write(*,*) write(*,*) " Initializing uL, uR..." uL = uL_real uR = uR_real write(*,*) " Initialization done." ! Seeding: set derivatives = 1 for uL. write(*,*) write(*,*) " Seeding uL to compute the Jacobian w.r.t. uL..." call ddt_seed(uL) write(*,*) " Seeding done." ! To compute the derivative with respect to uL: ! Note: 0.7_p2, 0.2_p2 are eigenvalue limiting coefficients (rondamly chosen). call inviscid_roe_n_ddt(uL, uR, njk, 0.7_p2, 0.2_p2, flux) write(*,*) " Flux:" do i = 1, 5 write(*,'(es26.15)') flux(i)%f * mag_njk end do write(*,*) do i = 1, 5 do j = 1, 5 jacobian(i,j) = flux(i)%df(j) * mag_njk end do end do write(*,*) write(*,*) " Left Jacobian computed by AD" write(*,'(5es25.15)') jacobian(1,:) write(*,'(5es25.15)') jacobian(2,:) write(*,'(5es25.15)') jacobian(3,:) write(*,'(5es25.15)') jacobian(4,:) write(*,'(5es25.15)') jacobian(5,:) write(*,*) ! To compute Jacobians by the hand-coded subroutine. ! Note: 0.7_p2, 0.2_p2 are eigenvalue limiting coefficients (rondamly chosen). call inviscid_roe_n_flux_jacobian(uL_real, uR_real, njk, 0.7_p2, 0.2_p2, jacobianL, jacobianR) jacobianL = jacobianL * mag_njk jacobianR = jacobianR * mag_njk write(*,*) write(*,*) " Left Jacobian computed by a hand-coded subroutine" write(*,'(5es25.15)') jacobianL(1,:) write(*,'(5es25.15)') jacobianL(2,:) write(*,'(5es25.15)') jacobianL(3,:) write(*,'(5es25.15)') jacobianL(4,:) write(*,'(5es25.15)') jacobianL(5,:) write(*,*) write(*,*) write(*,*) " Difference between the two" write(*,'(5es25.15)') jacobian(1,:)-jacobianL(1,:) write(*,'(5es25.15)') jacobian(2,:)-jacobianL(2,:) write(*,'(5es25.15)') jacobian(3,:)-jacobianL(3,:) write(*,'(5es25.15)') jacobian(4,:)-jacobianL(4,:) write(*,'(5es25.15)') jacobian(5,:)-jacobianL(5,:) write(*,*) do i = 1, 5 do j = 1, 5 if ( abs(jacobian(i,j)-jacobianL(i,j)) > epsilon(1.0)*max(1.0_p2,abs(jacobian(i,j))) ) then write(*,*) " AD test failed... Stop." stop endif end do end do write(*,*) " Zero, zero, zero! Congratulations!" write(*,*) !-------------------------------------------------------------------- !-------------------------------------------------------------------- !-------------------------------------------------------------------- write(*,*) write(*,*) "--- Check the right Jacobian ---" write(*,*) write(*,*) " Right Jacobian: Jacobian for the Roe flux with respect to the right state (uR)" ! Initialization write(*,*) write(*,*) " Initializing uL, uR..." uL = uL_real uR = uR_real write(*,*) " Initialization done." ! Seeding: set derivatives = 1 for uR. write(*,*) write(*,*) " Seeding uR to compute the Jacobian w.r.t. uR..." call ddt_seed(uR) write(*,*) " Seeding done." ! To compute the derivative wirh respect to uR: ! Note: 0.7_p2, 0.2_p2 are eigenvalue limiting coefficients (rondamly chosen). call inviscid_roe_n_ddt(uL, uR, njk, 0.7_p2,0.2_p2, flux) write(*,*) " Flux:" do i = 1, 5 write(*,'(es26.15)') flux(i)%f * mag_njk end do write(*,*) do i = 1, 5 do j = 1, 5 jacobian(i,j) = flux(i)%df(j) * mag_njk end do end do write(*,*) write(*,*) " Right Jacobian computed by AD" write(*,'(5es25.15)') jacobian(1,:) write(*,'(5es25.15)') jacobian(2,:) write(*,'(5es25.15)') jacobian(3,:) write(*,'(5es25.15)') jacobian(4,:) write(*,'(5es25.15)') jacobian(5,:) write(*,*) ! To compute Jacobians by the hand-coded subroutine. ! Note: 0.7_p2, 0.2_p2 are eigenvalue limiting coefficients (rondamly chosen). call inviscid_roe_n_flux_jacobian(uL_real, uR_real, njk, 0.7_p2,0.2_p2, jacobianL, jacobianR) jacobianL = jacobianL * mag_njk jacobianR = jacobianR * mag_njk write(*,*) write(*,*) " Right Jacobian computed by a hand-coded subroutine" write(*,'(5es25.15)') jacobianR(1,:) write(*,'(5es25.15)') jacobianR(2,:) write(*,'(5es25.15)') jacobianR(3,:) write(*,'(5es25.15)') jacobianR(4,:) write(*,'(5es25.15)') jacobianR(5,:) write(*,*) write(*,*) write(*,*) " Difference between the two" write(*,'(5es25.15)') jacobian(1,:)-jacobianR(1,:) write(*,'(5es25.15)') jacobian(2,:)-jacobianR(2,:) write(*,'(5es25.15)') jacobian(3,:)-jacobianR(3,:) write(*,'(5es25.15)') jacobian(4,:)-jacobianR(4,:) write(*,'(5es25.15)') jacobian(5,:)-jacobianR(5,:) write(*,*) do i = 1, 5 do j = 1, 5 if ( abs(jacobian(i,j)-jacobianR(i,j)) > epsilon(1.0)*max(1.0_p2,abs(jacobian(i,j))) ) then write(*,*) " AD test failed... Stop." stop endif end do end do write(*,*) " Zero, zero, zero! Congratulations!" write(*,*) write(*,*) " End of Roe flux Jacobian test." write(*,*) end subroutine test_inviscid_roe_flux_jacobian !******************************************************************************** !******************************************************************************** !******************************************************************************** !******************************************************************************** !* Example of automatic differentiation (AD) and verifying hand-coded Jacobian. !* !* Computation of Jacobians for the alpha-viscous flux for the 3D Navier-Stokes. !* !******************************************************************************** !******************************************************************************** !******************************************************************************** !******************************************************************************** subroutine test_viscous_alpha_flux_jacobian(uL_real, uR_real, gradwL, gradwR, & njk, ejk, mag_njk, mag_ejk ) use derivative_data_df5 implicit none integer , parameter :: p2 = selected_real_kind(15) real(p2) , dimension(5) , intent(in) :: uL_real, uR_real real(p2) , dimension(5,3), intent(in) :: gradwL, gradwR real(p2) , dimension(3) , intent(in) :: njk real(p2) , dimension(3) , intent(in) :: ejk real(p2) , intent(in) :: mag_njk real(p2) , intent(in) :: mag_ejk type(derivative_data_type_df5), dimension(5) :: flux, uL, uR real(p2) , dimension(5,5) :: jacobian real(p2) , dimension(5,5) :: jacobianL, jacobianR integer :: i, j write(*,*) write(*,'(a,5es25.15)') " uL = ", uL_real(:) write(*,'(a,5es25.15)') " uR = ", uR_real(:) write(*,'(a,3es25.15)') " njk = ", njk write(*,'(a,3es25.15)') " ejk = ", ejk write(*,*) write(*,'(a,5es25.15)') " gradwLx = ", gradwL(:,1) write(*,'(a,5es25.15)') " gradwRx = ", gradwR(:,1) write(*,'(a,5es25.15)') " gradwLy = ", gradwL(:,2) write(*,'(a,5es25.15)') " gradwRy = ", gradwR(:,2) write(*,'(a,5es25.15)') " gradwLz = ", gradwL(:,3) write(*,'(a,5es25.15)') " gradwRz = ", gradwR(:,3) write(*,*) write(*,*) write(*,*) "--- Check the left Jacobian ---" write(*,*) write(*,*) " Left Jacobian: Jacobian for the alpha-damping flux w.r.t. uL" ! Initialization write(*,*) write(*,*) " Initializing uL, uR..." uL = uL_real uR = uR_real write(*,*) " Initialization done." ! Seeding: set derivatives = 1 for uL. write(*,*) write(*,*) " Seeding uL to compute the Jacobian w.r.t. uL..." call ddt_seed(uL) write(*,*) " Seeding done." ! To compute the derivative wirh respect to uL: call viscous_alpha_ddt(uL, uR, gradwL, gradwR, njk,ejk,mag_ejk, flux) do i = 1, 5 do j = 1, 5 jacobian(i,j) = flux(i)%df(j) * mag_njk end do end do write(*,*) write(*,*) " Left Jacobian computed by AD" write(*,'(5es25.15)') jacobian(1,:) write(*,'(5es25.15)') jacobian(2,:) write(*,'(5es25.15)') jacobian(3,:) write(*,'(5es25.15)') jacobian(4,:) write(*,'(5es25.15)') jacobian(5,:) write(*,*) call viscous_alpha_flux_jacobian(uL_real, uR_real, njk,ejk,mag_ejk, jacobianL, jacobianR) jacobianL = jacobianL * mag_njk jacobianR = jacobianR * mag_njk write(*,*) write(*,*) " Left Jacobian computed by a hand-coded subroutine" write(*,'(5es25.15)') jacobianL(1,:) write(*,'(5es25.15)') jacobianL(2,:) write(*,'(5es25.15)') jacobianL(3,:) write(*,'(5es25.15)') jacobianL(4,:) write(*,'(5es25.15)') jacobianL(5,:) write(*,*) write(*,*) write(*,*) " Difference between the two" write(*,'(5es25.15)') jacobian(1,:)-jacobianL(1,:) write(*,'(5es25.15)') jacobian(2,:)-jacobianL(2,:) write(*,'(5es25.15)') jacobian(3,:)-jacobianL(3,:) write(*,'(5es25.15)') jacobian(4,:)-jacobianL(4,:) write(*,'(5es25.15)') jacobian(5,:)-jacobianL(5,:) write(*,*) do i = 1, 5 do j = 1, 5 if ( abs(jacobian(i,j)-jacobianL(i,j)) > epsilon(1.0)*max(1.0_p2,abs(jacobian(i,j))) ) then write(*,*) " AD test failed... Stop." stop endif end do end do write(*,*) " Zero, zero, zero! Congratulations!" write(*,*) !-------------------------------------------------------------------- !-------------------------------------------------------------------- !-------------------------------------------------------------------- write(*,*) write(*,*) "--- Check the right Jacobian ---" write(*,*) write(*,*) " Right Jacobian: Jacobian for the alpha-damping flux w.r.t. uR" ! Initialization write(*,*) write(*,*) " Initializing uL, uR..." uL = uL_real uR = uR_real write(*,*) " Initialization done." ! Seeding: set derivatives = 1 for uR. write(*,*) write(*,*) " Seeding uR to compute the Jacobian w.r.t. uR..." call ddt_seed(uR) write(*,*) " Seeding done." ! To compute the derivative wirh respect to uR: call viscous_alpha_ddt(uL, uR, gradwL, gradwR, njk,ejk,mag_ejk, flux) do i = 1, 5 do j = 1, 5 jacobian(i,j) = flux(i)%df(j) * mag_njk end do end do write(*,*) write(*,*) " Right Jacobian computed by AD" write(*,'(5es25.15)') jacobian(1,:) write(*,'(5es25.15)') jacobian(2,:) write(*,'(5es25.15)') jacobian(3,:) write(*,'(5es25.15)') jacobian(4,:) write(*,'(5es25.15)') jacobian(5,:) write(*,*) call viscous_alpha_flux_jacobian(uL_real, uR_real, njk,ejk,mag_ejk, jacobianL, jacobianR) jacobianL = jacobianL * mag_njk jacobianR = jacobianR * mag_njk write(*,*) write(*,*) " Right Jacobian computed by a hand-coded subroutine" write(*,'(5es25.15)') jacobianR(1,:) write(*,'(5es25.15)') jacobianR(2,:) write(*,'(5es25.15)') jacobianR(3,:) write(*,'(5es25.15)') jacobianR(4,:) write(*,'(5es25.15)') jacobianR(5,:) write(*,*) write(*,*) write(*,*) " Difference between the two" write(*,'(5es25.15)') jacobian(1,:)-jacobianR(1,:) write(*,'(5es25.15)') jacobian(2,:)-jacobianR(2,:) write(*,'(5es25.15)') jacobian(3,:)-jacobianR(3,:) write(*,'(5es25.15)') jacobian(4,:)-jacobianR(4,:) write(*,'(5es25.15)') jacobian(5,:)-jacobianR(5,:) write(*,*) do i = 1, 5 do j = 1, 5 if ( abs(jacobian(i,j)-jacobianR(i,j)) > epsilon(1.0)*max(1.0_p2,abs(jacobian(i,j))) ) then write(*,*) " AD test failed... Stop." stop endif end do end do write(*,*) " Zero, zero, zero! Congratulations!" write(*,*) write(*,*) " End of viscous alpha flux Jacobian test." write(*,*) end subroutine test_viscous_alpha_flux_jacobian !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Below are the flux subroutines written for automatic differentiation ! and the hand-coded Jacobian subroutines (long ones). ! ! inviscid_roe_n_ddt !<- Roe flux with derivative data type. ! inviscid_roe_n_flux_jacobian !<- Hand-coded exact Jacobian for Roe flux. ! ! viscous_alpha_ddt !<- Alpha-damping viscous flux flux with derivative data type. ! viscous_alpha_flux_jacobian !<- Hand-coded exact Jacobian for Alpha-damping viscous flux ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !******************************************************************************** !* -- 3D Roe's Flux Function and Jacobian -- !* !* NOTE: This version does not use any tangent vector. !* See "I do like CFD, VOL.1" about how tangent vectors are eliminated. !* !* This subroutine computes the Roe flux for the Euler equations !* in the direction, njk=[nx,ny,nz]. !* !* P. L. Roe, Approximate Riemann Solvers, Parameter Vectors and Difference !* Schemes, Journal of Computational Physics, 43, pp. 357-372. !* !* Conservative form of the Euler equations: !* !* dU/dt + dF/dx + dG/dy + dH/dz = 0 !* !* This subroutine computes the numerical flux for the flux in the direction, !* njk=[nx,ny,nz]: !* !* Fn = F*nx + G*ny + H*nz = | rho*qn | !* | rho*qn*u + p*nx | !* | rho*qn*v + p*ny | !* | rho*qn*w + p*nz | !* | rho*qn*H | (qn = u*nx + v*ny + w*nz) !* !* The Roe flux is implemented in the following form: !* !* Numerical flux = 1/2 [ Fn(UR) + Fn(UL) - |An|dU ], !* !* where !* !* An = dFn/dU, |An| = R|Lambda|L, dU = UR - UL. !* !* The dissipation term, |An|dU, is actually computed as !* !* sum_{k=1,4} |lambda_k| * (LdU)_k * r_k, !* !* where lambda_k is the k-th eigenvalue, (LdU)_k is the k-th wave strength, !* and r_k is the k-th right-eigenvector evaluated at the Roe-average state. !* !* Note: The 4th component is a combined contribution from two shear waves. !* They are combined to eliminate the tangent vectors. !* So, (LdU)_4 is not really a wave strength, and !* r_4 is not really an eigenvector. !* See "I do like CFD, VOL.1" about how tangent vectors are eliminated. !* !* Note: In the code, the vector of conserative variables are denoted by uc. !* !* ------------------------------------------------------------------------------ !* Input: ucL(1:5) = Left state (rhoL, rhoL*uL, rhoL*vL, rhoL*wR, rhoL*EL) !* ucR(1:5) = Right state (rhoR, rhoL*uR, rhoL*vR, rhoL*wR, rhoL*ER) !* njk(1:3) = unit face normal vector (nx, ny, nz), pointing from Left to Right. !* !* njk !* Face normal ^ o Right data point !* | . !* | . !* |. !* -------x-------- Face !* . Left and right states are !* . 1. Values at data points for 1st-order accuracy !* . 2. Extrapolated values at the face midpoint 'x' !* o Left data point for 2nd/higher-order accuracy. !* !* !* Output: numerical_flux(1:5) = The Roe flux with an entropy fix !* ------------------------------------------------------------------------------ !* !* Note: This function is a ddt-version, which means that each variable carries !* its derivatives, and that the resulting flux "numerical_flux" will have !* its derivatives in "numerical_flux%df". !* !* Note: This subroutine has been prepared for an educational purpose. !* It is not at all efficient. Think about how you can optimize it. !* One way to make it efficient is to reduce the number of local variables, !* by re-using temporary variables as many times as possible. !* !* Note: Please let me know if you find bugs. I'll greatly appreciate it and !* fix the bugs. !* !* Katate Masatsuka, November 2012. http://www.cfdbooks.com !******************************************************************************** subroutine inviscid_roe_n_ddt(ucL, ucR, njk, elimc_nonlinear,elimc_linear, numerical_flux) use derivative_data_df5 implicit none integer , parameter :: p2 = selected_real_kind(15) ! Double precision !Input type(derivative_data_type_df5), dimension(5), intent( in) :: ucL type(derivative_data_type_df5), dimension(5), intent( in) :: ucR real(p2) , dimension(3), intent( in) :: njk real(p2) , intent( in) :: elimc_nonlinear real(p2) , intent( in) :: elimc_linear !Output type(derivative_data_type_df5), dimension(5), intent(out) :: numerical_flux !Some constants real(p2) :: zero = 0.0_p2 real(p2) :: one = 1.0_p2 real(p2) :: two = 2.0_p2 real(p2) :: half = 0.5_p2 real(p2) :: gamma = 1.4_p2 ! Ratio of specific heats !Local variables ! ! L = Left ! R = Right ! No subscript = Roe average real(p2) :: nx, ny, nz ! Normal vector components type(derivative_data_type_df5) :: uL, uR, vL, vR, wL, wR ! Velocity components. type(derivative_data_type_df5) :: rhoL, rhoR, pL, pR ! Primitive variables. type(derivative_data_type_df5) :: qnL, qnR ! Normal velocities type(derivative_data_type_df5) :: aL, aR, HL, HR ! Speed of sound, Total enthalpy type(derivative_data_type_df5), dimension(5) :: fL ! Physical flux evaluated at ucL type(derivative_data_type_df5), dimension(5) :: fR ! Physical flux evaluated at ucR type(derivative_data_type_df5) :: RT ! RT = sqrt(rhoR/rhoL) type(derivative_data_type_df5) :: rho,u,v,w,H,a,qn ! Roe-averages type(derivative_data_type_df5) :: drho, dqn, dp ! Differences in rho, qn, p, e.g., dp=pR-pL type(derivative_data_type_df5), dimension(4) :: LdU ! Wave strengths = L*(UR-UL) type(derivative_data_type_df5) :: du, dv, dw ! Velocity differences type(derivative_data_type_df5), dimension(4) :: ws ! Wave speeds type(derivative_data_type_df5), dimension(4) :: dws ! Width of a parabolic fit for entropy fix type(derivative_data_type_df5), dimension(5,4) :: R ! Right-eigenvector matrix type(derivative_data_type_df5), dimension(5) :: diss ! Dissipation term integer :: k ! Face normal vector (unit vector) nx = njk(1) ny = njk(2) nz = njk(3) !Primitive and other variables. ! Left state rhoL = ucL(1) uL = ucL(2)/ucL(1) vL = ucL(3)/ucL(1) wL = ucL(4)/ucL(1) qnL = uL*nx + vL*ny + wL*nz pL = (gamma-one)*( ucL(5) - half*rhoL*(uL*uL+vL*vL+wL*wL) ) aL = ddt_sqrt(gamma*pL/rhoL) HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL) ! gamma*pL/rhoL/(gamma-one) = gamma*(Rgas*TL)/(gamma-one) = cp*TL ! Right state rhoR = ucR(1) uR = ucR(2)/ucR(1) vR = ucR(3)/ucR(1) wR = ucR(4)/ucR(1) qnR = uR*nx + vR*ny + wR*nz pR = (gamma-one)*( ucR(5) - half*rhoR*(uR*uR+vR*vR+wR*wR) ) aR = ddt_sqrt(gamma*pR/rhoR) HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR) !First compute the Roe-averaged quantities ! NOTE: See http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html for ! the Roe-averaged density. RT = ddt_sqrt(rhoR/rhoL) rho = RT*rhoL !Roe-averaged density u = (uL + RT*uR)/(one + RT) !Roe-averaged x-velocity v = (vL + RT*vR)/(one + RT) !Roe-averaged y-velocity w = (wL + RT*wR)/(one + RT) !Roe-averaged z-velocity H = (HL + RT*HR)/(one + RT) !Roe-averaged total enthalpy a = ddt_sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) ) !Roe-averaged speed of sound qn = u*nx + v*ny + w*nz !Roe-averaged face-normal velocity !Wave Strengths drho = rhoR - rhoL !Density difference dp = pR - pL !Pressure difference dqn = qnR - qnL !Normal velocity difference LdU(1) = (dp - rho*a*dqn )/(two*a*a) !Left-moving acoustic wave strength LdU(2) = drho - dp/(a*a) !Entropy wave strength LdU(3) = (dp + rho*a*dqn )/(two*a*a) !Right-moving acoustic wave strength LdU(4) = rho !Shear wave strength (not really, just a factor) !Absolute values of the wave Speeds ws(1) = ddt_abs(qn-a) !Left-moving acoustic wave ws(2) = ddt_abs(qn) !Entropy wave ws(3) = ddt_abs(qn+a) !Right-moving acoustic wave ws(4) = ddt_abs(qn) !Shear waves !Harten's entropy fix JCP(1983), 49, pp357-393 is applied to avoid vanishing !wave speeds by making a parabolic fit near ws = 0. It is an entropy fix for !nonlinear waves (avoids expnsion shock), but applied here as an eigenvalue !limiting, for the pusrpose of avoiding zero eigenvalues (wave speeds). !Nonlinear fields k = 1 dws(k) = elimc_nonlinear *a if ( ws(k) < dws(k) ) then ws(k) = half * ( ws(k)*ws(k)/dws(k)+dws(k) ) endif k = 3 dws(k) = elimc_nonlinear *a if ( ws(k) < dws(k) ) then ws(k) = half * ( ws(k)*ws(k)/dws(k)+dws(k) ) endif !Linear fields: elimc_linear is better to be small to avoid large dissiaption. k = 2 dws(k) = elimc_linear *a if ( ws(k) < dws(k) ) then ws(k) = half * ( ws(k)*ws(k)/dws(k)+dws(k) ) endif k = 4 dws(k) = elimc_linear *a if ( ws(k) < dws(k) ) then ws(k) = half * ( ws(k)*ws(k)/dws(k)+dws(k) ) endif !Right Eigenvectors !Note: Two shear wave components are combined into one, so that tangent vectors ! are not required. And that's why there are only 4 vectors here. ! See "I do like CFD, VOL.1" about how tangent vectors are eliminated. ! Left-moving acoustic wave R(1,1) = one R(2,1) = u - a*nx R(3,1) = v - a*ny R(4,1) = w - a*nz R(5,1) = H - a*qn ! Entropy wave R(1,2) = one R(2,2) = u R(3,2) = v R(4,2) = w R(5,2) = half*(u*u + v*v + w*w) ! Right-moving acoustic wave R(1,3) = one R(2,3) = u + a*nx R(3,3) = v + a*ny R(4,3) = w + a*nz R(5,3) = H + a*qn ! Two shear wave components combined into one (wave strength incorporated). du = uR - uL dv = vR - vL dw = wR - wL R(1,4) = zero R(2,4) = du - dqn*nx R(3,4) = dv - dqn*ny R(4,4) = dw - dqn*nz R(5,4) = u*du + v*dv + w*dw - qn*dqn !Dissipation Term: |An|(UR-UL) = R|Lambda|L*dU = sum_k of [ ws(k) * R(:,k) * L*dU(k) ] diss(:) = ws(1)*LdU(1)*R(:,1) + ws(2)*LdU(2)*R(:,2) & + ws(3)*LdU(3)*R(:,3) + ws(4)*LdU(4)*R(:,4) !Compute the physical flux: fL = Fn(UL) and fR = Fn(UR) fL(1) = rhoL*qnL fL(2) = rhoL*qnL * uL + pL*nx fL(3) = rhoL*qnL * vL + pL*ny fL(4) = rhoL*qnL * wL + pL*nz fL(5) = rhoL*qnL * HL fR(1) = rhoR*qnR fR(2) = rhoR*qnR * uR + pR*nx fR(3) = rhoR*qnR * vR + pR*ny fR(4) = rhoR*qnR * wR + pR*nz fR(5) = rhoR*qnR * HR ! This is the numerical flux: Roe flux = 1/2 *[ Fn(UL)+Fn(UR) - |An|(UR-UL) ] numerical_flux = half * (fL + fR - diss) !Normal max wave speed in the normal direction. ! wsn = abs(qn) + a end subroutine inviscid_roe_n_ddt !-------------------------------------------------------------------------------- !******************************************************************************** !* -- 3D Roe Flux Jacobian -- !* !* This subroutine computes the left and right Jacobians for the Roe flux !* for the Euler equations in the direction, njk=[nx,ny,nz]. !* !* Conservative form of the Euler equations: !* !* dU/dt + dF/dx + dG/dy + dH/dz = 0 !* !* The normal flux is defined by !* !* Fn = F*nx + G*ny + H*nz = | rho*qn | !* | rho*qn*u + p*nx | !* | rho*qn*v + p*ny | !* | rho*qn*w + p*nz | !* | rho*qn*H | (qn = u*nx + v*ny + w*nz) !* !* The Roe flux is given by !* !* Fn_Roe = 1/2 [ Fn(UR) + Fn(UL) - |An|dU ], !* !* where !* !* An = dFn/dU, |An| = R|Lambda|L, dU = UR - UL. !* !* The dissipation term, |An|dU, is actually computed as !* !* sum_{k=1,4} |lambda_k| * (LdU)_k * r_k, !* !* where lambda_k is the k-th eigenvalue, (LdU)_k is the k-th wave strength, !* and r_k is the k-th right-eigenvector evaluated at the Roe-average state. !* !* Note: The 4th component is a combined contribution from two shear waves. !* They are combined to eliminate the tangent vectors. !* So, (LdU)_4 is not really a wave strength, and !* r_4 is not really an eigenvector. !* See "I do like CFD, VOL.1" about how tangent vectors are eliminated. !* !* Note: In the code, the vector of conserative variables are denoted by uc. !* !* ------------------------------------------------------------------------------ !* Input: ucL(1:5) = Left state (rhoL, rhoL*uL, rhoL*vL, rhoL*wL, rhoL*EL) !* ucR(1:5) = Right state (rhoR, rhoR*uR, rhoR*vR, rhoR*wR, rhoR*ER) !* njk(1:3) = unit face normal vector (nx, ny, nz), pointing from Left to Right. !* !* njk !* Face normal ^ o Right data point !* | . !* | . !* |. !* -------x-------- Face !* . Left and right states are !* . 1. Values at data points for 1st-order accuracy !* . 2. Extrapolated values at the face midpoint 'x' !* o Left data point for 2nd/higher-order accuracy. !* !* !* Output: dFnducL: Derivative of the Roe flux w.r.t. the left state ucL !* dFnducR: Derivative of the Roe flux w.r.t. the left state ucR !* ------------------------------------------------------------------------------ !* !* Katate Masatsuka, November 2012. http://www.cfdbooks.com !******************************************************************************** subroutine inviscid_roe_n_flux_jacobian(ucL, ucR, njk, elimc_nonlinear, & elimc_linear, dFnducL, dFnducR ) implicit none integer , parameter :: p2 = selected_real_kind(15) ! Double precision !Input real(p2), dimension(5) , intent( in) :: ucL ! Left state (conservative variables) real(p2), dimension(5) , intent( in) :: ucR ! Right state (conservative variables) real(p2), dimension(3) , intent( in) :: njk ! Normal vector real(p2) , intent( in) :: elimc_nonlinear real(p2) , intent( in) :: elimc_linear !Output real(p2), dimension(5,5), intent(out) :: dFnducL ! Left Jacobian, d(Fn_Roe)/d(ucL) real(p2), dimension(5,5), intent(out) :: dFnducR ! Right Jacobian, d(Fn_Roe)/d(ucR) !Some constants real(p2) :: zero = 0.0_p2 real(p2) :: one = 1.0_p2 real(p2) :: two = 2.0_p2 real(p2) :: half = 0.5_p2 real(p2) :: gamma = 1.4_p2 ! Ratio of specific heats !Local variables: ! ! L = Left ! R = Right ! No subscript = Roe average real(p2) :: nx, ny, nz ! Normal vector components real(p2) :: uL, uR, vL, vR, wL, wR ! Velocity components. real(p2) :: rhoL, rhoR, pL, pR ! Primitive variables. real(p2) :: qnL, qnR ! Normal velocities real(p2) :: q2L, q2R ! Magnitude of the velocities real(p2) :: aL, aR, HL, HR ! Speed of sound, Total enthalpy real(p2) :: RT ! RT = sqrt(rhoR/rhoL) real(p2) :: rho,u,v,w,H,a,qn ! Roe-averages real(p2) :: q2 ! Squared Roe-averaged speed real(p2) :: drho,dqn,dp ! Differences in rho, qn, p, e.g., dp=pR-pL real(p2) :: du, dv, dw ! Velocity differences real(p2), dimension(4) :: LdU ! Wave strengths = L*(UR-UL) real(p2), dimension(4) :: ws_orig ! Wave speeds real(p2), dimension(4) :: ws ! Wave speeds real(p2), dimension(4) :: dws ! Width of a parabolic fit for entropy fix real(p2), dimension(5,4) :: R ! Right-eigenvector matrix ! Derivatives w.r.t. ucL real(p2), dimension(5) :: drho_ducL ! drho/ducL real(p2), dimension(5) :: dqn_ducL ! dqn/ducL real(p2), dimension(5) :: du_ducL ! du/ducL real(p2), dimension(5) :: dv_ducL ! dv/ducL real(p2), dimension(5) :: dw_ducL ! dw/ducL real(p2), dimension(5) :: dH_ducL ! dH/ducL real(p2), dimension(5) :: da_ducL ! da/ducL real(p2), dimension(5) :: dabs_qn_ducL ! d|qn|/ducL real(p2), dimension(5) :: ddrho_ducL ! d(drho)/ducL real(p2), dimension(5) :: ddp_ducL ! d(dp) /ducL real(p2), dimension(5) :: ddu_ducL ! d(du) /ducL real(p2), dimension(5) :: ddv_ducL ! d(dv) /ducL real(p2), dimension(5) :: ddw_ducL ! d(dw) /ducL real(p2), dimension(5) :: ddqn_ducL ! d(dqn) /ducL real(p2), dimension(5) :: dws1_ducL ! dws1/ducL real(p2), dimension(5) :: dws2_ducL ! dws2/ducL real(p2), dimension(5) :: dws3_ducL ! dws3/ducL real(p2), dimension(5) :: dws4_ducL ! dws4/ducL real(p2), dimension(5) :: ddws1_ducL ! ddws1/ducL real(p2), dimension(5) :: ddws2_ducL ! ddws2/ducL real(p2), dimension(5) :: ddws3_ducL ! ddws3/ducL real(p2), dimension(5) :: ddws4_ducL ! ddws4/ducL real(p2), dimension(5) :: ddws1_ducR ! ddws1/ducL real(p2), dimension(5) :: ddws2_ducR ! ddws2/ducL real(p2), dimension(5) :: ddws3_ducR ! ddws3/ducL real(p2), dimension(5) :: ddws4_ducR ! ddws4/ducL real(p2), dimension(5) :: dLdU1_ducL ! dLdU1/ducL real(p2), dimension(5) :: dLdU2_ducL ! dLdU2/ducL real(p2), dimension(5) :: dLdU3_ducL ! dLdU3/ducL real(p2), dimension(5) :: dLdU4_ducL ! dLdU4/ducL real(p2), dimension(5,5) :: dR1_ducL ! dr1/ducL real(p2), dimension(5,5) :: dR2_ducL ! dr2/ducL real(p2), dimension(5,5) :: dR3_ducL ! dr3/ducL real(p2), dimension(5,5) :: dR4_ducL ! dr4/ducL ! Derivatives w.r.t. ucR real(p2), dimension(5) :: drho_ducR ! drho/ducR real(p2), dimension(5) :: dqn_ducR ! dqn/ducR real(p2), dimension(5) :: du_ducR ! du/ducR real(p2), dimension(5) :: dv_ducR ! dv/ducR real(p2), dimension(5) :: dw_ducR ! dw/ducR real(p2), dimension(5) :: dH_ducR ! dH/ducR real(p2), dimension(5) :: da_ducR ! da/ducR real(p2), dimension(5) :: dabs_qn_ducR ! d|qn|/ducR real(p2), dimension(5) :: ddrho_ducR ! d(drho)/ducR real(p2), dimension(5) :: ddp_ducR ! d(dp) /ducR real(p2), dimension(5) :: ddu_ducR ! d(du) /ducR real(p2), dimension(5) :: ddv_ducR ! d(dv) /ducR real(p2), dimension(5) :: ddw_ducR ! d(dw) /ducR real(p2), dimension(5) :: ddqn_ducR ! d(dqn) /ducR real(p2), dimension(5) :: dws1_ducR ! dws1/ducR real(p2), dimension(5) :: dws2_ducR ! dws2/ducR real(p2), dimension(5) :: dws3_ducR ! dws3/ducR real(p2), dimension(5) :: dws4_ducR ! dws4/ducR real(p2), dimension(5) :: dLdU1_ducR ! dLdU1/ducR real(p2), dimension(5) :: dLdU2_ducR ! dLdU2/ducR real(p2), dimension(5) :: dLdU3_ducR ! dLdU3/ducR real(p2), dimension(5) :: dLdU4_ducR ! dLdU4/ducR real(p2), dimension(5,5) :: dR1_ducR ! dr1/ducR real(p2), dimension(5,5) :: dR2_ducR ! dr2/ducR real(p2), dimension(5,5) :: dR3_ducR ! dr3/ducR real(p2), dimension(5,5) :: dR4_ducR ! dr4/ducR integer :: k ! Face normal vector (unit vector) nx = njk(1) ny = njk(2) nz = njk(3) !Primitive and other variables. ! Left state rhoL = ucL(1) uL = ucL(2)/ucL(1) vL = ucL(3)/ucL(1) wL = ucL(4)/ucL(1) qnL = uL*nx + vL*ny + wL*nz q2L = uL*uL+vL*vL+wL*wL pL = (gamma-one)*( ucL(5) - half*rhoL*(uL*uL+vL*vL+wL*wL) ) aL = sqrt(gamma*pL/rhoL) HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL) ! Right state rhoR = ucR(1) uR = ucR(2)/ucR(1) vR = ucR(3)/ucR(1) wR = ucR(4)/ucR(1) qnR = uR*nx + vR*ny + wR*nz q2R = uR*uR+vR*vR+wR*wR pR = (gamma-one)*( ucR(5) - half*rhoR*(uR*uR+vR*vR+wR*wR) ) aR = sqrt(gamma*pR/rhoR) HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR) ! Compute the Roe-averaged quantities ! NOTE: See http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html for ! the Roe-averaged density. RT = sqrt(rhoR/rhoL) rho = RT*rhoL !Roe-averaged density u = (uL + RT*uR)/(one + RT) !Roe-averaged x-velocity v = (vL + RT*vR)/(one + RT) !Roe-averaged y-velocity w = (wL + RT*wR)/(one + RT) !Roe-averaged z-velocity H = (HL + RT*HR)/(one + RT) !Roe-averaged total enthalpy q2 = u*u + v*v + w*w !Square of Roe-averaged speed a = sqrt( (gamma-one)*(H-half*q2) ) !Roe-averaged speed of sound qn = u*nx + v*ny + w*nz !Roe-averaged face-normal velocity !Wave Strengths drho = rhoR - rhoL !Density difference dp = pR - pL !Pressure difference dqn = qnR - qnL !Normal velocity difference LdU(1) = (dp - rho*a*dqn )/(two*a*a) !Left-moving acoustic wave strength LdU(2) = drho - dp/(a*a) !Entropy wave strength LdU(3) = (dp + rho*a*dqn )/(two*a*a) !Right-moving acoustic wave strength LdU(4) = rho !Shear wave strength (not really, just a factor) !Absolute values of the wave Speeds ws(1) = abs(qn-a) !Left-moving acoustic wave ws(2) = abs(qn) !Entropy wave ws(3) = abs(qn+a) !Right-moving acoustic wave ws(4) = abs(qn) !Shear waves ws_orig = ws !Harten's entropy fix JCP(1983), 49, pp357-393 is applied to avoid vanishing !wave speeds by making a parabolic fit near ws = 0. It is an entropy fix for !nonlinear waves (avoids expnsion shock), but applied here as an eigenvalue !limiting, for the pusrpose of avoiding zero eigenvalues (wave speeds). !Nonlinear fields k = 1 dws(k) = elimc_nonlinear *a if ( ws(k) < dws(k) ) then ws(k) = half * ( ws(k)*ws(k)/dws(k)+dws(k) ) endif k = 3 dws(k) = elimc_nonlinear*a if ( ws(k) < dws(k) ) then ws(k) = half * ( ws(k)*ws(k)/dws(k)+dws(k) ) endif !Linear fields: elimc_linear is better to be small to avoid large dissiaption. k = 2 dws(k) = elimc_linear*a if ( ws(k) < dws(k) ) then ws(k) = half * ( ws(k)*ws(k)/dws(k)+dws(k) ) endif k = 4 dws(k) = elimc_linear*a if ( ws(k) < dws(k) ) then ws(k) = half * ( ws(k)*ws(k)/dws(k)+dws(k) ) endif !Right Eigenvectors !Note: Two shear wave components are combined into one, so that tangent vectors ! are not required. And that's why there are only 4 vectors here. ! See "I do like CFD, VOL.1" about how tangent vectors are eliminated. ! Left-moving acoustic wave R(1,1) = one R(2,1) = u - a*nx R(3,1) = v - a*ny R(4,1) = w - a*nz R(5,1) = H - a*qn ! Entropy wave R(1,2) = one R(2,2) = u R(3,2) = v R(4,2) = w R(5,2) = half*(u*u + v*v + w*w) ! Right-moving acoustic wave R(1,3) = one R(2,3) = u + a*nx R(3,3) = v + a*ny R(4,3) = w + a*nz R(5,3) = H + a*qn ! Two shear wave components combined into one (wave strength incorporated). du = uR - uL dv = vR - vL dw = wR - wL R(1,4) = zero R(2,4) = du - dqn*nx R(3,4) = dv - dqn*ny R(4,4) = dw - dqn*nz R(5,4) = u*du + v*dv + w*dw - qn*dqn ! We are now ready to compute the Jacobians for the Roe flux: ! ! Roe flux function -> Fn_Roe = 0.5*[Fn(ucL)+Fn(ucR)] - 0.5*sum_{k=1,4}|lambda_k|*(LdU)_k*r_k !-------------------------------------------------------------------------------- ! Part 1. Compute dFn_Roe/ducL ! ! dFn_Roe/ducL = d(0.5*Fn(ucL))/duL ! - 0.5*sum_{k=1,4} [d(|lambda_k|)/ducL]*(LdU)_k*r_k ! - 0.5*sum_{k=1,4} |lambda_k|*[d(LdU)_k/ducL]*r_k ! - 0.5*sum_{k=1,4} |lambda_k|*(LdU)_k*[dr_k/ducL] ! ! So, we proceed as follows: ! ! 1.1 Compute d(0.5*Fn(ucL))/duL ! 1.2 Compute various deriavives that will be used in the following steps. ! 1.3 Add the second term, - 0.5*sum_{k=1,4} [d(|lambda_k|)/ducL]*(LdU)_k*r_k ! 1.4 Add the third term, - 0.5*sum_{k=1,4} |lambda_k|*[d(LdU)_k/ducL]*r_k ! 1.5 Add the fourth term, - 0.5*sum_{k=1,4} |lambda_k|*(LdU)_k*[dr_k/ducL] ! !-------------------------------------- ! 1.1 Compute "d(0.5*Fn(ucL))/ducL" ! ! (See "I Do Like CFD, VOL.1", page 55, for the analytical Jacobian, dFn(u)/du) ! 1st column dFnducL(1,1) = zero dFnducL(2,1) = half*(gamma-one)*q2L*nx - uL*qnL dFnducL(3,1) = half*(gamma-one)*q2L*ny - vL*qnL dFnducL(4,1) = half*(gamma-one)*q2L*nz - wL*qnL dFnducL(5,1) = half*(gamma-one)*q2L*qnL - HL*qnL ! 2nd column dFnducL(1,2) = nx dFnducL(2,2) = uL*nx - (gamma-one)*uL*nx + qnL dFnducL(3,2) = vL*nx - (gamma-one)*uL*ny dFnducL(4,2) = wL*nx - (gamma-one)*uL*nz dFnducL(5,2) = HL*nx - (gamma-one)*uL*qnL ! 3rd column dFnducL(1,3) = ny dFnducL(2,3) = uL*ny - (gamma-one)*vL*nx dFnducL(3,3) = vL*ny - (gamma-one)*vL*ny + qnL dFnducL(4,3) = wL*ny - (gamma-one)*vL*nz dFnducL(5,3) = HL*ny - (gamma-one)*vL*qnL ! 4th column dFnducL(1,4) = nz dFnducL(2,4) = uL*nz - (gamma-one)*wL*nx dFnducL(3,4) = vL*nz - (gamma-one)*wL*ny dFnducL(4,4) = wL*nz - (gamma-one)*wL*nz + qnL dFnducL(5,4) = HL*nz - (gamma-one)*wL*qnL ! 5th column dFnducL(1,5) = zero dFnducL(2,5) = (gamma-one)*nx dFnducL(3,5) = (gamma-one)*ny dFnducL(4,5) = (gamma-one)*nz dFnducL(5,5) = gamma*qnL ! Factor 1/2 dFnducL = half * dFnducL !-------------------------------------- ! 1.2 Compute various deriavives that will be used in the following steps. ! (See "I Do Like CFD, VOL.1" for details) ! dqn/ducL dqn_ducL(1) = -half*(qnL+qn) / (rhoL+rho) dqn_ducL(2) = nx / (rhoL+rho) dqn_ducL(3) = ny / (rhoL+rho) dqn_ducL(4) = nz / (rhoL+rho) dqn_ducL(5) = zero ! d(|qn|)/ducL dabs_qn_ducL(1) = -half*sign(one,qn)*(qnL+qn) / (rhoL+rho) dabs_qn_ducL(2) = sign(one,qn)* nx / (rhoL+rho) dabs_qn_ducL(3) = sign(one,qn)* ny / (rhoL+rho) dabs_qn_ducL(4) = sign(one,qn)* nz / (rhoL+rho) dabs_qn_ducL(5) = zero ! da/ducL da_ducL(1) = half*(gamma-one)/a*( half*( uL*u+vL*v+wL*w + q2 ) & + half*(HL-H) - aL**2/(gamma-one) + half*(gamma-two)*q2L )/(rhoL+rho) da_ducL(2) = -half*(gamma-one)*(u+(gamma-one)*uL)/a / (rhoL+rho) da_ducL(3) = -half*(gamma-one)*(v+(gamma-one)*vL)/a / (rhoL+rho) da_ducL(4) = -half*(gamma-one)*(w+(gamma-one)*wL)/a / (rhoL+rho) da_ducL(5) = half*gamma*(gamma-one)/a / (rhoL+rho) ! drho/ducL drho_ducL(1) = half*rho/rhoL drho_ducL(2) = zero drho_ducL(3) = zero drho_ducL(4) = zero drho_ducL(5) = zero ! du/ducL du_ducL(1) = -half*(uL+u) / (rhoL+rho) du_ducL(2) = one / (rhoL+rho) du_ducL(3) = zero du_ducL(4) = zero du_ducL(5) = zero ! dv/ducL dv_ducL(1) = -half*(vL+v) / (rhoL+rho) dv_ducL(2) = zero dv_ducL(3) = one / (rhoL+rho) dv_ducL(4) = zero dv_ducL(5) = zero ! dw/ducL dw_ducL(1) = -half*(wL+w) / (rhoL+rho) dw_ducL(2) = zero dw_ducL(3) = zero dw_ducL(4) = one / (rhoL+rho) dw_ducL(5) = zero ! dH/ducL dH_ducL(1) = ( half*(HL-H) - aL**2/(gamma-one) + half*(gamma-two)*q2L ) / (rhoL+rho) dH_ducL(2) = ( one - gamma )*uL / (rhoL+rho) dH_ducL(3) = ( one - gamma )*vL / (rhoL+rho) dH_ducL(4) = ( one - gamma )*wL / (rhoL+rho) dH_ducL(5) = gamma / (rhoL+rho) ! d(rhoR-rhoL)/ducL = - drhoL/ducL = - (drhoL/dWL)*(dWL/ducL) = -(1,0,0,0,0)*dW/dU ddrho_ducL(1) = - ( one ) ddrho_ducL(2) = - ( zero ) ddrho_ducL(3) = - ( zero ) ddrho_ducL(4) = - ( zero ) ddrho_ducL(5) = - ( zero ) ! d(pR-pL)/ducL = - dpL/ducL = - (dpL/dWL)*(dWL/ducL) = -(0,0,0,0,1)*dW/dU ddp_ducL(1) = - ( half*(gamma-one)*q2L ) ddp_ducL(2) = - ( - (gamma-one)*uL ) ddp_ducL(3) = - ( - (gamma-one)*vL ) ddp_ducL(4) = - ( - (gamma-one)*wL ) ddp_ducL(5) = - ( gamma-one ) ! d(qnR-qnL)/ducL = - dqnL/ducL = - (dqnL/dWL)*(dWL/ducL) = -(0,nx,ny,nz,0)*dW/dU ddqn_ducL(1) = - (-qnL/rhoL) ddqn_ducL(2) = - ( nx/rhoL) ddqn_ducL(3) = - ( ny/rhoL) ddqn_ducL(4) = - ( nz/rhoL) ddqn_ducL(5) = - ( zero ) ! d(uR-uL)/ducL = - duL/ducL = - (duL/dWL)*(dWL/ducL) = -(0,1,0,0,0)*dW/dU ddu_ducL(1) = - ( -uL/rhoL) ddu_ducL(2) = - ( one/rhoL) ddu_ducL(3) = - ( zero ) ddu_ducL(4) = - ( zero ) ddu_ducL(5) = - ( zero ) ! d(vR-vL)/ducL = - dvL/ducL = - (dvL/dWL)*(dWL/ducL) = -(0,0,1,0,0)*dW/dU ddv_ducL(1) = - ( -vL/rhoL) ddv_ducL(2) = - ( zero ) ddv_ducL(3) = - ( one/rhoL) ddv_ducL(4) = - ( zero ) ddv_ducL(5) = - ( zero ) ! d(wR-wL)/ducL = - dwL/ducL = - (dwL/dWL)*(dWL/ducL) = -(0,0,0,1,0)*dW/dU ddw_ducL(1) = - ( -wL/rhoL) ddw_ducL(2) = - ( zero ) ddw_ducL(3) = - ( zero ) ddw_ducL(4) = - ( one/rhoL) ddw_ducL(5) = - ( zero ) ddws1_ducL = elimc_nonlinear*da_ducL ddws3_ducL = elimc_nonlinear*da_ducL ddws2_ducL = elimc_linear *da_ducL ddws4_ducL = elimc_linear *da_ducL !-------------------------------------- ! 1.3 Differentiate the absolute values of the wave speeds, and ! add the second term, - 0.5*sum_{k=1,4} [d(|lambda_k|)/ducL]*(LdU)_k*r_k ! dws1_ducL = d(|qn-a|)/dUcL ! ! Note on entropy fix: !  Let dws1' = half * ( ws(1)*ws(1)/dws(1)+dws(1) ) <- Entropy fix !  Then, dws1'/dUcL = (dws1'/dws1) * dws1_ducL ! = ws(1)/dws(1) * dws1_ducL ! Absolute value if (qn-a > zero) then dws1_ducL = dqn_ducL - da_ducL else dws1_ducL = - ( dqn_ducL - da_ducL ) endif ! Entropy fix/Eigenvalue limiting if ( ws_orig(1) < dws(1) ) then dws1_ducL = ws_orig(1)/dws(1) * dws1_ducL + half * (-ws_orig(1)**2/dws(1)**2 + one)*ddws1_ducL endif dFnducL(:,1) = dFnducL(:,1) - half * dws1_ducL(1)*LdU(1)*R(:,1) dFnducL(:,2) = dFnducL(:,2) - half * dws1_ducL(2)*LdU(1)*R(:,1) dFnducL(:,3) = dFnducL(:,3) - half * dws1_ducL(3)*LdU(1)*R(:,1) dFnducL(:,4) = dFnducL(:,4) - half * dws1_ducL(4)*LdU(1)*R(:,1) dFnducL(:,5) = dFnducL(:,5) - half * dws1_ducL(5)*LdU(1)*R(:,1) ! dws2_ducL = d(|qn|)/dUcL dws2_ducL = dabs_qn_ducL ! Entropy fix/Eigenvalue limiting if ( ws_orig(2) < dws(2) ) then dws2_ducL = ws_orig(2)/dws(2) * dws2_ducL + half * (-ws_orig(2)**2/dws(2)**2 + one)*ddws2_ducL endif dFnducL(:,1) = dFnducL(:,1) - half * dws2_ducL(1)*LdU(2)*R(:,2) dFnducL(:,2) = dFnducL(:,2) - half * dws2_ducL(2)*LdU(2)*R(:,2) dFnducL(:,3) = dFnducL(:,3) - half * dws2_ducL(3)*LdU(2)*R(:,2) dFnducL(:,4) = dFnducL(:,4) - half * dws2_ducL(4)*LdU(2)*R(:,2) dFnducL(:,5) = dFnducL(:,5) - half * dws2_ducL(5)*LdU(2)*R(:,2) ! dws3_ducL = d(|qn+a|)/dUcL ! ! Note on entropy fix: !  Let dws3' = half * ( ws(3)*ws(3)/dws(3)+dws(3) ) <- Entropy fix !  Then, dws3'/dUcL = (dws3'/dws3) * dws3_ducL ! = ws(3)/dws(3) * dws3_ducL ! Absolute value if (qn+a > zero) then dws3_ducL = dqn_ducL + da_ducL else dws3_ducL = - ( dqn_ducL + da_ducL ) endif ! Entropy fix/Eigenvalue limiting if ( ws_orig(3) < dws(3) ) then dws3_ducL = ws_orig(3)/dws(3) * dws3_ducL + half * (-ws_orig(3)**2/dws(3)**2 + one)*ddws3_ducL endif dFnducL(:,1) = dFnducL(:,1) - half * dws3_ducL(1)*LdU(3)*R(:,3) dFnducL(:,2) = dFnducL(:,2) - half * dws3_ducL(2)*LdU(3)*R(:,3) dFnducL(:,3) = dFnducL(:,3) - half * dws3_ducL(3)*LdU(3)*R(:,3) dFnducL(:,4) = dFnducL(:,4) - half * dws3_ducL(4)*LdU(3)*R(:,3) dFnducL(:,5) = dFnducL(:,5) - half * dws3_ducL(5)*LdU(3)*R(:,3) ! dws4_ducL = d(|qn|)/dUcL = dws1_ducL dws4_ducL = dabs_qn_ducL ! Entropy fix/Eigenvalue limiting if ( ws_orig(4) < dws(4) ) then dws4_ducL = ws_orig(4)/dws(4) * dws4_ducL + half * (-ws_orig(4)**2/dws(4)**2 + one)*ddws4_ducL endif dFnducL(:,1) = dFnducL(:,1) - half * dws4_ducL(1)*LdU(4)*R(:,4) dFnducL(:,2) = dFnducL(:,2) - half * dws4_ducL(2)*LdU(4)*R(:,4) dFnducL(:,3) = dFnducL(:,3) - half * dws4_ducL(3)*LdU(4)*R(:,4) dFnducL(:,4) = dFnducL(:,4) - half * dws4_ducL(4)*LdU(4)*R(:,4) dFnducL(:,5) = dFnducL(:,5) - half * dws4_ducL(5)*LdU(4)*R(:,4) !-------------------------------------- ! 1.4 Differentiate the wave strength, and ! add the third term, - 0.5*sum_{k=1,4} |lambda_k|*[d(LdU)_k/ducL]*r_k. ! dLdU1_ducL = d( dp/(2a^2) - rho/(2a)*dqn )/dUcL ! ! = dp*[d(1/(2a^2))/dUcL] - dqn*[d(rho/(2a))/dUcL] ! + [d(dp)/dUcL]/(2a^2) - [d(dqn)/dUcL]*rho/(2a) ! ! = -a^(-3)*dp*[da/dUcL] - dqn*( [drho/dUcL]/(2a) - rho/(2*a^2)*[da/dUcL] ) ! + [d(dp)/dUcL]/(2a^2) - [d(dqn)/dUcL]*rho/(2a) ! ! = ( -2*dp + rho*a*dqn )/(2a^3)*[da/dUcL] - dqn*[drho/dUcL]/(2a) ! + [d(dp)/dUcL]/(2a^2) - [d(dqn)/dUcL]*rho/(2a) dLdU1_ducL = half*(-two*dp+rho*a*dqn )/a**3 * (da_ducL) & - half*dqn/a * (drho_ducL) & + half*(ddp_ducL)/a**2 & - half*rho*(ddqn_ducL)/a dFnducL(:,1) = dFnducL(:,1) - half * ws(1)*dLdU1_ducL(1)*R(:,1) dFnducL(:,2) = dFnducL(:,2) - half * ws(1)*dLdU1_ducL(2)*R(:,1) dFnducL(:,3) = dFnducL(:,3) - half * ws(1)*dLdU1_ducL(3)*R(:,1) dFnducL(:,4) = dFnducL(:,4) - half * ws(1)*dLdU1_ducL(4)*R(:,1) dFnducL(:,5) = dFnducL(:,5) - half * ws(1)*dLdU1_ducL(5)*R(:,1) ! dLdU2_ducL = d( drho - dp/a^2 )/dUcL ! = [d(drho)/dUcL] - [d(dp)/dUcL]/a^2 + 2*dp/a^3*[da/dUcL] dLdU2_ducL = (ddrho_ducL) - (ddp_ducL)/a**2 + two*dp/a**3*(da_ducL) dFnducL(:,1) = dFnducL(:,1) - half * ws(2)*dLdU2_ducL(1)*R(:,2) dFnducL(:,2) = dFnducL(:,2) - half * ws(2)*dLdU2_ducL(2)*R(:,2) dFnducL(:,3) = dFnducL(:,3) - half * ws(2)*dLdU2_ducL(3)*R(:,2) dFnducL(:,4) = dFnducL(:,4) - half * ws(2)*dLdU2_ducL(4)*R(:,2) dFnducL(:,5) = dFnducL(:,5) - half * ws(2)*dLdU2_ducL(5)*R(:,2) ! dLdU3_ducL = d( dp/(2a^2) + rho/(2a)*dqn )/dUcL ! ! = dp*[d(1/(2a^2))/dUcL] + dqn*[d(rho/(2a))/dUcL] ! + [d(dp)/dUcL]/(2a^2) + [d(dqn)/dUcL]*rho/(2a) ! ! = -a^(-3)*dp*[da/dUcL] + dqn*( [drho/dUcL]/(2a) - rho/(2*a^2)*[da/dUcL] ) ! + [d(dp)/dUcL]/(2a^2) + [d(dqn)/dUcL]*rho/(2a) ! ! = ( -2*dp - rho*a*dqn )/(2a^3)*[da/dUcL] + dqn*[drho/dUcL]/(2a) ! + [d(dp)/dUcL]/(2a^2) + [d(dqn)/dUcL]*rho/(2a) dLdU3_ducL = half*(-two*dp-rho*a*dqn )/a**3 * (da_ducL) & + half*dqn/a * (drho_ducL) & + half*(ddp_ducL)/a**2 & + half*rho*(ddqn_ducL)/a dFnducL(:,1) = dFnducL(:,1) - half * ws(3)*dLdU3_ducL(1)*R(:,3) dFnducL(:,2) = dFnducL(:,2) - half * ws(3)*dLdU3_ducL(2)*R(:,3) dFnducL(:,3) = dFnducL(:,3) - half * ws(3)*dLdU3_ducL(3)*R(:,3) dFnducL(:,4) = dFnducL(:,4) - half * ws(3)*dLdU3_ducL(4)*R(:,3) dFnducL(:,5) = dFnducL(:,5) - half * ws(3)*dLdU3_ducL(5)*R(:,3) ! dLdU4_ducL = d(rho)/dUcL dLdU4_ducL = drho_ducL dFnducL(:,1) = dFnducL(:,1) - half * ws(4)*dLdU4_ducL(1)*R(:,4) dFnducL(:,2) = dFnducL(:,2) - half * ws(4)*dLdU4_ducL(2)*R(:,4) dFnducL(:,3) = dFnducL(:,3) - half * ws(4)*dLdU4_ducL(3)*R(:,4) dFnducL(:,4) = dFnducL(:,4) - half * ws(4)*dLdU4_ducL(4)*R(:,4) dFnducL(:,5) = dFnducL(:,5) - half * ws(4)*dLdU4_ducL(5)*R(:,4) !-------------------------------------- ! 1.5 Differentiate the right-eigenvectors, and ! add the fourth term, - 0.5*sum_{k=1,4} |lambda_k|*(LdU)_k*[dr_k/ducL] ! dR1_ducL = dR(:,1)/dUcL ! ! Left-moving acoustic wave ! ! Eigenvector -> Differentiated ! ! R(1,1) = one -> 0 ! R(2,1) = u - a*nx -> du/dUcL - da/dUcL*nx ! R(3,1) = v - a*ny -> dv/dUcL - da/dUcL*ny ! R(4,1) = w - a*nz -> dw/dUcL - da/dUcL*nz ! R(5,1) = H - a*qn -> dH/dUcL - da/dUcL*qn - dqn/dUcL*a dR1_ducL(1,:) = zero dR1_ducL(2,:) = (du_ducL) - (da_ducL) * nx dR1_ducL(3,:) = (dv_ducL) - (da_ducL) * ny dR1_ducL(4,:) = (dw_ducL) - (da_ducL) * nz dR1_ducL(5,:) = (dH_ducL) - (da_ducL) * qn - (dqn_ducL) * a dFnducL = dFnducL - half * ws(1)*LdU(1)*dR1_ducL ! dR2_ducL = dR(:,2)/dUcL ! ! Entropy wave ! ! Eigenvector -> Differentiated ! ! R(1,2) = one -> 0 ! R(2,2) = u -> du/dUcL ! R(3,2) = v -> dv/dUcL ! R(4,2) = w -> dw/dUcL ! R(5,2) = half*(u*u+v*v+w*w) -> u*du/dUcL + v*dv/dUcL + w*dw/dUcL dR2_ducL(1,:) = zero dR2_ducL(2,:) = (du_ducL) dR2_ducL(3,:) = (dv_ducL) dR2_ducL(4,:) = (dw_ducL) dR2_ducL(5,:) = u*(du_ducL) + v*(dv_ducL) + w*(dw_ducL) dFnducL = dFnducL - half * ws(2)*LdU(2)*dR2_ducL ! dR3_ducL = dR(:,3)/dUcL ! ! Right-moving acoustic wave ! ! Eigenvector -> Differentiated ! ! R(1,3) = one -> 0 ! R(2,3) = u + a*nx -> du/dUcL + da/dUcL*nx ! R(3,3) = v + a*ny -> dv/dUcL + da/dUcL*ny ! R(4,3) = w + a*nz -> dw/dUcL + da/dUcL*nz ! R(5,3) = H + a*qn -> dH/dUcL + da/dUcL*qn + dqn/dUcL*a dR3_ducL(1,:) = zero dR3_ducL(2,:) = (du_ducL) + (da_ducL) * nx dR3_ducL(3,:) = (dv_ducL) + (da_ducL) * ny dR3_ducL(4,:) = (dw_ducL) + (da_ducL) * nz dR3_ducL(5,:) = (dH_ducL) + (da_ducL) * qn + (dqn_ducL) * a dFnducL = dFnducL - half * ws(3)*LdU(3)*dR3_ducL ! dR4_ducL = dR(:,4)/dUcL ! Two shear wave components combined into one (wave strength incorporated). ! So, it is not really an eigenvector. ! ! Combined vector -> Differentiated ! ! R(1,4) = zero -> 0 ! R(2,4) = du - dqn*nx -> d(du)/dUcL - d(dqn)/dUcL*nx ! R(3,4) = dv - dqn*ny -> d(dv)/dUcL - d(dqn)/dUcL*ny ! R(4,4) = dw - dqn*nz -> d(dw)/dUcL - d(dqn)/dUcL*nz ! R(5,4) = u*du+v*dv+w*dw-qn*dqn -> du/dUcL*du + d(du)/dUcL*u ! + dv/dUcL*dv + d(dv)/dUcL*v ! + dw/dUcL*dw + d(dw)/dUcL*w ! - d(qn)/dUcL*dqn - d(dqn)/dUcL*qn dR4_ducL(1,:) = zero dR4_ducL(2,:) = (ddu_ducL) - (ddqn_ducL)*nx dR4_ducL(3,:) = (ddv_ducL) - (ddqn_ducL)*ny dR4_ducL(4,:) = (ddw_ducL) - (ddqn_ducL)*nz dR4_ducL(5,:) = du*( du_ducL) + dv*( dv_ducL) + dw*( dw_ducL) - dqn*( dqn_ducL) & + u*(ddu_ducL) + v*(ddv_ducL) + w*(ddw_ducL) - qn*(ddqn_ducL) dFnducL = dFnducL - half * ws(4)*LdU(4)*dR4_ducL !-------------------------------------------------------------------------------- ! Part 2. Compute dFn_Roe/ducR ! ! dFn_Roe/ducR = d(0.5*Fn(ucR))/duR ! - 0.5*sum_{k=1,4} [d(|lambda_k|)/ducR]*(LdU)_k*r_k ! - 0.5*sum_{k=1,4} |lambda_k|*[d(LdU)_k/ducR]*r_k ! - 0.5*sum_{k=1,4} |lambda_k|*(LdU)_k*[dr_k/ducR] ! ! So, we proceed as follows: ! ! 1.1 Compute d(0.5*Fn(ucR))/duR ! 1.2 Compute various deriavives that will be used in the following steps. ! 1.3 Add the second term, - 0.5*sum_{k=1,4} [d(|lambda_k|)/ducR]*(LdU)_k*r_k ! 1.4 Add the third term, - 0.5*sum_{k=1,4} |lambda_k|*[d(LdU)_k/ducR]*r_k ! 1.5 Add the fourth term, - 0.5*sum_{k=1,4} |lambda_k|*(LdU)_k*[dr_k/ducR] ! !-------------------------------------- ! 2.1 Compute "d(0.5*Fn(ucR))/ducR" ! ! (See "I Do Like CFD, VOL.1", page 55, for the analytical Jacobian, dFn(u)/du) ! 1st column dFnducR(1,1) = zero dFnducR(2,1) = half*(gamma-one)*q2R*nx - uR*qnR dFnducR(3,1) = half*(gamma-one)*q2R*ny - vR*qnR dFnducR(4,1) = half*(gamma-one)*q2R*nz - wR*qnR dFnducR(5,1) = half*(gamma-one)*q2R*qnR - HR*qnR ! 2nd column dFnducR(1,2) = nx dFnducR(2,2) = uR*nx - (gamma-one)*uR*nx + qnR dFnducR(3,2) = vR*nx - (gamma-one)*uR*ny dFnducR(4,2) = wR*nx - (gamma-one)*uR*nz dFnducR(5,2) = HR*nx - (gamma-one)*uR*qnR ! 3rd column dFnducR(1,3) = ny dFnducR(2,3) = uR*ny - (gamma-one)*vR*nx dFnducR(3,3) = vR*ny - (gamma-one)*vR*ny + qnR dFnducR(4,3) = wR*ny - (gamma-one)*vR*nz dFnducR(5,3) = HR*ny - (gamma-one)*vR*qnR ! 4th column dFnducR(1,4) = nz dFnducR(2,4) = uR*nz - (gamma-one)*wR*nx dFnducR(3,4) = vR*nz - (gamma-one)*wR*ny dFnducR(4,4) = wR*nz - (gamma-one)*wR*nz + qnR dFnducR(5,4) = HR*nz - (gamma-one)*wR*qnR ! 5th column dFnducR(1,5) = zero dFnducR(2,5) = (gamma-one)*nx dFnducR(3,5) = (gamma-one)*ny dFnducR(4,5) = (gamma-one)*nz dFnducR(5,5) = gamma*qnR ! Factor 1/2 dFnducR = half * dFnducR !-------------------------------------- ! 2.2 Compute various deriavives that will be used in the following steps. ! (See "I Do Like CFD, VOL.1" for details) ! dqn/ducR dqn_ducR(1) = -half*(qnR+qn) / (rhoR+rho) dqn_ducR(2) = nx / (rhoR+rho) dqn_ducR(3) = ny / (rhoR+rho) dqn_ducR(4) = nz / (rhoR+rho) dqn_ducR(5) = zero ! d(|qn|)/ducR dabs_qn_ducR(1) = -half*sign(one,qn)*(qnR+qn) / (rhoR+rho) dabs_qn_ducR(2) = sign(one,qn)* nx / (rhoR+rho) dabs_qn_ducR(3) = sign(one,qn)* ny / (rhoR+rho) dabs_qn_ducR(4) = sign(one,qn)* nz / (rhoR+rho) dabs_qn_ducR(5) = zero ! da/ducR da_ducR(1) = half*(gamma-one)/a*( half*( uR*u+vR*v+wR*w + q2 ) & + half*(HR-H) - aR**2/(gamma-one) + half*(gamma-two)*q2R )/(rhoR+rho) da_ducR(2) = -half*(gamma-one)*(u+(gamma-one)*uR)/a / (rhoR+rho) da_ducR(3) = -half*(gamma-one)*(v+(gamma-one)*vR)/a / (rhoR+rho) da_ducR(4) = -half*(gamma-one)*(w+(gamma-one)*wR)/a / (rhoR+rho) da_ducR(5) = half*gamma*(gamma-one)/a / (rhoR+rho) ! drho/ducR drho_ducR(1) = half*rho/rhoR drho_ducR(2) = zero drho_ducR(3) = zero drho_ducR(4) = zero drho_ducR(5) = zero ! du/ducR du_ducR(1) = -half*(uR+u) / (rhoR+rho) du_ducR(2) = one / (rhoR+rho) du_ducR(3) = zero du_ducR(4) = zero du_ducR(5) = zero ! dv/ducR dv_ducR(1) = -half*(vR+v) / (rhoR+rho) dv_ducR(2) = zero dv_ducR(3) = one / (rhoR+rho) dv_ducR(4) = zero dv_ducR(5) = zero ! dw/ducR dw_ducR(1) = -half*(wR+w) / (rhoR+rho) dw_ducR(2) = zero dw_ducR(3) = zero dw_ducR(4) = one / (rhoR+rho) dw_ducR(5) = zero ! dH/ducR dH_ducR(1) = ( half*(HR-H) - aR**2/(gamma-one) + half*(gamma-two)*q2R ) / (rhoR+rho) dH_ducR(2) = ( one - gamma )*uR / (rhoR+rho) dH_ducR(3) = ( one - gamma )*vR / (rhoR+rho) dH_ducR(4) = ( one - gamma )*wR / (rhoR+rho) dH_ducR(5) = gamma / (rhoR+rho) ! d(rhoR-rhoR)/ducR = drhoR/ducR = (drhoR/dWR)*(dWR/ducR) = (1,0,0,0,0)*dW/dU ddrho_ducR(1) = ( one ) ddrho_ducR(2) = ( zero ) ddrho_ducR(3) = ( zero ) ddrho_ducR(4) = ( zero ) ddrho_ducR(5) = ( zero ) ! d(pR-pR)/ducR = dpR/ducR = (dpR/dWR)*(dWR/ducR) = (0,0,0,0,1)*dW/dU ddp_ducR(1) = ( half*(gamma-one)*q2R ) ddp_ducR(2) = ( - (gamma-one)*uR ) ddp_ducR(3) = ( - (gamma-one)*vR ) ddp_ducR(4) = ( - (gamma-one)*wR ) ddp_ducR(5) = ( gamma-one ) ! d(qnR-qnR)/ducR = dqnR/ducR = (dqnR/dWR)*(dWR/ducR) = (0,nx,ny,nz,0)*dW/dU ddqn_ducR(1) = (-qnR/rhoR) ddqn_ducR(2) = ( nx/rhoR) ddqn_ducR(3) = ( ny/rhoR) ddqn_ducR(4) = ( nz/rhoR) ddqn_ducR(5) = ( zero ) ! d(uR-uR)/ducR = duR/ducR = (duR/dWR)*(dWR/ducR) = -(0,1,0,0,0)*dW/dU ddu_ducR(1) = ( -uR/rhoR) ddu_ducR(2) = ( one/rhoR) ddu_ducR(3) = ( zero ) ddu_ducR(4) = ( zero ) ddu_ducR(5) = ( zero ) ! d(vR-vR)/ducR = dvR/ducR = (dvR/dWR)*(dWR/ducR) = (0,0,1,0,0)*dW/dU ddv_ducR(1) = ( -vR/rhoR) ddv_ducR(2) = ( zero ) ddv_ducR(3) = ( one/rhoR) ddv_ducR(4) = ( zero ) ddv_ducR(5) = ( zero ) ! d(wR-wR)/ducR = dwR/ducR = (dwR/dWR)*(dWR/ducR) = (0,0,0,1,0)*dW/dU ddw_ducR(1) = ( -wR/rhoR) ddw_ducR(2) = ( zero ) ddw_ducR(3) = ( zero ) ddw_ducR(4) = ( one/rhoR) ddw_ducR(5) = ( zero ) ddws1_ducR = elimc_nonlinear*da_ducR ddws3_ducR = elimc_nonlinear*da_ducR ddws2_ducR = elimc_linear*da_ducR ddws4_ducR = elimc_linear*da_ducR !-------------------------------------- ! 2.3 Differentiate the absolute values of the wave speeds, and ! add the second term, - 0.5*sum_{k=1,4} [d(|lambda_k|)/ducR]*(LdU)_k*r_k ! dws1_ducR = d(|qn-a|)/dUcR ! ! Note on entropy fix: !  Let dws1' = half * ( ws(1)*ws(1)/dws(1)+dws(1) ) <- Entropy fix !  Then, dws1'/dUcR = (dws1'/dws1) * dws1_ducR + ( -ws(1)*ws(1)/dws(1)**2*ddws(1) + ddws(1) ) ! = ws(1)/dws(1) * dws1_ducR ! Absolute value if (qn-a > zero) then dws1_ducR = dqn_ducR - da_ducR else dws1_ducR = - ( dqn_ducR - da_ducR ) endif ! Entropy fix/Eigenvalue limiting if ( ws_orig(1) < dws(1) ) then dws1_ducR = ws_orig(1)/dws(1) * dws1_ducR + half * (-ws_orig(1)**2/dws(1)**2 + one)*ddws1_ducR endif dFnducR(:,1) = dFnducR(:,1) - half * dws1_ducR(1)*LdU(1)*R(:,1) dFnducR(:,2) = dFnducR(:,2) - half * dws1_ducR(2)*LdU(1)*R(:,1) dFnducR(:,3) = dFnducR(:,3) - half * dws1_ducR(3)*LdU(1)*R(:,1) dFnducR(:,4) = dFnducR(:,4) - half * dws1_ducR(4)*LdU(1)*R(:,1) dFnducR(:,5) = dFnducR(:,5) - half * dws1_ducR(5)*LdU(1)*R(:,1) ! dws2_ducR = d(|qn|)/dUcR dws2_ducR = dabs_qn_ducR ! Entropy fix/Eigenvalue limiting if ( ws_orig(2) < dws(2) ) then dws2_ducR = ws_orig(2)/dws(2) * dws2_ducR + half * (-ws_orig(2)**2/dws(2)**2 + one)*ddws2_ducR endif dFnducR(:,1) = dFnducR(:,1) - half * dws2_ducR(1)*LdU(2)*R(:,2) dFnducR(:,2) = dFnducR(:,2) - half * dws2_ducR(2)*LdU(2)*R(:,2) dFnducR(:,3) = dFnducR(:,3) - half * dws2_ducR(3)*LdU(2)*R(:,2) dFnducR(:,4) = dFnducR(:,4) - half * dws2_ducR(4)*LdU(2)*R(:,2) dFnducR(:,5) = dFnducR(:,5) - half * dws2_ducR(5)*LdU(2)*R(:,2) ! dws3_ducR = d(|qn+a|)/dUcR ! ! Note on entropy fix: !  Let dws3' = half * ( ws(3)*ws(3)/dws(3)+dws(3) ) <- Entropy fix !  Then, dws3'/dUcR = (dws3'/dws3) * dws3_ducR ! = ws(3)/dws(3) * dws3_ducR ! Absolute value if (qn+a > zero) then dws3_ducR = dqn_ducR + da_ducR else dws3_ducR = - ( dqn_ducR + da_ducR ) endif ! Entropy fix/Eigenvalue limiting if ( ws_orig(3) < dws(3) ) then dws3_ducR = ws_orig(3)/dws(3) * dws3_ducR + half * (-ws_orig(3)**2/dws(3)**2 + one)*ddws3_ducR endif dFnducR(:,1) = dFnducR(:,1) - half * dws3_ducR(1)*LdU(3)*R(:,3) dFnducR(:,2) = dFnducR(:,2) - half * dws3_ducR(2)*LdU(3)*R(:,3) dFnducR(:,3) = dFnducR(:,3) - half * dws3_ducR(3)*LdU(3)*R(:,3) dFnducR(:,4) = dFnducR(:,4) - half * dws3_ducR(4)*LdU(3)*R(:,3) dFnducR(:,5) = dFnducR(:,5) - half * dws3_ducR(5)*LdU(3)*R(:,3) ! dws4_ducR = d(|qn|)/dUcR = dws1_ducR dws4_ducR = dabs_qn_ducR ! Entropy fix/Eigenvalue limiting if ( ws_orig(4) < dws(4) ) then dws4_ducR = ws_orig(4)/dws(4) * dws4_ducR + half * (-ws_orig(4)**2/dws(4)**2 + one)*ddws4_ducR endif dFnducR(:,1) = dFnducR(:,1) - half * dws4_ducR(1)*LdU(4)*R(:,4) dFnducR(:,2) = dFnducR(:,2) - half * dws4_ducR(2)*LdU(4)*R(:,4) dFnducR(:,3) = dFnducR(:,3) - half * dws4_ducR(3)*LdU(4)*R(:,4) dFnducR(:,4) = dFnducR(:,4) - half * dws4_ducR(4)*LdU(4)*R(:,4) dFnducR(:,5) = dFnducR(:,5) - half * dws4_ducR(5)*LdU(4)*R(:,4) !-------------------------------------- ! 2.4 Differentiate the wave strength, and ! add the third term, - 0.5*sum_{k=1,4} |lambda_k|*[d(LdU)_k/ducR]*r_k. ! dLdU1_ducR = d( dp/(2a^2) - rho/(2a)*dqn )/dUcR ! ! = dp*[d(1/(2a^2))/dUcR] - dqn*[d(rho/(2a))/dUcR] ! + [d(dp)/dUcR]/(2a^2) - [d(dqn)/dUcR]*rho/(2a) ! ! = -a^(-3)*dp*[da/dUcR] - dqn*( [drho/dUcR]/(2a) - rho/(2*a^2)*[da/dUcR] ) ! + [d(dp)/dUcR]/(2a^2) - [d(dqn)/dUcR]*rho/(2a) ! ! = ( -2*dp + rho*a*dqn )/(2a^3)*[da/dUcR] - dqn*[drho/dUcR]/(2a) ! + [d(dp)/dUcR]/(2a^2) - [d(dqn)/dUcR]*rho/(2a) dLdU1_ducR = half*(-two*dp+rho*a*dqn )/a**3 * (da_ducR) & - half*dqn/a * (drho_ducR) & + half*(ddp_ducR)/a**2 & - half*rho*(ddqn_ducR)/a dFnducR(:,1) = dFnducR(:,1) - half * ws(1)*dLdU1_ducR(1)*R(:,1) dFnducR(:,2) = dFnducR(:,2) - half * ws(1)*dLdU1_ducR(2)*R(:,1) dFnducR(:,3) = dFnducR(:,3) - half * ws(1)*dLdU1_ducR(3)*R(:,1) dFnducR(:,4) = dFnducR(:,4) - half * ws(1)*dLdU1_ducR(4)*R(:,1) dFnducR(:,5) = dFnducR(:,5) - half * ws(1)*dLdU1_ducR(5)*R(:,1) ! dLdU2_ducR = d( drho - dp/a^2 )/dUcR ! = [d(drho)/dUcR] - [d(dp)/dUcR]/a^2 + 2*dp/a^3*[da/dUcR] dLdU2_ducR = (ddrho_ducR) - (ddp_ducR)/a**2 + two*dp/a**3*(da_ducR) dFnducR(:,1) = dFnducR(:,1) - half * ws(2)*dLdU2_ducR(1)*R(:,2) dFnducR(:,2) = dFnducR(:,2) - half * ws(2)*dLdU2_ducR(2)*R(:,2) dFnducR(:,3) = dFnducR(:,3) - half * ws(2)*dLdU2_ducR(3)*R(:,2) dFnducR(:,4) = dFnducR(:,4) - half * ws(2)*dLdU2_ducR(4)*R(:,2) dFnducR(:,5) = dFnducR(:,5) - half * ws(2)*dLdU2_ducR(5)*R(:,2) ! dLdU3_ducR = d( dp/(2a^2) + rho/(2a)*dqn )/dUcR ! ! = dp*[d(1/(2a^2))/dUcR] + dqn*[d(rho/(2a))/dUcR] ! + [d(dp)/dUcR]/(2a^2) + [d(dqn)/dUcR]*rho/(2a) ! ! = -a^(-3)*dp*[da/dUcR] + dqn*( [drho/dUcR]/(2a) - rho/(2*a^2)*[da/dUcR] ) ! + [d(dp)/dUcR]/(2a^2) + [d(dqn)/dUcR]*rho/(2a) ! ! = ( -2*dp - rho*a*dqn )/(2a^3)*[da/dUcR] + dqn*[drho/dUcR]/(2a) ! + [d(dp)/dUcR]/(2a^2) + [d(dqn)/dUcR]*rho/(2a) dLdU3_ducR = half*(-two*dp-rho*a*dqn )/a**3 * (da_ducR) & + half*dqn/a * (drho_ducR) & + half*(ddp_ducR)/a**2 & + half*rho*(ddqn_ducR)/a dFnducR(:,1) = dFnducR(:,1) - half * ws(3)*dLdU3_ducR(1)*R(:,3) dFnducR(:,2) = dFnducR(:,2) - half * ws(3)*dLdU3_ducR(2)*R(:,3) dFnducR(:,3) = dFnducR(:,3) - half * ws(3)*dLdU3_ducR(3)*R(:,3) dFnducR(:,4) = dFnducR(:,4) - half * ws(3)*dLdU3_ducR(4)*R(:,3) dFnducR(:,5) = dFnducR(:,5) - half * ws(3)*dLdU3_ducR(5)*R(:,3) ! dLdU4_ducR = d(rho)/dUcR dLdU4_ducR = drho_ducR dFnducR(:,1) = dFnducR(:,1) - half * ws(4)*dLdU4_ducR(1)*R(:,4) dFnducR(:,2) = dFnducR(:,2) - half * ws(4)*dLdU4_ducR(2)*R(:,4) dFnducR(:,3) = dFnducR(:,3) - half * ws(4)*dLdU4_ducR(3)*R(:,4) dFnducR(:,4) = dFnducR(:,4) - half * ws(4)*dLdU4_ducR(4)*R(:,4) dFnducR(:,5) = dFnducR(:,5) - half * ws(4)*dLdU4_ducR(5)*R(:,4) !-------------------------------------- ! 2.5 Differentiate the right-eigenvectors, and ! add the fourth term, - 0.5*sum_{k=1,4} |lambda_k|*(LdU)_k*[dr_k/ducL] ! dR1_ducR = dR(:,1)/dUcR ! ! Reft-moving acoustic wave ! ! Eigenvector -> Differentiated ! ! R(1,1) = one -> 0 ! R(2,1) = u - a*nx -> du/dUcR - da/dUcR*nx ! R(3,1) = v - a*ny -> dv/dUcR - da/dUcR*ny ! R(4,1) = w - a*nz -> dw/dUcR - da/dUcR*nz ! R(5,1) = H - a*qn -> dH/dUcR - da/dUcR*qn - dqn/dUcR*a dR1_ducR(1,:) = zero dR1_ducR(2,:) = (du_ducR) - (da_ducR) * nx dR1_ducR(3,:) = (dv_ducR) - (da_ducR) * ny dR1_ducR(4,:) = (dw_ducR) - (da_ducR) * nz dR1_ducR(5,:) = (dH_ducR) - (da_ducR) * qn - (dqn_ducR) * a dFnducR = dFnducR - half * ws(1)*LdU(1)*dR1_ducR ! dR2_ducR = dR(:,2)/dUcR ! ! Entropy wave ! ! Eigenvector -> Differentiated ! ! R(1,2) = one -> 0 ! R(2,2) = u -> du/dUcR ! R(3,2) = v -> dv/dUcR ! R(4,2) = w -> dw/dUcR ! R(5,2) = half*(u*u+v*v+w*w) -> u*du/dUcR + v*dv/dUcR + w*dw/dUcR dR2_ducR(1,:) = zero dR2_ducR(2,:) = (du_ducR) dR2_ducR(3,:) = (dv_ducR) dR2_ducR(4,:) = (dw_ducR) dR2_ducR(5,:) = u*(du_ducR) + v*(dv_ducR) + w*(dw_ducR) dFnducR = dFnducR - half * ws(2)*LdU(2)*dR2_ducR ! dR3_ducR = dR(:,3)/dUcR ! ! Right-moving acoustic wave ! ! Eigenvector -> Differentiated ! ! R(1,3) = one -> 0 ! R(2,3) = u + a*nx -> du/dUcR + da/dUcR*nx ! R(3,3) = v + a*ny -> dv/dUcR + da/dUcR*ny ! R(4,3) = w + a*nz -> dw/dUcR + da/dUcR*nz ! R(5,3) = H + a*qn -> dH/dUcR + da/dUcR*qn + dqn/dUcR*a dR3_ducR(1,:) = zero dR3_ducR(2,:) = (du_ducR) + (da_ducR) * nx dR3_ducR(3,:) = (dv_ducR) + (da_ducR) * ny dR3_ducR(4,:) = (dw_ducR) + (da_ducR) * nz dR3_ducR(5,:) = (dH_ducR) + (da_ducR) * qn + (dqn_ducR) * a dFnducR = dFnducR - half * ws(3)*LdU(3)*dR3_ducR ! dR4_ducR = dR(:,4)/dUcR ! Two shear wave components combined into one (wave strength incorporated). ! So, it is not really an eigenvector. ! ! Combined vector -> Differentiated ! ! R(1,4) = zero -> 0 ! R(2,4) = du - dqn*nx -> d(du)/dUcR - d(dqn)/dUcR*nx ! R(3,4) = dv - dqn*ny -> d(dv)/dUcR - d(dqn)/dUcR*ny ! R(4,4) = dw - dqn*nz -> d(dw)/dUcR - d(dqn)/dUcR*nz ! R(5,4) = u*du+v*dv+w*dw-qn*dqn -> du/dUcR*du + d(du)/dUcR*u ! + dv/dUcR*dv + d(dv)/dUcR*v ! + dw/dUcR*dw + d(dw)/dUcR*w ! - d(qn)/dUcR*dqn - d(dqn)/dUcR*qn dR4_ducR(1,:) = zero dR4_ducR(2,:) = (ddu_ducR) - (ddqn_ducR)*nx dR4_ducR(3,:) = (ddv_ducR) - (ddqn_ducR)*ny dR4_ducR(4,:) = (ddw_ducR) - (ddqn_ducR)*nz dR4_ducR(5,:) = du*( du_ducR) + dv*( dv_ducR) + dw*( dw_ducR) - dqn*( dqn_ducR) & + u*(ddu_ducR) + v*(ddv_ducR) + w*(ddw_ducR) - qn*(ddqn_ducR) dFnducR = dFnducR - half * ws(4)*LdU(4)*dR4_ducR end subroutine inviscid_roe_n_flux_jacobian !-------------------------------------------------------------------------------- !******************************************************************************** !* -- 3D Viscous Flux Function (alpha-damping scheme) and Jacobian --- !* !* This subroutine computes an approximate numerical flux for the viscous part of the !* projected Navier-Stokes flux in the direction njk=[nx,ny,nz]. !* !* Conservative form of the Navier-Stokes equations: !* !* dU/dt + dF/dx + dG/dy = 0 !* !* where !* !* F = Fi + Fv (Fi = inviscid part, Fv = viscous part) !* G = Gi + Gv (Gi = inviscid part, Gv = viscous part) !* !* Viscous flux is given by !* !* Fvn = Fv*nx + Gv*ny = | 0 | !* | - tauxn | !* | - tauyn | !* | - tauzn | !* | - taunv + qn | !* !* where 'qn' is the normal heat flux, and !* !* tauxn = tauxx*nx + tauxy*ny + tauxz*nz, !* tauyn = tauyx*nx + tauyy*ny + tauyz*nz, !* tauzn = tauzx*nx + tauzy*ny + tauzz*nz, !* taunv = tauxn*u + tauyn*v + tauzn*w. !* !* Traditional approach is to individually compute all the quantities in !* the viscous flux above at the interface, and then directly evaluate the !* viscous flux with them. The most critical part is the computation of the !* interface gradients which are required to compute the viscous stresses !* and heat fluxes. Gradients must be defined to ensure consistency, design !* accuracy, robustness, and h-ellipticity. !* !* Accurate and robust formulas for the interface gradient are given in !* Nishikawa-AIAA2010-5093: !* !* http://ossanworld.com/hiroakinishikawa/My_papers/nishikawa_AIAA-2010-5093.pdf !* !* We employ the one proposed for finite-volume schemes, which has been proved !* to be robust for unstructured-grid simulations. The value of the damping !* coefficient, alpha, is taken to be 4/3 !* which corresponds to 4-th order accurate diffusion scheme in 1D. !* The scheme is a generalized average-least-squares scheme, which includes !* widely-used schemes (edge-normal and Mathur-Murthy/face-tangent schemes). !* !* ------------------------------------------------------------------------------ !* Input: ucL(1:5) = Left state (rhoL, rhoL*uL, rhoL*vL, rhoL*wL, rhoL*EL) !* ucR(1:5) = Right state (rhoR, rhoR*uR, rhoR*vR, rhoR*wR, rhoR*ER) !* njk(1:3) = Unit face vector !* ejk(1:3) = Unit edge-vector !* mag_ejk = Magnitude of the edge vector !* !* njk !* Face normal ^ o Right data point !* | . !* | . ejk (edge vector from left to right) !* |. !* -------x-------- Face !* . !* . !* . !* o Left data point !* !* - mag_ejk is the length between the two data points. !* !* !* Output: numerical_flux(1:5) = Numerical viscous flux (alpha-damping scheme), !* and its derivatives as "numerical_flux%df". !* ------------------------------------------------------------------------------ !* !* NOTE: This subsroutine computes an approximate version of the alpha-damping !* viscous flux (edge-terms-only scheme) which ignores average LSQ gradients. !* The numerical flux is therefore not necessarily consistent. !* It is consistent only when the skewness parameter = 1.0. !* It is OK because this subroutine is designed for computing approximate !* derivative of the alpha-damping scheme. !* !* Note: You can input the temperature gradient instead of density and pressure !* gradients if you have it. I think many practical codes do so. !* !* Note: This function is a ddt-version, which means that each variable carries !* its derivatives, and that the resulting flux "numerical_flux" will have !* its derivatives in "numerical_flux%df". You can convert it to a real-value !* version by type(derivative_data_type_df5) -> real(p2). !* !* Note: This subroutine has been prepared for an educational purpose. !* It is not at all efficient. Think about how you can optimize it. !* One way to make it efficient is to reduce the number of local variables, !* by re-using temporary variables as many times as possible. !* !* Note: Please let me know if you find bugs. I'll greatly appreciate it and !* fix the bugs. !* !* Katate Masatsuka, December 2012. http://www.cfdbooks.com !******************************************************************************** subroutine viscous_alpha_ddt(ucL,ucR,gradwL,gradwR, njk,ejk,mag_ejk, numerical_flux) use derivative_data_df5 implicit none integer , parameter :: p2 = selected_real_kind(15) !Input type(derivative_data_type_df5), dimension(5) , intent( in) :: ucL !Left state (conservative) type(derivative_data_type_df5), dimension(5) , intent( in) :: ucR !Right state (conservative) real(p2) , dimension(5,3), intent( in) :: gradwL !Left gradient (primitive) real(p2) , dimension(5,3), intent( in) :: gradwR !Right gradient (primitive) real(p2) , dimension(3) , intent( in) :: njk !Unit directed area vector real(p2) , dimension(3) , intent( in) :: ejk !Unit edge vector real(p2) , intent( in) :: mag_ejk !Magnitude of the edge vector !Output type(derivative_data_type_df5), dimension(5), intent(out) :: numerical_flux !Numerical viscous flux !Some constants real(p2) :: zero = 0.0_p2 real(p2) :: one = 1.0_p2 real(p2) :: three = 3.0_p2 real(p2) :: half = 0.5_p2 real(p2) :: two_third = 2.0_p2/3.0_p2 real(p2) :: four_third = 4.0_p2/3.0_p2 real(p2) :: C = 110.5_p2 !Parameter for Sutherland's law real(p2) :: gamma = 1.4_p2 !Ratio of specific heats real(p2) :: Prandtl = 0.72_p2 !Prandtl number real(p2) :: Re_inf = 100.0_p2 !Free stream Reynolds number real(p2) :: T_inf = 300.0_p2 !Free stream temperature real(p2) :: M_inf = 2.0_p2 !Free stream Mach number !Local variables real(p2) :: alpha !Damping coefficient (see Nishikawa AIAA2010-5093) real(p2) :: Lr !Length scale (see Nishikawa AIAA2010-5093) type(derivative_data_type_df5) :: uL, uR ! x-velocity (Left and Right states) type(derivative_data_type_df5) :: vL, vR ! y-velocity (Left and Right states) type(derivative_data_type_df5) :: wL, wR ! z-velocity (Left and Right states) type(derivative_data_type_df5) :: rhoL, rhoR ! Density (Left and Right states) type(derivative_data_type_df5) :: pL, pR ! Pressure (Left and Right states) type(derivative_data_type_df5) :: TL, TR ! Temperature (Left and Right states) type(derivative_data_type_df5) :: tauxx, tauyy, tauzz !Viscous stresses: diagonal compontens type(derivative_data_type_df5) :: tauxy, tauyz, tauzx !Viscous stresses: off-diagonal components type(derivative_data_type_df5) :: tauyx, tauzy, tauxz !Viscous stresses: same as above by symmetry type(derivative_data_type_df5) :: qx, qy, qz !Heat flux components type(derivative_data_type_df5) :: tauxn, tauyn, tauzn !Normal stresses type(derivative_data_type_df5) :: qn !Normal heat flux type(derivative_data_type_df5) :: u, v, w, T, mu !Interface quantities type(derivative_data_type_df5), dimension(3) :: grad_u, grad_v, grad_w !Interface gradients of velocities type(derivative_data_type_df5), dimension(3) :: grad_rho, grad_p, grad_T !Interface gradients of rho, p, and T real(p2), dimension(3) :: grad_uL, grad_vL, grad_wL, grad_rL, grad_pL real(p2), dimension(3) :: grad_uR, grad_vR, grad_wR, grad_rR, grad_pR type(derivative_data_type_df5) :: rho, a2 !Interface values for density and (speed of sound)^2 integer :: ix, iy, iz !Index to extract the x-, y-, and z-component. alpha = one ix = 1 iy = 2 iz = 3 ! Left and right states and gradients: rhoL = ucL(1) uL = ucL(2)/ucL(1) vL = ucL(3)/ucL(1) wL = ucL(4)/ucL(1) pL = (gamma-one)*( ucL(5) - half*rhoL*(uL*uL+vL*vL+wL*wL) ) grad_rL = gradwL(1,:) grad_uL = gradwL(2,:) grad_vL = gradwL(3,:) grad_wL = gradwL(4,:) grad_pL = gradwL(5,:) rhoR = ucR(1) uR = ucR(2)/ucR(1) vR = ucR(3)/ucR(1) wR = ucR(4)/ucR(1) pR = (gamma-one)*( ucR(5) - half*rhoR*(uR*uR+vR*vR+wR*wR) ) grad_rR = gradwR(1,:) grad_uR = gradwR(2,:) grad_vR = gradwR(3,:) grad_wR = gradwR(4,:) grad_pR = gradwR(5,:) ! Temperature is identical to the speed of sound due to ! the nondimensionalization. TL = gamma*pL/rhoL TR = gamma*pR/rhoR ! Arithmetic averages of velocities and temperature. u = half*(uL + uR) v = half*(vL + vR) w = half*(wL + wR) T = half*(TL + TR) ! Sutherland's law in the nondimensional form. ! Note: The factor, M_inf/Re_inf, comes from nondimensionalization. mu = (one+C/T_inf)/(T+C/T_inf)*T**(three*half) * M_inf/Re_inf ! Damping coefficient, alpha: ! (1)alpha=1 gives the central-difference formula in 1D. ! (2)alpha=4/3 corresponds to the 4th-order diffusion scheme in 1D. alpha = four_third ! Lr = Length scale involving the skewness measure (see Nishikawa AIAA2010-5093) ! This is the key quantity for robust and accurate computations on skewed grids. Lr = half*abs( njk(ix)*ejk(ix) + njk(iy)*ejk(iy) + njk(iz)*ejk(iz) ) * mag_ejk ! Interface gradients from the derived diffusion scheme (Nishikawa-AIAA2010-5093). ! The second term is the damping term. grad_u = half*( (grad_uR + grad_uL) + alpha/Lr*(uR-uL)*njk ) grad_v = half*( (grad_vR + grad_vL) + alpha/Lr*(vR-vL)*njk ) grad_w = half*( (grad_wR + grad_wL) + alpha/Lr*(wR-wL)*njk ) ! The temperature gradient is computed from the interface density and the pressure ! gradients. ! Note: T = gamma*p/rho -> grad(T) = gamma*grad(p)/rho - (gamma*p/rho^2)*grad(rho) rho = half*(rhoR + rhoL) a2 = gamma*half*(pR + pL)/rho grad_rho = half*( (grad_rR + grad_rL) + alpha/Lr*(rhoR-rhoL)*njk ) grad_p = half*( (grad_pR + grad_pL) + alpha/Lr*( pR-pL )*njk ) grad_T = ( gamma*grad_p - a2*grad_rho) /rho ! Interface gradients have been computed: grad(u), grad(v), grad(w), grad(T). ! We now evaluate the physical viscous flux with them. ! Viscous stresses (Stokes' hypothesis is assumed) tauxx = mu*(four_third*grad_u(ix) - two_third*grad_v(iy) - two_third*grad_w(iz)) tauyy = mu*(four_third*grad_v(iy) - two_third*grad_u(ix) - two_third*grad_w(iz)) tauzz = mu*(four_third*grad_w(iz) - two_third*grad_u(ix) - two_third*grad_v(iy)) tauxy = mu*(grad_u(iy) + grad_v(ix)) tauxz = mu*(grad_u(iz) + grad_w(ix)) tauyz = mu*(grad_v(iz) + grad_w(iy)) tauyx = tauxy tauzx = tauxz tauzy = tauyz ! Heat fluxes: q = - mu*grad(T)/(Prandtl*(gamma-1)) qx = - mu*grad_T(ix)/(Prandtl*(gamma-one)) qy = - mu*grad_T(iy)/(Prandtl*(gamma-one)) qz = - mu*grad_T(iz)/(Prandtl*(gamma-one)) ! Normal components tauxn = tauxx*njk(ix) + tauxy*njk(iy) + tauxz*njk(iz) tauyn = tauyx*njk(ix) + tauyy*njk(iy) + tauyz*njk(iz) tauzn = tauzx*njk(ix) + tauzy*njk(iy) + tauzz*njk(iz) qn = qx*njk(ix) + qy*njk(iy) + qz*njk(iz) ! Evaluate the viscous flux at the interface numerical_flux(1) = zero numerical_flux(2) = - tauxn numerical_flux(3) = - tauyn numerical_flux(4) = - tauzn numerical_flux(5) = - (tauxn*u + tauyn*v + tauzn*w) + qn ! Normal max wave speed ! wsn = alpha*(mu/rho*gamma/Prandtl)/Lr end subroutine viscous_alpha_ddt !-------------------------------------------------------------------------------- !******************************************************************************** !* -- 3D Approximate Viscous Flux Jacobians for alpha-damping scheme --- !* !* This subroutine computes the left and right Jacobians for the alpha viscous !* flux in the direction, njk=[nx,ny,nz]. !* !* Conservative form of the Navier-Stokes equations: !* !* dU/dt + dF/dx + dG/dy = 0 !* !* where !* !* F = Fi + Fv (Fi = inviscid part, Fv = viscous part) !* G = Gi + Gv (Gi = inviscid part, Gv = viscous part) !* !* Viscous flux is given by !* !* Fvn = Fv*nx + Gv*ny = | 0 | !* | - tauxn | !* | - tauyn | !* | - tauzn | !* | - taunv + qn | !* !* where 'qn' is the normal heat flux, and !* !* tauxn = tauxx*nx + tauxy*ny + tauxz*nz, !* tauyn = tauyx*nx + tauyy*ny + tauyz*nz, !* tauzn = tauzx*nx + tauzy*ny + tauzz*nz, !* taunv = tauxn*u + tauyn*v + tauzn*w. !* !* Traditional approach is to individually compute all the quantities in !* the viscous flux above at the interface, and then directly evaluate the !* viscous flux with them. The most critical part is the computation of the !* interface gradients which are required to compute the viscous stresses !* and heat fluxes. Gradients must be defined to ensure consistency, design !* accuracy, robustness, and h-ellipticity. !* !* Accurate and robust formulas for the interface gradient are given in !* Nishikawa-AIAA2010-5093: !* !* http://ossanworld.com/hiroakinishikawa/My_papers/nishikawa_AIAA-2010-5093.pdf !* !* We employ the one proposed for finite-volume schemes, which has been proved !* to be robust for unstructured-grid simulations. The value of the damping !* coefficient, alpha, is taken to be 4/3 !* which corresponds to 4th order accurate diffusion scheme in 1D. !* The scheme is a generalized average-least-squares scheme, which includes !* widely-used schemes (edge-normal and Mathur-Murthy/face-tangent schemes). !* !* Note: In the code, the vector of conserative variables are denoted by uc. !* !* ------------------------------------------------------------------------------ !* Input: ucL(1:5) = Left state (rhoL, rhoL*uL, rhoL*vL, rhoL*wL, rhoL*EL) !* ucR(1:5) = Right state (rhoR, rhoR*uR, rhoR*vR, rhoR*wR, rhoR*ER) !* njk(1:3) = unit face normal vector (nx, ny, nz), pointing from Left to Right. !* ejk(1:3) = Unit edge-vector !* mag_ejk = Magnitude of the edge vector !* !* njk !* Face normal ^ o Right data point !* | . !* | . !* |. !* -------x-------- Face !* . Left and right states are values at data points. !* . !* o Left data point !* !* - mag_ejk is the length between the two data points. !* !* !* Output: dFnducL: Derivative of the alpha flux w.r.t. the left state ucL !* dFnducR: Derivative of the alpha flux w.r.t. the left state ucR !* ------------------------------------------------------------------------------ !* !* NOTE: This subsroutine computes an approximate version of the alpha-damping !* viscous flux (damping-terms-only scheme) which ignores LSQ gradients. !* !* !* Note: This subroutine has been prepared for an educational purpose. !* It is not at all efficient. Think about how you can optimize it. !* One way to make it efficient is to reduce the number of local variables, !* by re-using temporary variables as many times as possible. !* !* Note: Please let me (Hiro) know at sunmasen@hotmail.com if you find bugs. !* I'll greatly appreciate it and fix the bugs. !* !* Katate Masatsuka, December 2012. http://www.cfdbooks.com !******************************************************************************** subroutine viscous_alpha_flux_jacobian(ucL, ucR, njk,ejk,mag_ejk, dFnducL, dFnducR) implicit none integer , parameter :: p2 = selected_real_kind(15) ! Double precision !Input real(p2), dimension(5), intent( in) :: ucL ! Left state (conservative variables) real(p2), dimension(5), intent( in) :: ucR ! Right state (conservative variables) real(p2), dimension(3), intent( in) :: njk ! Normal vector real(p2), dimension(3), intent( in) :: ejk ! Unit edge vector real(p2), intent( in) :: mag_ejk ! Magnitude of the edge vector !Output real(p2), dimension(5,5), intent(out) :: dFnducL ! Left Jacobian, d(Fn_Roe)/d(ucL) real(p2), dimension(5,5), intent(out) :: dFnducR ! Right Jacobian, d(Fn_Roe)/d(ucR) !Some constants real(p2) :: zero = 0.0_p2 real(p2) :: one = 1.0_p2 real(p2) :: two = 2.0_p2 real(p2) :: three = 3.0_p2 real(p2) :: half = 0.5_p2 real(p2) :: two_third = 2.0_p2/3.0_p2 real(p2) ::four_third = 4.0_p2/3.0_p2 real(p2) :: C = 110.5_p2 !Parameter for Sutherland's law real(p2) :: gamma = 1.4_p2 !Ratio of specific heats real(p2) :: Prandtl = 0.72_p2 !Prandtl number real(p2) :: Re_inf = 100.0_p2 !Free stream Reynolds number real(p2) :: T_inf = 300.0_p2 !Free stream temperature real(p2) :: M_inf = 2.0_p2 !Free stream Mach number !Local variables real(p2) :: alpha !Damping coefficient (see Nishikawa AIAA2010-5093) real(p2) :: Lr !Length scale (see Nishikawa AIAA2010-5093) real(p2) :: rhoL, rhoR ! Density (Left and Right states) real(p2) :: pL, pR ! Pressure (Left and Right states) real(p2) :: uL, uR ! x-velocity (Left and Right states) real(p2) :: vL, vR ! y-velocity (Left and Right states) real(p2) :: wL, wR ! z-velocity (Left and Right states) real(p2) :: q2L, q2R ! Speed squared (Left and Right states) real(p2) :: TL, TR ! Temperature (Left and Right states) real(p2) :: sxx, syy, szz ! Deformation tensor real(p2) :: sxy, syz, sxz ! Deformation tensor real(p2) :: snx, sny, snz, snv ! Projected deformation tensor real(p2) :: u, v, w, T, mu !Interface quantities real(p2) :: grad_u(3) , grad_v(3), grad_w(3) !Interface Gradients of velocities real(p2) :: grad_rho(3), grad_p(3), grad_T(3) !Interface Gradients of rho, p, and T real(p2) :: rho, a2 !Interface values for density and (speed of sound)^2 integer :: ix, iy, iz !Index to extract the x-, y-, and z-component. ! Derivatives w.r.t. ucL real(p2) :: dmu_dTL real(p2) :: dmu_drhoL real(p2) :: dmu_dpL real(p2) :: dux_duL, duy_duL, duz_duL real(p2) :: dvx_dvL, dvy_dvL, dvz_dvL real(p2) :: dwx_dwL, dwy_dwL, dwz_dwL real(p2) :: dsxx_duL, dsyy_duL, dszz_duL real(p2) :: dsxy_duL, dsxz_duL, dsyz_duL real(p2) :: dsxx_dvL, dsyy_dvL, dszz_dvL real(p2) :: dsxy_dvL, dsxz_dvL, dsyz_dvL real(p2) :: dsxx_dwL, dsyy_dwL, dszz_dwL real(p2) :: dsxy_dwL, dsxz_dwL, dsyz_dwL real(p2) :: dsnx_duL, dsny_duL, dsnz_duL real(p2) :: dsnx_dvL, dsny_dvL, dsnz_dvL real(p2) :: dsnx_dwL, dsny_dwL, dsnz_dwL real(p2) :: dtauxx_duL, dtauxx_dvL, dtauxx_dwL real(p2) :: dtauyy_duL, dtauyy_dvL, dtauyy_dwL real(p2) :: dtauzz_duL, dtauzz_dvL, dtauzz_dwL real(p2) :: dtauxy_duL, dtauxy_dvL, dtauxy_dwL real(p2) :: dtauxz_duL, dtauxz_dvL, dtauxz_dwL real(p2) :: dtauyz_duL, dtauyz_dvL, dtauyz_dwL real(p2) :: dtaunx_drhoL, dtaunx_duL, dtaunx_dvL, dtaunx_dwL, dtaunx_dpL real(p2) :: dtauny_drhoL, dtauny_duL, dtauny_dvL, dtauny_dwL, dtauny_dpL real(p2) :: dtaunz_drhoL, dtaunz_duL, dtaunz_dvL, dtaunz_dwL, dtaunz_dpL real(p2) :: dtaunv_drhoL, dtaunv_duL, dtaunv_dvL, dtaunv_dwL, dtaunv_dpL real(p2) :: drhox_drhoL, drhoy_drhoL, drhoz_drhoL real(p2) :: dpx_dpL, dpy_dpL, dpz_dpL real(p2) :: dTx_drhoL, dTy_drhoL, dTz_drhoL real(p2) :: dTx_dpL, dTy_dpL, dTz_dpL real(p2) :: dqx_drhoL, dqy_drhoL, dqz_drhoL real(p2) :: dqx_dpL, dqy_dpL, dqz_dpL real(p2) :: dqn_drhoL, dqn_dpL real(p2), dimension(5,5) :: dFndwL ! Derivatives w.r.t. ucR real(p2) :: dmu_dTR real(p2) :: dmu_drhoR real(p2) :: dmu_dpR real(p2) :: dux_duR, duy_duR, duz_duR real(p2) :: dvx_dvR, dvy_dvR, dvz_dvR real(p2) :: dwx_dwR, dwy_dwR, dwz_dwR real(p2) :: dsxx_duR, dsyy_duR, dszz_duR real(p2) :: dsxy_duR, dsxz_duR, dsyz_duR real(p2) :: dsxx_dvR, dsyy_dvR, dszz_dvR real(p2) :: dsxy_dvR, dsxz_dvR, dsyz_dvR real(p2) :: dsxx_dwR, dsyy_dwR, dszz_dwR real(p2) :: dsxy_dwR, dsxz_dwR, dsyz_dwR real(p2) :: dsnx_duR, dsny_duR, dsnz_duR real(p2) :: dsnx_dvR, dsny_dvR, dsnz_dvR real(p2) :: dsnx_dwR, dsny_dwR, dsnz_dwR real(p2) :: dtauxx_duR, dtauxx_dvR, dtauxx_dwR real(p2) :: dtauyy_duR, dtauyy_dvR, dtauyy_dwR real(p2) :: dtauzz_duR, dtauzz_dvR, dtauzz_dwR real(p2) :: dtauxy_duR, dtauxy_dvR, dtauxy_dwR real(p2) :: dtauxz_duR, dtauxz_dvR, dtauxz_dwR real(p2) :: dtauyz_duR, dtauyz_dvR, dtauyz_dwR real(p2) :: dtaunx_drhoR, dtaunx_duR, dtaunx_dvR, dtaunx_dwR, dtaunx_dpR real(p2) :: dtauny_drhoR, dtauny_duR, dtauny_dvR, dtauny_dwR, dtauny_dpR real(p2) :: dtaunz_drhoR, dtaunz_duR, dtaunz_dvR, dtaunz_dwR, dtaunz_dpR real(p2) :: dtaunv_drhoR, dtaunv_duR, dtaunv_dvR, dtaunv_dwR, dtaunv_dpR real(p2) :: drhox_drhoR, drhoy_drhoR, drhoz_drhoR real(p2) :: dpx_dpR, dpy_dpR, dpz_dpR real(p2) :: dTx_drhoR, dTy_drhoR, dTz_drhoR real(p2) :: dTx_dpR, dTy_dpR, dTz_dpR real(p2) :: dqx_drhoR, dqy_drhoR, dqz_drhoR real(p2) :: dqx_dpR, dqy_dpR, dqz_dpR real(p2) :: dqn_drhoR, dqn_dpR real(p2), dimension(5,5) :: dFndwR ix = 1 iy = 2 iz = 3 ! Left and right states: rhoL = ucL(1) uL = ucL(2)/ucL(1) vL = ucL(3)/ucL(1) wL = ucL(4)/ucL(1) q2L = uL*uL+vL*vL+wL*wL pL = (gamma-one)*( ucL(5) - half*rhoL*q2L ) rhoR = ucR(1) uR = ucR(2)/ucR(1) vR = ucR(3)/ucR(1) wR = ucR(4)/ucR(1) q2R = uR*uR+vR*vR+wR*wR pR = (gamma-one)*( ucR(5) - half*rhoR*q2R ) ! Temperature is identical to the speed of sound due to ! the nondimensionalization. TL = gamma*pL/rhoL TR = gamma*pR/rhoR ! Arithmetic averages of velocities and temperature. u = half*(uL + uR) v = half*(vL + vR) w = half*(wL + wR) T = half*(TL + TR) ! Sutherland's law in the nondimensional form. ! Note: The factor, M_inf/Re_inf, arises from nondimensionalization. mu = (one+C/T_inf)/(T+C/T_inf)*T**(three*half) * M_inf/Re_inf ! Damping coefficient, alpha: ! (1)alpha=1 gives the central-difference formula in 1D. ! (2)alpha=4/3 corresponds to the 4th-order diffusion scheme in 1D. alpha = four_third ! Lr = Length scale involving the skewness measure (see Nishikawa AIAA2010-5093) ! This is the key quantity for robust and accurate computations on skewed grids. Lr = abs( njk(ix)*ejk(ix) + njk(iy)*ejk(iy) + njk(iz)*ejk(iz) ) * mag_ejk ! Interface gradients from the derived diffusion scheme (Nishikawa-AIAA2010-5093). ! The second term is the damping term. grad_u = alpha/Lr*(uR - uL)*njk ! Average LSQ gradient ignored grad_v = alpha/Lr*(vR - vL)*njk ! Average LSQ gradient ignored grad_w = alpha/Lr*(wR - wL)*njk ! Average LSQ gradient ignored ! The temperature gradient is computed from the interface density and the pressure ! gradients. ! Note: T = gamma*p/rho -> grad(T) = gamma*grad(p)/rho - (gamma*p/rho^2)*grad(rho) rho = half*(rhoR + rhoL) a2 = gamma*half*(pR + pL)/rho grad_rho = alpha/Lr*(rhoR - rhoL)*njk ! Average LSQ gradient ignored grad_p = alpha/Lr*( pR - pL )*njk ! Average LSQ gradient ignored grad_T = (gamma*grad_p - a2*grad_rho) /rho ! Viscous stresses (Stokes' hypothesis is assumed) ! tauxx = mu*sxx ! tauyy = mu*sxy ! tauzz = mu*szz sxx = four_third*grad_u(ix) - two_third*grad_v(iy) - two_third*grad_w(iz) syy = four_third*grad_v(iy) - two_third*grad_u(ix) - two_third*grad_w(iz) szz = four_third*grad_w(iz) - two_third*grad_u(ix) - two_third*grad_v(iy) ! tauxy = mu*sxy ! tauxz = mu*sxz ! tauyz = mu*syz sxy = grad_u(iy) + grad_v(ix) sxz = grad_u(iz) + grad_w(ix) syz = grad_v(iz) + grad_w(iy) ! Normal components ! taunx = mu*snx ! tauny = mu*sny ! taunz = mu*snz snx = sxx*njk(ix) + sxy*njk(iy) + sxz*njk(iz) sny = sxy*njk(ix) + syy*njk(iy) + syz*njk(iz) snz = sxz*njk(ix) + syz*njk(iy) + szz*njk(iz) ! Viscous work ! taunv = mu*snv = mu* ( snx*u + sny*v + snz*w ) snv = snx*u + sny*v + snz*w ! We are now ready to compute the approximate Jacobians for the alpha-damping flux: !-------------------------------------------------------------------------------- ! Part 1. Compute dFn_alpha/ducL ! ! dFn_alpha/ducL = [ zero, -d(taunx)/ducL, -d(tauny)/ducL, -d(taunz)/ducL, d(-taunv+qn)/ducL] ! = (dFn_alpha/dwL)*(dwL/ducL) ! ! where wL = [rhoL, uL, vL, wL, pL] - primitive variables. ! ! So, we proceed as follows: ! ! 1.1 Compute the derivatives of the viscosity: dmu/drhoL, dmu/dpL ! 1.2 Compute d(taunx)/dwL, -d(tauny)/dwL, -d(taunz)/dwL, -d(taunv)/dwL ! 1.3 Compute d(qn)/dwL ! 1.4 Compute the flux Jacbian with respect to the primitive variables: ! -> dFndwL ! 1.5 Compute the flux Jacbian with respect to the conservative variables: ! -> dFnducL = dFndwL*(dWL/ducL) !---------------------------------------------------------------------- ! 1.1 Compute the derivatives of the viscosity: dmu/drhoL, dmu/dpL ! dmu/dTL = d( (one+C/T_inf)/(T+C/T_inf)*T**(three*half) * M_inf/Re_inf )/dTL dmu_dTL = half*mu/(T+C/T_inf)*(one + three*(C/T_inf)/T) * ( half ) ! Note: dmu/dW = [dmu/drho, dmu/du, dmu/dv, dmu/dw, dmu/dp ] = [dmu/drho, 0, 0, 0, dmu/dp ] ! So, we only need dmu/drho and dmu/dp. ! dmu/drhoL = dmu/dTL*(dTL/drhoL); dT/drho = -gamma*p/rho**2 because T = gamma*p/rho. dmu_drhoL = dmu_dTL * (-gamma*pL/rhoL**2) ! dmu/dpL = dmu/dTL*(dTL/dpL); dT/dp = gamma/rho because T = gamma*p/rho. dmu_dpL = dmu_dTL * ( gamma/rhoL ) !---------------------------------------------------------------------- ! 1.2 Compute d(taunx)/dwL, -d(tauny)/dwL, -d(taunz)/dwL, -d(taunv)/dwL ! ! Derivatives of the interface velocity gradients: ! ! Interface gradients are given by ! ! grad_u = (Average LSQ gradient of u) + alpha/Lr*(uR - uL)*njk ! grad_v = (Average LSQ gradient of v) + alpha/Lr*(vR - vL)*njk ! grad_w = (Average LSQ gradient of w) + alpha/Lr*(wR - wL)*njk ! ! but for simplicity and compantness, we ignore the average LSQ gradients. ! Note: If the average LSQ gradients are computed on a compact stencil, ! we may add their contributions. dux_duL = alpha/Lr*(zero-one)*njk(ix) duy_duL = alpha/Lr*(zero-one)*njk(iy) duz_duL = alpha/Lr*(zero-one)*njk(iz) dvx_dvL = alpha/Lr*(zero-one)*njk(ix) dvy_dvL = alpha/Lr*(zero-one)*njk(iy) dvz_dvL = alpha/Lr*(zero-one)*njk(iz) dwx_dwL = alpha/Lr*(zero-one)*njk(ix) dwy_dwL = alpha/Lr*(zero-one)*njk(iy) dwz_dwL = alpha/Lr*(zero-one)*njk(iz) ! Derivatives of the stress tensor. ! Note: many terms are zero, e.g., dux/dvL=0, dwz/dvL=0, etc. ! dsxx_duL = four_third*dux_duL - two_third*( dvy_duL + dwz_duL ) ! dsyy_duL = four_third*dvy_duL - two_third*( dux_duL + dwz_duL ) ! dszz_duL = four_third*dwz_duL - two_third*( dux_duL + dvy_duL ) ! Ignoring terms that are zero, we obtain the following. dsxx_duL = four_third*dux_duL dsyy_duL = - two_third*( dux_duL ) dszz_duL = - two_third*( dux_duL ) ! dsxy_duL = duy_duL + dvx_duL ! dsxz_duL = duz_duL + dwx_duL ! dsyz_duL = dvz_duL + dwy_duL ! Ignoring terms that are zero, we obtain the following. dsxy_duL = duy_duL dsxz_duL = duz_duL dsyz_duL = zero ! dsxx_dvL = four_third*dux_dvL - two_third*( dvy_dvL + dwz_dvL ) ! dsyy_dvL = four_third*dvy_dvL - two_third*( dux_dvL + dwz_dvL ) ! dszz_dvL = four_third*dwz_dvL - two_third*( dux_dvL + dvy_dvL ) ! Ignoring terms that are zero, we obtain the following. dsxx_dvL = - two_third*( dvy_dvL ) dsyy_dvL = four_third*dvy_dvL dszz_dvL = - two_third*( dvy_dvL ) ! dsxy_dvL = duy_dvL + dvx_dvL ! dsxz_dvL = duz_dvL + dwx_dvL ! dsyz_dvL = dvz_dvL + dwy_dvL ! Ignoring terms that are zero, we obtain the following. dsxy_dvL = dvx_dvL dsxz_dvL = zero dsyz_dvL = dvz_dvL ! dsxx_dwL = four_third*dux_dwL - two_third*( dvy_dwL + dwz_dwL ) ! dsyy_dwL = four_third*dvy_dwL - two_third*( dux_dwL + dwz_dwL ) ! dszz_dwL = four_third*dwz_dwL - two_third*( dux_dwL + dvy_dwL ) ! Ignoring terms that are zero, we obtain the following. dsxx_dwL = - two_third*( dwz_dwL ) dsyy_dwL = - two_third*( dwz_dwL ) dszz_dwL = four_third*dwz_dwL ! dsxy_dwL = duy_dwL + dvx_dwL ! dsxz_dwL = duz_dwL + dwx_dwL ! dsyz_dwL = dvz_dwL + dwy_dwL ! Ignoring terms that are zero, we obtain the following. dsxy_dwL = zero dsxz_dwL = dwx_dwL dsyz_dwL = dwy_dwL dsnx_duL = dsxx_duL*njk(ix) + dsxy_duL*njk(iy) + dsxz_duL*njk(iz) dsny_duL = dsxy_duL*njk(ix) + dsyy_duL*njk(iy) + dsyz_duL*njk(iz) dsnz_duL = dsxz_duL*njk(ix) + dsyz_duL*njk(iy) + dszz_duL*njk(iz) dsnx_dvL = dsxx_dvL*njk(ix) + dsxy_dvL*njk(iy) + dsxz_dvL*njk(iz) dsny_dvL = dsxy_dvL*njk(ix) + dsyy_dvL*njk(iy) + dsyz_dvL*njk(iz) dsnz_dvL = dsxz_dvL*njk(ix) + dsyz_dvL*njk(iy) + dszz_dvL*njk(iz) dsnx_dwL = dsxx_dwL*njk(ix) + dsxy_dwL*njk(iy) + dsxz_dwL*njk(iz) dsny_dwL = dsxy_dwL*njk(ix) + dsyy_dwL*njk(iy) + dsyz_dwL*njk(iz) dsnz_dwL = dsxz_dwL*njk(ix) + dsyz_dwL*njk(iy) + dszz_dwL*njk(iz) ! Derivatives of the viscous stresses. ! Note: mu does not depend on the velocity. dtauxx_duL = mu*dsxx_duL dtauxx_dvL = mu*dsxx_dvL dtauxx_dwL = mu*dsxx_dwL dtauyy_duL = mu*dsyy_duL dtauyy_dvL = mu*dsyy_dvL dtauyy_dwL = mu*dsyy_dwL dtauzz_duL = mu*dszz_duL dtauzz_dvL = mu*dszz_dvL dtauzz_dwL = mu*dszz_dwL dtauxy_duL = mu*dsxy_duL dtauxy_dvL = mu*dsxy_dvL dtauxy_dwL = mu*dsxy_dwL dtauxz_duL = mu*dsxz_duL dtauxz_dvL = mu*dsxz_dvL dtauxz_dwL = mu*dsxz_dwL dtauyz_duL = mu*dsyz_duL dtauyz_dvL = mu*dsyz_dvL dtauyz_dwL = mu*dsyz_dwL ! Derivatives of the viscous stresses in the normal direction. ! Derivatives of taunx dtaunx_drhoL = dmu_drhoL * snx dtaunx_duL = dtauxx_duL*njk(ix) + dtauxy_duL*njk(iy) + dtauxz_duL*njk(iz) dtaunx_dvL = dtauxx_dvL*njk(ix) + dtauxy_dvL*njk(iy) + dtauxz_dvL*njk(iz) dtaunx_dwL = dtauxx_dwL*njk(ix) + dtauxy_dwL*njk(iy) + dtauxz_dwL*njk(iz) dtaunx_dpL = dmu_dpL * snx ! Derivatives of tauny dtauny_drhoL = dmu_drhoL * sny dtauny_duL = dtauxy_duL*njk(ix) + dtauyy_duL*njk(iy) + dtauyz_duL*njk(iz) dtauny_dvL = dtauxy_dvL*njk(ix) + dtauyy_dvL*njk(iy) + dtauyz_dvL*njk(iz) dtauny_dwL = dtauxy_dwL*njk(ix) + dtauyy_dwL*njk(iy) + dtauyz_dwL*njk(iz) dtauny_dpL = dmu_dpL * sny ! Derivatives of taunz dtaunz_drhoL = dmu_drhoL * snz dtaunz_duL = dtauxz_duL*njk(ix) + dtauyz_duL*njk(iy) + dtauzz_duL*njk(iz) dtaunz_dvL = dtauxz_dvL*njk(ix) + dtauyz_dvL*njk(iy) + dtauzz_dvL*njk(iz) dtaunz_dwL = dtauxz_dwL*njk(ix) + dtauyz_dwL*njk(iy) + dtauzz_dwL*njk(iz) dtaunz_dpL = dmu_dpL * snz ! Derivatives of taunv = mu* ( snx*u + sny*v + snz*w ) dtaunv_drhoL = dmu_drhoL * snv dtaunv_duL = mu * ( snx*(half) + dsnx_duL*u + dsny_duL*v + dsnz_duL*w ) dtaunv_dvL = mu * ( sny*(half) + dsnx_dvL*u + dsny_dvL*v + dsnz_dvL*w ) dtaunv_dwL = mu * ( snz*(half) + dsnx_dwL*u + dsny_dwL*v + dsnz_dwL*w ) dtaunv_dpL = dmu_dpL * snv !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! 1.3 Compute d(qn)/dwL ! First compute the derivatives of the interface gradients: ! ! rho = half*(rhoR + rhoL) ! a2 = gamma*half*(pR + pL)/rho ! grad_rho = (Average LSQ gradient of rho) + alpha/Lr*(rhoR - rhoL)*njk ! grad_p = (Average LSQ gradient of p) + alpha/Lr*( pR - pL )*njk ! grad_T = ( gamma*grad_p - a2*grad_rho) /rho ! ! Again, we ignore the average LSQ gradients. drhox_drhoL = alpha/Lr*(zero-one)*njk(ix) ! Average LSQ gradient ignored drhoy_drhoL = alpha/Lr*(zero-one)*njk(iy) ! Average LSQ gradient ignored drhoz_drhoL = alpha/Lr*(zero-one)*njk(iz) ! Average LSQ gradient ignored dTx_drhoL = (- gamma*grad_p(ix)/rho**2 & + two* gamma*half*(pR+pL)/rho**3*grad_rho(ix) ) * half & - a2/rho*drhox_drhoL dTy_drhoL = (- gamma*grad_p(iy)/rho**2 & + two* gamma*half*(pR+pL)/rho**3*grad_rho(iy) ) * half & - a2/rho*drhoy_drhoL dTz_drhoL = (- gamma*grad_p(iz)/rho**2 & + two* gamma*half*(pR+pL)/rho**3*grad_rho(iz) ) * half & - a2/rho*drhoz_drhoL dpx_dpL = alpha/Lr*(zero-one)*njk(ix) ! Average LSQ gradient ignored dpy_dpL = alpha/Lr*(zero-one)*njk(iy) ! Average LSQ gradient ignored dpz_dpL = alpha/Lr*(zero-one)*njk(iz) ! Average LSQ gradient ignored dTx_dpL = gamma*dpx_dpL/rho - ( gamma/rho**2*grad_rho(ix) ) * half dTy_dpL = gamma*dpy_dpL/rho - ( gamma/rho**2*grad_rho(iy) ) * half dTz_dpL = gamma*dpz_dpL/rho - ( gamma/rho**2*grad_rho(iz) ) * half ! Derivatives of the heat fluxes: ! ! qx = - mu*grad_T(ix)/(Prandtl*(gamma-one)) ! qy = - mu*grad_T(iy)/(Prandtl*(gamma-one)) ! qz = - mu*grad_T(iz)/(Prandtl*(gamma-one)) dqx_drhoL = - ( dmu_drhoL*grad_T(ix) + mu*dTx_drhoL )/(Prandtl*(gamma-one)) dqx_dpL = - ( dmu_dpL *grad_T(ix) + mu*dTx_dpL )/(Prandtl*(gamma-one)) dqy_drhoL = - ( dmu_drhoL*grad_T(iy) + mu*dTy_drhoL )/(Prandtl*(gamma-one)) dqy_dpL = - ( dmu_dpL *grad_T(iy) + mu*dTy_dpL )/(Prandtl*(gamma-one)) dqz_drhoL = - ( dmu_drhoL*grad_T(iz) + mu*dTz_drhoL )/(Prandtl*(gamma-one)) dqz_dpL = - ( dmu_dpL *grad_T(iz) + mu*dTz_dpL )/(Prandtl*(gamma-one)) ! Derivatives of the normal heat flux: qn = qx*njk(ix) + qy*njk(iy) + qz*njk(iz) dqn_drhoL = dqx_drhoL*njk(ix) + dqy_drhoL*njk(iy) + dqz_drhoL*njk(iz) dqn_dpL = dqx_dpL*njk(ix) + dqy_dpL*njk(iy) + dqz_dpL*njk(iz) !---------------------------------------------------------------------- ! 1.4 Compute the flux Jacbian with respect to the primitive variables: ! -> dFndwL ! 1st row: d(zero)/dW=zero (no viscous term in the continuity equation) dFndwL(1,1) = zero dFndwL(1,2) = zero dFndwL(1,3) = zero dFndwL(1,4) = zero dFndwL(1,5) = zero ! 2nd row: d(-taunx)/dW; taunx = tauxx*nx + tauxy*ny + tauxz*nz dFndwL(2,1) = - dtaunx_drhoL dFndwL(2,2) = - dtaunx_duL dFndwL(2,3) = - dtaunx_dvL dFndwL(2,4) = - dtaunx_dwL dFndwL(2,5) = - dtaunx_dpL ! 3rd row: d(-tauny)/dW; tauny = tauyx*nx + tauyy*ny + tauyz*nz dFndwL(3,1) = - dtauny_drhoL dFndwL(3,2) = - dtauny_duL dFndwL(3,3) = - dtauny_dvL dFndwL(3,4) = - dtauny_dwL dFndwL(3,5) = - dtauny_dpL ! 4th row: d(-taunz)/dW; taunz = tauzx*nx + tauzy*ny + tauzz*nz dFndwL(4,1) = - dtaunz_drhoL dFndwL(4,2) = - dtaunz_duL dFndwL(4,3) = - dtaunz_dvL dFndwL(4,4) = - dtaunz_dwL dFndwL(4,5) = - dtaunz_dpL ! 5th row: d(-taunv + qn)/dW; taunv = taunx*u + tauny*v + taunz*nz; qn=qx*nx+qy*ny+qz*nz dFndwL(5,1) = - dtaunv_drhoL + dqn_drhoL dFndwL(5,2) = - dtaunv_duL dFndwL(5,3) = - dtaunv_dvL dFndwL(5,4) = - dtaunv_dwL dFndwL(5,5) = - dtaunv_dpL + dqn_dpL !---------------------------------------------------------------------- ! 1.5 Compute the flux Jacbian with respect to the conservative variables: ! -> dFnducL = dFndwL*(dWL/ducL) ! ! dW/duc(1,:) = ( 1, 0, 0, 0, 0 ) ! dW/duc(2,:) = (-u/rho, 1/rho, 0, 0, 0 ) ! dW/duc(3,:) = (-v/rho, 0, 1/rho, 0, 0 ) ! dW/duc(4,:) = (-w/rho, 0, 0, 1/rho, 0 ) ! dW/duc(5,:) = ( q^2/2, -u, -v, -w, 1 )*(gamma-1) ! 1st row dFnducL(1,1) = zero dFnducL(1,2) = zero dFnducL(1,3) = zero dFnducL(1,4) = zero dFnducL(1,5) = zero ! 2nd row dFnducL(2,1) = dFndwL(2,1) & + dFndwL(2,2) * (-uL/rhoL) & + dFndwL(2,3) * (-vL/rhoL) & + dFndwL(2,4) * (-wL/rhoL) & + dFndwL(2,5) * half*(gamma-one)*q2L dFnducL(2,2) = dFndwL(2,2) * (one/rhoL) + dFndwL(2,5) * (-(gamma-one)*uL) dFnducL(2,3) = dFndwL(2,3) * (one/rhoL) + dFndwL(2,5) * (-(gamma-one)*vL) dFnducL(2,4) = dFndwL(2,4) * (one/rhoL) + dFndwL(2,5) * (-(gamma-one)*wL) dFnducL(2,5) = dFndwL(2,5) * (gamma-one) ! 3rd row dFnducL(3,1) = dFndwL(3,1) & + dFndwL(3,2) * (-uL/rhoL) & + dFndwL(3,3) * (-vL/rhoL) & + dFndwL(3,4) * (-wL/rhoL) & + dFndwL(3,5) * half*(gamma-one)*q2L dFnducL(3,2) = dFndwL(3,2) * (one/rhoL) + dFndwL(3,5) * (-(gamma-one)*uL) dFnducL(3,3) = dFndwL(3,3) * (one/rhoL) + dFndwL(3,5) * (-(gamma-one)*vL) dFnducL(3,4) = dFndwL(3,4) * (one/rhoL) + dFndwL(3,5) * (-(gamma-one)*wL) dFnducL(3,5) = dFndwL(3,5) * (gamma-one) ! 4th row dFnducL(4,1) = dFndwL(4,1) & + dFndwL(4,2) * (-uL/rhoL) & + dFndwL(4,3) * (-vL/rhoL) & + dFndwL(4,4) * (-wL/rhoL) & + dFndwL(4,5) * half*(gamma-one)*q2L dFnducL(4,2) = dFndwL(4,2) * (one/rhoL) + dFndwL(4,5) * (-(gamma-one)*uL) dFnducL(4,3) = dFndwL(4,3) * (one/rhoL) + dFndwL(4,5) * (-(gamma-one)*vL) dFnducL(4,4) = dFndwL(4,4) * (one/rhoL) + dFndwL(4,5) * (-(gamma-one)*wL) dFnducL(4,5) = dFndwL(4,5) * (gamma-one) ! 5th row dFnducL(5,1) = dFndwL(5,1) & + dFndwL(5,2) * (-uL/rhoL) & + dFndwL(5,3) * (-vL/rhoL) & + dFndwL(5,4) * (-wL/rhoL) & + dFndwL(5,5) * half*(gamma-one)*q2L dFnducL(5,2) = dFndwL(5,2) * (one/rhoL) + dFndwL(5,5) * (-(gamma-one)*uL) dFnducL(5,3) = dFndwL(5,3) * (one/rhoL) + dFndwL(5,5) * (-(gamma-one)*vL) dFnducL(5,4) = dFndwL(5,4) * (one/rhoL) + dFndwL(5,5) * (-(gamma-one)*wL) dFnducL(5,5) = dFndwL(5,5) * (gamma-one) !-------------------------------------------------------------------------------- ! Part 2. Compute dFn_alpha/ducR ! ! dFn_alpha/ducR = [ zero, -d(taunx)/ducR, -d(tauny)/ducR, -d(taunz)/ducR, d(-taunv+qn)/ducR] ! = (dFn_alpha/dwR)*(dwR/ducR) ! ! where wR = [rhoR, uR, vR, wR, pR] - primitive variables. ! ! So, we proceed as follows: ! ! 2.1 Compute the derivatives of the viscosity: dmu/drhoR, dmu/dpR ! 2.2 Compute d(taunx)/dwR, -d(tauny)/dwR, -d(taunz)/dwR, -d(taunv)/dwR ! 2.3 Compute d(qn)/dwR ! 2.4 Compute the flux Jacbian with respect to the primitive variables: ! -> dFndwR ! 2.5 Compute the flux Jacbian with respect to the conservative variables: ! -> dFnducR = dFndwR*(dWR/ducR) !---------------------------------------------------------------------- ! 2.1 Compute the derivatives of the viscosity: dmu/drhoR, dmu/dpR ! dmu/dTR = d( (one+C/T_inf)/(T+C/T_inf)*T**(three*half) * M_inf/Re_inf )/dTR dmu_dTR = half*mu/(T+C/T_inf)*(one + three*(C/T_inf)/T) * ( half ) ! Note: dmu/dW = [dmu/drho, dmu/du, dmu/dv, dmu/dw, dmu/dp ] = [dmu/drho, 0, 0, 0, dmu/dp ] ! So, we only need dmu/drho and dmu/dp. ! dmu/drhoR = dmu/dTR*(dTR/drhoR); dT/drho = -gamma*p/rho**2 because T = gamma*p/rho. dmu_drhoR = dmu_dTR * (-gamma*pR/rhoR**2) ! dmu/dpR = dmu/dTR*(dTR/dpR); dT/dp = gamma/rho because T = gamma*p/rho. dmu_dpR = dmu_dTR * ( gamma/rhoR ) !---------------------------------------------------------------------- ! 2.2 Compute d(taunx)/dwR, -d(tauny)/dwR, -d(taunz)/dwR, -d(taunv)/dwR ! ! Derivatives of the interface velocity gradients: ! ! Interface gradients are given by ! ! grad_u = (Average LSQ gradient of u) + alpha/Lr*(uR - uL)*njk ! grad_v = (Average LSQ gradient of v) + alpha/Lr*(vR - vL)*njk ! grad_w = (Average LSQ gradient of w) + alpha/Lr*(wR - wL)*njk ! ! but for simplicity and compantness, we ignore the average LSQ gradients. ! Note: If the average LSQ gradients are computed on a compact stencil, ! we may add their contributions. dux_duR = alpha/Lr*(one-zero)*njk(ix) duy_duR = alpha/Lr*(one-zero)*njk(iy) duz_duR = alpha/Lr*(one-zero)*njk(iz) dvx_dvR = alpha/Lr*(one-zero)*njk(ix) dvy_dvR = alpha/Lr*(one-zero)*njk(iy) dvz_dvR = alpha/Lr*(one-zero)*njk(iz) dwx_dwR = alpha/Lr*(one-zero)*njk(ix) dwy_dwR = alpha/Lr*(one-zero)*njk(iy) dwz_dwR = alpha/Lr*(one-zero)*njk(iz) ! Derivatives of the stress tensor. ! Note: many terms are zero, e.g., dux/dvR=0, dwz/dvR=0, etc. ! dsxx_duR = four_third*dux_duR - two_third*( dvy_duR + dwz_duR ) ! dsyy_duR = four_third*dvy_duR - two_third*( dux_duR + dwz_duR ) ! dszz_duR = four_third*dwz_duR - two_third*( dux_duR + dvy_duR ) ! Ignoring terms that are zero, we obtain the following. dsxx_duR = four_third*dux_duR dsyy_duR = - two_third*( dux_duR ) dszz_duR = - two_third*( dux_duR ) ! dsxy_duR = duy_duR + dvx_duR ! dsxz_duR = duz_duR + dwx_duR ! dsyz_duR = dvz_duR + dwy_duR ! Ignoring terms that are zero, we obtain the following. dsxy_duR = duy_duR dsxz_duR = duz_duR dsyz_duR = zero ! dsxx_dvR = four_third*dux_dvR - two_third*( dvy_dvR + dwz_dvR ) ! dsyy_dvR = four_third*dvy_dvR - two_third*( dux_dvR + dwz_dvR ) ! dszz_dvR = four_third*dwz_dvR - two_third*( dux_dvR + dvy_dvR ) ! Ignoring terms that are zero, we obtain the following. dsxx_dvR = - two_third*( dvy_dvR ) dsyy_dvR = four_third*dvy_dvR dszz_dvR = - two_third*( dvy_dvR ) ! dsxy_dvR = duy_dvR + dvx_dvR ! dsxz_dvR = duz_dvR + dwx_dvR ! dsyz_dvR = dvz_dvR + dwy_dvR ! Ignoring terms that are zero, we obtain the following. dsxy_dvR = dvx_dvR dsxz_dvR = zero dsyz_dvR = dvz_dvR ! dsxx_dwR = four_third*dux_dwR - two_third*( dvy_dwR + dwz_dwR ) ! dsyy_dwR = four_third*dvy_dwR - two_third*( dux_dwR + dwz_dwR ) ! dszz_dwR = four_third*dwz_dwR - two_third*( dux_dwR + dvy_dwR ) ! Ignoring terms that are zero, we obtain the following. dsxx_dwR = - two_third*( dwz_dwR ) dsyy_dwR = - two_third*( dwz_dwR ) dszz_dwR = four_third*dwz_dwR ! dsxy_dwR = duy_dwR + dvx_dwR ! dsxz_dwR = duz_dwR + dwx_dwR ! dsyz_dwR = dvz_dwR + dwy_dwR ! Ignoring terms that are zero, we obtain the following. dsxy_dwR = zero dsxz_dwR = dwx_dwR dsyz_dwR = dwy_dwR dsnx_duR = dsxx_duR*njk(ix) + dsxy_duR*njk(iy) + dsxz_duR*njk(iz) dsny_duR = dsxy_duR*njk(ix) + dsyy_duR*njk(iy) + dsyz_duR*njk(iz) dsnz_duR = dsxz_duR*njk(ix) + dsyz_duR*njk(iy) + dszz_duR*njk(iz) dsnx_dvR = dsxx_dvR*njk(ix) + dsxy_dvR*njk(iy) + dsxz_dvR*njk(iz) dsny_dvR = dsxy_dvR*njk(ix) + dsyy_dvR*njk(iy) + dsyz_dvR*njk(iz) dsnz_dvR = dsxz_dvR*njk(ix) + dsyz_dvR*njk(iy) + dszz_dvR*njk(iz) dsnx_dwR = dsxx_dwR*njk(ix) + dsxy_dwR*njk(iy) + dsxz_dwR*njk(iz) dsny_dwR = dsxy_dwR*njk(ix) + dsyy_dwR*njk(iy) + dsyz_dwR*njk(iz) dsnz_dwR = dsxz_dwR*njk(ix) + dsyz_dwR*njk(iy) + dszz_dwR*njk(iz) ! Derivatives of the viscous stresses. ! Note: mu does not depend on the velocity. dtauxx_duR = mu*dsxx_duR dtauxx_dvR = mu*dsxx_dvR dtauxx_dwR = mu*dsxx_dwR dtauyy_duR = mu*dsyy_duR dtauyy_dvR = mu*dsyy_dvR dtauyy_dwR = mu*dsyy_dwR dtauzz_duR = mu*dszz_duR dtauzz_dvR = mu*dszz_dvR dtauzz_dwR = mu*dszz_dwR dtauxy_duR = mu*dsxy_duR dtauxy_dvR = mu*dsxy_dvR dtauxy_dwR = mu*dsxy_dwR dtauxz_duR = mu*dsxz_duR dtauxz_dvR = mu*dsxz_dvR dtauxz_dwR = mu*dsxz_dwR dtauyz_duR = mu*dsyz_duR dtauyz_dvR = mu*dsyz_dvR dtauyz_dwR = mu*dsyz_dwR ! Derivatives of the viscous stresses in the normal direction. ! Derivatives of taunx dtaunx_drhoR = dmu_drhoR * snx dtaunx_duR = dtauxx_duR*njk(ix) + dtauxy_duR*njk(iy) + dtauxz_duR*njk(iz) dtaunx_dvR = dtauxx_dvR*njk(ix) + dtauxy_dvR*njk(iy) + dtauxz_dvR*njk(iz) dtaunx_dwR = dtauxx_dwR*njk(ix) + dtauxy_dwR*njk(iy) + dtauxz_dwR*njk(iz) dtaunx_dpR = dmu_dpR * snx ! Derivatives of tauny dtauny_drhoR = dmu_drhoR * sny dtauny_duR = dtauxy_duR*njk(ix) + dtauyy_duR*njk(iy) + dtauyz_duR*njk(iz) dtauny_dvR = dtauxy_dvR*njk(ix) + dtauyy_dvR*njk(iy) + dtauyz_dvR*njk(iz) dtauny_dwR = dtauxy_dwR*njk(ix) + dtauyy_dwR*njk(iy) + dtauyz_dwR*njk(iz) dtauny_dpR = dmu_dpR * sny ! Derivatives of taunz dtaunz_drhoR = dmu_drhoR * snz dtaunz_duR = dtauxz_duR*njk(ix) + dtauyz_duR*njk(iy) + dtauzz_duR*njk(iz) dtaunz_dvR = dtauxz_dvR*njk(ix) + dtauyz_dvR*njk(iy) + dtauzz_dvR*njk(iz) dtaunz_dwR = dtauxz_dwR*njk(ix) + dtauyz_dwR*njk(iy) + dtauzz_dwR*njk(iz) dtaunz_dpR = dmu_dpR * snz ! Derivatives of taunv = mu* ( snx*u + sny*v + snz*w ) dtaunv_drhoR = dmu_drhoR * snv dtaunv_duR = mu * ( snx*(half) + dsnx_duR*u + dsny_duR*v + dsnz_duR*w ) dtaunv_dvR = mu * ( sny*(half) + dsnx_dvR*u + dsny_dvR*v + dsnz_dvR*w ) dtaunv_dwR = mu * ( snz*(half) + dsnx_dwR*u + dsny_dwR*v + dsnz_dwR*w ) dtaunv_dpR = dmu_dpR * snv !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! 2.3 Compute d(qn)/dwR ! First compute the derivatives of the interface gradients: ! ! rho = half*(rhoR + rhoR) ! a2 = gamma*half*(pR + pR)/rho ! grad_rho = (Average LSQ gradient of rho) + alpha/Lr*(rhoR - rhoL)*njk ! grad_p = (Average LSQ gradient of p) + alpha/Lr*( pR - pL )*njk ! grad_T = ( gamma*grad_p - a2*grad_rho) /rho ! ! Again, we ignore the average LSQ gradients. drhox_drhoR = alpha/Lr*(one-zero)*njk(ix) ! Average LSQ gradient ignored drhoy_drhoR = alpha/Lr*(one-zero)*njk(iy) ! Average LSQ gradient ignored drhoz_drhoR = alpha/Lr*(one-zero)*njk(iz) ! Average LSQ gradient ignored dTx_drhoR = (- gamma*grad_p(ix)/rho**2 & + two* gamma*half*(pR+pL)/rho**3*grad_rho(ix) ) * half & - a2/rho*drhox_drhoR dTy_drhoR = (- gamma*grad_p(iy)/rho**2 & + two* gamma*half*(pR+pL)/rho**3*grad_rho(iy) ) * half & - a2/rho*drhoy_drhoR dTz_drhoR = (- gamma*grad_p(iz)/rho**2 & + two* gamma*half*(pR+pL)/rho**3*grad_rho(iz) ) * half & - a2/rho*drhoz_drhoR dpx_dpR = alpha/Lr*(one-zero)*njk(ix) ! Average LSQ gradient ignored dpy_dpR = alpha/Lr*(one-zero)*njk(iy) ! Average LSQ gradient ignored dpz_dpR = alpha/Lr*(one-zero)*njk(iz) ! Average LSQ gradient ignored dTx_dpR = gamma*dpx_dpR/rho - ( gamma/rho**2*grad_rho(ix) ) * half dTy_dpR = gamma*dpy_dpR/rho - ( gamma/rho**2*grad_rho(iy) ) * half dTz_dpR = gamma*dpz_dpR/rho - ( gamma/rho**2*grad_rho(iz) ) * half ! Derivatives of the heat fluxes: ! ! qx = - mu*grad_T(ix)/(Prandtl*(gamma-one)) ! qy = - mu*grad_T(iy)/(Prandtl*(gamma-one)) ! qz = - mu*grad_T(iz)/(Prandtl*(gamma-one)) dqx_drhoR = - ( dmu_drhoR*grad_T(ix) + mu*dTx_drhoR )/(Prandtl*(gamma-one)) dqx_dpR = - ( dmu_dpR *grad_T(ix) + mu*dTx_dpR )/(Prandtl*(gamma-one)) dqy_drhoR = - ( dmu_drhoR*grad_T(iy) + mu*dTy_drhoR )/(Prandtl*(gamma-one)) dqy_dpR = - ( dmu_dpR *grad_T(iy) + mu*dTy_dpR )/(Prandtl*(gamma-one)) dqz_drhoR = - ( dmu_drhoR*grad_T(iz) + mu*dTz_drhoR )/(Prandtl*(gamma-one)) dqz_dpR = - ( dmu_dpR *grad_T(iz) + mu*dTz_dpR )/(Prandtl*(gamma-one)) ! Derivatives of the normal heat flux: qn = qx*njk(ix) + qy*njk(iy) + qz*njk(iz) dqn_drhoR = dqx_drhoR*njk(ix) + dqy_drhoR*njk(iy) + dqz_drhoR*njk(iz) dqn_dpR = dqx_dpR*njk(ix) + dqy_dpR*njk(iy) + dqz_dpR*njk(iz) !---------------------------------------------------------------------- ! 2.4 Compute the flux Jacbian with respect to the primitive variables: ! -> dFndwR ! 1st row: d(zero)/dW=zero (no viscous term in the continuity equation) dFndwR(1,1) = zero dFndwR(1,2) = zero dFndwR(1,3) = zero dFndwR(1,4) = zero dFndwR(1,5) = zero ! 2nd row: d(-taunx)/dW; taunx = tauxx*nx + tauxy*ny + tauxz*nz dFndwR(2,1) = - dtaunx_drhoR dFndwR(2,2) = - dtaunx_duR dFndwR(2,3) = - dtaunx_dvR dFndwR(2,4) = - dtaunx_dwR dFndwR(2,5) = - dtaunx_dpR ! 3rd row: d(-tauny)/dW; tauny = tauyx*nx + tauyy*ny + tauyz*nz dFndwR(3,1) = - dtauny_drhoR dFndwR(3,2) = - dtauny_duR dFndwR(3,3) = - dtauny_dvR dFndwR(3,4) = - dtauny_dwR dFndwR(3,5) = - dtauny_dpR ! 4th row: d(-taunz)/dW; taunz = tauzx*nx + tauzy*ny + tauzz*nz dFndwR(4,1) = - dtaunz_drhoR dFndwR(4,2) = - dtaunz_duR dFndwR(4,3) = - dtaunz_dvR dFndwR(4,4) = - dtaunz_dwR dFndwR(4,5) = - dtaunz_dpR ! 5th row: d(-taunv + qn)/dW; taunv = taunx*u + tauny*v + taunz*nz; qn=qx*nx+qy*ny+qz*nz dFndwR(5,1) = - dtaunv_drhoR + dqn_drhoR dFndwR(5,2) = - dtaunv_duR dFndwR(5,3) = - dtaunv_dvR dFndwR(5,4) = - dtaunv_dwR dFndwR(5,5) = - dtaunv_dpR + dqn_dpR !---------------------------------------------------------------------- ! 2.5 Compute the flux Jacbian with respect to the conservative variables: ! -> dFnducR = dFndwR*(dWR/ducR) ! ! dW/duc(1,:) = ( 1, 0, 0, 0, 0 ) ! dW/duc(2,:) = (-u/rho, 1/rho, 0, 0, 0 ) ! dW/duc(3,:) = (-v/rho, 0, 1/rho, 0, 0 ) ! dW/duc(4,:) = (-w/rho, 0, 0, 1/rho, 0 ) ! dW/duc(5,:) = ( q^2/2, -u, -v, -w, 1 )*(gamma-1) ! 1st row dFnducR(1,1) = zero dFnducR(1,2) = zero dFnducR(1,3) = zero dFnducR(1,4) = zero dFnducR(1,5) = zero ! 2nd row dFnducR(2,1) = dFndwR(2,1) & + dFndwR(2,2) * (-uR/rhoR) & + dFndwR(2,3) * (-vR/rhoR) & + dFndwR(2,4) * (-wR/rhoR) & + dFndwR(2,5) * half*(gamma-one)*q2R dFnducR(2,2) = dFndwR(2,2) * (one/rhoR) + dFndwR(2,5) * (-(gamma-one)*uR) dFnducR(2,3) = dFndwR(2,3) * (one/rhoR) + dFndwR(2,5) * (-(gamma-one)*vR) dFnducR(2,4) = dFndwR(2,4) * (one/rhoR) + dFndwR(2,5) * (-(gamma-one)*wR) dFnducR(2,5) = dFndwR(2,5) * (gamma-one) ! 3rd row dFnducR(3,1) = dFndwR(3,1) & + dFndwR(3,2) * (-uR/rhoR) & + dFndwR(3,3) * (-vR/rhoR) & + dFndwR(3,4) * (-wR/rhoR) & + dFndwR(3,5) * half*(gamma-one)*q2R dFnducR(3,2) = dFndwR(3,2) * (one/rhoR) + dFndwR(3,5) * (-(gamma-one)*uR) dFnducR(3,3) = dFndwR(3,3) * (one/rhoR) + dFndwR(3,5) * (-(gamma-one)*vR) dFnducR(3,4) = dFndwR(3,4) * (one/rhoR) + dFndwR(3,5) * (-(gamma-one)*wR) dFnducR(3,5) = dFndwR(3,5) * (gamma-one) ! 4th row dFnducR(4,1) = dFndwR(4,1) & + dFndwR(4,2) * (-uR/rhoR) & + dFndwR(4,3) * (-vR/rhoR) & + dFndwR(4,4) * (-wR/rhoR) & + dFndwR(4,5) * half*(gamma-one)*q2R dFnducR(4,2) = dFndwR(4,2) * (one/rhoR) + dFndwR(4,5) * (-(gamma-one)*uR) dFnducR(4,3) = dFndwR(4,3) * (one/rhoR) + dFndwR(4,5) * (-(gamma-one)*vR) dFnducR(4,4) = dFndwR(4,4) * (one/rhoR) + dFndwR(4,5) * (-(gamma-one)*wR) dFnducR(4,5) = dFndwR(4,5) * (gamma-one) ! 5th row dFnducR(5,1) = dFndwR(5,1) & + dFndwR(5,2) * (-uR/rhoR) & + dFndwR(5,3) * (-vR/rhoR) & + dFndwR(5,4) * (-wR/rhoR) & + dFndwR(5,5) * half*(gamma-one)*q2R dFnducR(5,2) = dFndwR(5,2) * (one/rhoR) + dFndwR(5,5) * (-(gamma-one)*uR) dFnducR(5,3) = dFndwR(5,3) * (one/rhoR) + dFndwR(5,5) * (-(gamma-one)*vR) dFnducR(5,4) = dFndwR(5,4) * (one/rhoR) + dFndwR(5,5) * (-(gamma-one)*wR) dFnducR(5,5) = dFndwR(5,5) * (gamma-one) end subroutine viscous_alpha_flux_jacobian !-------------------------------------------------------------------------------- end program test_flux_jacobians