Warning: This page is an HTML page, i.e. contains markup. Follow the download-link to obtain an unaltered copy of the file the content of which is displayed below.

1
C     ------------
2
CN    NAME: R I E M A N N _ V T 
3
C     ------------
4

        
5
CP    PURPOSE:
6
CP    THIS PROGRAM COMPUTES THE SOLUTION OF A 1D   
7
CP    RELATIVISTIC RIEMANN PROBLEM WITH ARBITRARY TANGENTIAL VELOCITIES
8
CP    (FOR CONSTANT-GAMMA IDEAL GASES) 
9
CP    WITH INITIAL DATA UL IF X<R0 AND UR IF X>R0   
10
CP    IN THE WHOLE SPATIAL DOMAIN [R0 - 0.5,R0 + 0.5]  
11
C
12

        
13
CC    COMMENTS:
14
CC    SEE PONS, MARTI AND MUELLER, JFM, 2000
15
CC
16
CC    WRITTEN BY:     Jose-Maria Marti
17
CC                    Departamento de Astronomia y Astrofisica 
18
CC                    Universidad de Valencia 
19
CC                    46100 Burjassot (Valencia), Spain
20
CC                    jose-maria.marti@uv.es
21
CC    AND
22
CC                    Ewald Mueller
23
CC                    Max-Planck-Institut fuer Astrophysik
24
CC                    Karl-Schwarzschild-Str. 1
25
CC                    85741 Garching, Germany
26
CC                    emueller@mpa-garching.mpg.de
27
C
28

        
29
      PROGRAM RIEMANN_VT
30

        
31
      IMPLICIT NONE
32

        
33
      INCLUDE 'npoints'
34

        
35
C     ------
36
C     COMMON BLOCKS
37
C     ------
38

        
39
      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, VELTL, WL, 
40
     &                 RHOR, PR, UR, HR, CSR, VELR, VELTR, WR 
41
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, VELTL, WL, 
42
     &                 RHOR, PR, UR, HR, CSR, VELR, VELTR, WR
43

        
44
      DOUBLE PRECISION RHOLS, ULS, HLS, CSLS, VELLS, VELTLS, VSHOCKL 
45
      COMMON /LS/      RHOLS, ULS, HLS, CSLS, VELLS, VELTLS, VSHOCKL
46

        
47
      DOUBLE PRECISION RHORS, URS, HRS, CSRS, VELRS, VELTRS, VSHOCKR 
48
      COMMON /RS/      RHORS, URS, HRS, CSRS, VELRS, VELTRS, VSHOCKR
49

        
50
      DOUBLE PRECISION GB
51
      COMMON /GB/      GB
52

        
53
      DOUBLE PRECISION QX(NG), QW(NG)
54
      COMMON /INT/     QX, QW
55
 
56
C     -----------
57
C     INTERNAL VARIABLES
58
C     -----------
59

        
60
      INTEGER          MN, N, I, ILOOP 
61
      PARAMETER        (MN = 400)
62

        
63
      DOUBLE PRECISION TOL, PMIN, PMAX, DVEL1, DVEL2, CHECK, R0
64
      PARAMETER        (R0=0.D0)
65
 
66
      DOUBLE PRECISION PS, VELS
67

        
68
      DOUBLE PRECISION RHOA(MN), PA(MN), VELA(MN), VELTA(MN),
69
     &                 UA(MN), HA(MN)
70
  
71
      DOUBLE PRECISION A
72
 
73
      DOUBLE PRECISION RAD(MN), X1, X2, X3, X4, X5, T, LAM, LAMS
74

        
75
      DOUBLE PRECISION KL, KR, CONST
76
 
77
C     ----------
78
C     INITIAL STATES 
79
C     ----------
80

        
81
      WRITE(*,*) ' ADIABATIC INDEX OF THE GAS: '
82
      READ (*,*)   GB
83
 
84
      WRITE(*,*) ' TIME FOR THE SOLUTION: '
85
      READ (*,*)   T
86
 
87
C     -----
88
C     LEFT STATE 
89
C     -----
90

        
91
      WRITE(*,*) ' ---- LEFT STATE ---- '
92
      WRITE(*,*) '      PRESSURE       : '
93
      READ (*,*) PL
94
      WRITE(*,*) '      DENSITY      : '
95
      READ (*,*) RHOL
96
      WRITE(*,*) '      FLOW VELOCITY-X: '
97
      READ (*,*) VELL
98
      WRITE(*,*) '      FLOW VELOCITY-Y: '
99
      READ (*,*) VELTL
100
 
101
C     ------
102
C     RIGHT STATE 
103
C     ------
104

        
105
      WRITE(*,*) ' ---- RIGHT STATE --- '
106
      WRITE(*,*) '      PRESSURE       : '
107
      READ (*,*) PR
108
      WRITE(*,*) '      DENSITY      : '
109
      READ (*,*) RHOR
110
      WRITE(*,*) '      FLOW VELOCITY-X: '
111
      READ (*,*) VELR
