1 : C COPMSP SOURCE CHAT 05/01/12 22:20:44 5004
2 : SUBROUTINE COPMSP(NSP,JPL,JLL,WVEC_L,WVEC_R,NVECT,TVECT,
3 : & mpyn,lrecp,lrecv,nlcg,nlcd)
4 : C************************************************************************
5 : C
6 : C PROJET : CASTEM 2000
7 : C
8 : C NOM : COPMSP ('convection for multispecies')
9 : C
10 : C DESCRIPTION : Voir CLI252 (appele par CLI252)
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 - MULTISPECIES GAS
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 nsp -- number of species (total);
40 : c
41 : c wvec_l -- vector of the primitive variables
42 : c (rho,u_x,u_y,p) at the left cell;
43 : c
44 : c wvec_r -- vector of the primitive variables
45 : c (rho,u_x,u_y,p) at the right cell;
46 : c
47 : c nvect -- normal vector to the interface (2 components in 2D);
48 : c
49 : c tvect -- tangential vector to the interface;
50 : c
51 : c mpyn -- pointer to the vectors of the primitive variables
52 : c (Y_1,Y_2,...Y_(nsp-1)) at the left and the right cells;
53 : c
54 : c lrecp -- pointer to the vector of specific heats at constant pressure
55 : c (size of the vector is equal to number of species (nsp));
56 : c
57 : c lrecv -- pointer to the vector of specific heats at constant volume
58 : c (size of the vector is equal to number of species (nsp));
59 : c
60 : c nlcg -- "local" number corresponding to the left cell;
61 : c
62 : c nlcd -- "local" number corresponding to the right cell;
63 : c----------------------------------------------------------------------
64 : c
65 : c OUTPUT:
66 : c
67 : c jll -- jakobian matrix (3+nsp) by (3+nsp) -
68 : c derivatives of the numerical
69 : c flux function with respect to the conservative variables
70 : c from the left cell;
71 : c
72 : c jpl -- jakobian matrix (3+nsp) by (3+nsp) -
73 : c derivatives of the numerical
74 : c flux function with respect to the primitive variables
75 : c from the left cell.
76 : c----------------------------------------------------------------------
77 : IMPLICIT INTEGER(I-N)
78 : integer nsp,jll,jpl,lrecp,lrecv,nlcg,nlcd
79 : real*8 wvec_l(4),wvec_r(4)
80 : real*8 nvect(2),tvect(2)
81 : real*8 alpha,beta
82 : real*8 gal,gar,gm1l,gm1r,temph
83 : real*8 n_x,n_y,t_x,t_y
84 : real*8 un_l, un_r, ut_l, ut_r
85 : real*8 ml,mr,Mplus,Mmin,mmid
86 : real*8 mpl_m, mmin_m,am
87 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l
88 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r
89 : real*8 Pplus,Pmin,pmid
90 : real*8 hr_l,hr_r,det,top,bot
91 : real*8 br1,br2,br3,br4,temp_l,temp_r,brac_l,brac_r
92 : real*8 aleft, arigh
93 : real*8 damr_l,damr_r,damu_l,damu_r
94 : real*8 damv_l,damv_r,damp_l,damp_r
95 : real*8 damg_l,damg_r
96 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
97 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
98 : real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
99 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
100 : real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
101 : real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
102 : real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
103 : real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
104 : real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
105 : real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
106 : real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
107 : real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
108 : real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
109 : real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
110 : real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
111 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
112 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
113 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
114 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
115 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
116 : integer i,j,k
117 : parameter(alpha = 0.1875D0, beta = 0.125D0)
118 : C------------------------------------------------------------
119 : -INC SMCHPOI
120 : POINTEUR MPYN.MPOVAL
121 : C-------------------------------------------------------------
122 : -INC SMLREEL
123 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
124 : C-------------------------------------------------------------
125 : C******* Les fractionines massiques **************************
126 : C-------------------------------------------------------------
127 : SEGMENT FRAMAS
128 : REAL*8 YET(NSP)
129 : ENDSEGMENT
130 : POINTEUR YL.FRAMAS, YR.FRAMAS
131 : C-------------------------------------------------------
132 : C********** Les CP's and CV's ***********************
133 : C-------------------------------------------------------
134 : SEGMENT GCONST
135 : REAL*8 GC(NSP)
136 : ENDSEGMENT
137 : POINTEUR CP.GCONST, CV.GCONST
138 : C-------------------------------------------------------------
139 : C******** Segments for the elementary matrixes *************
140 : C-------------------------------------------------------------
141 : SEGMENT JACEL
142 : REAL*8 JAC(3+NSP,3+NSP)
143 : ENDSEGMENT
144 : POINTEUR JTL.JACEL, JL.JACEL, JR.JACEL,
145 : & WL.JACEL
146 : c----------------------------------------
147 : SEGINI JTL
148 : SEGINI JL
149 : SEGINI JR
150 : SEGINI WL
151 : C-------------------------------------------------------------
152 : C********** Segments for the vectors ***********************
153 : C-------------------------------------------------------------
154 : SEGMENT VECEL
155 : REAL*8 VV(NSP)
156 : ENDSEGMENT
157 : POINTEUR DMLY_L.VECEL, DMLY_R.VECEL,
158 : & dmry_l.vecel, dmry_r.vecel,
159 : & dMpy_l.vecel, dMpy_r.vecel,
160 : & dMmy_l.vecel, dMmy_r.vecel,
161 : & dmiy_l.vecel, dmiy_r.vecel,
162 : & d3my_l.vecel, d3my_r.vecel,
163 : & d2my_l.vecel, d2my_r.vecel,
164 : & dPpy_l.vecel, dPpy_r.vecel,
165 : & dPmy_l.vecel, dPmy_r.vecel,
166 : & dpiy_l.vecel, dpiy_r.vecel,
167 : & dgdyl.vecel, dgdyr.vecel
168 : C----------------------------------------------
169 : SEGINI DMLY_L, DMLY_R,
170 : & dmry_l, dmry_r,
171 : & dMpy_l, dMpy_r,
172 : & dMmy_l, dMmy_r,
173 : & dmiy_l, dmiy_r,
174 : & d3my_l, d3my_r,
175 : & d2my_l, d2my_r,
176 : & dPpy_l, dPpy_r,
177 : & dPmy_l, dPmy_r,
178 : & dpiy_l, dpiy_r,
179 : & dgdyl, dgdyr
180 : C-------------------------------------------------------------
181 : C********** Segments for the flux-vector *******************
182 : C-------------------------------------------------------------
183 : SEGMENT FUNEL
184 : REAL*8 FU(3+NSP)
185 : ENDSEGMENT
186 : POINTEUR f.funel
187 : C-------------------------
188 : SEGINI f
189 : C------------------------------------------------------------
190 : SEGINI YL, YR
191 : SEGACT MPYN
192 : DO 100 I=1,(NSP-1)
193 : YL.YET(I)=MPYN.VPOCHA(NLCG,I)
194 : YR.YET(I)=MPYN.VPOCHA(NLCD,I)
195 : 100 CONTINUE
196 : C----------------------------------------
197 : SEGINI CP, CV
198 : MLRECP = LRECP
199 : MLRECV = LRECV
200 : SEGACT MLRECP, MLRECV
201 : DO 101 I=1,(NSP-1)
202 : CP.GC(I)=MLRECP.PROG(I)
203 : CV.GC(I)=MLRECV.PROG(I)
204 : 101 CONTINUE
205 : CP.GC(NSP)=MLRECP.PROG(NSP)
206 : CV.GC(NSP)=MLRECV.PROG(NSP)
207 : c-------------------------------------------------------------
208 : c Computing GAMMA at the left cell and its derivatives
209 : c with respect to the primitive variables Y_i
210 : c-------------------------------------------------------------
211 : top=0.0D0
212 : bot=0.0D0
213 : do 40 i=1,(nsp-1)
214 : top=top+yl.yet(i)*(cp.gc(i)-cp.gc(nsp))
215 : bot=bot+yl.yet(i)*(cv.gc(i)-cv.gc(nsp))
216 : 40 continue
217 : top=cp.gc(nsp)+top
218 : bot=cv.gc(nsp)+bot
219 : gal=top/bot
220 : gm1l=gal-1.0d0
221 : c-------------------------------------------------------------
222 : do 41 i=1,(nsp-1)
223 : dgdyl.vv(i)=(cp.gc(i)-cp.gc(nsp)-
224 : & gal*(cv.gc(i)-cv.gc(nsp)))/bot
225 : 41 continue
226 : c-------------------------------------------------------------
227 : c Computing GAMMA at the right cell and its derivatives
228 : c with respect to the primitive variables Y_i
229 : c-------------------------------------------------------------
230 : top=0.0D0
231 : bot=0.0D0
232 : do 42 i=1,(nsp-1)
233 : top=top+yr.yet(i)*(cp.gc(i)-cp.gc(nsp))
234 : bot=bot+yr.yet(i)*(cv.gc(i)-cv.gc(nsp))
235 : 42 continue
236 : top=cp.gc(nsp)+top
237 : bot=cv.gc(nsp)+bot
238 : gar=top/bot
239 : gm1r=gar-1.0d0
240 : c-------------------------------------------------------------
241 : do 43 i=1,(nsp-1)
242 : dgdyr.vv(i)=(cp.gc(i)-cp.gc(nsp)-
243 : & gar*(cv.gc(i)-cv.gc(nsp)))/bot
244 : 43 continue
245 : c-------------------------------------------------------------
246 : n_x=nvect(1)
247 : n_y=nvect(2)
248 : t_x=tvect(1)
249 : t_y=tvect(2)
250 : c----------------------------
251 : c----------------------------
252 : rold_l=wvec_l(1)
253 : uold_l=wvec_l(2)
254 : vold_l=wvec_l(3)
255 : pold_l=wvec_l(4)
256 : c-----------------------
257 : rold_r=wvec_r(1)
258 : uold_r=wvec_r(2)
259 : vold_r=wvec_r(3)
260 : pold_r=wvec_r(4)
261 : c------------------------------------------------------------------
262 : c Computation of the specific total energy on the left and right.
263 : c------------------------------------------------------------------
264 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
265 : eold_l=eold_l+pold_l/(gm1l*rold_l)
266 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
267 : eold_r=eold_r+pold_r/(gm1r*rold_r)
268 : c-------------------------------------------------------------------
269 : c Computation of the speed of sound and its derivatives;
270 : c numerical speed of sound at the interface is taken as an average
271 : c of the speeds of sounds of the neighbouring cells
272 : c-------------------------------------------------------------------
273 : aleft=sqrt(gal*pold_l/rold_l)
274 : arigh=sqrt(gar*pold_r/rold_r)
275 : am=0.5d0*(aleft+arigh)
276 : c-------------------------------------------------------------------
277 : damr_r=-arigh/(4.0d0*rold_r)
278 : damu_r=0.0d0
279 : damv_r=0.0d0
280 : damp_r=gar/(4.0d0*arigh*rold_r)
281 : damg_r=arigh/(4.0d0*gar)
282 : c-----------------------
283 : damr_l=-aleft/(4.0d0*rold_l)
284 : damu_l=0.0d0
285 : damv_l=0.0d0
286 : damp_l=gal/(4.0d0*aleft*rold_l)
287 : damg_l=aleft/(4.0d0*gal)
288 : c-------------------------------------------------------------------
289 : c Computing numerical Mach number and its derivatives,
290 : c see p.370, under (A1).
291 : c-------------------------------------------------------------------
292 : un_l=uold_l*n_x+vold_l*n_y
293 : un_r=uold_r*n_x+vold_r*n_y
294 : ut_l=uold_l*t_x+vold_l*t_y
295 : ut_r=uold_r*t_x+vold_r*t_y
296 : c-------------------------------------------------------------------
297 : ml=un_l/am
298 : mr=un_r/am
299 : c-------------------------------------------------------------------
299 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
300 : c see (19a) and (19b), p.367.
301 : c-------------------------------------------------------------------
302 : if(abs(ml) .ge. 1.0d0) then
303 : Mplus=(ml+abs(ml))/2.0d0
304 : else
305 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
306 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
307 : endif
308 : c-------------------------------------------------------------------
309 : if(abs(mr) .ge. 1.0d0) then
310 : Mmin=(mr-abs(mr))/2.0d0
311 : else
312 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
313 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
314 : endif
315 : c-------------------------------------------------------------------
316 : c Derivatives of ml and mr with respect to both sets of primitive
317 : c variables: left and right.
318 : c-------------------------------------------------------------------
319 : temp_l=-un_l/(am*am)
320 : temp_r=-un_r/(am*am)
321 : c--------
322 : dmlr_l=temp_l*damr_l
323 : dmlr_r=temp_l*damr_r
324 : dmrr_l=temp_r*damr_l
325 : dmrr_r=temp_r*damr_r
326 : c--------
327 : dmlu_l=n_x/am+temp_l*damu_l
328 : dmlu_r=temp_l*damu_r
329 : dmru_l=temp_r*damu_l
330 : dmru_r=n_x/am+temp_r*damu_r
331 : c--------
332 : dmlv_l=n_y/am+temp_l*damv_l
333 : dmlv_r=temp_l*damv_r
334 : dmrv_l=temp_r*damv_l
335 : dmrv_r=n_y/am+temp_r*damv_r
336 : c--------
337 : dmlp_l=temp_l*damp_l
338 : dmlp_r=temp_l*damp_r
339 : dmrp_l=temp_r*damp_l
340 : dmrp_r=temp_r*damp_r
341 : c--------
342 : do 44 i=1,(nsp-1)
343 : dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i)
344 : dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i)
345 : dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i)
346 : dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i)
347 : 44 continue
348 : c-----------------------------------------------------------
349 : c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
350 : c-----------------------------------------------------------
351 : mmid=Mplus+Mmin
352 : c-----------------------------------------------------------
353 : c Computing the derivatives of M+ and M-
354 : c-----------------------------------------------------------
355 : if(ml .ge. 1.0d0) then
356 : dMpr_l=dmlr_l
357 : dMpu_l=dmlu_l
358 : dMpv_l=dmlv_l
359 : dMpp_l=dmlp_l
360 : do 45 i=1,(nsp-1)
361 : dMpy_l.vv(i)=dmly_l.vv(i)
362 : 45 continue
363 : c--------------------
364 : dMpr_r=dmlr_r
365 : dMpu_r=dmlu_r
366 : dMpv_r=dmlv_r
367 : dMpp_r=dmlp_r
368 : do 46 i=1,(nsp-1)
369 : dMpy_r.vv(i)=dmly_r.vv(i)
370 : 46 continue
371 : else
372 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
373 : temph=(ml+1.0d0)/2.0d0
374 : dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
375 : dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
376 : dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
377 : dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
378 : do 47 i=1,(nsp-1)
379 : dMpy_l.vv(i)=(temph+4.0d0*beta*ml*
380 : & (ml*ml-1.0d0))*dmly_l.vv(i)
381 : 47 continue
382 : c--------------------
383 : dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
384 : dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
385 : dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
386 : dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
387 : do 48 i=1,(nsp-1)
388 : dMpy_r.vv(i)=(temph+4.0d0*beta*ml*
389 : & (ml*ml-1.0d0))*dmly_r.vv(i)
390 : 48 continue
391 : else
392 : dMpr_l=0.0d0
393 : dMpu_l=0.0d0
394 : dMpv_l=0.0d0
395 : dMpp_l=0.0d0
396 : do 49 i=1,(nsp-1)
397 : dMpy_l.vv(i)=0.0d0
398 : 49 continue
399 : c---------------------
400 : dMpr_r=0.0d0
401 : dMpu_r=0.0d0
402 : dMpv_r=0.0d0
403 : dMpp_r=0.0d0
404 : do 50 i=1,(nsp-1)
405 : dMpy_r.vv(i)=0.0d0
406 : 50 continue
407 : endif
408 : endif
409 : c-----------------------------------------------------------
410 : if(mr .ge. 1.0d0) then
411 : dMmr_l=0.0d0
412 : dMmu_l=0.0d0
413 : dMmv_l=0.0d0
414 : dMmp_l=0.0d0
415 : do 51 i=1,(nsp-1)
416 : dMmy_l.vv(i)=0.0d0
417 : 51 continue
418 : c---------------------
419 : dMmr_r=0.0d0
420 : dMmu_r=0.0d0
421 : dMmv_r=0.0d0
422 : dMmp_r=0.0d0
423 : do 52 i=1,(nsp-1)
424 : dMmy_r.vv(i)=0.0d0
425 : 52 continue
426 : else
427 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
428 : temph=(-mr+1.0d0)/2.0d0
429 : dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
430 : dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
431 : dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
432 : dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
433 : do 53 i=1,(nsp-1)
434 : dMmy_l.vv(i)=(temph-4.0d0*beta*mr*
435 : & (mr*mr-1.0d0))*dmry_l.vv(i)
436 : 53 continue
437 : c--------------------
438 : dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
439 : dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
440 : dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
441 : dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
442 : do 54 i=1,(nsp-1)
443 : dMmy_r.vv(i)=(temph-4.0d0*beta*mr*
444 : & (mr*mr-1.0d0))*dmry_r.vv(i)
445 : 54 continue
446 : else
447 : dMmr_l=dmrr_l
448 : dMmu_l=dmru_l
449 : dMmv_l=dmrv_l
450 : dMmp_l=dmrp_l
451 : do 55 i=1,(nsp-1)
452 : dMmy_l.vv(i)=dmry_l.vv(i)
453 : 55 continue
454 : c--------------------
455 : dMmr_r=dmrr_r
456 : dMmu_r=dmru_r
457 : dMmv_r=dmrv_r
458 : dMmp_r=dmrp_r
459 : do 56 i=1,(nsp-1)
460 : dMmy_r.vv(i)=dmry_r.vv(i)
461 : 56 continue
462 : endif
463 : endif
464 : c-----------------------------------------------------------------
465 : c computing the derivatives of m_{1/2} (notation as in the paper)
466 : c-----------------------------------------------------------------
467 : dmir_l=dMpr_l+dMmr_l
468 : dmir_r=dMpr_r+dMmr_r
469 : c-------------
470 : dmiu_l=dMpu_l+dMmu_l
471 : dmiu_r=dMpu_r+dMmu_r
472 : c-------------
473 : dmiv_l=dMpv_l+dMmv_l
474 : dmiv_r=dMpv_r+dMmv_r
475 : c-------------
476 : dmip_l=dMpp_l+dMmp_l
477 : dmip_r=dMpp_r+dMmp_r
478 : c-------------
479 : do 57 i=1,(nsp-1)
480 : dmiy_l.vv(i)=dMpy_l.vv(i)+dMmy_l.vv(i)
481 : dmiy_r.vv(i)=dMpy_r.vv(i)+dMmy_r.vv(i)
482 : 57 continue
483 : c----------------------------------------------------------------
484 : c computing the main convective variables and their derivatives
485 : c mpl_m is m^{+}_{1/2} (paper's notation) and
486 : c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
487 : c----------------------------------------------------------------
488 : if(mmid .ge. 0.0d0) then
489 : mpl_m = mmid
490 : d2mr_l=dmir_l
491 : d2mu_l=dmiu_l
492 : d2mv_l=dmiv_l
493 : d2mp_l=dmip_l
494 : do 58 i=1,(nsp-1)
495 : d2my_l.vv(i)=dmiy_l.vv(i)
496 : 58 continue
497 : c------------
498 : d2mr_r=dmir_r
499 : d2mu_r=dmiu_r
500 : d2mv_r=dmiv_r
501 : d2mp_r=dmip_r
502 : do 59 i=1,(nsp-1)
503 : d2my_r.vv(i)=dmiy_r.vv(i)
504 : 59 continue
505 : c------------
506 : else
507 : mpl_m = 0.0d0
508 : d2mr_l=0.0d0
509 : d2mu_l=0.0d0
510 : d2mv_l=0.0d0
511 : d2mp_l=0.0d0
512 : do 60 i=1,(nsp-1)
513 : d2my_l.vv(i)=0.0d0
514 : 60 continue
515 : c------------
516 : d2mr_r=0.0d0
517 : d2mu_r=0.0d0
518 : d2mv_r=0.0d0
519 : d2mp_r=0.0d0
520 : do 61 i=1,(nsp-1)
521 : d2my_r.vv(i)=0.0d0
522 : 61 continue
523 : endif
524 : c---------------------------------------------
525 : if(mmid .le. 0.0d0) then
526 : mmin_m = mmid
527 : d3mr_l=dmir_l
528 : d3mu_l=dmiu_l
529 : d3mv_l=dmiv_l
530 : d3mp_l=dmip_l
531 : do 62 i=1,(nsp-1)
532 : d3my_l.vv(i)=dmiy_l.vv(i)
533 : 62 continue
534 : c------------
535 : d3mr_r=dmir_r
536 : d3mu_r=dmiu_r
537 : d3mv_r=dmiv_r
538 : d3mp_r=dmip_r
539 : do 63 i=1,(nsp-1)
540 : d3my_r.vv(i)=dmiy_r.vv(i)
541 : 63 continue
542 : c------------
543 : else
544 : mmin_m = 0.0d0
545 : d3mr_l=0.0d0
546 : d3mu_l=0.0d0
547 : d3mv_l=0.0d0
548 : d3mp_l=0.0d0
549 : do 64 i=1,(nsp-1)
550 : d3my_l.vv(i)=0.0d0
551 : 64 continue
552 : c------------
553 : d3mr_r=0.0d0
554 : d3mu_r=0.0d0
555 : d3mv_r=0.0d0
556 : d3mp_r=0.0d0
557 : do 65 i=1,(nsp-1)
558 : d3my_r.vv(i)=0.0d0
559 : 65 continue
560 : endif
561 : c---------------------------------------------------------------
561 : c Computing the calligraphic P+ and P- with their derivatives,
562 : c see (21a) & (21b) on p.368.
563 : c---------------------------------------------------------------
564 : if(ml .ge. 1.0d0) then
565 : Pplus = 1.0d0
566 : else
567 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
568 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
569 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
570 : else
571 : Pplus = 0.0d0
572 : endif
573 : endif
574 : c---------------------------------------------------------------
575 : if(mr .ge. 1.0d0) then
576 : Pmin = 0.0d0
577 : else
578 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
579 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
580 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
581 : else
582 : Pmin = 1.0d0
583 : endif
584 : endif
585 : c---------------------------------------------------------------
586 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
587 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
588 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
589 : c--------------
590 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
591 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
592 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
593 : c---------------------------------------------------------------
594 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
595 : dPpr_l=brac_l*dmlr_l
596 : dPpr_r=brac_l*dmlr_r
597 : c------------
598 : dPpu_l=brac_l*dmlu_l
599 : dPpu_r=brac_l*dmlu_r
600 : c------------
601 : dPpv_l=brac_l*dmlv_l
602 : dPpv_r=brac_l*dmlv_r
603 : c------------
604 : dPpp_l=brac_l*dmlp_l
605 : dPpp_r=brac_l*dmlp_r
606 : c------------
607 : do 66 i=1,(nsp-1)
608 : dPpy_l.vv(i)=brac_l*dmly_l.vv(i)
609 : dPpy_r.vv(i)=brac_l*dmly_r.vv(i)
610 : 66 continue
611 : c------------
612 : else
613 : dPpr_l=0.0d0
614 : dPpr_r=0.0d0
615 : c-----------
616 : dPpu_l=0.0d0
617 : dPpu_r=0.0d0
618 : c-----------
619 : dPpv_l=0.0d0
620 : dPpv_r=0.0d0
621 : c-----------
622 : dPpp_l=0.0d0
623 : dPpp_r=0.0d0
624 : c-----------
625 : do 67 i=1,(nsp-1)
626 : dPpy_l.vv(i)=0.0d0
627 : dPpy_r.vv(i)=0.0d0
628 : 67 continue
629 : c-----------
630 : endif
631 : c---------------------------------------------------------------
632 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
633 : dPmr_l=brac_r*dmrr_l
634 : dPmr_r=brac_r*dmrr_r
635 : c------------
636 : dPmu_l=brac_r*dmru_l
637 : dPmu_r=brac_r*dmru_r
638 : c------------
639 : dPmv_l=brac_r*dmrv_l
640 : dPmv_r=brac_r*dmrv_r
641 : c------------
642 : dPmp_l=brac_r*dmrp_l
643 : dPmp_r=brac_r*dmrp_r
644 : c------------
645 : do 68 i=1,(nsp-1)
646 : dPmy_l.vv(i)=brac_r*dmry_l.vv(i)
647 : dPmy_r.vv(i)=brac_r*dmry_r.vv(i)
648 : 68 continue
649 : c------------
650 : else
651 : dPmr_l=0.0d0
652 : dPmr_r=0.0d0
653 : c-----------
654 : dPmu_l=0.0d0
655 : dPmu_r=0.0d0
656 : c-----------
657 : dPmv_l=0.0d0
658 : dPmv_r=0.0d0
659 : c-----------
660 : dPmp_l=0.0d0
661 : dPmp_r=0.0d0
662 : c-----------
663 : do 69 i=1,(nsp-1)
664 : dPmy_l.vv(i)=0.0d0
665 : dPmy_r.vv(i)=0.0d0
666 : 69 continue
667 : c-----------
668 : endif
669 : c-------------------------------------------------------------------
670 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
671 : c-------------------------------------------------------------------
672 : pmid=Pplus*pold_l + Pmin*pold_r
673 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
674 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
675 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
676 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
677 : do 70 i=1,(nsp-1)
678 : dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r
679 : 70 continue
680 : c----------------------------
681 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
682 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
683 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
684 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
685 : do 71 i=1,(nsp-1)
686 : dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r
687 : 71 continue
688 : c---------------------------------------------------------------------
689 : c computing JACOBIAN as a derivative of the numerical flux function
690 : c with respect to the primitive variables.
691 : c Notation: jl(i,j) --- is the derivative of the i-component of the
692 : c flux function with respect to the j-component of the
693 : c vector of primitive variables of the left state.
694 : c jr(i,j) --- is the derivative of the i-component of the
695 : c flux function with respect to the j-component of the
696 : c vector of primitive variables of the right state.
697 : c---------------------------------------------------------------------
698 : f.fu(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
699 : c---------------------------------------------------------------------
700 : jl.jac(1,1)=damr_l*f.fu(1)/am+am*(d2mr_l*rold_l+mpl_m)
701 : jl.jac(1,1)=jl.jac(1,1)+am*d3mr_l*rold_r
702 : jl.jac(1,2)=damu_l*f.fu(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
703 : jl.jac(1,3)=damv_l*f.fu(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
704 : jl.jac(1,4)=damp_l*f.fu(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
705 : do 72 i=1,(nsp-1)
706 : jl.jac(1,4+i)=damg_l*dgdyl.vv(i)*f.fu(1)/am+
707 : & am*(d2my_l.vv(i)*rold_l+d3my_l.vv(i)*rold_r)
708 : 72 continue
709 : c------------------------------------
710 : jr.jac(1,1)=damr_r*f.fu(1)/am+am*(d2mr_r*rold_l+mmin_m)
711 : jr.jac(1,1)=jr.jac(1,1)+am*d3mr_r*rold_r
712 : jr.jac(1,2)=damu_r*f.fu(1)/am+am*(d2mu_r*rold_l+d3mu_r*rold_r)
713 : jr.jac(1,3)=damv_r*f.fu(1)/am+am*(d2mv_r*rold_l+d3mv_r*rold_r)
714 : jr.jac(1,4)=damp_r*f.fu(1)/am+am*(d2mp_r*rold_l+d3mp_r*rold_r)
715 : do 73 i=1,(nsp-1)
716 : jr.jac(1,4+i)=damg_r*dgdyr.vv(i)*f.fu(1)/am+
717 : & am*(d2my_r.vv(i)*rold_l+d3my_r.vv(i)*rold_r)
718 : 73 continue
719 : c------------------------------------
720 : c------------------------------------
721 : br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
722 : br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
723 : f.fu(2)=n_x*(am*br1+pmid)+am*t_x*br2
724 : c------------------
725 : det=n_x*t_y-n_y*t_x
726 : c---------------------------------------------------------
727 : br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
728 : br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
729 : jl.jac(2,1)=n_x*(damr_l*br1+am*br3+dpir_l)
730 : jl.jac(2,1)=jl.jac(2,1)+t_x*damr_l*br2+am*t_x*br4
731 : c-------------------
732 : br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
733 : br3=br3+d3mu_l*rold_r*un_r
734 : br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y)/det
735 : br4=br4+d3mu_l*rold_r*ut_r
736 : jl.jac(2,2)=n_x*(damu_l*br1+am*br3+dpiu_l)
737 : jl.jac(2,2)=jl.jac(2,2)+t_x*damu_l*br2+am*t_x*br4
738 : c-------------------
739 : br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x)/det
740 : br3=br3+d3mv_l*rold_r*un_r
741 : br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
742 : br4=br4+d3mv_l*rold_r*ut_r
743 : jl.jac(2,3)=n_x*(damv_l*br1+am*br3+dpiv_l)
744 : jl.jac(2,3)=jl.jac(2,3)+t_x*damv_l*br2+am*t_x*br4
745 : c-------------------
746 : br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
747 : br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
748 : jl.jac(2,4)=n_x*(damp_l*br1+am*br3+dpip_l)
749 : jl.jac(2,4)=jl.jac(2,4)+t_x*damp_l*br2+am*t_x*br4
750 : c-------------------
751 : do 74 i=1,(nsp-1)
752 : br3=d2my_l.vv(i)*rold_l*un_l+d3my_l.vv(i)*rold_r*un_r
753 : br4=d2my_l.vv(i)*rold_l*ut_l+d3my_l.vv(i)*rold_r*ut_r
754 : jl.jac(2,4+i)=n_x*(damg_l*dgdyl.vv(i)*br1+
755 : & am*br3+dpiy_l.vv(i))
756 : jl.jac(2,4+i)=jl.jac(2,4+i)+
757 : & t_x*damg_l*dgdyl.vv(i)*br2+am*t_x*br4
758 : 74 continue
759 : c-------------------------------------------------------------
760 : br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
761 : br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
762 : jr.jac(2,1)=n_x*(damr_r*br1+am*br3+dpir_r)
763 : jr.jac(2,1)=jr.jac(2,1)+t_x*damr_r*br2+am*t_x*br4
764 : c-------------------
765 : br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
766 : br3=br3+d3mu_r*rold_r*un_r
767 : br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
768 : br4=br4+d3mu_r*rold_r*ut_r
769 : jr.jac(2,2)=n_x*(damu_r*br1+am*br3+dpiu_r)
770 : jr.jac(2,2)=jr.jac(2,2)+t_x*damu_r*br2+am*t_x*br4
771 : c-------------------
772 : br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
773 : br3=br3+d3mv_r*rold_r*un_r
774 : br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
775 : br4=br4+d3mv_r*rold_r*ut_r
776 : jr.jac(2,3)=n_x*(damv_r*br1+am*br3+dpiv_r)
777 : jr.jac(2,3)=jr.jac(2,3)+t_x*damv_r*br2+am*t_x*br4
778 : c-------------------
779 : br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
780 : br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
781 : jr.jac(2,4)=n_x*(damp_r*br1+am*br3+dpip_r)
782 : jr.jac(2,4)=jr.jac(2,4)+t_x*damp_r*br2+am*t_x*br4
783 : c-------------------
784 : do 75 i=1,(nsp-1)
785 : br3=d2my_r.vv(i)*rold_l*un_l+d3my_r.vv(i)*rold_r*un_r
786 : br4=d2my_r.vv(i)*rold_l*ut_l+d3my_r.vv(i)*rold_r*ut_r
787 : jr.jac(2,4+i)=n_x*(damg_r*dgdyr.vv(i)*br1+
788 : & am*br3+dpiy_r.vv(i))
789 : jr.jac(2,4+i)=jr.jac(2,4+i)+
790 : & t_x*damg_r*dgdyr.vv(i)*br2+am*t_x*br4
791 : 75 continue
792 : c-------------------------------------------------------------
793 : c-------------------------------------------------------------
794 : br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
795 : br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
796 : f.fu(3)=n_y*(am*br1+pmid)+am*t_y*br2
797 : br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
798 : br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
799 : jl.jac(3,1)=n_y*(damr_l*br1+am*br3+dpir_l)
800 : jl.jac(3,1)=jl.jac(3,1)+t_y*damr_l*br2+am*t_y*br4
801 : c-------------------
802 : br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
803 : br3=br3+d3mu_l*rold_r*un_r
804 : br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y/det)
805 : br4=br4+d3mu_l*rold_r*ut_r
806 : jl.jac(3,2)=n_y*(damu_l*br1+am*br3+dpiu_l)
807 : jl.jac(3,2)=jl.jac(3,2)+t_y*damu_l*br2+am*t_y*br4
808 : c-------------------
809 : br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x/det)
810 : br3=br3+d3mv_l*rold_r*un_r
811 : br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
812 : br4=br4+d3mv_l*rold_r*ut_r
813 : jl.jac(3,3)=n_y*(damv_l*br1+am*br3+dpiv_l)
814 : jl.jac(3,3)=jl.jac(3,3)+t_y*damv_l*br2+am*t_y*br4
815 : c-------------------
816 : br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
817 : br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
818 : jl.jac(3,4)=n_y*(damp_l*br1+am*br3+dpip_l)
819 : jl.jac(3,4)=jl.jac(3,4)+t_y*damp_l*br2+am*t_y*br4
820 : c-------------------
821 : do 76 i=1,(nsp-1)
822 : br3=d2my_l.vv(i)*rold_l*un_l+d3my_l.vv(i)*rold_r*un_r
823 : br4=d2my_l.vv(i)*rold_l*ut_l+d3my_l.vv(i)*rold_r*ut_r
824 : jl.jac(3,4+i)=n_y*(damg_l*dgdyl.vv(i)*br1+
825 : & am*br3+dpiy_l.vv(i))
826 : jl.jac(3,4+i)=jl.jac(3,4+i)+
827 : & t_y*damg_l*dgdyl.vv(i)*br2+am*t_y*br4
828 : 76 continue
829 : c-------------------------------------------------------------
830 : br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
831 : br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
832 : jr.jac(3,1)=n_y*(damr_r*br1+am*br3+dpir_r)
833 : jr.jac(3,1)=jr.jac(3,1)+t_y*damr_r*br2+am*t_y*br4
834 : c-------------------
835 : br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
836 : br3=br3+d3mu_r*rold_r*un_r
837 : br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
838 : br4=br4+d3mu_r*rold_r*ut_r
839 : jr.jac(3,2)=n_y*(damu_r*br1+am*br3+dpiu_r)
840 : jr.jac(3,2)=jr.jac(3,2)+t_y*damu_r*br2+am*t_y*br4
841 : c-------------------
842 : br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
843 : br3=br3+d3mv_r*rold_r*un_r
844 : br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
845 : br4=br4+d3mv_r*rold_r*ut_r
846 : jr.jac(3,3)=n_y*(damv_r*br1+am*br3+dpiv_r)
847 : jr.jac(3,3)=jr.jac(3,3)+t_y*damv_r*br2+am*t_y*br4
848 : c-------------------
849 : br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
850 : br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
851 : jr.jac(3,4)=n_y*(damp_r*br1+am*br3+dpip_r)
852 : jr.jac(3,4)=jr.jac(3,4)+t_y*damp_r*br2+am*t_y*br4
853 : c-------------------
854 : do 77 i=1,(nsp-1)
855 : br3=d2my_r.vv(i)*rold_l*un_l+d3my_r.vv(i)*rold_r*un_r
856 : br4=d2my_r.vv(i)*rold_l*ut_l+d3my_r.vv(i)*rold_r*ut_r
857 : jr.jac(3,4+i)=n_y*(damg_r*dgdyr.vv(i)*br1+
858 : & am*br3+dpiy_r.vv(i))
859 : jr.jac(3,4+i)=jr.jac(3,4+i)+
860 : & t_y*damg_r*dgdyr.vv(i)*br2+am*t_y*br4
861 : 77 continue
862 : c-------------------------------------------------------------
863 : c --------------------- f4 ----------------------------
864 : c-------------------------------------------------------------
865 : hr_l=rold_l*(uold_l*uold_l+vold_l*vold_l)/2.0d0+gal*pold_l/gm1l
866 : hr_r=rold_r*(uold_r*uold_r+vold_r*vold_r)/2.0d0+gar*pold_r/gm1r
867 : f.fu(4)=am*(mpl_m*hr_l+mmin_m*hr_r)
868 : c---------------------
869 : br1=d2mr_l*hr_l+mpl_m*(uold_l*uold_l+vold_l*vold_l)/2.0d0
870 : br1=br1+d3mr_l*hr_r
871 : jl.jac(4,1)=damr_l*f.fu(4)/am+am*br1
872 : c---------------------
873 : br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
874 : br1=br1+d3mu_l*hr_r
875 : jl.jac(4,2)=damu_l*f.fu(4)/am+am*br1
876 : c---------------------
877 : br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
878 : br1=br1+d3mv_l*hr_r
879 : jl.jac(4,3)=damv_l*f.fu(4)/am+am*br1
880 : c---------------------
881 : br1=d2mp_l*hr_l+mpl_m*gal/gm1l
882 : br1=br1+d3mp_l*hr_r
883 : jl.jac(4,4)=damp_l*f.fu(4)/am+am*br1
884 : c---------------------
885 : do 78 i=1,(nsp-1)
886 : br1=d2my_l.vv(i)*hr_l+mpl_m*(-pold_l)*dgdyl.vv(i)/(gm1l*gm1l)
887 : br1=br1+d3my_l.vv(i)*hr_r
888 : jl.jac(4,4+i)=damg_l*dgdyl.vv(i)*f.fu(4)/am+am*br1
889 : 78 continue
890 : c----------------------------------------------------------
891 : c----------------------------------------------------------
892 : br1=d2mr_r*hr_l+mmin_m*(uold_r*uold_r+vold_r*vold_r)/2.0d0
893 : br1=br1+d3mr_r*hr_r
894 : jr.jac(4,1)=damr_r*f.fu(4)/am+am*br1
895 : c---------------------
896 : br1=d2mu_r*hr_l+mmin_m*uold_r*rold_r
897 : br1=br1+d3mu_r*hr_r
898 : jr.jac(4,2)=damu_r*f.fu(4)/am+am*br1
899 : c---------------------
900 : br1=d2mv_r*hr_l+mmin_m*vold_r*rold_r
901 : br1=br1+d3mv_r*hr_r
902 : jr.jac(4,3)=damv_r*f.fu(4)/am+am*br1
903 : c---------------------
904 : br1=d2mp_r*hr_l+mmin_m*gar/gm1r
905 : br1=br1+d3mp_r*hr_r
906 : jr.jac(4,4)=damp_r*f.fu(4)/am+am*br1
907 : c---------------------
908 : do 79 i=1,(nsp-1)
909 : br1=d2my_r.vv(i)*hr_l+mmin_m*
910 : & (-pold_r)*dgdyr.vv(i)/(gm1r*gm1r)
911 : br1=br1+d3my_r.vv(i)*hr_r
912 : jr.jac(4,4+i)=damg_r*dgdyr.vv(i)*f.fu(4)/am+am*br1
913 : 79 continue
914 : c-------------------------------------------------------------
915 : c ------------------ f5++ ------------------------------------
916 : c-------------------------------------------------------------
917 : do 80 i=1,(nsp-1)
918 : f.fu(4+i)=am*(mpl_m*rold_l*yl.yet(i)+mmin_m*rold_r*yr.yet(i))
919 : c---------------------
920 : jl.jac(4+i,1)=damr_l*f.fu(4+i)/am+
921 : & am*(d2mr_l*rold_l*yl.yet(i)+mpl_m*yl.yet(i))
922 : jl.jac(4+i,1)=jl.jac(4+i,1)+am*d3mr_l*rold_r*yr.yet(i)
923 : jl.jac(4+i,2)=damu_l*f.fu(4+i)/am+am*(d2mu_l*rold_l*yl.yet(i)+
924 : & d3mu_l*rold_r*yr.yet(i))
925 : jl.jac(4+i,3)=damv_l*f.fu(4+i)/am+am*(d2mv_l*rold_l*yl.yet(i)+
926 : & d3mv_l*rold_r*yr.yet(i))
927 : jl.jac(4+i,4)=damp_l*f.fu(4+i)/am+am*(d2mp_l*rold_l*yl.yet(i)+
928 : & d3mp_l*rold_r*yr.yet(i))
929 : do 81 j=5,(4+nsp-1)
930 : if((4+i).eq.j) then
931 : jl.jac(4+i,j)=damg_l*dgdyl.vv(j-4)*f.fu(4+i)/am+
932 : & am*(d2my_l.vv(j-4)*rold_l*yl.yet(i)+mpl_m*rold_l+
933 : & d3my_l.vv(j-4)*rold_r*yr.yet(i))
934 : else
935 : jl.jac(4+i,j)=damg_l*dgdyl.vv(j-4)*f.fu(4+i)/am+
936 : & am*(d2my_l.vv(j-4)*rold_l*yl.yet(i)+
937 : & d3my_l.vv(j-4)*rold_r*yr.yet(i))
938 : endif
939 : 81 continue
940 : c------------------------------------
941 : jr.jac(4+i,1)=damr_r*f.fu(4+i)/am+
942 : & am*(d2mr_r*rold_l*yl.yet(i)+mmin_m*yr.yet(i))
943 : jr.jac(4+i,1)=jr.jac(4+i,1)+am*d3mr_r*rold_r*yr.yet(i)
944 : jr.jac(4+i,2)=damu_r*f.fu(4+i)/am+am*(d2mu_r*rold_l*yl.yet(i)+
945 : & d3mu_r*rold_r*yr.yet(i))
946 : jr.jac(4+i,3)=damv_r*f.fu(4+i)/am+am*(d2mv_r*rold_l*yl.yet(i)+
947 : & d3mv_r*rold_r*yr.yet(i))
948 : jr.jac(4+i,4)=damp_r*f.fu(4+i)/am+am*(d2mp_r*rold_l*yl.yet(i)+
949 : & d3mp_r*rold_r*yr.yet(i))
950 : do 82 j=1,(nsp-1)
951 : if(i.eq.j) then
952 : jr.jac(4+i,4+j)=damg_r*dgdyr.vv(j)*f.fu(4+i)/am+
953 : & am*(d2my_r.vv(j)*rold_l*yl.yet(i)+mmin_m*rold_r+
954 : & d3my_r.vv(j)*rold_r*yr.yet(i))
955 : else
956 : jr.jac(4+i,4+j)=damg_r*dgdyr.vv(j)*f.fu(4+i)/am+
957 : & am*(d2my_r.vv(j)*rold_l*yl.yet(i)+
958 : & d3my_r.vv(j)*rold_r*yr.yet(i))
959 : endif
960 : 82 continue
961 : 80 continue
962 : C-------------------------------------------------------------
963 : C This operation is made for the outlet boundary conditions
964 : C when the pressure at the outlet is given.
965 : C To find the flux it was assumed that rho_l = rho_r,
966 : C velocity_l=velocity_r and yl = yr
967 : C therefore we hace to modify the jacobian
968 : C-------------------------------------------------------------
969 : DO 105 I=1,(3+NSP)
970 : DO 106 J=1,(3+NSP)
971 : IF(J .EQ. 4) THEN
972 : JL.JAC(I,J) = JL.JAC(I,J)
973 : ELSE
974 : JL.JAC(I,J)=JL.JAC(I,J)+JR.JAC(I,J)
975 : ENDIF
976 : 106 CONTINUE
977 : 105 CONTINUE
978 : c-------------------------------------------------------------
979 : c matrix wl(i,j) represents the derivative of the i-component
980 : c of the vector of primitive variables of the left state with
981 : c respect to the j-component of the vector of the conservative
982 : c variables of the left state.
983 : c
984 : c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
985 : c vector of primitive variables;
986 : c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
987 : c vector of conservative variables.
988 : c-------------------------------------------------------------
989 : wl.jac(1,1)=1.0d0
990 : wl.jac(1,2)=0.0d0
991 : wl.jac(1,3)=0.0d0
992 : wl.jac(1,4)=0.0d0
993 : do 83 i=1,(nsp-1)
994 : wl.jac(1,4+i)=0.0d0
995 : 83 continue
996 : c------------------------------
997 : wl.jac(2,1)=-uold_l/rold_l
998 : wl.jac(2,2)=1.0d0/rold_l
999 : wl.jac(2,3)=0.0d0
1000 : wl.jac(2,4)=0.0d0
1001 : do 84 i=1,(nsp-1)
1002 : wl.jac(2,4+i)=0.0d0
1003 : 84 continue
1004 : c------------------------------
1005 : wl.jac(3,1)=-vold_l/rold_l
1006 : wl.jac(3,2)=0.0d0
1007 : wl.jac(3,3)=1.0d0/rold_l
1008 : wl.jac(3,4)=0.0d0
1009 : do 85 i=1,(nsp-1)
1010 : wl.jac(3,4+i)=0.0d0
1011 : 85 continue
1012 : c------------------------------
1013 : br1=0.0d0
1014 : do 86 i=1,(nsp-1)
1015 : br1=br1+dgdyl.vv(i)*yl.yet(i)
1016 : 86 continue
1017 : br1=br1*pold_l/(rold_l*gm1l)
1018 : wl.jac(4,1)=gm1l*(uold_l*uold_l+vold_l*vold_l)/2.0d0-br1
1019 : wl.jac(4,2)=-uold_l*gm1l
1020 : wl.jac(4,3)=-vold_l*gm1l
1021 : wl.jac(4,4)=gm1l
1022 : do 87 i=1,(nsp-1)
1023 : wl.jac(4,4+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l)
1024 : 87 continue
1025 : c------------------------------
1026 : do 88 i=1,(nsp-1)
1027 : do 89 j=1,4
1028 : wl.jac(4+i,j)=0.0d0
1029 : if(j.eq.1) wl.jac(4+i,j)=-yl.yet(i)/rold_l
1030 : 89 continue
1031 : c------------
1032 : do 890 j=5,(4+nsp-1)
1033 : wl.jac(4+i,j)=0.0d0
1034 : if(4+i.eq.j) then
1035 : wl.jac(4+i,j)=1.0d0/rold_l
1036 : endif
1037 : 890 continue
1038 : 88 continue
1039 : c------------------------------
1040 : c----------------------------------------------
1041 : do 1 i=1,(3+nsp)
1042 : do 2 j=1,(3+nsp)
1043 : jtl.jac(i,j)=0.0d0
1044 : do 3 k=1,(3+nsp)
1045 : jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j)
1046 : 3 continue
1047 : 2 continue
1048 : 1 continue
1049 : c----------------------------------------------------------------------
1050 : SEGSUP DMLY_L, DMLY_R,
1051 : & dmry_l, dmry_r,
1052 : & dMpy_l, dMpy_r,
1053 : & dMmy_l, dMmy_r,
1054 : & dmiy_l, dmiy_r,
1055 : & d3my_l, d3my_r,
1056 : & d2my_l, d2my_r,
1057 : & dPpy_l, dPpy_r,
1058 : & dPmy_l, dPmy_r,
1059 : & dpiy_l, dpiy_r,
1060 : & dgdyl, dgdyr
1061 : C--------------------------------------------
1062 : SEGSUP f
1063 : C--------------------------------------------
1064 : jpl = jl
1065 : SEGDES JL
1066 : SEGDES JR
1067 : SEGSUP WL
1068 : C--------------------------------------------
1069 : jll=jtl
1070 : SEGDES JTL
1071 : SEGDES YL
1072 : SEGDES YR
1073 : SEGDES CP
1074 : SEGDES CV
1075 : SEGDES MLRECP, MLRECV
1076 : C--------------------------------------------
1077 : return
1078 : end
1079 :
1080 :
1081 :
1082 :
1083 :
1084 :
1085 :
1086 :
1087 :
1088 :
1089 :
1090 :
1091 :
1092 :
1093 :
1094 :
1095 :
1096 :
1097 :
1098 :
1099 :
1100 :
1101 :
1102 :
1103 :
1104 :
1105 :
1106 :
1107 :
1108 :
1109 :
1110 :
1111 :
1112 :
© Cast3M 2003 - All rights reserved.
Disclaimer