1 : C CONJP5 SOURCE CHAT 05/01/12 22:17:49 5004
2 : c Voir konjp4.eso
3 : c----------------------------------------------------------------------
4 : c GENERAL DESCRIPTION:
5 : c This subroutine provides the jacobians which are the derivatives
6 : c of the numerical flux function defined at "wall" cell interface
7 : c with respect to the primitive variables of the left and right
8 : c cells (relative to the cell interface).
9 : c
10 : c EQUATIONS: 3D Euler equations of gas dynamics
11 : c
12 : c
13 : c REFERENCE: JCP, 129, 364-382 (1996)
14 : c " A Sequel to AUSM: AUSM+ ".
15 : c----------------------------------------------------------------------
16 : c INPUT:
17 : c
18 : c alpha -- parameter of the AUSM+ scheme in the Pressure function;
19 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
20 : c
21 : c beta -- parameter of the AUSM+ scheme in the Mach function;
22 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
23 : c
24 : c wvec_l -- vector of the primitive variables (rho,ux,uy,uz,p) at the
25 : c left cell;
26 : c
27 : c wvec_r -- vector of the primitive variables (rho,ux,uy,uz,p) at the
28 : c right cell;
29 : c
30 : c nvect -- normal vector to the interface (3 components in 3D);
31 : c
32 : c tvec1 -- tangential vector to the interface;
33 : c
34 : c tvec2 -- tangential vector to the interface;
35 : c
36 : c ga -- ratio of the specific heats (assumed constant)
37 : c----------------------------------------------------------------------
38 : c
39 : c OUTPUT:
40 : c
41 : c jl -- jakobian matrix 5 by 5 - derivatives of the numerical
42 : c flux function with respect to the primitive variables
43 : c from the left cell;
44 : c
45 : c jr -- jakobian matrix 5 by 5 - derivatives of the numerical
46 : c flux function with respect to the primitive variables
47 : c from the right cell.
48 : c----------------------------------------------------------------------
49 : SUBROUTINE CONJP5(JL,JR,WVEC_L,WVEC_R,NVECT,TVECT1,TVECT2,
50 : & ga)
51 : IMPLICIT INTEGER(I-N)
52 : real*8 wvec_l(5),wvec_r(5)
53 : real*8 nvect(3),tvect1(3),tvect2(3)
54 : real*8 jl(5,5),jr(5,5)
55 : real*8 alpha,beta
56 : real*8 ga, gm1
57 : real*8 n_x,n_y,n_z,t1_x,t1_y,t1_z,t2_x,t2_y,t2_z
58 : real*8 un_l,un_r
59 : real*8 ml,mr,Mplus,Mmin
60 : real*8 am
61 : real*8 rold_l,uold_l,vold_l,wold_l,pold_l,eold_l
62 : real*8 rold_r,uold_r,vold_r,wold_r,pold_r,eold_r
63 : real*8 Pplus,Pmin
64 : real*8 temp_l,temp_r,brac_l,brac_r
65 : real*8 aleft, arigh
66 : real*8 damr_l,damr_r,damu_l,damu_r
67 : real*8 damv_l,damv_r,damp_l,damp_r
68 : real*8 damw_l,damw_r
69 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
70 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
71 : real*8 dmlw_l,dmlw_r
72 : real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
73 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
74 : real*8 dmrw_l,dmrw_r
75 : real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
76 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
77 : real*8 dPpw_l,dPpw_r
78 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
79 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
80 : real*8 dPmw_l,dPmw_r
81 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
82 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
83 : real*8 dpiw_l,dpiw_r
84 : real*8 c11,c12,c13,c21,c22,c23,c31,c32,c33,det
85 : real*8 zc11,zc12,zc13,zc21,zc22,zc23,zc31,zc32,zc33
86 : integer i
87 : parameter(alpha = 0.1875D0, beta = 0.125D0)
88 : c-----------------------
89 : gm1=ga-1.0d0
90 : c-----------------------
91 : n_x=nvect(1)
92 : n_y=nvect(2)
93 : n_z=nvect(3)
94 : c-----------------------
95 : t1_x=tvect1(1)
96 : t1_y=tvect1(2)
97 : t1_z=tvect1(3)
98 : c-----------------------
99 : t2_x=tvect2(1)
100 : t2_y=tvect2(2)
101 : t2_z=tvect2(3)
102 : c--------------------------------------
103 : rold_l=wvec_l(1)
104 : uold_l=wvec_l(2)
105 : vold_l=wvec_l(3)
106 : wold_l=wvec_l(4)
107 : pold_l=wvec_l(5)
108 : c-----------------------
109 : rold_r=wvec_r(1)
110 : uold_r=wvec_r(2)
111 : vold_r=wvec_r(3)
112 : wold_r=wvec_r(4)
113 : pold_r=wvec_r(5)
114 : c------------------------------------------------------------------
115 : c Computation of the specific total energy on the left and right.
116 : c------------------------------------------------------------------
117 : eold_l=(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
118 : eold_l=eold_l+pold_l/(gm1*rold_l)
119 : eold_r=(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
120 : eold_r=eold_r+pold_r/(gm1*rold_r)
121 : c-------------------------------------------------------------------
122 : c Computation of the speed of sound and its derivatives;
123 : c numerical speed of sound at the interface is taken as an average
124 : c of the speeds of sounds of the neighbouring cells
125 : c---------------------------------------------------------------------
126 : aleft=dsqrt(ga*pold_l/rold_l)
127 : arigh=dsqrt(ga*pold_r/rold_r)
128 : am=0.5d0*(aleft+arigh)
129 : c--------------------------------------------------------------------
130 : damr_r=-arigh/(4.0d0*rold_r)
131 : damu_r=0.0d0
132 : damv_r=0.0d0
133 : damw_r=0.0d0
134 : damp_r=ga/(4.0d0*arigh*rold_r)
135 : c-----------------------
136 : damr_l=-aleft/(4.0d0*rold_l)
137 : damu_l=0.0d0
138 : damv_l=0.0d0
139 : damw_l=0.0d0
140 : damp_l=ga/(4.0d0*aleft*rold_l)
141 : c-------------------------------------------------------------------
142 : c Computing numerical Mach number and its derivatives,
143 : c see p.370, under (A1).
144 : c-------------------------------------------------------------------
145 : un_l=uold_l*n_x+vold_l*n_y+wold_l*n_z
146 : un_r=uold_r*n_x+vold_r*n_y+wold_r*n_z
147 : c----------------------------------------
148 : ml=un_l/am
149 : mr=un_r/am
150 : c-------------------------------------------------------------------
150 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
151 : c see (19a) and (19b), p.367.
152 : c-------------------------------------------------------------------
153 : if(dabs(ml) .ge. 1.0d0) then
154 : Mplus=(ml+dabs(ml))/2.0d0
155 : else
156 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
157 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
158 : endif
159 : c-----------------
160 : if(dabs(mr) .ge. 1.0d0) then
161 : Mmin=(mr-dabs(mr))/2.0d0
162 : else
163 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
164 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
165 : endif
166 : c-------------------------------------------------------------------
167 : c Derivatives of ml and mr with respect to both sets of primitive
168 : c variables: left and right.
169 : c-------------------------------------------------------------------
170 : temp_l=-un_l/(am*am)
171 : temp_r=-un_r/(am*am)
172 : c--------
173 : dmlr_l=temp_l*damr_l
174 : dmlr_r=temp_l*damr_r
175 : dmrr_l=temp_r*damr_l
176 : dmrr_r=temp_r*damr_r
177 : c--------
178 : dmlu_l=n_x/am+temp_l*damu_l
179 : dmlu_r=temp_l*damu_r
180 : dmru_l=temp_r*damu_l
181 : dmru_r=n_x/am+temp_r*damu_r
182 : c--------
183 : dmlv_l=n_y/am+temp_l*damv_l
184 : dmlv_r=temp_l*damv_r
185 : dmrv_l=temp_r*damv_l
186 : dmrv_r=n_y/am+temp_r*damv_r
187 : c--------
188 : dmlw_l=n_z/am+temp_l*damw_l
189 : dmlw_r=temp_l*damw_r
190 : dmrw_l=temp_r*damw_l
191 : dmrw_r=n_z/am+temp_r*damw_r
192 : c--------
193 : dmlp_l=temp_l*damp_l
194 : dmlp_r=temp_l*damp_r
195 : dmrp_l=temp_r*damp_l
196 : dmrp_r=temp_r*damp_r
197 : c---------------------------------------------------------------
197 : c Computing the calligraphic P+ and P- with their derivatives
198 : c see (21a) & (21b) on p.368.
199 : c---------------------------------------------------------------
200 : if(ml .ge. 1.0d0) then
201 : Pplus = 1.0d0
202 : else
203 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
204 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
205 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
206 : else
207 : Pplus = 0.0d0
208 : endif
209 : endif
210 : c---------------------------------------------------------------
211 : if(mr .ge. 1.0d0) then
212 : Pmin = 0.0d0
213 : else
214 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
215 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
216 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
217 : else
218 : Pmin = 1.0d0
219 : endif
220 : endif
221 : c---------------------------------------------------------------
222 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
223 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
224 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
225 : c--------------
226 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
227 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
228 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
229 : c---------------------------------------------------------------
230 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
231 : dPpr_l=brac_l*dmlr_l
232 : dPpr_r=brac_l*dmlr_r
233 : c------------
234 : dPpu_l=brac_l*dmlu_l
235 : dPpu_r=brac_l*dmlu_r
236 : c------------
237 : dPpv_l=brac_l*dmlv_l
238 : dPpv_r=brac_l*dmlv_r
239 : c------------
240 : dPpw_l=brac_l*dmlw_l
241 : dPpw_r=brac_l*dmlw_r
242 : c------------
243 : dPpp_l=brac_l*dmlp_l
244 : dPpp_r=brac_l*dmlp_r
245 : c------------
246 : else
247 : dPpr_l=0.0d0
248 : dPpr_r=0.0d0
249 : c-----------
250 : dPpu_l=0.0d0
251 : dPpu_r=0.0d0
252 : c-----------
253 : dPpv_l=0.0d0
254 : dPpv_r=0.0d0
255 : c-----------
256 : dPpw_l=0.0d0
257 : dPpw_r=0.0d0
258 : c-----------
259 : dPpp_l=0.0d0
260 : dPpp_r=0.0d0
261 : c-----------
262 : endif
263 : c---------------------------------------------------------------
264 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
265 : dPmr_l=brac_r*dmrr_l
266 : dPmr_r=brac_r*dmrr_r
267 : c------------
268 : dPmu_l=brac_r*dmru_l
269 : dPmu_r=brac_r*dmru_r
270 : c------------
271 : dPmv_l=brac_r*dmrv_l
272 : dPmv_r=brac_r*dmrv_r
273 : c------------
274 : dPmw_l=brac_r*dmrw_l
275 : dPmw_r=brac_r*dmrw_r
276 : c------------
277 : dPmp_l=brac_r*dmrp_l
278 : dPmp_r=brac_r*dmrp_r
279 : c------------
280 : else
281 : dPmr_l=0.0d0
282 : dPmr_r=0.0d0
283 : c-----------
284 : dPmu_l=0.0d0
285 : dPmu_r=0.0d0
286 : c-----------
287 : dPmv_l=0.0d0
288 : dPmv_r=0.0d0
289 : c-----------
290 : dPmw_l=0.0d0
291 : dPmw_r=0.0d0
292 : c-----------
293 : dPmp_l=0.0d0
294 : dPmp_r=0.0d0
295 : c-----------
296 : endif
297 : c---------------------------------------------------------------------
298 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
299 : c---------------------------------------------------------------------
300 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
301 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
302 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
303 : dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
304 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
305 : c----------------------------
306 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
307 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
308 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
309 : dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
310 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
311 : c---------------------------------------------------------------------
312 : c computing JACOBIAN as a derivative of the numerical flux function
313 : c with respect to the primitive variables.
314 : c Notation: jl(i,j) --- is the derivative of the i-component of the
315 : c flux function with respect to the j-component of the
316 : c vector of primitive variables of the left state.
317 : c jr(i,j) --- is the derivative of the i-component of the
318 : c flux function with respect to the j-component of the
319 : c vector of primitive variables of the right state.
320 : c---------------------------------------------------------------------
321 : jl(1,1)=0.0D0
322 : jl(1,2)=0.0D0
323 : jl(1,3)=0.0D0
324 : jl(1,4)=0.0D0
325 : jl(1,5)=0.0D0
326 : c------------------------------------
327 : jr(1,1)=0.0D0
328 : jr(1,2)=0.0D0
329 : jr(1,3)=0.0D0
330 : jr(1,4)=0.0D0
331 : jr(1,5)=0.0D0
332 : c------------------------------------
333 : c------------------------------------
334 : jl(2,1)=dpir_l*n_x
335 : c------------------------
336 : jl(2,2)=dpiu_l*n_x
337 : c------------------------
338 : jl(2,3)=dpiv_l*n_x
339 : c------------------------
340 : jl(2,4)=dpiw_l*n_x
341 : c------------------------
342 : jl(2,5)=dpip_l*n_x
343 : c-------------------------------------------------
344 : c-------------------------------------------------
345 : jl(3,1)=dpir_l*n_y
346 : c------------------------
347 : jl(3,2)=dpiu_l*n_y
348 : c------------------------
349 : jl(3,3)=dpiv_l*n_y
350 : c------------------------
351 : jl(3,4)=dpiw_l*n_y
352 : c------------------------
353 : jl(3,5)=dpip_l*n_y
354 : c-------------------------------------------------
355 : c-------------------------------------------------
356 : jl(4,1)=dpir_l*n_z
357 : c------------------------
358 : jl(4,2)=dpiu_l*n_z
359 : c------------------------
360 : jl(4,3)=dpiv_l*n_z
361 : c------------------------
362 : jl(4,4)=dpiw_l*n_z
363 : c------------------------
364 : jl(4,5)=dpip_l*n_z
365 : c-------------------------------------------------------
366 : c derivatives with respect to the right
367 : c set of primitive variables
368 : c-------------------------------------------------------
369 : jr(2,1)=dpir_r*n_x
370 : c------------------------
371 : jr(2,2)=dpiu_r*n_x
372 : c------------------------
373 : jr(2,3)=dpiv_r*n_x
374 : c------------------------
375 : jr(2,4)=dpiw_r*n_x
376 : c------------------------
377 : jr(2,5)=dpip_r*n_x
378 : c-------------------------------------------------------
379 : jr(3,1)=dpir_r*n_y
380 : c------------------------
381 : jr(3,2)=dpiu_r*n_y
382 : c------------------------
383 : jr(3,3)=dpiv_r*n_y
384 : c------------------------
385 : jr(3,4)=dpiw_r*n_y
386 : c------------------------
387 : jr(3,5)=dpip_r*n_y
388 : c--------------------------------------------------------
389 : jr(4,1)=dpir_r*n_z
390 : c------------------------
391 : jr(4,2)=dpiu_r*n_z
392 : c------------------------
393 : jr(4,3)=dpiv_r*n_z
394 : c------------------------
395 : jr(4,4)=dpiw_r*n_z
396 : c------------------------
397 : jr(4,5)=dpip_r*n_z
398 : c-------------------------------------------------------
399 : c-------------------------------------------------------
400 : jl(5,1)=0.0D0
401 : jl(5,2)=0.0D0
402 : jl(5,3)=0.0D0
403 : jl(5,4)=0.0D0
404 : jl(5,5)=0.0D0
405 : c-------------------------------------------------
406 : c-------------------------------------------------
407 : jr(5,1)=0.0D0
408 : jr(5,2)=0.0D0
409 : jr(5,3)=0.0D0
410 : jr(5,4)=0.0D0
411 : jr(5,5)=0.0D0
412 : c-------------------------------------------------------------
413 : c-----------------------------------------------
414 : c Taking into account the dependancy of variables
415 : c-----------------------------------------------
416 : c11=t1_y*t2_z - t1_z*t2_y
417 : c12=n_y*t2_z - n_z*t2_y
418 : c13=n_y*t1_z - n_z*t1_y
419 : c-------------------------------------
420 : c21=t1_x*t2_z - t1_z*t2_x
421 : c22=n_x*t2_z - n_z*t2_x
422 : c23=n_x*t1_z - n_z*t1_x
423 : c-------------------------------------
424 : c31=t1_x*t2_y - t1_y*t2_x
425 : c32=n_x*t2_y - n_y*t2_x
426 : c33=n_x*t1_y - n_y*t1_x
427 : det=n_x*c11 - n_y*c21 + n_z*c31
428 : c----------------------------------------------------------------------
429 : ZC11=-NVECT(1)*C11-TVECT1(1)*C12+TVECT2(1)*C13
430 : ZC12=-NVECT(2)*C11-TVECT1(2)*C12+TVECT2(2)*C13
431 : ZC13=-NVECT(3)*C11-TVECT1(3)*C12+TVECT2(3)*C13
432 : C---------------------------------
433 : ZC21=NVECT(1)*C21+TVECT1(1)*C22-TVECT2(1)*C23
434 : ZC22=NVECT(2)*C21+TVECT1(2)*C22-TVECT2(2)*C23
435 : ZC23=NVECT(3)*C21+TVECT1(3)*C22-TVECT2(3)*C23
436 : C---------------------------------
437 : ZC31=-NVECT(1)*C31-TVECT1(1)*C32+TVECT2(1)*C33
438 : ZC32=-NVECT(2)*C31-TVECT1(2)*C32+TVECT2(2)*C33
439 : ZC33=-NVECT(3)*C31-TVECT1(3)*C32+TVECT2(3)*C33
440 : c---------------------------------------------------------------------
441 : do 11 i=1,5
442 : jl(i,1)=jl(i,1)+jr(i,1)
443 : jl(i,2)=jl(i,2)+jr(i,2)*zc11/det+jr(i,3)*zc21/det+
444 : & jr(i,4)*zc31/det
445 : jl(i,3)=jl(i,3)+jr(i,2)*zc12/det+jr(i,3)*zc22/det+
446 : & jr(i,4)*zc32/det
447 : jl(i,4)=jl(i,4)+jr(i,2)*zc13/det+jr(i,3)*zc23/det+
448 : & jr(i,4)*zc33/det
449 : jl(i,5)=jl(i,5)+jr(i,5)
450 : 11 continue
451 : c----------------------------------------------------------------------
452 : return
453 : end
454 :
455 :
456 :
457 :
458 :
459 :
460 :
461 :
462 :
463 :
464 :
465 :
466 :
467 :
468 :
469 :
470 :
471 :
472 :
473 :
© Cast3M 2003 - All rights reserved.
Disclaimer