112
      WRITE(*,*) '      FLOW VELOCITY-Y: '
113
      READ (*,*) VELTR
114

        
115

        
116
C     ---------------------------------------
117
C     SPECIFIC INTERNAL ENERGY, SPECIFIC ENTHALPY, SOUND SPEED, 
118
C     ADIABATIC CONSTANT AND FLOW LORENTZ FACTORS IN THE INITIAL STATES 
119
C     ---------------------------------------
120

        
121
      UL = PL/(GB - 1.D0)/RHOL
122
      UR = PR/(GB - 1.D0)/RHOR
123

        
124
      HL = 1.D0 + GB*UL
125
      HR = 1.D0 + GB*UR
126

        
127
      CSL= DSQRT((GB - 1.D0)*(HL - 1.D0)/HL)
128
      CSR= DSQRT((GB - 1.D0)*(HR - 1.D0)/HR)
129

        
130
      KL = PL/RHOL**GB
131
      KR = PR/RHOR**GB
132

        
133
      WL = 1.D0/DSQRT(1.D0 - VELL*VELL - VELTL*VELTL)
134
      WR = 1.D0/DSQRT(1.D0 - VELR*VELR - VELTR*VELTR)
135

        
136
C     ------
137
C     COEFFICIENTS FOR NUMERICAL INTEGRATION IN RAREFACTIONS
138
C     ------
139

        
140
      CALL GAULEG(-1.D0,1.D0,QX,QW,NG)
141

        
142
C     -------- 
143
C     NUMBER OF POINTS 
144
C     --------
145

        
146
      N   = 400
147

        
148
C     ------------- 
149
C     TOLERANCE FOR THE SOLUTION 
150
C     -------------
151

        
152
      TOL = 0.D0
153

        
154
C
155

        
156
      ILOOP = 0
157

        
158
      IF ((PL.EQ.PR).AND.(VELL.EQ.VELR)) THEN
159
 
160
         PS      = PL
161
         VELS    = VELL
162

        
163
         VSHOCKL = VELL
164
         RHOLS   = (PS/KL)**(1.D0/GB)
165
         ULS     = PS/(GB - 1.D0)/RHOLS
166
         
167
         VSHOCKR = VELL
168
         RHORS   = (PS/KR)**(1.D0/GB)
169
         URS     = PS/(GB - 1.D0)/RHORS
170
                
171
      ELSE
172
        
173
        PMIN  = (PL + PR)/2.D0
174
        PMAX  = PMIN
175
 
176
5       ILOOP = ILOOP + 1
177

        
178
        PMIN  = 0.5D0*MAX(PMIN,0.D0)
179
        PMAX  = 2.D0*PMAX
180
 
181

        
182
        CALL GETDVEL2(PMIN, DVEL1)
183

        
184
        CALL GETDVEL2(PMAX, DVEL2)
185
        
186
        CHECK = DVEL1*DVEL2
187
        IF (CHECK.GT.0.D0) GOTO 5
188

        
189

        
190
C     --------------------------- 
191
C     PRESSURE AND FLOW VELOCITY IN THE INTERMEDIATE STATES 
192
C     ---------------------------
193

        
194
        CALL GETP2(PMIN, PMAX, TOL, PS)
195

        
196
        VELS = 0.5D0*(VELLS + VELRS)
197

        
198
        WRITE(*,*) 'VELS = ', VELS, 'PS = ', PS
199

        
200
      ENDIF
201
      
202
C     -------
203
C     SOLUTION ON THE NUMERICAL MESH 
204
C     -------
205

        
206
C     -----------
207
C     POSITIONS OF THE WAVES
208
C     -----------
209

        
210
      IF (PL.GE.PS) THEN
211

        
212
        CONST = HL*WL*VELTL
213

        
214
        CALL FLAMB(KL, CONST, PL, VELL, 'L', LAM )
215

        
216
        CALL FLAMB(KL, CONST, PS, VELS, 'L', LAMS)
217

        
218
        X1 = R0 + LAM *T
219
        X2 = R0 + LAMS*T 
220

        
221
      ELSE
222

        
223
        X1 = R0 + VSHOCKL*T 
224
        X2 = X1
225

        
226
      END IF
227
 
228
      X3 = R0 + VELS*T 
229
 
230
      IF (PR.GE.PS) THEN
231

        
232
        CONST = HR*WR*VELTR
233

        
234
        CALL FLAMB(KR, CONST, PS, VELS, 'R', LAMS)
235

        
236
        CALL FLAMB(KR, CONST, PR, VELR, 'R', LAM )
237

        
238
        X4 = R0 + LAMS*T 
239
        X5 = R0 + LAM *T
240
 
241
      ELSE
242

        
243
        X4 = R0 + VSHOCKR*T 
244
        X5 = X4
245

        
246
      END IF
247

        
248
C     ---------- 
249
C     SOLUTION ON THE MESH 
250
C     ----------
251

        
252
      DO 10 I=1,N
