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