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