253

        
254
        RAD(I) = R0 + DFLOAT(I)/DFLOAT(N) - 0.5D0
255

        
256
 10     CONTINUE
257
 
258
      DO 120 I=1,N
259

        
260
        IF (RAD(I).LE.X1) THEN
261

        
262
          PA(I)    = PL
263
          RHOA(I)  = RHOL
264
          VELA(I)  = VELL
265
          VELTA(I) = VELTL  
266
          UA(I)    = UL
267
          HA(I)    = 1.D0 + GB*UA(I)
268

        
269
        ELSE IF (RAD(I).LE.X2) THEN
270

        
271
          A = (RAD(I) - R0)/T 
272

        
273
          CALL RAREF2(A, PS, RHOL, PL, UL, CSL, VELL, VELTL,  
274
     &                'L', RHOA(I), PA(I), UA(I), VELA(I), VELTA(I))
275

        
276
        ELSE IF (RAD(I).LE.X3) THEN
277

        
278
          PA(I)    = PS
279
          RHOA(I)  = RHOLS
280
          VELA(I)  = VELS
281
          UA(I)    = ULS
282
          HA(I)    = 1.D0 + GB*UA(I)
283
          CONST    = HL*WL*VELTL
284
          VELTA(I) = CONST*DSQRT((1.D0 - VELS*VELS)/
285
     &               (CONST**2 + HA(I)**2))
286

        
287
        ELSE IF (RAD(I).LE.X4) THEN
288

        
289
          PA(I)    = PS
290
          RHOA(I)  = RHORS
291
          VELA(I)  = VELS
292
          UA(I)    = URS
293
          HA(I)    = 1.D0 + GB*UA(I)
294
          CONST    = HR*WR*VELTR
295
          VELTA(I) = CONST*DSQRT((1.D0 - VELS*VELS)/
296
     &               (CONST**2 + HA(I)**2))
297

        
298
        ELSE IF (RAD(I).LE.X5) THEN
299

        
300
          A = (RAD(I) - R0)/T 
301
 
302
          CALL RAREF2(A, PS, RHOR, PR, UR, CSR, VELR, VELTR,  
303
     &                'R', RHOA(I), PA(I), UA(I), VELA(I), VELTA(I))
304

        
305
        ELSE
306

        
307
          PA(I)    = PR
308
          RHOA(I)  = RHOR
309
          VELA(I)  = VELR
310
          VELTA(I) = VELTR
311
          UA(I)    = UR
312
          HA(I)    = 1.D0 + GB*UA(I)
313

        
314
        END IF
315

        
316
120     CONTINUE
317

        
318
      OPEN (3,FILE='solution.dat',FORM='FORMATTED',STATUS='NEW')
319

        
320
      WRITE(3,150) N, T 
321
 150  FORMAT(I5,1X,F10.5)
322

        
323
      DO 60 I=1,N
324
        WRITE(3,200) RAD(I),PA(I),RHOA(I),VELA(I),VELTA(I),UA(I)
325
60      CONTINUE
326

        
327
200   FORMAT(5(E14.8,1X)) 
328
 
329
      CLOSE(3)
330

        
331
      STOP
332
      END
333
 
334
C     ---------- 
335
CN    NAME: G E T D V E L 2
336
C     ----------
337

        
338
CP    PURPOSE: 
339
CP    COMPUTE THE DIFFERENCE IN FLOW SPEED BETWEEN LEFT AND RIGHT INTERMEDIATE 
340
CP    STATES FOR GIVEN LEFT AND RIGHT STATES AND PRESSURE 
341
C 
342

        
343
CC    COMMENTS:
344
CC    NONE
345
 
346
      SUBROUTINE GETDVEL2( P, DVEL )
347

        
348
      IMPLICIT NONE
349

        
350
      INCLUDE 'npoints'
351

        
352
C     ------
353
C     ARGUMENTS
354
C     ------
355
 
356
      DOUBLE PRECISION P, DVEL
357

        
358
C     ------
359
C     COMMON BLOCKS
360
C     ------
361

        
362
      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, VELTL, WL, 
363
     &                 RHOR, PR, UR, HR, CSR, VELR, VELTR, WR 
364
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, VELTL, WL, 
365
     &                 RHOR, PR, UR, HR, CSR, VELR, VELTR, WR
366
 
367
      DOUBLE PRECISION RHOLS, ULS, HLS, CSLS, VELLS, VELTLS, VSHOCKL 
368
      COMMON /LS/      RHOLS, ULS, HLS, CSLS, VELLS, VELTLS, VSHOCKL
369

        
370
      DOUBLE PRECISION RHORS, URS, HRS, CSRS, VELRS, VELTRS, VSHOCKR 
371
      COMMON /RS/      RHORS, URS, HRS, CSRS, VELRS, VELTRS, VSHOCKR
372

        
373
      DOUBLE PRECISION GB
374
      COMMON /GB/      GB
375

        
376
      DOUBLE PRECISION QX(NG), QW(NG)
377
      COMMON /INT/     QX, QW
