-
Notifications
You must be signed in to change notification settings - Fork 4
/
eos_initHelmholtz.F90
235 lines (198 loc) · 8.06 KB
/
eos_initHelmholtz.F90
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
!!****if* source/physics/Eos/EosMain/Helmholtz/SpeciesBased/eos_initHelmholtz
!!
!! NAME
!!
!! eos_initHelmholtz
!!
!! SYNOPSIS
!!
!! call eos_initHelmholtz()
!!
!! DESCRIPTION
!!
!! Initialize the Helmholtz EOS. The table data is read in on processor 0
!! and broadcast to all other processors at the start of FLASH execution.
!! This routine first checks for a binary copy of the table (helm_extended_table.bdat),
!! and then for the ASCII version (helm_extended_table.dat). If the binary table
!! file was not found, it is created by this routine for subsequent use.
!!
!! ARGUMENTS
!!
!! none
!!
!! PARAMETERS
!!
!! eos_coulombMult[Real, default 1.0] -- Coulomb correction multiplier. Set
!! to zero to ignore the Coloumb correction.
!! eos_tolerance[Real, 1.0e-8] -- Convergence tolerance for the Newton-Rhapson
!! iterations
!! eos_maxNewton[Integer, 50] -- Maximum number of Newton-Rhapson iterations
!! eos_forceConstantInput -- This switch forces the Eos implementation
!! to never modify EINT or PRES in MODE_DENS_EI
!! and MODE_DENS_PRES. If this is .false. (the
!! default), calls to Eos may slightly modify
!! these input variables in order to preserve
!! thermodynamic equilibrium.
!!
!! NOTES
!!
!! Helmholtz law Eos defines two mesh-based parameters GAMC_VAR and GAME_VAR in Flash.h
!!
!!***
#ifdef DEBUG_ALL
#define DEBUG_EOS
#endif
subroutine eos_initHelmholtz()
use Eos_data, ONLY : eos_type, eos_meshMe, &
eos_eintSwitch, eos_smallt
use eos_helmData
use Driver_interface, ONLY : Driver_abortFlash
use RuntimeParameters_interface, ONLY : RuntimeParameters_get
implicit none
! vector_eos.fh computes the vector length from nxb, nyb, nzb, so
! this information must be provided
#include "Flash.h"
#include "constants.h"
#include "Eos.h"
include 'Flash_mpi.h'
integer:: unitEos =2
integer :: i, j
real :: tstp, dstp
integer :: istat, ierr
!get the runtime parameters
call RuntimeParameters_get('smallt', eos_smallt)
call RuntimeParameters_get('eos_tolerance', eos_tol)
call RuntimeParameters_get('eos_maxNewton', eos_maxNewton)
call RuntimeParameters_get('eos_coulombMult', eos_coulombMult)
#ifdef DEBUG_EOS
print *, 'in eos_initHelmholtz'
#endif
call RuntimeParameters_get('eos_coulombAbort', eos_coulombAbort)
#ifdef DEBUG_EOS
print *, 'done with RuntimeParameters_get (eos_coulombAbort)'
#endif
#ifndef EINT_VAR
if (eos_eintSwitch > 0.0) then
call Driver_abortFlash("[Eos_init] eintSwitch is nonzero, but EINT_VAR not defined!")
end if
#endif
#ifdef USE_EOS_YE
write(*,*)"USE_EOS_YE should not be defined with Helmholtz/SpeciesBased EOS. Use Helmholtz/Ye instead!"
call Driver_abortFlash("[Eos_init] Use Helmholtz/Ye with USE_EOS_YE mode")
#endif
call RuntimeParameters_get("eos_forceConstantInput",eos_forceConstantInput)
if (eos_meshMe==MASTER_PE) then
open(unit=unitEos,file='helm_extended_table.bdat',status='old',iostat=istat)
close(unit=unitEos)
print*,'about to open file'
istat=1
if (istat.ne.0) then
write(*,*) '[Eos_init] Cannot open helm_extended_table.bdat!'
write(*,*) '[Eos_init] Trying old helm_extended_table.dat!'
open(unit=unitEos,file='helm_extended_table.dat',status='old',iostat=istat)
if (istat .ne. 0) then
write(*,*) '[Eos_init] ERROR: opening helm_extended_table.dat!'
call Driver_abortFlash("[Eos_init] ERROR: opening helm_extended_table.dat")
endif
!..read the helmholtz free energy table
do j=1,EOSJMAX
do i=1,EOSIMAX
read(unitEos,*) eos_f(i,j),eos_fd(i,j),eos_ft(i,j),&
eos_fdd(i,j),eos_ftt(i,j), &
eos_fdt(i,j), eos_fddt(i,j),eos_fdtt(i,j),eos_fddtt(i,j)
enddo
enddo
!..read the pressure derivative with density table
do j=1,EOSJMAX
do i=1,EOSIMAX
read(unitEos,*) eos_dpdf(i,j),eos_dpdfd(i,j),&
eos_dpdft(i,j),eos_dpdfdt(i,j)
enddo
enddo
!..read the electron chemical potential table
do j=1,EOSJMAX
do i=1,EOSIMAX
read(unitEos,*) eos_ef(i,j),eos_efd(i,j),&
eos_eft(i,j),eos_efdt(i,j)
enddo
enddo
!..read the number density table
do j=1,EOSJMAX
do i=1,EOSIMAX
read(unitEos,*) eos_xf(i,j),eos_xfd(i,j),&
eos_xft(i,j),eos_xfdt(i,j)
enddo
enddo
!..close up the data file and write a message
close(unitEos)
!..dump binary version of table for later use
istat = EOSIMAX*EOSJMAX
call eos_writeHfet(istat, &
eos_f,eos_fd,eos_ft,eos_fdd,&
eos_ftt,eos_fdt,eos_fddt,eos_fdtt,eos_fddtt, &
eos_dpdf,eos_dpdfd,eos_dpdft,eos_dpdfdt, &
eos_ef,eos_efd,eos_eft,eos_efdt, &
eos_xf,eos_xfd,eos_xft,eos_xfdt)
!..read binary version of table
else
istat = EOSIMAX*EOSJMAX
call eos_readHfet(istat, &
eos_f,eos_fd,eos_ft,eos_fdd,&
eos_ftt,eos_fdt,eos_fddt,eos_fdtt,eos_fddtt, &
eos_dpdf,eos_dpdfd,eos_dpdft,eos_dpdfdt, &
eos_ef,eos_efd,eos_eft,eos_efdt, &
eos_xf,eos_xfd,eos_xft,eos_xfdt)
endif
endif
!..broadcast to rest of processors
istat = EOSIMAX*EOSJMAX
call MPI_BCAST(eos_f, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_fd, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_ft, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_fdd, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_ftt, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_fdt, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_fddt, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_fdtt, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_fddtt, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_dpdf, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_dpdfd, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_dpdft, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_dpdfdt, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_ef, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_efd, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_eft, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_efdt, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_xf, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_xfd, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_xft, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
call MPI_BCAST(eos_xfdt, istat, FLASH_REAL, MASTER_PE, MPI_COMM_WORLD, ierr)
! Switched to using the extended table from cococubed.
eos_tlo = 3.0e0
tstp = (13.0e0 - eos_tlo)/float(EOSJMAX-1)
eos_tstpi = 1.0e0/tstp
eos_dlo = -12.0e0
dstp = (15.0e0 - eos_dlo)/float(EOSIMAX-1)
eos_dstpi = 1.0e0/dstp
do j=1,EOSJMAX
eos_t(j) = 10.0e0**(eos_tlo + (j-1)*tstp)
do i=1,EOSIMAX
eos_d(i) = 10.0e0**(eos_dlo + (i-1)*dstp)
enddo
enddo
!..store the temperature and density differences and their inverses
do j=1,EOSJMAX-1
eos_dt(j) = eos_t(j+1) - eos_t(j)
eos_dtSqr(j) = eos_dt(j)*eos_dt(j)
eos_dtInv(j) = 1.0e0/eos_dt(j)
eos_dtSqrInv(j) = 1.0e0/eos_dtSqr(j)
enddo
do i=1,EOSIMAX-1
eos_dd(i) = eos_d(i+1) - eos_d(i)
eos_ddSqr(i) = eos_dd(i)*eos_dd(i)
eos_ddInv(i) = 1.0e0/eos_dd(i)
eos_ddSqrInv(i) = 1.0e0/eos_ddSqr(i)
enddo
eos_type=EOS_HLM
return
end subroutine eos_initHelmholtz