C RACREA SOURCE BECC 09/12/07 21:16:02 6579 subroutine racrea(nordpo, Tmaxcv, & Rgas_l, acv_l, Rg_rss, acv_z, Rgas_r, acv_r, & q, k0, & rho_l, p_l, u_l, T_l, et_l, gam_l, & rho_r, p_r, u_r, T_r, et_r, gam_r, & rho_ls, p_ls, u_ls, T_ls, c_ls, d_ls, & rho_z, p_rss, u_rss, T_g_z, Dsw, & rho_rs, p_rs, u_rs, T_rs, c_rs, d_rs, & logan, lognc) C c************************************************************************ c c project : c c name : racrea c c description : computation of the shock-shock field-by-field c decomposition of the reactive Riemann problem for c thermally perfect gases. C C We have to distinguish between 5 different states C C l = left state; C burnt state (csi = csimax) C | C | Left GNL wave C | C ls = left star state; C connected to l by a GNL wave C (csi = csimax) C | C | Contact discontinuity C | C rss = right star star state; C connected to rs by a reactive wave C (csi = csimax) C | C | Reactive shock (weak defl.) C | C rss = right star state; C connected to r by a GNL wave C (csi = csimin) C | C | Right GNL wave C | C r = right state; C unburnt state (csi = csimin) c c language : fortran 77 c c author : a. beccantini den/dm2s/sfme/ltmf c c************************************************************************ c c called by : fcolre c c called program : funco1, funrea c c************************************************************************ c c c**** input c c nordpo = polynomials degree of cv(t) c tmaxcv = threshold temperature for cv(t) c rgas_l = gas constant (J/kg/K)in l (and in ls) c acv_l = cv_l = \sum_{i=0,nordpo} (acv_l{i} T^i) c rgas_rss = gas constant in rss c acv_z = cv_rss = \sum_{i=0,nordpo} (acv_z{i} T^i) c rgas_r = gas constant r (and in rs) c acv_r = cv_r = \sum_{i=0,nordpo} (acv_r{i} T^i) c c q = released heat per mass unit from rs to rss c k0 = fundamental flame speed c C rho_l = density in l C p_l = pressure in c u_l = velocity in l C T_l = temperature in l (redundant with rho_l and c p_l) c et_l = (sensible) internal energy in l (redundant c with rho_l and p_l) c gam_l = specific heat ratio in l (redundant c with rho_l and p_l) c C rho_r = density in r C p_r = pressure in r c u_r = velocity in r C T_r = temperature in r (redundant with rho_r and c p_r) c et_r = (sensible) internal energy in r (redundant c with rho_r and p_r) c gam_r = specific heat ratio in r (redundant c with rho_r and p_r) c c**** output c c C rho_ls = density in ls C p_ls = pressure in ls c u_ls = velocity in ls c T_ls = temperature in ls c c_ls = sound speed in ls c C rho_z = density in rss C p_rss = pressure in rss c u_rss = velocity in rss c T_rss = temperature in rss c Dsw = reactive shock velocity in the laboratory C frame C C rho_rs = density in rs C p_rs = pressure in rs c u_rs = velocity in rs c T_rs = temperature in rs c c_rs = sound speed in rs c c logan = an anomaly is detected c lognc = the apporach has some problems. c c************************************************************************ c c 07/12/2009 Created. c c implicit none integer nmax, nordpo, icon parameter (nmax=50) real*8 eps, opeps, zero parameter (eps=1.0d-3,opeps=(1.d0+eps),zero=1.0d-12) real*8 Tmaxcv, Rgas_l, Rg_rss, Rgas_r & , acv_l(0:nordpo), acv_r(0:nordpo), acv_z(0:nordpo) & , q, k0 & , rho_l, p_l, u_l, T_l, et_l, gam_l & , rho_r, p_r, u_r, T_r, et_r, gam_r & , rho_ls, p_ls, u_ls, et_ls, c_ls & , rho_z, p_rss, u_rss, et_rss, Dsw & , rho_rs, p_rs, u_rs, et_rs, c_rs real*8 Tref, pref, cref, uref, c_l, c_r, u_rm real*8 u_ls1, p_ls1, T_ls, T_ls2 C real*8 du_ls_ov_dT_ls, dp_ls_ov_dT_ls real*8 der3, der4 real*8 u_rss1, p_rss1, T_rs, T_rs2 real*8 cell1, cell2, cell3, cell4 & , cell5, cell6, cell7, cell8, cell9 & , cv_rs, cv_ls real*8 T_g_z, T_g_z1 C real*8 du_rss_ov_dT_rs, dp_rss_ov_dT_rs real*8 der1, der2 real*8 detj, dT_ls, dT_rs real*8 umx_ls, upx_ls, rc1_ls, d_ls real*8 umx_rs, upx_rs, rc1_rs, d_rs logical lognc, logan, logeq, loe_ls, loe_rs logical logdeb parameter (logdeb = .false.) CC CC**** We write the inputs CC C write(*,*) 'racrea.f' C write(*,'(I2,E12.4)') nordpo, tmaxcv C write(*,'(6E12.4)') rgas_l, (acv_l(icon), icon=0,nordpo) C write(*,'(6E12.4)') rg_rss, (acv_z(icon), icon=0,nordpo) C write(*,'(6E12.4)') rgas_r, (acv_r(icon), icon=0,nordpo) C write(*,'(2E12.4)') q, k0 C write(*,'(6E12.4)') rho_l, p_l, u_l, T_l, et_l, gam_l C write(*,'(6E12.4)') rho_r, p_r, u_r, T_r, et_r, gam_r C C**** Reference scales for a compressible solver. C For a low Mach number solver, the cut off speed is better than C the sound speed C Tref = 0.5d0 * ((T_l + T_r) + abs (T_l - T_r)) C WARNING: c_l and c_r are input of funco1 C and are not only used to compute the reference speed c_l = (gam_l * p_l / rho_l) ** 0.5d0 c_r = (gam_r * p_r / rho_r) ** 0.5d0 cref = 0.5d0 * ((c_l + c_r) + abs (c_l - c_r)) uref = cref C pref = 0.5d0 * ((p_l + p_r) + abs (p_l - p_r)) C C**** Redefinition of the right speed u_rm. C u_rm, used in the function funco1, which is designed to evaluate C a left GNL wave and not a right one. C Note that, after calling funco1 to compute the rs state, the C speed u_rs must be multiplied by -1. C u_rm = -1.0D0 * u_r C C**** Initialization of the states ls, rs, rss C C Starting point ls = l C u_ls = u_l p_ls = p_l T_ls = T_l C C Starting point rs = r C u_rs = u_r p_rs = p_r T_rs = T_r C We need to define *_rs since used in funrea C We initialize T_g_z = T_l since rss is a burnt state, like the C left one. rho_rs = rho_r et_rs = et_r T_g_z = T_l C C Starting point rss C NB Note that rss state is here considered function of T_rs only C (and of the guess temperature for T_rss, here called T_g_z, C in the case of thermally perfect gases). C C call funrea((nordpo+1), Tmaxcv, Rg_rss, acv_z, & rho_rs, p_rs, u_rs, et_rs, q, k0, T_g_z, & rho_z, p_rss, u_rss, Dsw, et_rss, logan) C if (logan) then write(*,'(/A8)') 'racrea.f' write(*,*) 'icon =', 0 C write(*,*) 'rho_l, p_l, u_l, T_l, et_l, gam_l' C write(*,'(6E16.6)') rho_l, p_l, u_l, T_l, et_l, gam_l C write(*,*) 'rho_r, p_r, u_r, T_r, et_r, gam_r' C write(*,'(6E16.6)') rho_r, p_r, u_r, T_r, et_r, gam_r write(*,*) 'ANOMALY DETECTED 1' goto 9999 endif T_g_z = p_rss / (Rg_rss * rho_z) C C**** Iterative procedure to compute the intermediate states C We use a quasi-Newton approach, with numerical evaluation of the C derivatives, involving the two variables T_rs and T_ls C (temperature of the states rs and ls). C The temperature for the rss state is C evaluated at the end of each iteration and is exact in the C particular case of calorically perfect gases. For the moment we C do not involve this variable in the quasi-Newton approach. C Summarizing we have to compute C C F = (f_1, f_2)^T = 0 C C i.e. C C f_1(T_ls, T_rs; T_rss) = u_ls(T_ls, T_rs; T_rss) - C u_rss(T_ls, T_rs; T_rss) = 0 C f_2(T_ls, T_rs; T_rss) = p_ls(T_ls, T_rs; T_rss) - C p_rss(T_ls, T_rs; T_rss) = 0 C C At each iteration we compute X = (T_ls, T_rs)^T by solving C C F(X_1;T_rss1) + F'(X_1;T_rss1) (X - X_1) = 0 C C (X - X_1) = - (F'(X_1;T_rss1))^{-1} F(X_1;T_rss1) C C Then C T_rss = fun(X,T_rss1) C do icon = 1, nmax, 1 C C Linearisation of the problem around ls1 = ls C rss1 = rss C rs1 = rs C We need to save the old values of u_rs and p_rs. C u_ls1 = u_ls p_ls1 = p_ls C C Computation of of the derivatives C dp_ls/dT_ls and du_ls/dT_ls C T_ls2 = T_ls * opeps call funco1((nordpo+1), Tmaxcv, Rgas_l, acv_l, & rho_l, T_l, (Rgas_l * T_l), et_l, u_l, gam_l, c_l, & T_ls2, cell1, cv_ls, rho_ls, p_ls, et_ls, u_ls, & cell2, cell3, cell4, & cell5, cell6, cell7, & logeq, cell8, cell9) C State _ls output of funco1 only used to evaluate derivatives! der3 = (u_ls - u_ls1) / (T_ls * eps) der4 = (p_ls - p_ls1) / (T_ls * eps) C write(*,*) 'T_ls , eps', T_ls, eps C write(*,*) 'der3, der4 = ', der3, der4 C C We do not need save the variables _rs in _rs1. C We need save the old values of u_rss1 and p_rss1 and the C guess temperature for rss. C T_g_z1 = T_g_z u_rss1 = u_rss p_rss1 = p_rss C C Computation of of the derivatives C dp_rss/dT_rs2 and du_rss/dT_rs2 C T_rs2 = T_rs * opeps call funco1((nordpo+1), Tmaxcv, Rgas_r, acv_r, & rho_r, T_r, (Rgas_r * T_r), et_r, u_rm, gam_r, c_r, & T_rs2, cell1, cv_rs, rho_rs, p_rs, et_rs, u_rs, & cell2, cell3, cell4, & cell5, cell6, cell7, & logeq, cell8, cell9) u_rs = -1.0D0 * u_rs C -1 * necessary since funco1 is designed for left GNL waves. C State _rs output of funco1, namely C rho_rs, p_rs, et_rs, u_rs C only used to evaluate derivatives as input of funrea! call funrea((nordpo+1), Tmaxcv, Rg_rss, acv_z, & rho_rs, p_rs, u_rs, et_rs, q, k0, T_g_z, & rho_z, p_rss, u_rss, Dsw, et_rss, logan) if (logan) then write(*,'(/A8)') 'racrea.f' write(*,*) 'icon =', icon write(*,*) 'ANOMALY DETECTED 2' goto 9999 endif C State _rss output of funrea only used to evaluate derivatives! der1 = (u_rss - u_rss1) / (T_rs * eps) der2 = (p_rss - p_rss1) / (T_rs * eps) C write(*,*) der1, der2 C C | || | | | C |der3 -der1||dT_ls| | u_ls1 - u_rss1 | C | || | = - | | C |der4 -der2||dT_rs| | p_ls1 - p_rss1 | C | || | | | C C | | | || | C |dT_ls| 1 |der2 -der1|| u_ls1 - u_rss1 | C | | = ------ | || | C |dT_rs| det(M) |der4 -der3 || p_ls1 - p_rss1 | C | | | || | C detj = (der4 * der1) & - (der3 * der2) C dT_ls = (der2 * (u_ls1 - u_rss1)) - & (der1 * (p_ls1 - p_rss1)) dT_ls = dT_ls / detj dT_rs = (der4 * (u_ls1 - u_rss1)) - & (der3 * (p_ls1 - p_rss1)) dT_rs = dT_rs / detj C We limited the jump if (abs (dT_ls) .gt. (0.5d0 * T_ls)) then dT_ls = sign(dT_ls, (0.5d0 * T_ls)) endif if (abs (dT_rs) .gt. (0.5d0 * T_rs)) then dT_rs = sign(dT_rs, (0.5d0 * T_rs)) endif C T_ls = T_ls + dT_ls T_rs = T_rs + dT_rs C C Computation of the states. C ls call funco1((nordpo+1), Tmaxcv, Rgas_l, acv_l, & rho_l, T_l, (Rgas_l * T_l), et_l, u_l, gam_l, c_l, & T_ls, cell1, cv_ls, rho_ls, p_ls, et_ls, u_ls, & cell2, umx_ls, upx_ls, & rc1_ls, cell6, cell7, & loe_ls, cell8, cell9) C rs call funco1((nordpo+1), Tmaxcv, Rgas_r, acv_r, & rho_r, T_r, (Rgas_r * T_r), et_r, u_rm, gam_r, c_r, & T_rs, cell1, cv_rs, rho_rs, p_rs, et_rs, u_rs, & cell2, umx_rs, upx_rs, & rc1_rs, cell6, cell7, & loe_rs, cell8, cell9) u_rs = -1.0D0 * u_rs C rss call funrea((nordpo+1), Tmaxcv, Rg_rss, acv_z, & rho_rs, p_rs, u_rs, et_rs, q, k0, T_g_z, & rho_z, p_rss, u_rss, Dsw, et_rss, logan) if (logan) then write(*,'(/A8)') 'racrea.f' write(*,*) 'icon =', icon write(*,*) 'ANOMALY DETECTED 3' goto 9999 endif C C Computation of T_g_z C (T_rss = fun(X,T_rss1)). C T_g_z = p_rss / (Rg_rss * rho_z) C C write(*,*) icon, T_ls, T_rs, T_g_z C lognc = abs (dT_rs) .lt. (zero * tref) lognc = lognc .and. (abs (dT_ls) .lt. (zero * tref)) lognc = lognc .and. (abs (T_g_z - T_g_z1) .lt. & (zero * tref)) lognc = lognc .and. (abs(p_rss - p_ls) .lt. (zero * pref)) lognc = lognc .and. (abs(u_rss - u_ls) .lt. (zero * uref)) if (lognc) goto 9998 enddo 9998 continue C write(*,*) icon C C**** Sound speed and shock speed C c_ls = ((Rgas_l + cv_ls) * p_ls) / (cv_ls * rho_ls) c_ls = c_ls ** 0.5d0 c_rs = ((Rgas_r + cv_rs) * p_rs) / (cv_rs * rho_rs) c_rs = c_rs ** 0.5d0 lognc = .not. lognc if (T_rs .le. (1.00001D0 * T_r)) then C In case of rarefaction wave we take d_rs = d_r d_rs = u_r + c_r else d_rs = u_r + sqrt(rc1_rs / umx_rs / upx_rs) endif if (T_ls .le. (1.00001D0 * T_l)) then C In case of rarefaction wave we take d_ls = d_l d_ls = u_l - c_l else d_ls = u_l - sqrt(rc1_ls / umx_ls / upx_ls) endif if (lognc) then write(*,'(/A8)') 'racrea.f' write(*,*) 'Quasi-Newton method does not converge' endif if (logdeb) then write(*,*) 'Intersection ' write(*,*) rho_l, p_l, u_l, T_l write(*,*) rho_ls, p_ls, u_ls, T_ls, d_ls, loe_rs write(*,*) rho_z, p_rss, u_rss, T_g_z write(*,*) rho_rs, p_rs, u_rs, T_rs, d_rs, loe_rs write(*,*) rho_r, p_r, u_r, T_r endif C 9999 continue return end