378

        
379
C     -----
380
C     LEFT WAVE 
381
C     -----
382

        
383
      CALL GETVEL2(P, RHOL, PL, UL,  HL,  CSL,  VELL,  VELTL,  WL, 'L',
384
     &                RHOLS,    ULS, HLS, CSLS, VELLS, VELTLS, VSHOCKL)
385

        
386
C     -----
387
C     RIGHT WAVE
388
C     -----
389

        
390
      CALL GETVEL2(P, RHOR, PR, UR,  HR,  CSR,  VELR,  VELTR,  WR, 'R',
391
     &                RHORS,    URS, HRS, CSRS, VELRS, VELTRS, VSHOCKR)
392

        
393
      DVEL = VELLS - VELRS
394

        
395
      RETURN
396
      END
397
C     ------- 
398
CN    NAME: G E T P 2
399
C     -------
400

        
401
CP    PURPOSE: 
402
CP    FIND THE PRESSURE IN THE INTERMEDIATE STATE OF A RIEMANN PROBLEM IN 
403
CP    RELATIVISTIC HYDRODYNAMICS 
404
C 
405

        
406
CC    COMMENTS: 
407
CC    THIS ROUTINE USES A COMBINATION OF INTERVAL BISECTION AND INVERSE 
408
CC    QUADRATIC INTERPOLATION TO FIND THE ROOT IN A SPECIFIED INTERVAL. 
409
CC    IT IS ASSUMED THAT DVEL(PMIN) AND DVEL(PMAX) HAVE OPPOSITE SIGNS WITHOUT 
410
CC    A CHECK. 
411
CC    ADAPTED FROM "COMPUTER METHODS FOR MATHEMATICAL COMPUTATION", 
412
CC    BY G. E. FORSYTHE, M. A. MALCOLM, AND C. B. MOLER, 
413
CC    PRENTICE-HALL, ENGLEWOOD CLIFFS N.J. 
414

        
415
      SUBROUTINE GETP2( PMIN, PMAX, TOL, PS )
416

        
417
      IMPLICIT NONE
418

        
419
C     ----- 
420
C     ARGUMENTS 
421
C     -----
422

        
423
      DOUBLEPRECISION PMIN, PMAX, TOL, PS
424

        
425
C     ------- 
426
C     COMMON BLOCKS 
427
C     -------
428

        
429
      DOUBLE PRECISION GB
430
      COMMON /GB/      GB
431

        
432
      DOUBLE PRECISION RHOL, PL, UL, HL, CSL, VELL, VELTL, WL, 
433
     &                 RHOR, PR, UR, HR, CSR, VELR, VELTR, WR 
434
      COMMON /STATES/  RHOL, PL, UL, HL, CSL, VELL, VELTL, WL, 
435
     &                 RHOR, PR, UR, HR, CSR, VELR, VELTR, WR
436

        
437
C     --------- 
438
C     INTERNAL VARIABLES 
439
C     ---------
440

        
441
      DOUBLEPRECISION A, B, C, D, E, EPS, FA, FB, FC, TOL1, 
442
     &                XM, P, Q, R, S
443

        
444
C     ------------- 
445
C     COMPUTE MACHINE PRECISION 
446
C     -------------
447

        
448
      EPS  = 1.D0 
449
10    EPS  = EPS/2.D0 
450
      TOL1 = 1.D0 + EPS 
451
      IF( TOL1 .GT. 1.D0 ) GO TO 10
452

        
453
C     ------- 
454
C     INITIALIZATION 
455
C     -------
456

        
457
      A  = PMIN 
458
      B  = PMAX 
459
      CALL GETDVEL2(A,FA) 
460
      CALL GETDVEL2(B,FB)
461

        
462
C     ----- 
463
C     BEGIN STEP 
464
C     -----
465

        
466
20    C  = A 
467
      FC = FA 
468
      D  = B - A 
469
      E  = D 
470
30    IF( DABS(FC) .GE. DABS(FB) )GO TO 40 
471
      A  = B 
472
      B  = C 
473
      C  = A 
474
      FA = FB 
475
      FB = FC 
476
      FC = FA
477

        
478
C     -------- 
479
C     CONVERGENCE TEST 
480
C     --------
481

        
482
40    TOL1 = 2.D0*EPS*DABS(B) + 0.5D0*TOL 
483
      XM   = 0.5D0*(C - B) 
484
      IF( DABS(XM) .LE. TOL1 ) GO TO 90 
485
      IF( FB .EQ. 0.D0 ) GO TO 90
486

        
487
C     ------------ 
488
C     IS BISECTION NECESSARY? 
489
C     ------------
490

        
491
      IF( DABS(E) .LT. TOL1 ) GO TO 70 
492
      IF( DABS(FA) .LE. DABS(FB) ) GO TO 70
493

        
494
C     ------------------ 
495
C     IS QUADRATIC INTERPOLATION POSSIBLE? 
496
C     ------------------
497

        
498
      IF( A .NE. C ) GO TO 50
