-
Notifications
You must be signed in to change notification settings - Fork 37
Expand file tree
/
Copy pathtemporal.f90
More file actions
182 lines (143 loc) · 5.81 KB
/
Copy pathtemporal.f90
File metadata and controls
182 lines (143 loc) · 5.81 KB
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
module temporal
!DO NOT FIX THESE WARNINGS - THEY ARE FOR UNUSED VARIABLES THAT MAY BE LEVERAGED
!IN LATER RELEASES
!/home/zettergm/zettergmdata/GEMINI/temporal/temporal.f90:65:0: warning: unused
!parameter ‘ns’ [-Wunused-parameter]
! pure subroutine
!dt_calc(tcfl,ns,Ts,vs1,vs2,vs3,B1,B2,B3,dx1i,dx2i,dx3i,potsolve,cour1,cour2,cour3,dt) ^
!/home/zettergm/zettergmdata/GEMINI/temporal/temporal.f90:65:0: warning: unused parameter ‘b1’ [-Wunused-parameter]
!/home/zettergm/zettergmdata/GEMINI/temporal/temporal.f90:65:0: warning: unused parameter ‘b2’ [-Wunused-parameter]
!/home/zettergm/zettergmdata/GEMINI/temporal/temporal.f90:65:0: warning: unused parameter ‘b3’ [-Wunused-parameter]
!/home/zettergm/zettergmdata/GEMINI/temporal/temporal.f90:65:0: warning: unused parameter ‘potsolve’ [-Wunused-parameter]
use mpi, only: MPI_COMM_WORLD,MPI_STATUS_IGNORE
use phys_consts, only: kB,mu0,ms,lsp,pi, wp, debug
use mpimod, only: mpi_realprec, tagdt, lid, myid
use mesh, only: curvmesh
implicit none
private
public :: dt_comm
contains
subroutine dt_comm(t,tout,tglowout,flagglow,tcfl,ns,Ts,vs1,vs2,vs3,B1,B2,B3,x,potsolve,dt)
real(wp), intent(in) :: t,tout,tglowout,tcfl
real(wp), dimension(-1:,-1:,-1:,:), intent(in) :: ns,Ts,vs1,vs2,vs3
real(wp), dimension(-1:,-1:,-1:), intent(in) :: B1,B2,B3
type(curvmesh), intent(in) :: x
integer, intent(in) :: flagglow,potsolve
real(wp), intent(out) :: dt
real(wp), dimension(lsp) :: cour1,cour2,cour3
integer :: iid,isp, ierr
real(wp) :: dttmp
call dt_calc(tcfl,ns,Ts,vs1,vs2,vs3,B1,B2,B3,x%dl1i,x%dl2i,x%dl3i,potsolve,cour1,cour2,cour3,dt)
if (myid/=0) then
call mpi_send(dt,1,mpi_realprec,0,tagdt,MPI_COMM_WORLD,ierr) !send what I think dt should be
call mpi_recv(dt,1,mpi_realprec,0,tagdt,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr) !receive roots decision
else
!FIGURE OUT GLOBAL DT REQUIRED FOR STABILITY
do iid=1,lid-1
call mpi_recv(dttmp,1,mpi_realprec,iid,tagdt,MPI_COMM_WORLD,MPI_STATUS_IGNORE,ierr)
if (dttmp < dt) dt=dttmp
end do
!CHECK WHETHER WE'D OVERSTEP OUR TARGET OUTPUT TIME
!GLOW OUTPUT HAS PRIORITY SINCE IT WILL OUTPUT MORE OFTEN
if ((flagglow/=0).and.(t+dt>tglowout)) then
dt=tglowout-t
print *, 'GLOW is throttling dt...'
end if
if (t+dt>tout) then
dt=tout-t
if (debug) print *, 'Slowing down for an output...'
end if
!DON'T ALLOW ZERO DT
dt = max(dt, 1e-6_wp)
!SEND GLOBAL DT TO ALL WORKERS
do iid=1,lid-1
call mpi_send(dt,1,mpi_realprec,iid,tagdt,MPI_COMM_WORLD,ierr)
end do
if (debug) then
print *, 'dt figured to be: ',dt
print *, 'x1,x2,x3 courant numbers (root process only!): '
do isp=1,lsp
print '(a,f4.2,a,f4.2,a,f4.2)', ' ',cour1(isp),', ',cour2(isp),', ',cour3(isp)
!! these are roots courant numbers
end do
print *, 'Min and max density: ',minval(pack(ns(:,:,:,7),.true.)),maxval(pack(ns(:,:,:,7),.true.))
endif
end if
end subroutine dt_comm
pure subroutine dt_calc(tcfl,ns,Ts,vs1,vs2,vs3,B1,B2,B3,dx1i,dx2i,dx3i,potsolve,cour1,cour2,cour3,dt)
!------------------------------------------------------------
!-------COMPUTE TIME STEP SUCH THAT STABILITY CONDITION IS
!-------SATISFIED. NOTE THAT THE DIFFERENTIALS ARE ASSUMED
!-------TO HAVE UNITS OF DISTANCE, SO THEY MUST IMPLICITLY
!-------INCLUDE THE METRIC FACTORS.
!------------------------------------------------------------
real(wp), intent(in) :: tcfl
real(wp), dimension(-1:,-1:,-1:,:), intent(in) :: ns,Ts,vs1,vs2,vs3
real(wp), dimension(-1:,-1:,-1:), intent(in) :: B1,B2,B3
real(wp), dimension(:,:,:), intent(in) :: dx1i
real(wp), dimension(:,:,:), intent(in) :: dx2i
real(wp), dimension(:,:,:), intent(in) :: dx3i
integer, intent(in) :: potsolve
real(wp), dimension(lsp), intent(out) :: cour1,cour2,cour3
real(wp), intent(out) :: dt
real(wp), dimension(lsp) :: gridrate1,gridrate2,gridrate3
real(wp) :: vsnd
real(wp) :: rhom,Bmag,vA
integer :: lx1,lx2,lx3,ix1,ix2,ix3,isp
lx1=size(Ts,1)-4
lx2=size(Ts,2)-4
lx3=size(Ts,3)-4
gridrate1=0._wp
gridrate2=0._wp
gridrate3=0._wp
!EVALUATE TIME STEP AGAINST LOCAL SOUND SPEED AND ADVECTION
do isp=1,lsp
do ix3=1,lx3
do ix2=1,lx2
do ix1=1,lx1
if (isp<lsp) then
vsnd=sqrt(kB*Ts(ix1,ix2,ix3,isp)/ms(isp)+5._wp/3._wp*kB*Ts(ix1,ix2,ix3,lsp)/ms(isp))
else
vsnd=0._wp
end if
gridrate1(isp)=max((vsnd+abs(vs1(ix1,ix2,ix3,isp)))/dx1i(ix1,ix2,ix3),gridrate1(isp))
gridrate2(isp)=max(abs(vs2(ix1,ix2,ix3,isp))/dx2i(ix1,ix2,ix3),gridrate2(isp))
gridrate3(isp)=max(abs(vs3(ix1,ix2,ix3,isp))/dx3i(ix1,ix2,ix3),gridrate3(isp))
end do
end do
end do
end do
! !CHECK GRIDRATE MAX AGAINST LOCAL ALFVEN SPEED (IF THIS SIMULATION IS INDUCTIVE)
! !NOTE THAT THIS SHOULD REALLY INCLUDE MAGNETOSONIC (FAST) MODES IN GRIDRATE DETERMINATION
! !AS OF 2/10/2016 THIS CODE IS NOT USED AT ALL, BUT IS KEPT FOR FUTURE DEVELOPMENT
! if (potsolve == 2) then
! do ix3=1,lx3
! do ix2=1,lx2
! do ix1=1,lx1
! rhom=0._wp
! do isp=1,lsp
! rhom=rhom+ms(isp)*ns(ix1,ix2,ix3,isp)
! end do
! Bmag=sqrt(B1(ix1,ix2,ix3)**2+B2(ix1,ix2,ix3)**2+B3(ix1,ix2,ix3)**2)
!
! vA=Bmag/sqrt(rhom*mu0)
!
! do isp=1,lsp
! gridrate1(isp)=max(vA/dx1i(ix1),gridrate1(isp))
! gridrate2(isp)=max(vA/dx2i(ix2),gridrate2(isp))
! gridrate3(isp)=max(vA/dx3i(ix3),gridrate3(isp))
! end do
! end do
! end do
! end do
! end if
!ENFORCE A MINIMUM VALUE FOR THE GRIDRATE (WHICH IS CFL/DT, GRID POINTS PER SECOND)
gridrate1=max(gridrate1, 1e-10_wp)
gridrate2=max(gridrate2, 1e-10_wp)
gridrate3=max(gridrate3, 1e-10_wp)
dt=tcfl*min(minval(1._wp/gridrate1),minval(1._wp/gridrate2),minval(1._wp/gridrate3))
cour1=dt*gridrate1
cour2=dt*gridrate2
cour3=dt*gridrate3
end subroutine dt_calc
end module temporal