-
Notifications
You must be signed in to change notification settings - Fork 35
/
pgas2d_dual.F
135 lines (134 loc) · 4.55 KB
/
pgas2d_dual.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#define FORTRAN
c=======================================================================
c////////////////////// SUBROUTINE PGAS2D_DUAL \\\\\\\\\\\\\\\\\\\\\\\
c
subroutine pgas2d_dual
& (dslice, eslice, geslice, pslice,
& uslice, vslice, wslice, eta1, eta2,
& idim, jdim, i1, i2, j1, j2, gamma, pmin, ierror)
c
c COMPUTES GAS PRESSURE ON A SLICE (DUAL ENERGY VERSION)
c
c written by: Greg Bryan
c date: March, 1994
c modified1:
c
c PURPOSE: Computes gas pressure needed for interpolations for PPM
c solver. EOS used is p=(gamma-1)*d*[E-(u**2+v**2+w**2)/2],
c (eqn 2.1) Only works for gamma law gas with gamma > 0.
c This particular version uses the dual energy formalism
c (see PPM gravity method paper) to determine the pressure
c from a combination of the total energy and the gas energy.
c
c INPUTS:
c dslice - density slice
c eslice - total specific energy slice
c eta1 - selection parameter for gas energy (typically ~0.001)
c eta2 - selection parameter for total energy (typically ~0.1)
c gamma - gamma law parameter
c geslice - specific gas energy slice
c i1,i2 - starting and ending addresses for dimension 1
c idim - declared leading dimension of slices
c j1,j2 - starting and ending addresses for dimension 2
c jdim - declared second dimension of slices
c pmin - minimum allowed pressure
c uslice - velocity-1 slice
c vslice - velocity-2 slice
c wslice - velocity-3 slice
c
c OUTPUT ARGUMENTS:
c pslice - pressure slice
c
c EXTERNALS:
c
c LOCALS:
c
c-----------------------------------------------------------------------
c
implicit NONE
#define FORTRAN
#include "fortran_types.h"
c
c-----------------------------------------------------------------------
c
c argument declarations
c
INTG_PREC i1, i2, idim, j1, j2, jdim
R_PREC eta1, eta2, gamma, pmin
R_PREC dslice(idim,jdim), eslice(idim,jdim), pslice(idim,jdim),
& uslice(idim,jdim), vslice(idim,jdim), wslice(idim,jdim),
& geslice(idim,jdim)
INTG_PREC ierror
c
c locals
c
INTG_PREC i, j, im1, ip1
R_PREC demax, ge1, ge2, ke
c
c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\///////////////////////////////
c=======================================================================
c
c compute the pressure
C print*, 'DEBUG_PRESSURE pgas2d_dual.F'
C print*, 'DEBUG_PRESSURE i1 i2 ',i1,i2
C print*, 'DEBUG_PRESSURE j1 j2 ',j1,j2
C print*, 'DEBUG_PRESSURE idim jdim ',idim,jdim
C print*, 'DEBUG_PRESSURE eta ',eta1,eta2
C print*, 'DEBUG_PRESSURE pmin ',pmin
do j=j1,j2
do i=i1,i2
c
c Compute the specific energy from the total energy
c
ke = 0.5_RKIND*(uslice(i,j)**2 + vslice(i,j)**2
& + wslice(i,j)**2)
ge1 = eslice(i,j) - ke
c
c Find the maximum nearby total energy (not specific)
c
im1 = max(i-1,i1)
ip1 = min(i+1,i2)
demax = max(dslice(i ,j)*eslice(i ,j),
& dslice(im1,j)*eslice(im1,j),
& dslice(ip1,j)*eslice(ip1,j))
c
c If the ratio of the gas energy to the max nearby total energy is
c is > eta2 then use the gas energy computed from the total energy
c to update the specific energy (total energy is better).
c
if (geslice(i,j) .le. 0._RKIND) write(6,*)
& 'pga:',ierror,i,j,geslice(i,j),demax
if (ge1*dslice(i,j)/demax .gt. eta2) geslice(i,j) = ge1
if (geslice(i,j) .le. 0._RKIND) write(6,*)
& 'pgb:',ierror,i,j,
& geslice(i,j),dslice(i,j),demax,eslice(i,j)
if (geslice(i,j) .le. 0._RKIND)
& ierror = ERROR_PGAS2D_DUAL_GE_LT_0
c
c If the ratio of the specific gas energy to specific total energy
c (in this cell) is < eta1 then use the gas energy to update
c the specific energy (gas energy is better).
c
if (ge1/eslice(i,j) .gt. eta1) then
ge2 = ge1
else
ge2 = geslice(i,j)
endif
c
c If pressure is below the minimum, set it to the minimum
c
ge2 = max(ge2, pmin/((gamma - 1._RKIND)*dslice(i,j)))
c
c Update the total energy
c
eslice(i,j) = eslice(i,j) - ge1 + ge2
c
c Compute the pressure with the total energy derived gas energy
c
pslice(i,j) = (gamma - 1._RKIND)*dslice(i,j)*ge2
c
enddo
enddo
c
return
end