499

        
500
C     ---------- 
501
C     LINEAR INTERPOLATION 
502
C     ----------
503

        
504
      S = FB/FA 
505
      P = 2.D0*XM*S 
506
      Q = 1.D0 - S 
507
      GO TO 60
508

        
509
C     ---------------- 
510
C     INVERSE QUADRATIC INTERPOLATION 
511
C     ----------------
512

        
513
50    Q = FA/FC 
514
      R = FB/FC 
515
      S = FB/FA 
516
      P = S*(2.D0*XM*Q*(Q - R) - (B - A)*(R - 1.D0)) 
517
      Q = (Q - 1.D0)*(R - 1.D0)*(S - 1.D0)
518

        
519
C     ------ 
520
C     ADJUST SIGNS 
521
C     ------
522

        
523
60    IF( P .GT. 0.D0 ) Q = -Q 
524
      P = DABS(P)
525

        
526
C     -------------- 
527
C     IS INTERPOLATION ACCEPTABLE? 
528
C     --------------
529

        
530
      IF( (2.D0*P) .GE. (3.D0*XM*Q-DABS(TOL1*Q)) ) GO TO 70 
531
      IF( P .GE. DABS(0.5D0*E*Q) ) GO TO 70 
532
      E = D 
533
      D = P/Q 
534
      GO TO 80
535

        
536
C     ----- 
537
C     BISECTION 
538
C     -----
539

        
540
70    D = XM 
541
      E = D
542

        
543
C     ------- 
544
C     COMPLETE STEP 
545
C     -------
546

        
547
80    A  = B 
548
      FA = FB 
549
      IF( DABS(D) .GT. TOL1 ) B = B+D 
550
      IF( DABS(D) .LE. TOL1 ) B = B+DSIGN(TOL1,XM) 
551
      CALL GETDVEL2(B,FB) 
552
      IF( (FB*(FC/DABS(FC))) .GT. 0.D0) GO TO 20 
553
      GO TO 30
554

        
555
C     -- 
556
C     DONE 
557
C     --
558

        
559
90    PS = B
560

        
561
      RETURN 
562
      END
563

        
564
C     ------
565
CN    NAME: F L A M B
566
C     ------
567

        
568
CP    PURPOSE:
569
CP    COMPUTE THE VALUE OF THE SELF-SIMILARITY VARIABLE INSIDE A RAREFACTION 
570
CP    CONNECTED TO A SPECIFIED LEFT / RIGHT STATE
571
C
572

        
573
CC    COMMENTS:
574
CC    NONE
575

        
576
      SUBROUTINE FLAMB(K, A, P, VEL, S, XI)
577

        
578
      IMPLICIT NONE
579

        
580
C     --------
581
C     ARGUMENTS
582
C     --------
583

        
584
      DOUBLE PRECISION K, A, P, VEL
585

        
586
      CHARACTER*1      S
587

        
588
      DOUBLE PRECISION XI
589

        
590
C     -------
591
C     COMMON BLOCKS
592
C     -------
593

        
594
      DOUBLE PRECISION G
595
      COMMON /GB/      G
596

        
597
C     --------------
598
C     INTERNAL VARIABLES
599
C     --------------
600
      
601
      DOUBLE PRECISION SIGN
602
      DOUBLE PRECISION RHO, H, CS2, VELT2, V2, BETA, DISC
603

        
604
      IF (S.EQ.'L') SIGN = -1.D0
605

        
606
      IF (S.EQ.'R') SIGN =  1.D0
607
  
608
      RHO   = (P/K)**(1.D0/G)
609
      CS2   = G*(G - 1.D0)*P/(G*P +(G - 1.D0)*RHO)
610
      H     = 1.D0/(1.D0 - CS2/(G - 1.D0))
611

        
612
      VELT2 = (1.D0 - VEL*VEL)*A*A/(H*H + A*A)
613
      V2    = VELT2 + VEL*VEL
614
     
615
      BETA  = (1.D0 - V2)*CS2/(1.D0 - CS2)
616
      DISC  = DSQRT(BETA*(1.D0 + BETA - VEL*VEL))
617
      
618
      XI = (VEL + SIGN*DISC)/(1.D0 + BETA)
619

        
620
      RETURN
621
      END
622

        
623
C     -------- 
624
CN    NAME: R A R E F 2
625
C     --------
626

        
627
CP    PURPOSE: 
628
CP    COMPUTE THE FLOW STATE IN A RAREFACTION FOR GIVEN PRE-WAVE STATE 
629
C
630
 
631
CC    COMMENTS: 
632
CC    THE VELOCITY IN THE RAREFACTION IS WRITTEN IN TERMS OF THE PRESCRIBED 
633
CC    LEFT / RIGHT STATE  AND PRESSURE ACCORDING TO EXPRESSIONS (3.25) AND 
634
CC    (3.26) OF REZZOLLA, ZANOTTI AND PONS, JFM, 2002.
635
CC    THE INTEGRAL IN THE VELOCITY EXPRESSION IS COMPUTED THROUGH A GAUSSIAN 
636
CC    QUADRATURE 
637

        
638
      SUBROUTINE RAREF2(XI, PS, RHOA, PA, UA, CSA, VELA, VELTA, S, 
639
     &                      RHO,  P,  U,       VEL,  VELT)
