1 : C CONJ3D SOURCE CHAT 05/01/12 22:17:23 5004
2 : c----------------------------------------------------------------------
3 : c GENERAL DESCRIPTION:
4 : c This subroutine provides the jacobians which are the derivatives
5 : c of the numerical flux function defined at the cell interface
6 : c with respect to the conservative variables of the left and right
7 : c cells (relative to the cell interface).
8 : c
9 : c EQUATIONS: 3D Euler equations of gas dynamics
10 : c
11 : c
12 : c REFERENCE: JCP, 129, 364-382 (1996)
13 : c " A Sequel to AUSM: AUSM+ ".
14 : c----------------------------------------------------------------------
15 : c INPUT:
16 : c
17 : c alpha -- parameter of the AUSM+ scheme in the Pressure function;
18 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
19 : c
20 : c beta -- parameter of the AUSM+ scheme in the Mach function;
21 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
22 : c
23 : c wvec_l -- vector of the primitive variables (rho,ux,uy,uz,p) at the
24 : c left cell;
25 : c
26 : c wvec_r -- vector of the primitive variables (rho,ux,uy,uz,p) at the
27 : c right cell;
28 : c
29 : c nvect -- normal vector to the interface (3 components in 3D);
30 : c
31 : c tvec1 -- tangential vector to the interface;
32 : c
33 : c tvec2 -- tangential vector to the interface;
34 : c
35 : c ga -- ratio of the specific heats (assumed constant)
36 : c----------------------------------------------------------------------
37 : c
38 : c OUTPUT:
39 : c
40 : c jtl -- jakobian matrix 5 by 5 - derivatives of the numerical
41 : c flux function with respect to the conservative variables
42 : c from the left cell;
43 : c
44 : c jtr -- jakobian matrix 5 by 5 - derivatives of the numerical
45 : c flux function with respect to the conservative variables
46 : c from the right cell.
47 : c----------------------------------------------------------------------
48 : SUBROUTINE CONJ3D(JTL,JTR,WVEC_L,WVEC_R,NVECT,TVEC1,TVEC2,
49 : & ga)
50 : IMPLICIT INTEGER(I-N)
51 : real*8 wvec_l(5),wvec_r(5)
52 : real*8 nvect(3),tvec1(3),tvec2(3)
53 : real*8 jl(5,5),jr(5,5),f(5)
54 : real*8 wl(5,5),wr(5,5)
55 : real*8 jtl(5,5),jtr(5,5)
56 : real*8 alpha,beta
57 : real*8 ga, gm1,temph
58 : real*8 n_x,n_y,n_z,t1_x,t1_y,t1_z,t2_x,t2_y,t2_z
59 : real*8 un_l,un_r,ut1_l,ut1_r,ut2_l,ut2_r
60 : real*8 ml,mr,Mplus,Mmin,mmid
61 : real*8 mpl_m, mmin_m,am
62 : real*8 rold_l,uold_l,vold_l,wold_l,pold_l,eold_l
63 : real*8 rold_r,uold_r,vold_r,wold_r,pold_r,eold_r
64 : real*8 Pplus,Pmin,pmid
65 : real*8 hr_l,hr_r,det,c11,c12,c13,c21,c22,c23,c31,c32,c33
66 : real*8 br1,temp_l,temp_r,brac_l,brac_r
67 : real*8 aleft, arigh
68 : real*8 damr_l,damr_r,damu_l,damu_r
69 : real*8 damv_l,damv_r,damp_l,damp_r
70 : real*8 damw_l,damw_r
71 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
72 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
73 : real*8 dmlw_l,dmlw_r
74 : real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
75 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
76 : real*8 dmrw_l,dmrw_r
77 : real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
78 : real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
79 : real*8 dMpw_l,dMpw_r
80 : real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
81 : real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
82 : real*8 dMmw_l,dMmw_r
83 : real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
84 : real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
85 : real*8 dmiw_l,dmiw_r
86 : real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
87 : real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
88 : real*8 d3mw_l,d3mw_r
89 : real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
90 : real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
91 : real*8 d2mw_l,d2mw_r
92 : real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
93 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
94 : real*8 dPpw_l,dPpw_r
95 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
96 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
97 : real*8 dPmw_l,dPmw_r
98 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
99 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
100 : real*8 dpiw_l,dpiw_r
101 : real*8 ff,fs,ft,dffr,dffu,dffv,dffw,dffp
102 : real*8 dfsr,dfsu,dfsv,dfsw,dfsp,dftr,dftu,dftv,dftw,dftp
103 : integer i,j,k
104 : parameter(alpha = 0.1875D0, beta = 0.125D0)
105 : c-----------------------
106 : gm1=ga-1.0d0
107 : c-----------------------
108 : n_x=nvect(1)
109 : n_y=nvect(2)
110 : n_z=nvect(3)
111 : c-------------------
112 : t1_x=tvec1(1)
113 : t1_y=tvec1(2)
114 : t1_z=tvec1(3)
115 : c-------------------
116 : t2_x=tvec2(1)
117 : t2_y=tvec2(2)
118 : t2_z=tvec2(3)
119 : c-------------------------------------------
120 : c These are coefficients for the expressions
121 : c
122 : c un = (c11*ux + c12*uy + c13*uz)/det
123 : c ut1 = (c21*ux + c22*uy + c23*uz)/det
124 : c ut2 = (c31*ux + c32*uy + c33*uz)/det
125 : c-------------------------------------------
126 : c11=t1_y*t2_z - t1_z*t2_y
127 : c12=t1_z*t2_x - t1_x*t2_z
128 : c13=t1_x*t2_y - t1_y*t2_x
129 : det=n_x*c11 + n_y*c12 + n_z*c13
130 : c--------------------------------------
131 : c21=n_z*t2_y - n_y*t2_z
132 : c22=n_x*t2_z - n_z*t2_x
133 : c23=n_y*t2_x - n_x*t2_y
134 : c--------------------------------------
135 : c31=n_y*t1_z - n_z*t1_y
136 : c32=n_z*t1_x - n_x*t1_z
137 : c33=n_x*t1_y - n_y*t1_x
138 : c--------------------------------------
139 : rold_l=wvec_l(1)
140 : uold_l=wvec_l(2)
141 : vold_l=wvec_l(3)
142 : wold_l=wvec_l(4)
143 : pold_l=wvec_l(5)
144 : c-----------------------
145 : rold_r=wvec_r(1)
146 : uold_r=wvec_r(2)
147 : vold_r=wvec_r(3)
148 : wold_r=wvec_r(4)
149 : pold_r=wvec_r(5)
150 : c------------------------------------------------------------------
151 : c Computation of the specific total energy on the left and right.
152 : c------------------------------------------------------------------
153 : eold_l=(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
154 : eold_l=eold_l+pold_l/(gm1*rold_l)
155 : eold_r=(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
156 : eold_r=eold_r+pold_r/(gm1*rold_r)
157 : c-------------------------------------------------------------------
158 : c Computation of the speed of sound and its derivatives;
159 : c numerical speed of sound at the interface is taken as an average
160 : c of the speeds of sounds of the neighbouring cells
161 : c---------------------------------------------------------------------
162 : aleft=dsqrt(ga*pold_l/rold_l)
163 : arigh=dsqrt(ga*pold_r/rold_r)
164 : am=0.5d0*(aleft+arigh)
165 : c--------------------------------------------------------------------
166 : damr_r=-arigh/(4.0d0*rold_r)
167 : damu_r=0.0d0
168 : damv_r=0.0d0
169 : damw_r=0.0d0
170 : damp_r=ga/(4.0d0*arigh*rold_r)
171 : c-----------------------
172 : damr_l=-aleft/(4.0d0*rold_l)
173 : damu_l=0.0d0
174 : damv_l=0.0d0
175 : damw_l=0.0d0
176 : damp_l=ga/(4.0d0*aleft*rold_l)
177 : c-------------------------------------------------------------------
178 : c Computing numerical Mach number and its derivatives,
179 : c see p.370, under (A1).
180 : c-------------------------------------------------------------------
181 : un_l=uold_l*n_x+vold_l*n_y+wold_l*n_z
182 : un_r=uold_r*n_x+vold_r*n_y+wold_r*n_z
183 : c----------
184 : ut1_l=uold_l*t1_x+vold_l*t1_y+wold_l*t1_z
185 : ut1_r=uold_r*t1_x+vold_r*t1_y+wold_r*t1_z
186 : c----------
187 : ut2_l=uold_l*t2_x+vold_l*t2_y+wold_l*t2_z
188 : ut2_r=uold_r*t2_x+vold_r*t2_y+wold_r*t2_z
189 : c----------------------------------------
190 : ml=un_l/am
191 : mr=un_r/am
192 : c-------------------------------------------------------------------
192 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
193 : c see (19a) and (19b), p.367.
194 : c-------------------------------------------------------------------
195 : if(dabs(ml) .ge. 1.0d0) then
196 : Mplus=(ml+dabs(ml))/2.0d0
197 : else
198 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
199 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
200 : endif
201 : c-----------------
202 : if(dabs(mr) .ge. 1.0d0) then
203 : Mmin=(mr-dabs(mr))/2.0d0
204 : else
205 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
206 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
207 : endif
208 : c-------------------------------------------------------------------
209 : c Derivatives of ml and mr with respect to both sets of primitive
210 : c variables: left and right.
211 : c-------------------------------------------------------------------
212 : temp_l=-un_l/(am*am)
213 : temp_r=-un_r/(am*am)
214 : c--------
215 : dmlr_l=temp_l*damr_l
216 : dmlr_r=temp_l*damr_r
217 : dmrr_l=temp_r*damr_l
218 : dmrr_r=temp_r*damr_r
219 : c--------
220 : dmlu_l=n_x/am+temp_l*damu_l
221 : dmlu_r=temp_l*damu_r
222 : dmru_l=temp_r*damu_l
223 : dmru_r=n_x/am+temp_r*damu_r
224 : c--------
225 : dmlv_l=n_y/am+temp_l*damv_l
226 : dmlv_r=temp_l*damv_r
227 : dmrv_l=temp_r*damv_l
228 : dmrv_r=n_y/am+temp_r*damv_r
229 : c--------
230 : dmlw_l=n_z/am+temp_l*damw_l
231 : dmlw_r=temp_l*damw_r
232 : dmrw_l=temp_r*damw_l
233 : dmrw_r=n_z/am+temp_r*damw_r
234 : c--------
235 : dmlp_l=temp_l*damp_l
236 : dmlp_r=temp_l*damp_r
237 : dmrp_l=temp_r*damp_l
238 : dmrp_r=temp_r*damp_r
239 : c-----------------------------------------------------------
240 : c mmid is m_{1/2} (notation as in the paper)
241 : c----------------------------------------------------------
242 : mmid=Mplus+Mmin
243 : c-----------------------------------------------------------
244 : c Computing the derivatives of M+ and M-
245 : c-----------------------------------------------------------
246 : if(ml .ge. 1.0d0) then
247 : dMpr_l=dmlr_l
248 : dMpu_l=dmlu_l
249 : dMpv_l=dmlv_l
250 : dMpw_l=dmlw_l
251 : dMpp_l=dmlp_l
252 : c--------------------
253 : dMpr_r=dmlr_r
254 : dMpu_r=dmlu_r
255 : dMpv_r=dmlv_r
256 : dMpw_r=dmlw_r
257 : dMpp_r=dmlp_r
258 : else
259 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
260 : temph=(ml+1.0d0)/2.0d0
261 : dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
262 : dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
263 : dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
264 : dMpw_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_l
265 : dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
266 : c--------------------
267 : dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
268 : dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
269 : dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
270 : dMpw_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_r
271 : dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
272 : else
273 : dMpr_l=0.0d0
274 : dMpu_l=0.0d0
275 : dMpv_l=0.0d0
276 : dMpw_l=0.0d0
277 : dMpp_l=0.0d0
278 : c---------------------
279 : dMpr_r=0.0d0
280 : dMpu_r=0.0d0
281 : dMpv_r=0.0d0
282 : dMpw_r=0.0d0
283 : dMpp_r=0.0d0
284 : endif
285 : endif
286 : c-----------------------------------------------------------
287 : if(mr .ge. 1.0d0) then
288 : dMmr_l=0.0d0
289 : dMmu_l=0.0d0
290 : dMmv_l=0.0d0
291 : dMmw_l=0.0d0
292 : dMmp_l=0.0d0
293 : c---------------------
294 : dMmr_r=0.0d0
295 : dMmu_r=0.0d0
296 : dMmv_r=0.0d0
297 : dMmw_r=0.0d0
298 : dMmp_r=0.0d0
299 : else
300 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
301 : temph=(-mr+1.0d0)/2.0d0
302 : dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
303 : dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
304 : dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
305 : dMmw_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_l
306 : dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
307 : c--------------------
308 : dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
309 : dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
310 : dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
311 : dMmw_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_r
312 : dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
313 : else
314 : dMmr_l=dmrr_l
315 : dMmu_l=dmru_l
316 : dMmv_l=dmrv_l
317 : dMmw_l=dmrw_l
318 : dMmp_l=dmrp_l
319 : c--------------------
320 : dMmr_r=dmrr_r
321 : dMmu_r=dmru_r
322 : dMmv_r=dmrv_r
323 : dMmw_r=dmrw_r
324 : dMmp_r=dmrp_r
325 : endif
326 : endif
327 : c-----------------------------------------------------------------
328 : c computing the derivatives of m_{1/2} (notation as in the paper)
329 : c-----------------------------------------------------------------
330 : dmir_l=dMpr_l+dMmr_l
331 : dmir_r=dMpr_r+dMmr_r
332 : c-------------
333 : dmiu_l=dMpu_l+dMmu_l
334 : dmiu_r=dMpu_r+dMmu_r
335 : c-------------
336 : dmiv_l=dMpv_l+dMmv_l
337 : dmiv_r=dMpv_r+dMmv_r
338 : c-------------
339 : dmiw_l=dMpw_l+dMmw_l
340 : dmiw_r=dMpw_r+dMmw_r
341 : c-------------
342 : dmip_l=dMpp_l+dMmp_l
343 : dmip_r=dMpp_r+dMmp_r
344 : c-------------------------------------------------------------
345 : c computing the main convective variables and their derivatives
346 : c mpl_m is m^{+}_{1/2} (paper's notation) and
347 : c mmin_m is m^{-}_{1/2} (paper's notation)
348 : c-------------------------------------------------------------
349 : if(mmid .gt. 0.0d0) then
350 : mpl_m = mmid
351 : d2mr_l=dmir_l
352 : d2mu_l=dmiu_l
353 : d2mv_l=dmiv_l
354 : d2mw_l=dmiw_l
355 : d2mp_l=dmip_l
356 : c------------
357 : d2mr_r=dmir_r
358 : d2mu_r=dmiu_r
359 : d2mv_r=dmiv_r
360 : d2mw_r=dmiw_r
361 : d2mp_r=dmip_r
362 : c------------
363 : else
364 : mpl_m = 0.0d0
365 : d2mr_l=0.0d0
366 : d2mu_l=0.0d0
367 : d2mv_l=0.0d0
368 : d2mw_l=0.0d0
369 : d2mp_l=0.0d0
370 : c------------
371 : d2mr_r=0.0d0
372 : d2mu_r=0.0d0
373 : d2mv_r=0.0d0
374 : d2mw_r=0.0d0
375 : d2mp_r=0.0d0
376 : endif
377 : c---------------------------------------------
378 : if(mmid .lt. 0.0d0) then
379 : mmin_m = mmid
380 : d3mr_l=dmir_l
381 : d3mu_l=dmiu_l
382 : d3mv_l=dmiv_l
383 : d3mw_l=dmiw_l
384 : d3mp_l=dmip_l
385 : c------------
386 : d3mr_r=dmir_r
387 : d3mu_r=dmiu_r
388 : d3mv_r=dmiv_r
389 : d3mw_r=dmiw_r
390 : d3mp_r=dmip_r
391 : c------------
392 : else
393 : mmin_m = 0.0d0
394 : d3mr_l=0.0d0
395 : d3mu_l=0.0d0
396 : d3mv_l=0.0d0
397 : d3mw_l=0.0d0
398 : d3mp_l=0.0d0
399 : c------------
400 : d3mr_r=0.0d0
401 : d3mu_r=0.0d0
402 : d3mv_r=0.0d0
403 : d3mw_r=0.0d0
404 : d3mp_r=0.0d0
405 : endif
406 : c---------------------------------------------------------------
406 : c Computing the calligraphic P+ and P- with their derivatives
407 : c see (21a) & (21b) on p.368.
408 : c---------------------------------------------------------------
409 : if(ml .ge. 1.0d0) then
410 : Pplus = 1.0d0
411 : else
412 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
413 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
414 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
415 : else
416 : Pplus = 0.0d0
417 : endif
418 : endif
419 : c---------------------------------------------------------------
420 : if(mr .ge. 1.0d0) then
421 : Pmin = 0.0d0
422 : else
423 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
424 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
425 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
426 : else
427 : Pmin = 1.0d0
428 : endif
429 : endif
430 : c---------------------------------------------------------------
431 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
432 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
433 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
434 : c--------------
435 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
436 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
437 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
438 : c---------------------------------------------------------------
439 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
440 : dPpr_l=brac_l*dmlr_l
441 : dPpr_r=brac_l*dmlr_r
442 : c------------
443 : dPpu_l=brac_l*dmlu_l
444 : dPpu_r=brac_l*dmlu_r
445 : c------------
446 : dPpv_l=brac_l*dmlv_l
447 : dPpv_r=brac_l*dmlv_r
448 : c------------
449 : dPpw_l=brac_l*dmlw_l
450 : dPpw_r=brac_l*dmlw_r
451 : c------------
452 : dPpp_l=brac_l*dmlp_l
453 : dPpp_r=brac_l*dmlp_r
454 : c------------
455 : else
456 : dPpr_l=0.0d0
457 : dPpr_r=0.0d0
458 : c-----------
459 : dPpu_l=0.0d0
460 : dPpu_r=0.0d0
461 : c-----------
462 : dPpv_l=0.0d0
463 : dPpv_r=0.0d0
464 : c-----------
465 : dPpw_l=0.0d0
466 : dPpw_r=0.0d0
467 : c-----------
468 : dPpp_l=0.0d0
469 : dPpp_r=0.0d0
470 : c-----------
471 : endif
472 : c---------------------------------------------------------------
473 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
474 : dPmr_l=brac_r*dmrr_l
475 : dPmr_r=brac_r*dmrr_r
476 : c------------
477 : dPmu_l=brac_r*dmru_l
478 : dPmu_r=brac_r*dmru_r
479 : c------------
480 : dPmv_l=brac_r*dmrv_l
481 : dPmv_r=brac_r*dmrv_r
482 : c------------
483 : dPmw_l=brac_r*dmrw_l
484 : dPmw_r=brac_r*dmrw_r
485 : c------------
486 : dPmp_l=brac_r*dmrp_l
487 : dPmp_r=brac_r*dmrp_r
488 : c------------
489 : else
490 : dPmr_l=0.0d0
491 : dPmr_r=0.0d0
492 : c-----------
493 : dPmu_l=0.0d0
494 : dPmu_r=0.0d0
495 : c-----------
496 : dPmv_l=0.0d0
497 : dPmv_r=0.0d0
498 : c-----------
499 : dPmw_l=0.0d0
500 : dPmw_r=0.0d0
501 : c-----------
502 : dPmp_l=0.0d0
503 : dPmp_r=0.0d0
504 : c-----------
505 : endif
506 : c---------------------------------------------------------------------
507 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
508 : c---------------------------------------------------------------------
509 : pmid=Pplus*pold_l + Pmin*pold_r
510 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
511 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
512 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
513 : dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
514 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
515 : c----------------------------
516 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
517 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
518 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
519 : dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
520 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
521 : c---------------------------------------------------------------------
522 : c computing JACOBIAN as a derivative of the numerical flux function
523 : c with respect to the primitive variables.
524 : c Notation: jl(i,j) --- is the derivative of the i-component of the
525 : c flux function with respect to the j-component of the
526 : c vector of primitive variables of the left state.
527 : c jr(i,j) --- is the derivative of the i-component of the
528 : c flux function with respect to the j-component of the
529 : c vector of primitive variables of the right state.
530 : c---------------------------------------------------------------------
531 : f(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
532 : c---------------------------------------------------------------------
533 : jl(1,1)=damr_l*f(1)/am+am*(d2mr_l*rold_l+mpl_m)
534 : jl(1,1)=jl(1,1)+am*d3mr_l*rold_r
535 : jl(1,2)=damu_l*f(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
536 : jl(1,3)=damv_l*f(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
537 : jl(1,4)=damw_l*f(1)/am+am*(d2mw_l*rold_l+d3mw_l*rold_r)
538 : jl(1,5)=damp_l*f(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
539 : c------------------------------------
540 : jr(1,1)=damr_r*f(1)/am+am*(d2mr_r*rold_l+mmin_m)
541 : jr(1,1)=jr(1,1)+am*d3mr_r*rold_r
542 : jr(1,2)=damu_r*f(1)/am+am*(d2mu_r*rold_l+d3mu_r*rold_r)
543 : jr(1,3)=damv_r*f(1)/am+am*(d2mv_r*rold_l+d3mv_r*rold_r)
544 : jr(1,4)=damw_r*f(1)/am+am*(d2mw_r*rold_l+d3mw_r*rold_r)
545 : jr(1,5)=damp_r*f(1)/am+am*(d2mp_r*rold_l+d3mp_r*rold_r)
546 : c------------------------------------
547 : c------------------------------------
548 : ff=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
549 : dffr=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
550 : dffu=d2mu_l*rold_l*un_l+mpl_m*rold_l*c11/det
551 : dffu=dffu+d3mu_l*rold_r*un_r
552 : dffv=d2mv_l*rold_l*un_l+mpl_m*rold_l*c12/det
553 : dffv=dffv+d3mv_l*rold_r*un_r
554 : dffw=d2mw_l*rold_l*un_l+mpl_m*rold_l*c13/det
555 : dffw=dffw+d3mw_l*rold_r*un_r
556 : dffp=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
557 : c------------------------------------------------
558 : fs=mpl_m*rold_l*ut1_l+mmin_m*rold_r*ut1_r
559 : dfsr=d2mr_l*rold_l*ut1_l+mpl_m*ut1_l
560 : dfsr=dfsr+d3mr_l*rold_r*ut1_r
561 : dfsu=d2mu_l*rold_l*ut1_l+mpl_m*rold_l*c21/det
562 : dfsu=dfsu+d3mu_l*rold_r*ut1_r
563 : dfsv=d2mv_l*rold_l*ut1_l+mpl_m*rold_l*c22/det
564 : dfsv=dfsv+d3mv_l*rold_r*ut1_r
565 : dfsw=d2mw_l*rold_l*ut1_l+mpl_m*rold_l*c23/det
566 : dfsw=dfsw+d3mw_l*rold_r*ut1_r
567 : dfsp=d2mp_l*rold_l*ut1_l+d3mp_l*rold_r*ut1_r
568 : c-------------------------------------------------
569 : ft=mpl_m*rold_l*ut2_l+mmin_m*rold_r*ut2_r
570 : dftr=d2mr_l*rold_l*ut2_l+mpl_m*ut2_l
571 : dftr=dftr+d3mr_l*rold_r*ut2_r
572 : dftu=d2mu_l*rold_l*ut2_l+mpl_m*rold_l*c31/det
573 : dftu=dftu+d3mu_l*rold_r*ut2_r
574 : dftv=d2mv_l*rold_l*ut2_l+mpl_m*rold_l*c32/det
575 : dftv=dftv+d3mv_l*rold_r*ut2_r
576 : dftw=d2mw_l*rold_l*ut2_l+mpl_m*rold_l*c33/det
577 : dftw=dftw+d3mw_l*rold_r*ut2_r
578 : dftp=d2mp_l*rold_l*ut2_l+d3mp_l*rold_r*ut2_r
579 : c-------------------------------------------------
580 : f(2)=am*(ff*n_x+fs*t1_x+ft*t2_x)+pmid*n_x
581 : c-------------------------------------------------
582 : jl(2,1)=damr_l*(f(2)-pmid*n_x)/am+dpir_l*n_x
583 : jl(2,1)=jl(2,1)+am*(dffr*n_x+dfsr*t1_x+dftr*t2_x)
584 : jl(2,2)=damu_l*(f(2)-pmid*n_x)/am+dpiu_l*n_x
585 : jl(2,2)=jl(2,2)+am*(dffu*n_x+dfsu*t1_x+dftu*t2_x)
586 : jl(2,3)=damv_l*(f(2)-pmid*n_x)/am+dpiv_l*n_x
587 : jl(2,3)=jl(2,3)+am*(dffv*n_x+dfsv*t1_x+dftv*t2_x)
588 : jl(2,4)=damw_l*(f(2)-pmid*n_x)/am+dpiw_l*n_x
589 : jl(2,4)=jl(2,4)+am*(dffw*n_x+dfsw*t1_x+dftw*t2_x)
590 : jl(2,5)=damp_l*(f(2)-pmid*n_x)/am+dpip_l*n_x
591 : jl(2,5)=jl(2,5)+am*(dffp*n_x+dfsp*t1_x+dftp*t2_x)
592 : c-------------------------------------------------
593 : f(3)=am*(ff*n_y+fs*t1_y+ft*t2_y)+pmid*n_y
594 : c-------------------------------------------------
595 : jl(3,1)=damr_l*(f(3)-pmid*n_y)/am+dpir_l*n_y
596 : jl(3,1)=jl(3,1)+am*(dffr*n_y+dfsr*t1_y+dftr*t2_y)
597 : jl(3,2)=damu_l*(f(3)-pmid*n_y)/am+dpiu_l*n_y
598 : jl(3,2)=jl(3,2)+am*(dffu*n_y+dfsu*t1_y+dftu*t2_y)
599 : jl(3,3)=damv_l*(f(3)-pmid*n_y)/am+dpiv_l*n_y
600 : jl(3,3)=jl(3,3)+am*(dffv*n_y+dfsv*t1_y+dftv*t2_y)
601 : jl(3,4)=damw_l*(f(3)-pmid*n_y)/am+dpiw_l*n_y
602 : jl(3,4)=jl(3,4)+am*(dffw*n_y+dfsw*t1_y+dftw*t2_y)
603 : jl(3,5)=damp_l*(f(3)-pmid*n_y)/am+dpip_l*n_y
604 : jl(3,5)=jl(3,5)+am*(dffp*n_y+dfsp*t1_y+dftp*t2_y)
605 : c-------------------------------------------------
606 : f(4)=am*(ff*n_z+fs*t1_z+ft*t2_z)+pmid*n_z
607 : c-------------------------------------------------
608 : jl(4,1)=damr_l*(f(4)-pmid*n_z)/am+dpir_l*n_z
609 : jl(4,1)=jl(4,1)+am*(dffr*n_z+dfsr*t1_z+dftr*t2_z)
610 : jl(4,2)=damu_l*(f(4)-pmid*n_z)/am+dpiu_l*n_z
611 : jl(4,2)=jl(4,2)+am*(dffu*n_z+dfsu*t1_z+dftu*t2_z)
612 : jl(4,3)=damv_l*(f(4)-pmid*n_z)/am+dpiv_l*n_z
613 : jl(4,3)=jl(4,3)+am*(dffv*n_z+dfsv*t1_z+dftv*t2_z)
614 : jl(4,4)=damw_l*(f(4)-pmid*n_z)/am+dpiw_l*n_z
615 : jl(4,4)=jl(4,4)+am*(dffw*n_z+dfsw*t1_z+dftw*t2_z)
616 : jl(4,5)=damp_l*(f(4)-pmid*n_z)/am+dpip_l*n_z
617 : jl(4,5)=jl(4,5)+am*(dffp*n_z+dfsp*t1_z+dftp*t2_z)
618 : c-------------------------------------------------------
619 : c derivatives with respect to the right
620 : c set of primitive variables
621 : c-------------------------------------------------------
622 : dffr=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
623 : dffu=d2mu_r*rold_l*un_l+mmin_m*rold_r*c11/det
624 : dffu=dffu+d3mu_r*rold_r*un_r
625 : dffv=d2mv_r*rold_l*un_l+mmin_m*rold_r*c12/det
626 : dffv=dffv+d3mv_r*rold_r*un_r
627 : dffw=d2mw_r*rold_l*un_l+mmin_m*rold_r*c13/det
628 : dffw=dffw+d3mw_r*rold_r*un_r
629 : dffp=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
630 : c------------------------------------------------------
631 : dfsr=d2mr_r*rold_l*ut1_l+mmin_m*ut1_r
632 : dfsr=dfsr+d3mr_r*rold_r*ut1_r
633 : dfsu=d2mu_r*rold_l*ut1_l+mmin_m*rold_r*c21/det
634 : dfsu=dfsu+d3mu_r*rold_r*ut1_r
635 : dfsv=d2mv_r*rold_l*ut1_l+mmin_m*rold_r*c22/det
636 : dfsv=dfsv+d3mv_r*rold_r*ut1_r
637 : dfsw=d2mw_r*rold_l*ut1_l+mmin_m*rold_r*c23/det
638 : dfsw=dfsw+d3mw_r*rold_r*ut1_r
639 : dfsp=d2mp_r*rold_l*ut1_l+d3mp_r*rold_r*ut1_r
640 : c------------------------------------------------------
641 : dftr=d2mr_r*rold_l*ut2_l+mmin_m*ut2_r
642 : dftr=dftr+d3mr_r*rold_r*ut2_r
643 : dftu=d2mu_r*rold_l*ut2_l+mmin_m*rold_r*c31/det
644 : dftu=dftu+d3mu_r*rold_r*ut2_r
645 : dftv=d2mv_r*rold_l*ut2_l+mmin_m*rold_r*c32/det
646 : dftv=dftv+d3mv_r*rold_r*ut2_r
647 : dftw=d2mw_r*rold_l*ut2_l+mmin_m*rold_r*c33/det
648 : dftw=dftw+d3mw_r*rold_r*ut2_r
649 : dftp=d2mp_r*rold_l*ut2_l+d3mp_r*rold_r*ut2_r
650 : c-------------------------------------------------------
651 : jr(2,1)=damr_r*(f(2)-pmid*n_x)/am+dpir_r*n_x
652 : jr(2,1)=jr(2,1)+am*(dffr*n_x+dfsr*t1_x+dftr*t2_x)
653 : jr(2,2)=damu_r*(f(2)-pmid*n_x)/am+dpiu_r*n_x
654 : jr(2,2)=jr(2,2)+am*(dffu*n_x+dfsu*t1_x+dftu*t2_x)
655 : jr(2,3)=damv_r*(f(2)-pmid*n_x)/am+dpiv_r*n_x
656 : jr(2,3)=jr(2,3)+am*(dffv*n_x+dfsv*t1_x+dftv*t2_x)
657 : jr(2,4)=damw_r*(f(2)-pmid*n_x)/am+dpiw_r*n_x
658 : jr(2,4)=jr(2,4)+am*(dffw*n_x+dfsw*t1_x+dftw*t2_x)
659 : jr(2,5)=damp_r*(f(2)-pmid*n_x)/am+dpip_r*n_x
660 : jr(2,5)=jr(2,5)+am*(dffp*n_x+dfsp*t1_x+dftp*t2_x)
661 : c-------------------------------------------------------
662 : jr(3,1)=damr_r*(f(3)-pmid*n_y)/am+dpir_r*n_y
663 : jr(3,1)=jr(3,1)+am*(dffr*n_y+dfsr*t1_y+dftr*t2_y)
664 : jr(3,2)=damu_r*(f(3)-pmid*n_y)/am+dpiu_r*n_y
665 : jr(3,2)=jr(3,2)+am*(dffu*n_y+dfsu*t1_y+dftu*t2_y)
666 : jr(3,3)=damv_r*(f(3)-pmid*n_y)/am+dpiv_r*n_y
667 : jr(3,3)=jr(3,3)+am*(dffv*n_y+dfsv*t1_y+dftv*t2_y)
668 : jr(3,4)=damw_r*(f(3)-pmid*n_y)/am+dpiw_r*n_y
669 : jr(3,4)=jr(3,4)+am*(dffw*n_y+dfsw*t1_y+dftw*t2_y)
670 : jr(3,5)=damp_r*(f(3)-pmid*n_y)/am+dpip_r*n_y
671 : jr(3,5)=jr(3,5)+am*(dffp*n_y+dfsp*t1_y+dftp*t2_y)
672 : c--------------------------------------------------------
673 : jr(4,1)=damr_r*(f(4)-pmid*n_z)/am+dpir_r*n_z
674 : jr(4,1)=jr(4,1)+am*(dffr*n_z+dfsr*t1_z+dftr*t2_z)
675 : jr(4,2)=damu_r*(f(4)-pmid*n_z)/am+dpiu_r*n_z
676 : jr(4,2)=jr(4,2)+am*(dffu*n_z+dfsu*t1_z+dftu*t2_z)
677 : jr(4,3)=damv_r*(f(4)-pmid*n_z)/am+dpiv_r*n_z
678 : jr(4,3)=jr(4,3)+am*(dffv*n_z+dfsv*t1_z+dftv*t2_z)
679 : jr(4,4)=damw_r*(f(4)-pmid*n_z)/am+dpiw_r*n_z
680 : jr(4,4)=jr(4,4)+am*(dffw*n_z+dfsw*t1_z+dftw*t2_z)
681 : jr(4,5)=damp_r*(f(4)-pmid*n_z)/am+dpip_r*n_z
682 : jr(4,5)=jr(4,5)+am*(dffp*n_z+dfsp*t1_z+dftp*t2_z)
683 : c-------------------------------------------------------
684 : c-------------------------------------------------------
685 : c-------------------------------------------------------
686 : hr_l=(eold_l+pold_l/rold_l)*rold_l
687 : hr_r=(eold_r+pold_r/rold_r)*rold_r
688 : f(5)=am*(mpl_m*hr_l+mmin_m*hr_r)
689 : c-------------------------------------------------
690 : br1=d2mr_l*hr_l+d3mr_l*hr_r
691 : br1=br1+mpl_m*(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
692 : jl(5,1)=damr_l*f(5)/am+am*br1
693 : c-------------------------------------------------
694 : br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
695 : br1=br1+d3mu_l*hr_r
696 : jl(5,2)=damu_l*f(5)/am+am*br1
697 : c-------------------------------------------------
698 : br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
699 : br1=br1+d3mv_l*hr_r
700 : jl(5,3)=damv_l*f(5)/am+am*br1
701 : c-------------------------------------------------
702 : br1=d2mw_l*hr_l+mpl_m*wold_l*rold_l
703 : br1=br1+d3mw_l*hr_r
704 : jl(5,4)=damw_l*f(5)/am+am*br1
705 : c-------------------------------------------------
706 : br1=d2mp_l*hr_l+mpl_m*ga/gm1
707 : br1=br1+d3mp_l*hr_r
708 : jl(5,5)=damp_l*f(5)/am+am*br1
709 : c-------------------------------------------------
710 : c-------------------------------------------------
711 : br1=d2mr_r*hr_l+d3mr_r*hr_r
712 : br1=br1+mmin_m*(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
713 : jr(5,1)=damr_r*f(5)/am+am*br1
714 : c---------------------
715 : br1=d2mu_r*hr_l+mmin_m*uold_r*rold_r
716 : br1=br1+d3mu_r*hr_r
717 : jr(5,2)=damu_r*f(5)/am+am*br1
718 : c---------------------
719 : br1=d2mv_r*hr_l+mmin_m*vold_r*rold_r
720 : br1=br1+d3mv_r*hr_r
721 : jr(5,3)=damv_r*f(5)/am+am*br1
722 : c---------------------
723 : br1=d2mw_r*hr_l+mmin_m*wold_r*rold_r
724 : br1=br1+d3mw_r*hr_r
725 : jr(5,4)=damw_r*f(5)/am+am*br1
726 : c---------------------
727 : br1=d2mp_r*hr_l+mmin_m*ga/gm1
728 : br1=br1+d3mp_r*hr_r
729 : jr(5,5)=damp_r*f(5)/am+am*br1
730 : c-------------------------------------------------------------
731 : c matrix wl(i,j) represents the derivative of the i-component
732 : c of the vector of primitive variables of the left state with
733 : c respect to the j-component of the vector of the conservative
734 : c variables of the left state.
735 : c
736 : c Here: (rho, u, v, w, p) - vector of primitive variables;
737 : c (rho, rho u, rh o v, rho w, rho e) - conservative variables.
738 : c-------------------------------------------------------------
739 : wl(1,1)=1.0d0
740 : wl(1,2)=0.0d0
741 : wl(1,3)=0.0d0
742 : wl(1,4)=0.0d0
743 : wl(1,5)=0.0d0
744 : c------------------------------
745 : wl(2,1)=-uold_l/rold_l
746 : wl(2,2)=1.0d0/rold_l
747 : wl(2,3)=0.0d0
748 : wl(2,4)=0.0d0
749 : wl(2,5)=0.0d0
750 : c------------------------------
751 : wl(3,1)=-vold_l/rold_l
752 : wl(3,2)=0.0d0
753 : wl(3,3)=1.0d0/rold_l
754 : wl(3,4)=0.0d0
755 : wl(3,5)=0.0d0
756 : c------------------------------
757 : wl(4,1)=-wold_l/rold_l
758 : wl(4,2)=0.0d0
759 : wl(4,3)=0.0d0
760 : wl(4,4)=1.0d0/rold_l
761 : wl(4,5)=0.0d0
762 : c------------------------------
763 : wl(5,1)=gm1*(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
764 : wl(5,2)=-uold_l*gm1
765 : wl(5,3)=-vold_l*gm1
766 : wl(5,4)=-wold_l*gm1
767 : wl(5,5)=gm1
768 : c------------------------------
769 : c------------------------------
770 : wr(1,1)=1.0d0
771 : wr(1,2)=0.0d0
772 : wr(1,3)=0.0d0
773 : wr(1,4)=0.0d0
774 : wr(1,5)=0.0d0
775 : c------------------------------
776 : wr(2,1)=-uold_r/rold_r
777 : wr(2,2)=1.0d0/rold_r
778 : wr(2,3)=0.0d0
779 : wr(2,4)=0.0d0
780 : wr(2,5)=0.0d0
781 : c------------------------------
782 : wr(3,1)=-vold_r/rold_r
783 : wr(3,2)=0.0d0
784 : wr(3,3)=1.0d0/rold_r
785 : wr(3,4)=0.0d0
786 : wr(3,5)=0.0d0
787 : c------------------------------
788 : wr(4,1)=-wold_r/rold_r
789 : wr(4,2)=0.0d0
790 : wr(4,3)=0.0d0
791 : wr(4,4)=1.0d0/rold_r
792 : wr(4,5)=0.0d0
793 : c------------------------------
794 : wr(5,1)=gm1*(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
795 : wr(5,2)=-uold_r*gm1
796 : wr(5,3)=-vold_r*gm1
797 : wr(5,4)=-wold_r*gm1
798 : wr(5,5)=gm1
799 : c----------------------------------------------
800 : c----------------------------------------------
801 : do 1 i=1,5
802 : do 2 j=1,5
803 : jtl(i,j)=0.0d0
804 : jtr(i,j)=0.0d0
805 : do 3 k=1,5
806 : jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j)
807 : jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j)
808 : 3 continue
809 : 2 continue
810 : 1 continue
811 : c----------------------------------------------------------------------
812 : return
813 : end
814 :
815 :
816 :
817 :
818 :
819 :
820 :
821 :
822 :
823 :
824 :
825 :
826 :
827 :
828 :
829 :
830 :
831 :
832 :
833 :
© Cast3M 2003 - All rights reserved.
Disclaimer