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