640

        
641
      IMPLICIT NONE
642

        
643
      INCLUDE 'npoints'
644

        
645
C     ------
646
C     ARGUMENTS
647
C     ------
648

        
649
      DOUBLE PRECISION XI, PS, RHOA, PA, UA, CSA, VELA, VELTA
650

        
651
      CHARACTER*1      S
652

        
653
      DOUBLE PRECISION RHO, P, U, VEL, VELT
654

        
655
C     ------
656
C     COMMON BLOCKS
657
C     ------
658

        
659
      DOUBLE PRECISION GB
660
      COMMON /GB/      GB
661

        
662
      DOUBLE PRECISION QX(NG), QW(NG)
663
      COMMON /INT/     QX, QW
664

        
665
C     --------
666
C     INTERNAL VARIABLES
667
C     --------
668

        
669
      INTEGER          I
670

        
671
      DOUBLE PRECISION HA, WA, SIGN
672

        
673
      DOUBLE PRECISION CONST, K, XIO, XIP, PO, H
674

        
675
      DOUBLE PRECISION SUMW, DIFW, INTEGRAL, XX, RRHO, CCS2, HH, 
676
     &                 FUNR, A, FP, DFDP
677

        
678
      HA = 1.D0 + UA + PA/RHOA
679

        
680
      WA = 1.D0/DSQRT(1.D0 - VELA*VELA - VELTA*VELTA)
681

        
682
      CONST = HA*WA*VELTA
683

        
684
      K     = PA/RHOA**GB
685

        
686
      IF (S.EQ.'L') SIGN = -1.D0
687

        
688
      IF (S.EQ.'R') SIGN =  1.D0
689

        
690
      CALL FLAMB(K, CONST, PA, VELA, S, XIO)
691

        
692
      PO = PA
693
      P  = 0.95D0*PA
694

        
695
 20   CONTINUE
696

        
697
      SUMW = 0.5D0*(P + PA)
698
      DIFW = 0.5D0*(P - PA)
699
      INTEGRAL = 0.D0
700

        
701
      DO 10 I = 1, NG
702
         
703
         XX   = DIFW*QX(I) + SUMW
704
         RRHO = (XX/K)**(1.D0/GB)
705
         CCS2 = GB*(GB - 1.D0)*XX/(GB*XX + (GB - 1.0)*RRHO) 
706
         HH   = 1.D0/(1.D0 - CCS2/(GB - 1.D0)) 
707

        
708
         FUNR = DSQRT(HH*HH + CONST*CONST*(1.D0 - CCS2))/
709
     &         (HH*HH + CONST*CONST)/(RRHO*DSQRT(CCS2))
710
 
711
         INTEGRAL = INTEGRAL + DIFW*QW(I)*FUNR
712
        
713
 10      CONTINUE
714

        
715
      A    = SIGN*INTEGRAL + 0.5D0*DLOG((1.D0 + VELA)/(1.D0 - VELA))
716
      VEL  = DTANH(A)
717
 
718
      CALL FLAMB(K, CONST, P, VEL, S, XIP)
719

        
720
      FP   = XIP - XI
721

        
722
      DFDP = (XIP - XIO)/(P - PO)
723

        
724
      PO  = P
725
      XIO = XIP 
726

        
727
      P = P - FP/DFDP
728
      P = DMAX1(P,PS)
729

        
730
      IF (DABS(FP).GT.1.D-10) GOTO 20
731

        
732
      RHO = (P/K)**(1.D0/GB)
733

        
734
      U   = P/RHO/(GB -1.D0)
735
      
736
      H   = 1.D0 + U + P/RHO
737

        
738
      VELT = CONST*DSQRT((1.D0 - VEL*VEL)/(CONST*CONST + H*H))
739

        
740
      RETURN
741
      END
742

        
743
C     --------- 
744
CN    NAME: G E T V E L 2
745
C     ---------
746

        
747
CP    PURPOSE: 
748
CP    COMPUTE THE FLOW VELOCITY BEHIND A RAREFACTION OR SHOCK IN TERMS OF THE 
749
CP    POST-WAVE PRESSURE FOR A GIVEN STATE AHEAD THE WAVE IN A RELATIVISTIC 
750
CP    FLOW 
751
C 
752

        
753
CC    COMMENTS: 
754
CC    THIS ROUTINE CLOSELY FOLLOWS THE EXPRESSIONS IN PONS, MARTI AND MUELLER, 
755
CC    JFM, 2002
756

        
757
      SUBROUTINE GETVEL2(P, RHOA, PA, UA, HA, CSA, VELA, VELTA, WA, S,
758
     &                      RHO,      U,  H,  CS,  VEL,  VELT,  VSHOCK)
