1 : C CWMSP SOURCE CHAT 05/01/12 22:33:42 5004
2 : SUBROUTINE CWMSP(NSP,JLL,WVEC_L,WVEC_R,NVECT,TVECT,
3 : & mpyn,lrecp,lrecv,nlcg)
4 : C************************************************************************
5 : C
6 : C PROJET : CASTEM 2000
7 : C
8 : C NOM : CWMSP ('convection at wall for multispicies')
9 : C
10 : C DESCRIPTION : Voir KONMSP (appele par KONMSP.ESO)
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 jacobian which is the derivatives
21 : c of the numerical flux function defined at the wall cell interface
22 : c with respect to the conservative variables of the left (pre-wall)
23 : c cell.
24 : c
25 : c EQUATIONS: 2D Euler equations of gas dynamics - MULTICPECIES 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 (rho,u_x,u_y,p) at the
42 : c left cell;
43 : c
44 : c wvec_r -- vector of the primitive variables (rho,u_x,u_y,p) at the
45 : c 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 cell at wall.
61 : c----------------------------------------------------------------------
62 : c
63 : c OUTPUT:
64 : c
65 : c jll -- jakobian matrix (3+nsp) by (3+nsp) -
66 : c derivatives of the numerical
67 : c flux function with respect to the conservative variables
68 : c from the left cell;
69 : c
70 : c----------------------------------------------------------------------
71 : IMPLICIT INTEGER(I-N)
72 : integer nsp,jll,lrecp,lrecv,nlcg
73 : real*8 wvec_l(4),wvec_r(4)
74 : real*8 nvect(2),tvect(2)
75 : real*8 alpha,beta
76 : real*8 gal,gar,gm1l,gm1r
77 : real*8 n_x,n_y
78 : real*8 un_l, un_r
79 : real*8 ml,mr,Mplus,Mmin,am
80 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l
81 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r
82 : real*8 Pplus,Pmin
83 : real*8 top,bot
84 : real*8 br1,temp_l,temp_r,brac_l,brac_r
85 : real*8 aleft, arigh
86 : real*8 damr_l,damr_r,damu_l,damu_r
87 : real*8 damv_l,damv_r,damp_l,damp_r
88 : real*8 damg_l,damg_r
89 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
90 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
91 : real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
92 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
93 : real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
94 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
95 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
96 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
97 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
98 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
99 : real*8 tcoef, bcoef
100 : integer i,j,k
101 : parameter(alpha = 0.1875D0, beta = 0.125D0)
102 : C------------------------------------------------------------
103 : -INC SMCHPOI
104 : POINTEUR MPYN.MPOVAL
105 : C-------------------------------------------------------------
106 : -INC SMLREEL
107 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
108 : C-------------------------------------------------------------
109 : C******* Les fractionines messiques **************************
110 : C-------------------------------------------------------------
111 : SEGMENT FRAMAS
112 : REAL*8 YET(NSP)
113 : ENDSEGMENT
114 : POINTEUR YL.FRAMAS, YR.FRAMAS
115 : C-------------------------------------------------------
116 : C********** Les CP's and CV's ***********************
117 : C-------------------------------------------------------
118 : SEGMENT GCONST
119 : REAL*8 GC(NSP)
120 : ENDSEGMENT
121 : POINTEUR CP.GCONST, CV.GCONST
122 : C-------------------------------------------------------------
123 : C******** Segments for the elementary matrixes *************
124 : C-------------------------------------------------------------
125 : SEGMENT JACEL
126 : REAL*8 JAC(3+NSP,3+NSP)
127 : ENDSEGMENT
128 : POINTEUR JTL.JACEL, JTR.JACEL, JL.JACEL, JR.JACEL,
129 : & WL.JACEL, WR.JACEL
130 : C-------------------------------------------------------------
131 : C********** Segments for the vectors ***********************
132 : C-------------------------------------------------------------
133 : SEGMENT VECEL
134 : REAL*8 VV(NSP)
135 : ENDSEGMENT
136 : POINTEUR dmly_l.vecel, dmly_r.vecel,
137 : & dmry_l.vecel, dmry_r.vecel,
138 : & dPpy_l.vecel, dPpy_r.vecel,
139 : & dPmy_l.vecel, dPmy_r.vecel,
140 : & dpiy_l.vecel, dpiy_r.vecel,
141 : & dgdyl.vecel, dgdyr.vecel
142 : C-------------------------------------------
143 : SEGINI dmly_l, dmly_r,
144 : & dmry_l, dmry_r,
145 : & dPpy_l, dPpy_r,
146 : & dPmy_l, dPmy_r,
147 : & dpiy_l, dpiy_r,
148 : & dgdyl, dgdyr
149 : C------------------------------------------------------------
150 : SEGINI YL, YR
151 : SEGACT MPYN
152 : DO 100 I=1,(NSP-1)
153 : YL.YET(I)=MPYN.VPOCHA(NLCG,I)
154 : YR.YET(I)=MPYN.VPOCHA(NLCG,I)
155 : 100 CONTINUE
156 : C----------------------------------------
157 : SEGINI CP, CV
158 : MLRECP = LRECP
159 : MLRECV = LRECV
160 : SEGACT MLRECP, MLRECV
161 : DO 101 I=1,(NSP-1)
162 : CP.GC(I)=MLRECP.PROG(I)
163 : CV.GC(I)=MLRECV.PROG(I)
164 : 101 CONTINUE
165 : CP.GC(NSP)=MLRECP.PROG(NSP)
166 : CV.GC(NSP)=MLRECV.PROG(NSP)
167 : c-------------------------------------------------------------
168 : c Computing GAMMA at the left cell and its derivatives
169 : c with respect to the primitive variables Y_i
170 : c-------------------------------------------------------------
171 : top=0.0D0
172 : bot=0.0D0
173 : do 40 i=1,(nsp-1)
174 : top=top+yl.yet(i)*(cp.gc(i)-cp.gc(nsp))
175 : bot=bot+yl.yet(i)*(cv.gc(i)-cv.gc(nsp))
176 : 40 continue
177 : top=cp.gc(nsp)+top
178 : bot=cv.gc(nsp)+bot
179 : gal=top/bot
180 : gm1l=gal-1.0d0
181 : c-------------------------------------------------------------
182 : do 41 i=1,(nsp-1)
183 : dgdyl.vv(i)=(cp.gc(i)-cp.gc(nsp)-
184 : & gal*(cv.gc(i)-cv.gc(nsp)))/bot
185 : 41 continue
186 : c-------------------------------------------------------------
187 : c Computing GAMMA at the right cell and its derivatives
188 : c with respect to the primitive variables Y_i
189 : c-------------------------------------------------------------
190 : top=0.0D0
191 : bot=0.0D0
192 : do 42 i=1,(nsp-1)
193 : top=top+yr.yet(i)*(cp.gc(i)-cp.gc(nsp))
194 : bot=bot+yr.yet(i)*(cv.gc(i)-cv.gc(nsp))
195 : 42 continue
196 : top=cp.gc(nsp)+top
197 : bot=cv.gc(nsp)+bot
198 : gar=top/bot
199 : gm1r=gar-1.0d0
200 : c-------------------------------------------------------------
201 : do 43 i=1,(nsp-1)
202 : dgdyr.vv(i)=(cp.gc(i)-cp.gc(nsp)-
203 : & gar*(cv.gc(i)-cv.gc(nsp)))/bot
204 : 43 continue
205 : c-------------------------------------------------------------
206 : n_x=nvect(1)
207 : n_y=nvect(2)
208 : c----------------------------
209 : c----------------------------
210 : rold_l=wvec_l(1)
211 : uold_l=wvec_l(2)
212 : vold_l=wvec_l(3)
213 : pold_l=wvec_l(4)
214 : c-----------------------
215 : rold_r=wvec_r(1)
216 : uold_r=wvec_r(2)
217 : vold_r=wvec_r(3)
218 : pold_r=wvec_r(4)
219 : c------------------------------------------------------------------
220 : c Computation of the specific total energy on the left and right.
221 : c------------------------------------------------------------------
222 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
223 : eold_l=eold_l+pold_l/(gm1l*rold_l)
224 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
225 : eold_r=eold_r+pold_r/(gm1r*rold_r)
226 : c-------------------------------------------------------------------
227 : c Computation of the speed of sound and its derivatives;
228 : c numerical speed of sound at the interface is taken as an average
229 : c of the speeds of sounds of the neighbouring cells
230 : c-------------------------------------------------------------------
231 : aleft=sqrt(gal*pold_l/rold_l)
232 : arigh=sqrt(gar*pold_r/rold_r)
233 : am=0.5d0*(aleft+arigh)
234 : c-------------------------------------------------------------------
235 : damr_r=-arigh/(4.0d0*rold_r)
236 : damu_r=0.0d0
237 : damv_r=0.0d0
238 : damp_r=gar/(4.0d0*arigh*rold_r)
239 : damg_r=arigh/(4.0d0*gar)
240 : c-----------------------
241 : damr_l=-aleft/(4.0d0*rold_l)
242 : damu_l=0.0d0
243 : damv_l=0.0d0
244 : damp_l=gal/(4.0d0*aleft*rold_l)
245 : damg_l=aleft/(4.0d0*gal)
246 : c-------------------------------------------------------------------
247 : c Computing numerical Mach number and its derivatives,
248 : c see p.370, under (A1).
249 : c-------------------------------------------------------------------
250 : un_l=uold_l*n_x+vold_l*n_y
251 : un_r=uold_r*n_x+vold_r*n_y
252 : c-------------------------------------------------------------------
253 : ml=un_l/am
254 : mr=un_r/am
255 : c-------------------------------------------------------------------
255 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
256 : c see (19a) and (19b), p.367.
257 : c-------------------------------------------------------------------
258 : if(abs(ml) .ge. 1.0d0) then
259 : Mplus=(ml+abs(ml))/2.0d0
260 : else
261 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
262 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
263 : endif
264 : c-------------------------------------------------------------------
265 : if(abs(mr) .ge. 1.0d0) then
266 : Mmin=(mr-abs(mr))/2.0d0
267 : else
268 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
269 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
270 : endif
271 : c-------------------------------------------------------------------
272 : c Derivatives of ml and mr with respect to both sets of primitive
273 : c variables: left and right.
274 : c-------------------------------------------------------------------
275 : temp_l=-un_l/(am*am)
276 : temp_r=-un_r/(am*am)
277 : c--------
278 : dmlr_l=temp_l*damr_l
279 : dmlr_r=temp_l*damr_r
280 : dmrr_l=temp_r*damr_l
281 : dmrr_r=temp_r*damr_r
282 : c--------
283 : dmlu_l=n_x/am+temp_l*damu_l
284 : dmlu_r=temp_l*damu_r
285 : dmru_l=temp_r*damu_l
286 : dmru_r=n_x/am+temp_r*damu_r
287 : c--------
288 : dmlv_l=n_y/am+temp_l*damv_l
289 : dmlv_r=temp_l*damv_r
290 : dmrv_l=temp_r*damv_l
291 : dmrv_r=n_y/am+temp_r*damv_r
292 : c--------
293 : dmlp_l=temp_l*damp_l
294 : dmlp_r=temp_l*damp_r
295 : dmrp_l=temp_r*damp_l
296 : dmrp_r=temp_r*damp_r
297 : c---------------------------------
298 : do 44 i=1,(nsp-1)
299 : dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i)
300 : dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i)
301 : dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i)
302 : dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i)
303 : 44 continue
304 : c---------------------------------------------------------------
304 : c Computing the calligraphic P+ and P- with their derivatives,
305 : c see (21a) & (21b) on p.368.
306 : c---------------------------------------------------------------
307 : if(ml .ge. 1.0d0) then
308 : Pplus = 1.0d0
309 : else
310 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
311 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
312 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
313 : else
314 : Pplus = 0.0d0
315 : endif
316 : endif
317 : c---------------------------------------------------------------
318 : if(mr .ge. 1.0d0) then
319 : Pmin = 0.0d0
320 : else
321 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
322 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
323 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
324 : else
325 : Pmin = 1.0d0
326 : endif
327 : endif
328 : c---------------------------------------------------------------
329 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
330 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
331 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
332 : c--------------
333 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
334 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
335 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
336 : c---------------------------------------------------------------
337 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
338 : dPpr_l=brac_l*dmlr_l
339 : dPpr_r=brac_l*dmlr_r
340 : c------------
341 : dPpu_l=brac_l*dmlu_l
342 : dPpu_r=brac_l*dmlu_r
343 : c------------
344 : dPpv_l=brac_l*dmlv_l
345 : dPpv_r=brac_l*dmlv_r
346 : c------------
347 : dPpp_l=brac_l*dmlp_l
348 : dPpp_r=brac_l*dmlp_r
349 : c----------------------------------------------
350 : do 66 i=1,(nsp-1)
351 : dPpy_l.vv(i)=brac_l*dmly_l.vv(i)
352 : dPpy_r.vv(i)=brac_l*dmly_r.vv(i)
353 : 66 continue
354 : c------------
355 : else
356 : dPpr_l=0.0d0
357 : dPpr_r=0.0d0
358 : c-----------
359 : dPpu_l=0.0d0
360 : dPpu_r=0.0d0
361 : c-----------
362 : dPpv_l=0.0d0
363 : dPpv_r=0.0d0
364 : c-----------
365 : dPpp_l=0.0d0
366 : dPpp_r=0.0d0
367 : c-----------
368 : do 67 i=1,(nsp-1)
369 : dPpy_l.vv(i)=0.0d0
370 : dPpy_r.vv(i)=0.0d0
371 : 67 continue
372 : c-----------
373 : endif
374 : c---------------------------------------------------------------
375 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
376 : dPmr_l=brac_r*dmrr_l
377 : dPmr_r=brac_r*dmrr_r
378 : c------------
379 : dPmu_l=brac_r*dmru_l
380 : dPmu_r=brac_r*dmru_r
381 : c------------
382 : dPmv_l=brac_r*dmrv_l
383 : dPmv_r=brac_r*dmrv_r
384 : c------------
385 : dPmp_l=brac_r*dmrp_l
386 : dPmp_r=brac_r*dmrp_r
387 : c------------
388 : do 68 i=1,(nsp-1)
389 : dPmy_l.vv(i)=brac_r*dmry_l.vv(i)
390 : dPmy_r.vv(i)=brac_r*dmry_r.vv(i)
391 : 68 continue
392 : c------------
393 : else
394 : dPmr_l=0.0d0
395 : dPmr_r=0.0d0
396 : c-----------
397 : dPmu_l=0.0d0
398 : dPmu_r=0.0d0
399 : c-----------
400 : dPmv_l=0.0d0
401 : dPmv_r=0.0d0
402 : c-----------
403 : dPmp_l=0.0d0
404 : dPmp_r=0.0d0
405 : c-----------
406 : do 69 i=1,(nsp-1)
407 : dPmy_l.vv(i)=0.0d0
408 : dPmy_r.vv(i)=0.0d0
409 : 69 continue
410 : c-----------
411 : endif
412 : c-------------------------------------------------------------------
413 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
414 : c-------------------------------------------------------------------
415 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
416 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
417 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
418 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
419 : do 70 i=1,(nsp-1)
420 : dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r
421 : 70 continue
422 : c----------------------------
423 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
424 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
425 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
426 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
427 : do 71 i=1,(nsp-1)
428 : dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r
429 : 71 continue
430 : c---------------------------------------------------------------------
431 : c computing JACOBIAN as a derivative of the numerical flux function
432 : c with respect to the primitive variables.
433 : c Notation: jl(i,j) --- is the derivative of the i-component of the
434 : c flux function with respect to the j-component of the
435 : c vector of primitive variables of the left state.
436 : c jr(i,j) --- is the derivative of the i-component of the
437 : c flux function with respect to the j-component of the
438 : c vector of primitive variables of the right state.
439 : c---------------------------------------------------------------------
440 : SEGINI JL, JR
441 : c---------------------------------------------------------------------
442 : jl.jac(1,1)=0.0d0
443 : jl.jac(1,2)=0.0d0
444 : jl.jac(1,3)=0.0d0
445 : jl.jac(1,4)=0.0d0
446 : do 72 i=1,(nsp-1)
447 : jl.jac(1,4+i)=0.0d0
448 : 72 continue
449 : c------------------------------------
450 : jr.jac(1,1)=0.0d0
451 : jr.jac(1,2)=0.0d0
452 : jr.jac(1,3)=0.0d0
453 : jr.jac(1,4)=0.0d0
454 : do 73 i=1,(nsp-1)
455 : jr.jac(1,4+i)=0.0d0
456 : 73 continue
457 : c------------------------------------
458 : c------------------------------------
459 : c---------------------------------------------------------
460 : jl.jac(2,1)=n_x*dpir_l
461 : c-------------------
462 : jl.jac(2,2)=n_x*dpiu_l
463 : c-------------------
464 : jl.jac(2,3)=n_x*dpiv_l
465 : c-------------------
466 : jl.jac(2,4)=n_x*dpip_l
467 : c-------------------
468 : do 74 i=1,(nsp-1)
469 : jl.jac(2,4+i)=n_x*dpiy_l.vv(i)
470 : 74 continue
471 : c-------------------------------------------------------------
472 : jr.jac(2,1)=n_x*dpir_r
473 : c-------------------
474 : jr.jac(2,2)=n_x*dpiu_r
475 : c-------------------
476 : jr.jac(2,3)=n_x*dpiv_r
477 : c-------------------
478 : jr.jac(2,4)=n_x*dpip_r
479 : c-------------------
480 : do 75 i=1,(nsp-1)
481 : jr.jac(2,4+i)=n_x*dpiy_r.vv(i)
482 : 75 continue
483 : c-------------------------------------------------------------
484 : c-------------------------------------------------------------
485 : jl.jac(3,1)=n_y*dpir_l
486 : c-------------------
487 : jl.jac(3,2)=n_y*dpiu_l
488 : c-------------------
489 : jl.jac(3,3)=n_y*dpiv_l
490 : c-------------------
491 : jl.jac(3,4)=n_y*dpip_l
492 : c-------------------
493 : do 76 i=1,(nsp-1)
494 : jl.jac(3,4+i)=n_y*dpiy_l.vv(i)
495 : 76 continue
496 : c-------------------------------------------------------------
497 : jr.jac(3,1)=n_y*dpir_r
498 : c-------------------
499 : jr.jac(3,2)=n_y*dpiu_r
500 : c-------------------
501 : jr.jac(3,3)=n_y*dpiv_r
502 : c-------------------
503 : jr.jac(3,4)=n_y*dpip_r
504 : c-------------------
505 : do 77 i=1,(nsp-1)
506 : jr.jac(3,4+i)=n_y*dpiy_r.vv(i)
507 : 77 continue
508 : c-------------------------------------------------------------
509 : c ---------------------- f4 -------------------------------
510 : c-------------------------------------------------------------
511 : c---------------------
512 : jl.jac(4,1)=0.0d0
513 : c---------------------
514 : jl.jac(4,2)=0.0d0
515 : c---------------------
516 : jl.jac(4,3)=0.0d0
517 : c---------------------
518 : jl.jac(4,4)=0.0d0
519 : c---------------------
520 : do 78 i=1,(nsp-1)
521 : jl.jac(4,4+i)=0.0d0
522 : 78 continue
523 : c----------------------------------------------------------
524 : c----------------------------------------------------------
525 : jr.jac(4,1)=0.0d0
526 : c---------------------
527 : jr.jac(4,2)=0.0d0
528 : c---------------------
529 : jr.jac(4,3)=0.0d0
530 : c---------------------
531 : jr.jac(4,4)=0.0d0
532 : c---------------------
533 : do 79 i=1,(nsp-1)
534 : jr.jac(4,4+i)=0.0d0
535 : 79 continue
536 : c-------------------------------------------------------------
537 : c ---------------------- f5 --------------------------------
538 : c-------------------------------------------------------------
539 : do 80 i=1,(nsp-1)
540 : c---------------------
541 : jl.jac(4+i,1)=0.0d0
542 : jl.jac(4+i,2)=0.0d0
543 : jl.jac(4+i,3)=0.0d0
544 : jl.jac(4+i,4)=0.0d0
545 : do 81 j=5,(4+nsp-1)
546 : jl.jac(4+i,j)=0.0d0
547 : 81 continue
548 : c------------------------------------
549 : jr.jac(4+i,1)=0.0d0
550 : jr.jac(4+i,2)=0.0d0
551 : jr.jac(4+i,3)=0.0d0
552 : jr.jac(4+i,4)=0.0d0
553 : do 82 j=1,(nsp-1)
554 : jr.jac(4+i,4+j)=0.0d0
555 : 82 continue
556 : 80 continue
557 : c-------------------------------------------------------------
558 : c matrix wl(i,j) represents the derivative of the i-component
559 : c of the vector of primitive variables of the left state with
560 : c respect to the j-component of the vector of the conservative
561 : c variables of the left state.
562 : c
563 : c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
564 : c vector of primitive variables;
565 : c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
566 : c vector of conservative variables.
567 : c-------------------------------------------------------------
568 : SEGINI WL, WR
569 : c-------------------------------------------------------------
570 : wl.jac(1,1)=1.0d0
571 : wl.jac(1,2)=0.0d0
572 : wl.jac(1,3)=0.0d0
573 : wl.jac(1,4)=0.0d0
574 : do 83 i=1,(nsp-1)
575 : wl.jac(1,4+i)=0.0d0
576 : 83 continue
577 : c------------------------------
578 : wl.jac(2,1)=-uold_l/rold_l
579 : wl.jac(2,2)=1.0d0/rold_l
580 : wl.jac(2,3)=0.0d0
581 : wl.jac(2,4)=0.0d0
582 : do 84 i=1,(nsp-1)
583 : wl.jac(2,4+i)=0.0d0
584 : 84 continue
585 : c------------------------------
586 : wl.jac(3,1)=-vold_l/rold_l
587 : wl.jac(3,2)=0.0d0
588 : wl.jac(3,3)=1.0d0/rold_l
589 : wl.jac(3,4)=0.0d0
590 : do 85 i=1,(nsp-1)
591 : wl.jac(3,4+i)=0.0d0
592 : 85 continue
593 : c------------------------------
594 : br1=0.0d0
595 : do 86 i=1,(nsp-1)
596 : br1=br1+dgdyl.vv(i)*yl.yet(i)
597 : 86 continue
598 : br1=br1*pold_l/(rold_l*gm1l)
599 : wl.jac(4,1)=gm1l*(uold_l*uold_l+vold_l*vold_l)/2.0d0-br1
600 : wl.jac(4,2)=-uold_l*gm1l
601 : wl.jac(4,3)=-vold_l*gm1l
602 : wl.jac(4,4)=gm1l
603 : do 87 i=1,(nsp-1)
604 : wl.jac(4,4+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l)
605 : 87 continue
606 : c------------------------------
607 : do 88 i=1,(nsp-1)
608 : do 89 j=1,4
609 : wl.jac(4+i,j)=0.0d0
610 : if(j.eq.1) wl.jac(4+i,j)=-yl.yet(i)/rold_l
611 : 89 continue
612 : c------------
613 : do 890 j=5,(4+nsp-1)
614 : wl.jac(4+i,j)=0.0d0
615 : if(4+i.eq.j) then
616 : wl.jac(4+i,j)=1.0d0/rold_l
617 : endif
618 : 890 continue
619 : 88 continue
620 : c------------------------------
621 : c------------------------------
622 : wr.jac(1,1)=1.0d0
623 : wr.jac(1,2)=0.0d0
624 : wr.jac(1,3)=0.0d0
625 : wr.jac(1,4)=0.0d0
626 : do 90 i=1,(nsp-1)
627 : wr.jac(1,4+i)=0.0d0
628 : 90 continue
629 : c------------------------------
630 : wr.jac(2,1)=-uold_r/rold_r
631 : wr.jac(2,2)=1.0d0/rold_r
632 : wr.jac(2,3)=0.0d0
633 : wr.jac(2,4)=0.0d0
634 : do 91 i=1,(nsp-1)
635 : wr.jac(2,4+i)=0.0d0
636 : 91 continue
637 : c------------------------------
638 : wr.jac(3,1)=-vold_r/rold_r
639 : wr.jac(3,2)=0.0d0
640 : wr.jac(3,3)=1.0d0/rold_r
641 : wr.jac(3,4)=0.0d0
642 : do 92 i=1,(nsp-1)
643 : wr.jac(3,4+i)=0.0d0
644 : 92 continue
645 : c------------------------------
646 : br1=0.0d0
647 : do 860 i=1,(nsp-1)
648 : br1=br1+dgdyr.vv(i)*yr.yet(i)
649 : 860 continue
650 : br1=br1*pold_r/(rold_r*gm1r)
651 : wr.jac(4,1)=gm1r*(uold_r*uold_r+vold_r*vold_r)/2.0d0-br1
652 : wr.jac(4,2)=-uold_r*gm1r
653 : wr.jac(4,3)=-vold_r*gm1r
654 : wr.jac(4,4)=gm1r
655 : do 93 i=1,(nsp-1)
656 : wr.jac(4,4+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r)
657 : 93 continue
658 : c------------------------------
659 : do 94 i=1,(nsp-1)
660 : do 95 j=1,4
661 : wr.jac(4+i,j)=0.0d0
662 : if(j.eq.1) wr.jac(4+i,j)=-yr.yet(i)/rold_r
663 : 95 continue
664 : c----------------------
665 : do 950 j=5,(4+nsp-1)
666 : wr.jac(4+i,j)=0.0d0
667 : if(4+i.eq.j) wr.jac(4+i,j)=1.0d0/rold_r
668 : 950 continue
669 : 94 continue
670 : c----------------------------------------------
671 : SEGINI JTL, JTR
672 : c----------------------------------------------
673 : do 1 i=1,(3+nsp)
674 : do 2 j=1,(3+nsp)
675 : jtl.jac(i,j)=0.0d0
676 : jtr.jac(i,j)=0.0d0
677 : do 3 k=1,(3+nsp)
678 : jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j)
679 : jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j)
680 : 3 continue
681 : 2 continue
682 : 1 continue
683 : c----------------------------------------------------------------------
684 : tcoef=nvect(1)*tvect(2)+tvect(1)*nvect(2)
685 : bcoef=nvect(1)*tvect(2)-tvect(1)*nvect(2)
686 : c----------------------------------------------------------------------
687 : do 11 i=1,(3+nsp)
688 : jtl.jac(i,1)=jtl.jac(i,1)+jtr.jac(i,1)
689 : jtl.jac(i,2)=jtl.jac(i,2)+jtr.jac(i,2)*(-tcoef/bcoef)+
690 : & jtr.jac(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef
691 : jtl.jac(i,3)=jtl.jac(i,3)+jtr.jac(i,2)*
692 : & (-2.0d0*nvect(2)*tvect(2)/bcoef)+
693 : & jtr.jac(i,3)*tcoef/bcoef
694 : jtl.jac(i,4)=jtl.jac(i,4)+jtr.jac(i,4)
695 : do 777 j=5,(3+nsp)
696 : jtl.jac(i,j)=jtl.jac(i,j)+jtr.jac(i,j)
697 : 777 continue
698 : 11 continue
699 : c----------------------------------------------------------------------
700 : SEGDES JL
701 : SEGDES JR
702 : SEGDES WL
703 : SEGDES WR
704 : SEGDES JTR
705 : jll=jtl
706 : SEGDES JTL
707 : SEGDES YL
708 : SEGDES YR
709 : SEGDES CP
710 : SEGDES CV
711 : SEGDES MLRECP, MLRECV
712 : C----------------------------------------------------------------------
713 : SEGDES dmly_l, dmly_r,
714 : & dmry_l, dmry_r,
715 : & dPpy_l, dPpy_r,
716 : & dPmy_l, dPmy_r,
717 : & dpiy_l, dpiy_r,
718 : & dgdyl, dgdyr
719 : C----------------------------------------------------------------------
720 : return
721 : end
722 :
723 :
724 :
725 :
726 :
727 :
728 :
729 :
730 :
731 :
732 :
733 :
734 :
735 :
736 :
737 :
738 :
739 :
740 :
741 :
742 :
743 :
744 :
745 :
746 :
747 :
748 :
749 :
750 :
751 :
752 :
753 :
754 :
© Cast3M 2003 - All rights reserved.
Disclaimer