C CONJP3 SOURCE CB215821 16/04/21 21:15:58 8920 SUBROUTINE CONJP3(JL,JR,WVEC_L,WVEC_R,NVECT,TVECT, & ga) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : CONJP3 C C DESCRIPTION : Voir KONJP2 C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF C C************************************************************************ c---------------------------------------------------------------------- c GENERAL DESCRIPTION: c This subroutine provides jacobian jtl at the interface c next to wall. c For full descriptions please see file 'conjp2.eso'. c---------------------------------------------------------------------- c INPUT: c c alpha -- parameter of the AUSM+ scheme in the Pressure function; c ( -3/4 <= alpha <= 3/16 ) c c beta -- parameter of the AUSM+ scheme in the Mach function; c ( -1/16 <= beta <= 1/2 ) c c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the c left cell; c c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the c right cell; c c nvect -- normal vector to the interface (2 components in 2D); c c tvect -- tangential vector to the interface; c c ga -- ratio of the specific heats (assumed constant) c---------------------------------------------------------------------- c c OUTPUT: c c jl -- jakobian matrix 4 by 4 - derivatives of the numerical c flux function with respect to the primitive variables c from the left cell; c c jr -- jakobian matrix 4 by 4 - derivatives of the numerical c flux function with respect to the conservative variables c from the right cell. c---------------------------------------------------------------------- IMPLICIT INTEGER(I-N) real*8 wvec_l(4),wvec_r(4) real*8 nvect(2),tvect(2) real*8 jl(4,4),jr(4,4) real*8 alpha,beta real*8 ga, gm1 real*8 n_x,n_y real*8 un_l, un_r real*8 ml,mr,Mplus,Mmin real*8 rold_l,uold_l,vold_l,pold_l,eold_l real*8 rold_r,uold_r,vold_r,pold_r,eold_r real*8 Pplus,Pmin real*8 temp_l,temp_r,brac_l,brac_r real*8 aleft, arigh, am real*8 damr_l,damr_r,damu_l,damu_r real*8 damv_l,damv_r,damp_l,damp_r real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r real*8 dmrr_l,dmrr_r,dmru_l,dmru_r real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r real*8 dpir_l,dpir_r,dpiu_l,dpiu_r real*8 dpiv_l,dpiv_r,dpip_l,dpip_r real*8 tcoef,bcoef integer i parameter(alpha = 0.1875D0, beta = 0.125D0) c----------------------- gm1=ga-1.0d0 c----------------------- n_x=nvect(1) n_y=nvect(2) c---------------------------- c----------------------- rold_l=wvec_l(1) uold_l=wvec_l(2) vold_l=wvec_l(3) pold_l=wvec_l(4) c----------------------- rold_r=wvec_r(1) uold_r=wvec_r(2) vold_r=wvec_r(3) pold_r=wvec_r(4) c----------------------- eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0 eold_l=eold_l+pold_l/(gm1*rold_l) eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0 eold_r=eold_r+pold_r/(gm1*rold_r) c------------------------------------------------------------------- c Computation of the speed of sound and its derivatives c--------------------------------------------------------------------- aleft=SQRT(ga*pold_l/rold_l) arigh=SQRT(ga*pold_r/rold_r) am=0.5d0*(aleft+arigh) c-------------------------------------------------------------------- damr_r=-arigh/(4.0d0*rold_r) damu_r=0.0d0 damv_r=0.0d0 damp_r=ga/(4.0d0*arigh*rold_r) c----------------------- damr_l=-aleft/(4.0d0*rold_l) damu_l=0.0d0 damv_l=0.0d0 damp_l=ga/(4.0d0*aleft*rold_l) c---------------------------------------------------------------------- c computing numerical Mach number and its derivatives c---------------------------------------------------------------------- un_l=uold_l*n_x+vold_l*n_y un_r=uold_r*n_x+vold_r*n_y c---------------------------------------- ml=un_l/am mr=un_r/am c------------------------------------- c Mplus and Mmin are calligraphic c lettes M+ and M- from the paper c------------------------------------- if(ABS(ml) .ge. 1.0d0) then Mplus=(ml+ABS(ml))/2.0d0 else Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0 Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0) endif c----------------- if(ABS(mr) .ge. 1.0d0) then Mmin=(mr-ABS(mr))/2.0d0 else Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0 Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0) endif c--------------------------------------- c Derivatives of ml and mr with respect c to both sets of primitive variables: c left and right c-------------------------------------- temp_l=-un_l/(am*am) temp_r=-un_r/(am*am) c-------- dmlr_l=temp_l*damr_l dmlr_r=temp_l*damr_r dmrr_l=temp_r*damr_l dmrr_r=temp_r*damr_r c-------- dmlu_l=n_x/am+temp_l*damu_l dmlu_r=temp_l*damu_r dmru_l=temp_r*damu_l dmru_r=n_x/am+temp_r*damu_r c-------- dmlv_l=n_y/am+temp_l*damv_l dmlv_r=temp_l*damv_r dmrv_l=temp_r*damv_l dmrv_r=n_y/am+temp_r*damv_r c-------- dmlp_l=temp_l*damp_l dmlp_r=temp_l*damp_r dmrp_l=temp_r*damp_l dmrp_r=temp_r*damp_r c--------------------------------------------------------------- c Computing the calligraphic P+ and P- with their derivatives c--------------------------------------------------------------- if(ml .ge. 1.0d0) then Pplus = 1.0d0 else if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0 Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0) else Pplus = 0.0d0 endif endif c--------------------------------------------------------------- if(mr .ge. 1.0d0) then Pmin = 0.0d0 else if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0 Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0) else Pmin = 1.0d0 endif endif c--------------------------------------------------------------- brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0 brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0) brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0) c-------------- brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0 brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0) brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0) c--------------------------------------------------------------- if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then dPpr_l=brac_l*dmlr_l dPpr_r=brac_l*dmlr_r c------------ dPpu_l=brac_l*dmlu_l dPpu_r=brac_l*dmlu_r c------------ dPpv_l=brac_l*dmlv_l dPpv_r=brac_l*dmlv_r c------------ dPpp_l=brac_l*dmlp_l dPpp_r=brac_l*dmlp_r c------------ else dPpr_l=0.0d0 dPpr_r=0.0d0 c----------- dPpu_l=0.0d0 dPpu_r=0.0d0 c----------- dPpv_l=0.0d0 dPpv_r=0.0d0 c----------- dPpp_l=0.0d0 dPpp_r=0.0d0 c----------- endif c--------------------------------------------------------------- if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then dPmr_l=brac_r*dmrr_l dPmr_r=brac_r*dmrr_r c------------ dPmu_l=brac_r*dmru_l dPmu_r=brac_r*dmru_r c------------ dPmv_l=brac_r*dmrv_l dPmv_r=brac_r*dmrv_r c------------ dPmp_l=brac_r*dmrp_l dPmp_r=brac_r*dmrp_r c------------ else dPmr_l=0.0d0 dPmr_r=0.0d0 c----------- dPmu_l=0.0d0 dPmu_r=0.0d0 c----------- dPmv_l=0.0d0 dPmv_r=0.0d0 c----------- dPmp_l=0.0d0 dPmp_r=0.0d0 c----------- endif c--------------------------------------------------------------------- c computing pmid -- p_{1/2} and its derivatives c--------------------------------------------------------------------- dpir_l=dPpr_l*pold_l+dPmr_l*pold_r dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r c---------------------------- dpir_r=dPpr_r*pold_l+dPmr_r*pold_r dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r c--------------------------------------------------------------------- c !!!!!!!!!!!!!!!!!! JACOBIAN !!!!!!!!!!!!!!!!!!!!!!!!!!!! c--------------------------------------------------------------------- jl(1,1)=0.0D0 jl(1,2)=0.0D0 jl(1,3)=0.0D0 jl(1,4)=0.0D0 c------------------------------------ jr(1,1)=0.0D0 jr(1,2)=0.0D0 jr(1,3)=0.0D0 jr(1,4)=0.0D0 c------------------------------------ ccccc f222222222222 ------------- c------------------------------------ c--------------------------------------------------------- jl(2,1)=n_x*dpir_l c------------------- jl(2,2)=n_x*dpiu_l c------------------- jl(2,3)=n_x*dpiv_l c------------------- jl(2,4)=n_x*dpip_l c------------------------------------------------------------- jr(2,1)=n_x*dpir_r c------------------- jr(2,2)=n_x*dpiu_r c------------------- jr(2,3)=n_x*dpiv_r c------------------- jr(2,4)=n_x*dpip_r c------------------------------------------------------------- c------------ f33333333333333333333 --------------------- c------------------------------------------------------------- jl(3,1)=n_y*dpir_l c------------------- jl(3,2)=n_y*dpiu_l c------------------- jl(3,3)=n_y*dpiv_l c------------------- jl(3,4)=n_y*dpip_l c------------------------------------------------------------- jr(3,1)=n_y*dpir_r c------------------- jr(3,2)=n_y*dpiu_r c------------------- jr(3,3)=n_y*dpiv_r c------------------- jr(3,4)=n_y*dpip_r c------------------------------------------------------------- c ------ f44444444444444444444444444444 --------- c------------------------------------------------------------- jl(4,1)=0.0D0 c--------------------- jl(4,2)=0.0D0 c--------------------- jl(4,3)=0.0D0 c--------------------- jl(4,4)=0.0D0 c---------------------------------------------------------- c---------------------------------------------------------- jr(4,1)=0.0D0 c--------------------- jr(4,2)=0.0D0 c--------------------- jr(4,3)=0.0D0 c--------------------- jr(4,4)=0.0D0 c---------------------------------------------------------- c----------------------------- tcoef=nvect(1)*tvect(2)+tvect(1)*nvect(2) bcoef=nvect(1)*tvect(2)-tvect(1)*nvect(2) c----------------------------- do 11 i=1,4 jl(i,1)=jl(i,1)+jr(i,1) jl(i,2)=jl(i,2)+jr(i,2)*(-tcoef/bcoef)+ & jr(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef jl(i,3)=jl(i,3)+jr(i,2)*(-2.0d0*nvect(2)*tvect(2)/bcoef)+ & jr(i,3)*tcoef/bcoef jl(i,4)=jl(i,4)+jr(i,4) 11 continue c---------------------------------------------------------------------- return end