759

        
760
      IMPLICIT NONE
761

        
762
      INCLUDE 'npoints'
763

        
764
C     --------
765
C     ARGUMENTS
766
C     --------
767

        
768
      DOUBLE PRECISION P, RHOA, PA, UA, HA, CSA, VELA, VELTA, WA
769
      CHARACTER*1      S
770
      DOUBLE PRECISION RHO, U, H, CS, VEL, VELT, VSHOCK
771

        
772
C     -----
773
C     COMMON BLOCKS
774
C     -----
775

        
776
      DOUBLE PRECISION GB
777
      COMMON /GB/      GB
778

        
779
      DOUBLE PRECISION QX(NG), QW(NG)
780
      COMMON /INT/     QX, QW
781

        
782
C     --------
783
C     INTERNAL VARIABLES
784
C     --------
785

        
786
      INTEGER          I
787

        
788
      DOUBLE PRECISION A, B, C, SIGN
789
      DOUBLE PRECISION J, WSHOCK
790
      DOUBLE PRECISION K, CONST, SUMW, DIFW, INTEGRAL, XX, RRHO,
791
     &                 HH, CCS2, FUNR
792

        
793
C     ------------
794
C     LEFT OR RIGHT PROPAGATING WAVE
795
C     ------------
796

        
797
      IF (S.EQ.'L') SIGN = -1.D0
798

        
799
      IF (S.EQ.'R') SIGN =  1.D0
800

        
801
      IF (P.GE.PA) THEN
802

        
803
C       ---
804
C       SHOCK 
805
C       ---
806

        
807
        A  = 1.D0 - (GB - 1.D0)*(P - PA)/GB/P
808
        B  = 1.D0 - A
809
        C  = HA*(PA - P)/RHOA - HA**2
810
 
811
C       ---------------- 
812
C       CHECK FOR UNPHYSICAL ENTHALPIES 
813
C       ----------------
814

        
815
        IF (C.GT.(B**2/4.D0/A)) STOP 
816
     & 'GETVEL2: UNPHYSICAL SPECIFIC ENTHALPY IN INTERMEDIATE STATE'
817
 
818
C       ---------------------------------
819
C       SPECIFIC ENTHALPY AT THE LEFT OF THE CONTACT DISCONTINUITY 
820
C       (OBTAINED FROM THE EQUATION OF STATE AND THE TAUB ADIABAT)
821
C       ---------------------------------
822

        
823
        H = (-B + DSQRT(B**2 - 4.D0*A*C))/2.D0/A
824
 
825
C       -----------------------------------
826
C       DENSITY AT THE LEFT OF THE CONTACT DISCONTINUITY
827
C       (OBTAINED FROM SPECIFIC ENTHALPY AND THE EQUATION OF STATE)
828
C       -----------------------------------
829

        
830
        RHO = GB*P/(GB - 1.D0)/(H - 1.D0)
831

        
832
C       ----------------------------------
833
C       SPECIFIC INT. ENERGY AT THE LEFT OF THE CONTACT DISCONTINUITY 
834
C       (OBTAINED FROM THE EQUATION OF STATE)
835
C       ----------------------------------
836

        
837
        U = P/(GB - 1.D0)/RHO
838

        
839
C       -----------------------------
840
C       MASS FLUX ACROSS LEFT WAVE
841
C       (OBTAINED FROM THE RANKINE-HUGONIOT RELATIONS)
842
C       -----------------------------
843

        
844
        J = SIGN*DSQRT((P - PA)/(HA/RHOA - H/RHO))
845
 
846
C       ------------------------------
847
C       SHOCK VELOCITY 
848
C       (OBTAINED FROM THE DEFINITION OF MASS FLUX)
849
C       ------------------------------
850

        
851
        A      = J**2 + (RHOA*WA)**2
852
        B      = -2.D0*VELA*RHOA**2*WA**2
853
        C      = (RHOA*WA*VELA)**2 - J**2
854

        
855
        VSHOCK = (-B + SIGN*DSQRT(B*B - 4.D0*A*C))/(2.D0*A)
856
        WSHOCK = 1.D0/DSQRT(1.D0 - VSHOCK**2)
857
 
858
C       --------------------------
859
C       VELOCITY AT THE LEFT OF THE CONTACT DISCONTINUITY 
860
C       --------------------------
861

        
862
        A    = WSHOCK*(P - PA)/J + HA*WA*VELA
863
        B    = HA*WA + (P - PA)*(WSHOCK*VELA/J + 1.D0/RHOA/WA)
864
 
865
        VEL  = A/B
866
 
867
        A    = HA*WA*VELTA
868

        
869
        VELT = A*DSQRT((1.D0 - VEL*VEL)/(H*H + A*A))
870
 
871
C       -----------------------------
872
C       SOUND SPEED AT THE LEFT OF THE CONTACT DISCONTINUITY
873
C       -----------------------------
874

        
875
        CS = DSQRT(GB*P/RHO/H)
876
 
877
      ELSE
