-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathstepND_mhd_gr.f90
153 lines (146 loc) · 4.6 KB
/
stepND_mhd_gr.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
!!--------------------------------------------------------------------
!! Computes one timestep
!! Change this subroutine to change the timestepping algorithm
!!--------------------------------------------------------------------
subroutine step
use dimen_mhd
use debug
use loguns
use bound
use eos
use hterms
use options
use part
use part_in
use rates
use timestep
use setup_params
use xsph
use cons2prim, only:conservative2primitive,getv_from_pmom
!
!--define local variables
!
implicit none
integer :: i
real :: hdt
!
!--allow for tracing flow
!
if (trace) write(iprint,*) ' Entering subroutine step'
!
!--set initial quantities
!
hdt = 0.5*dt
do i=1,npart
xin(:,i) = x(:,i)
pmomin(:,i) = pmom(:,i)
Bevolin(:,i) = Bevol(:,i)
rhoin(i) = rho(i)
hhin(i) = hh(i)
enin(i) = en(i)
alphain(:,i) = alpha(:,i)
enddo
!
!--Mid-point Predictor step
!
do i=1,npart
if (itype(i).EQ.itypebnd .or. itype(i).EQ.itypebnd2) then ! fixed particles
pmom(:,i) = pmomin(:,i)
rho(i) = rhoin(i)
Bevol(:,i) = Bevolin(:,i)
if (iener.NE.0) en(i) = enin(i)
hh(i) = hhin(i)
x(:,i) = xin(:,i) + hdt*(vel(1:ndim,i) + xsphfac*xsphterm(1:ndim,i))
alpha(:,i) = alphain(:,i)
else
pmom(:,i) = pmomin(:,i) + hdt*force(:,i)
if (iener.NE.0) en(i) = enin(i) + hdt*dendt(i)
!--use updated v to change position
!print*,i,'predictor'
call getv_from_pmom(xin(:,i),pmom(:,i),vel(:,i),en(i),pr(i),rho(i),dens(i),uu(i))
x(:,i) = xin(:,i) + hdt*(vel(1:ndim,i) + xsphfac*xsphterm(1:ndim,i))
!--do an 'iteration' of this to be on the safe side
!call getv_from_pmom(x(:,i),pmom(:,i),vel(:,i),en(i),pr(i),rhoin(i),dens(i),uu(i))
x(:,i) = xin(:,i) + hdt*(vel(1:ndim,i) + xsphfac*xsphterm(1:ndim,i))
if (imhd.NE.0) Bevol(:,i) = Bevolin(:,i) + hdt*dBevoldt(:,i)
if (ihvar.EQ.1) then
! hh(i) = hfact(pmass(i)/rho(i))**dndim ! my version
hh(i) = hhin(i)*(rhoin(i)/rho(i))**dndim ! Joe's
elseif (ihvar.EQ.2 .OR. ihvar.EQ.3) then
hh(i) = hhin(i) + hdt*dhdt(i)
endif
if (icty.GE.1) rho(i) = rhoin(i) + hdt*drhodt(i)
if (ANY(iavlim.NE.0)) alpha(:,i) = min(alphain(:,i) + hdt*daldt(:,i),1.0)
endif
enddo
!
!--calculate all derivatives
!
call derivs
!print*,'corrector'
!
!--Mid-point Corrector step
!
do i=1,npart
if (itype(i).EQ.itypebnd .or. itype(i).eq.itypebnd2) then
pmom(:,i) = pmomin(:,i)
rho(i) = rhoin(i)
Bevol(:,i) = Bevolin(:,i)
if (iener.NE.0) en(i) = enin(i)
x(:,i) = xin(:,i) + hdt*(vel(1:ndim,i) + xsphfac*xsphterm(1:ndim,i))
alpha(:,i) = alphain(:,i)
hh(i) = hhin(i)
else
pmom(:,i) = pmomin(:,i) + dt*force(:,i)
if (iener.NE.0) en(i) = enin(i) + dt*dendt(i)
!--use updated v to change position
call getv_from_pmom(xin(:,i),pmom(:,i),vel(:,i),en(i),pr(i),rho(i),dens(i),uu(i))
x(:,i) = xin(:,i) + dt*(vel(1:ndim,i) + xsphfac*xsphterm(1:ndim,i))
!--do an 'iteration' of this (with updated x) to be on the safe side
call getv_from_pmom(x(:,i),pmom(:,i),vel(:,i),en(i),pr(i),rho(i),dens(i),uu(i))
x(:,i) = xin(:,i) + dt*(vel(1:ndim,i) + xsphfac*xsphterm(1:ndim,i))
if (ihvar.EQ.2) then
hh(i) = hhin(i) + dt*dhdt(i)
if (hh(i).LE.0.) then
write(iprint,*) 'step: hh -ve ',i,hh(i)
call quit
endif
endif
if (icty.GE.1) then
rho(i) = rhoin(i) + dt*drhodt(i)
if (rho(i).LE.0.) then
write(iprint,*) 'step: rho -ve ',i,rho(i)
call quit
endif
endif
if (ANY(iavlim.NE.0)) alpha(:,i) = min(alphain(:,i) + dt*daldt(:,i),1.0)
if (imhd.NE.0) Bevol(:,i) = Bevolin(:,i) + dt*dBevoldt(:,i)
endif
enddo
!
!--update density using a full summation every so often
!
if (MOD(nsteps,ndirect).EQ.0) then
call iterate_density
do i=1,npart
rhoin(i) = rho(i)
enddo
endif
!
!--if doing divergence correction then do correction to magnetic field
!
! if (idivBzero.NE.0) call divBcorrect
!
if (ANY(ibound.NE.0)) call boundary ! inflow/outflow/periodic boundary conditions
!
!--set new timestep from courant/forces condition
!
dt = min(C_force*dtforce,C_cour*dtcourant)
!
!--need to calculate velocity from particle momentum for next predictor step
! (no need to call this from output if called here)
!
call conservative2primitive
if (trace) write (iprint,*) ' Exiting subroutine step'
return
end subroutine step