878

        
879
C       ------
880
C       RAREFACTION 
881
C       ------
882

        
883
        CONST = HA*WA*VELTA
884

        
885
        K     = PA/RHOA**GB
886

        
887
C       --------------------------
888
C       DENSITY AT THE LEFT SIDE OF THE CONTACT DISCONTINUITY
889
C       --------------------------
890

        
891
        RHO = (P/K)**(1.D0/GB)
892
 
893
C       ----------------------------------
894
C       SPECIFIC INT. ENERGY AT THE LEFT OF THE CONTACT DISCONTINUITY 
895
C       (OBTAINED FROM THE EQUATION OF STATE)
896
C       ----------------------------------
897

        
898
        U = P/(GB - 1.D0)/RHO
899
        H = 1.D0 + GB*U
900
 
901
C       ----------------------------------
902
C       SOUND SPEED AT THE LEFT OF THE CONTACT DISCONTINUITY 
903
C       ----------------------------------
904

        
905
        CS = DSQRT(GB*P/(RHO + GB*P/(GB - 1.D0)))
906
 
907
C       ----------------------------------
908
C       VELOCITY AT THE LEFT OF THE CONTACT DISCONTINUITY 
909
C       ----------------------------------
910

        
911
C       ------
912
C       INTEGRAL 
913
C       ------
914

        
915
        SUMW = 0.5D0*(P + PA)
916

        
917
        DIFW = 0.5D0*(P - PA)
918

        
919
        INTEGRAL = 0.D0
920

        
921
        DO 110 I = 1, NG
922

        
923
           XX  = DIFW*QX(I) + SUMW
924
           RRHO = (XX/K)**(1.D0/GB)
925
           CCS2 = GB*(GB - 1.D0)*XX/(GB*XX + (GB - 1.0)*RRHO) 
926
           HH   = 1.D0/(1.D0 - CCS2/(GB - 1.D0)) 
927

        
928
           FUNR = DSQRT(HH*HH + CONST*CONST*(1.D0 - CCS2))/
929
     &           (HH*HH + CONST*CONST)/(RRHO*DSQRT(CCS2))
930
 
931
           INTEGRAL = INTEGRAL + DIFW*QW(I)*FUNR
932
        
933
 110       CONTINUE
934

        
935
        A    = SIGN*INTEGRAL + 0.5D0*DLOG((1.D0 + VELA)/(1.D0 - VELA))
936
        VEL  = DTANH(A)
937
        VELT = CONST*DSQRT((1.D0 - VEL*VEL)/(CONST**2 + H**2))
938
        
939
      END IF
940
 
941
      RETURN
942
      END
943

        
944
C     ------
945
CN    NAME: G A U L E G 
946
C     ------
947

        
948
CP    PURPOSE:
949
CP    COMPUTE ABCISSAS AND WEIGHTS FOR GAUSS-LEGENDRE QUADRATURE INTEGRATION
950
C
951

        
952
CC    COMMENTS:
953
CC    ADAPTED FROM PRESS ET AL., "NUMERICAL RECIPES", CAMBRIDGE, 1988
954

        
955
      SUBROUTINE GAULEG(X1,X2,X,W,N)
956

        
957
      IMPLICIT NONE
958

        
959
C     --------
960
C     ARGUMENTS
961
C     --------
962

        
963
      INTEGER          N
964
      DOUBLE PRECISION X1,X2,X(N),W(N)
965

        
966
C     ---------
967
C     INTERNAL VARIABLES
968
C     ---------
969

        
970
      DOUBLE PRECISION EPS
971
      PARAMETER        (EPS=3.D-14)
972
      INTEGER          I, J, M
973
      DOUBLE PRECISION P1, P2, P3, PP, XL, XM, Z, Z1
974

        
975
      M  = (N + 1)/2
976
      XM = 0.5D0*(X2 + X1)
977
      XL = 0.5D0*(X2 - X1)
978

        
979
      DO 12 I = 1, M
980

        
981
        Z = COS(3.141592654D0*(I - .25D0)/(N + .5D0))
982
1       CONTINUE
983

        
984
        P1 = 1.D0
985
        P2 = 0.D0
986
          
987
        DO 11 J = 1, N
988

        
989
          P3 = P2
990
          P2 = P1
991
          P1 = ((2.D0*J - 1.D0)*Z*P2 - (J - 1.D0)*P3)/J
992
11        CONTINUE
993

        
994
        PP = N*(Z*P1 - P2)/(Z*Z - 1.D0)
995
        Z1 = Z
996
        Z  = Z1 - P1/PP
997
        
998
        IF (ABS(Z -Z1).GT.EPS) GOTO 1
999

        
1000
        X(I)     = XM - XL*Z
1001
        X(N+1-I) = XM + XL*Z
1002
        W(I)     = 2.D0*XL/((1.D0 - Z*Z)*PP*PP)
1003
        W(N+1-I) = W(I)
1004
12      CONTINUE
1005

        
1006
      RETURN
1007
      END