-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathndspmhd.f90
386 lines (373 loc) · 17.9 KB
/
ndspmhd.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
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
!------------------------------------------------------------------------------!
! NDSPMHD: A Smoothed Particle (Magneto)Hydrodynamics code for (astrophysical) !
! fluid dynamics simulations in 1, 2 and 3 spatial dimensions. !
! !
! (c) 2002-2015 Daniel Price !
! !
! http://users.monash.edu.au/~dprice/ndspmhd !
! daniel.price@monash.edu -or- dprice@cantab.net (forwards to current address) !
! !
! NDSPMHD comes with ABSOLUTELY NO WARRANTY. !
! This is free software; and you are welcome to redistribute !
! it under the terms of the GNU General Public License !
! (see LICENSE file for details) and the provision that !
! this notice remains intact. If you modify this file, please !
! note section 2a) of the GPLv2 states that: !
! !
! a) You must cause the modified files to carry prominent notices !
! stating that you changed the files and the date of any change. !
! !
! ChangeLog: see below !
!------------------------------------------------------------------------------!
!!---------------------------------------------------------------------!!
!! _ _ _ !!
!! _ __ __| |___ _ __ _ __ ___ | |__ __| | !!
!! | '_ \ / _` / __| '_ \| '_ ` _ \| '_ \ / _` | !!
!! | | | | (_| \__ \ |_) | | | | | | | | | (_| | !!
!! |_| |_|\__,_|___/ .__/|_| |_| |_|_| |_|\__,_| !!
!! |_| !!
!! _ _ _ _ _ _ _ _ _ _ _ _ _ !!
!! / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ !!
!! ( B | y ) ( D | a | n | i | e | l ) ( P | r | i | c | e ) !!
!! \_/ \_/ \_/ \_/ \_/ \_/ \_/ \_/ \_/ \_/ \_/ \_/ \_/ !!
!! !!
!!---------------------------------------------------------------------!!
!! An N-D SPH code to handle compressible gas dynamics with MHD !!
!! !!
!! Written in Fortran 90 !!
!! !!
!! By Daniel Price, Institute of Astronomy, Cambridge, UK, 2002-2004 !!
!! University of Exeter, UK 2004-2008 !!
!! Monash University, Melbourne, Australia 2008- !!
!! !!
!! Email: daniel.price@monash.edu !!
!! !!
!! This version is designed to be as modular (and thus as adaptable) !!
!! as possible, as a testbed for SPH algorithms !!
!! !!
!! Specific features include: !!
!! !!
!! * options to choose continuity equation or density by !!
!! direct summation !!
!! !!
!! * options to choose different energy eqn implementations !!
!! !!
!! * different artificial viscosity prescriptions !!
!! !!
!! * choice of whether to use average h, average gradient !!
!! of kernel or Springel/Hernquist (2001) type correction !!
!! terms for varying h !!
!! !!
!! * Morris and Monaghan (1997) artificial viscosity switch !!
!! (turns off artificial viscosity away from shocks) !!
!! !!
!! * ghost particle boundaries (reflective/periodic) !!
!! !!
!! Although this code is original, I learnt my SPH from !!
!! other SPH codes written by Joe Monaghan and Matthew Bate, and so !!
!! some parts bear similarities to these codes. !!
!! !!
!!---------------------------------------------------------------------!!
program ndspmhd
use debug
use loguns
use versn
implicit none
integer, parameter :: maxruns = 200
integer :: i,iprev,irun, nruns
character(len=120), dimension(maxruns) :: runname
runname(1) = 'a'
!
!--version number
!
version = 'NDSPMHD-3D-SR-v5-4'
! * special relativity
! * a LOT of other changes over the years
! version = 'NDSPMHD-3D-v5-3'
! * multiple runnames off command line (does them in order)
! * dissipation switches for resistivity and conductivity
! * anticlumping term implemented as a modified kernel gradient
! * gradh terms calculated for anticlumping kernel (**removed)
! * particle info explicitly passed to density, density_partial
! * external_forces subroutine cleaned up + potential calculation
! version = 'NDSPMHD-3D-v5-2_18_05_2004' (saved 2:50pm)
! *** this version used to obtain swave and hydro shocks for thesis 10/5/04 ***
! *** this version used to obtain 1D MHD shock tube results for thesis 13/5/04
! *** also used for Bx peak test in thesis/paper III 17/5/04
! * 1D turbulence setup revised/set_powerspec
! * fixed particle boundaries work in >1D, initialise cleaned up
! * hyperbolic/parabolic divergence cleaning
! * div B = 0 by projection method in 2D - poisson eq. by direct sum in 2D,3D
! * accretion disc setups
! * step (predictor-corrector) could be bad for disks - look out!
! * lots of crap to do with MHD AV/switches in rates
! * lots of things for GR code (grutils, conservative2primitive_gr etc)
! * tested on shear flows
! * preliminary riemann solver routine
! * compiler warnings fixed in lots of subroutines (mostly unused variables)
! -> link renamed so doesn't conflict with internal function
! * kernel plotting utility, more kernels added
! * quite a few changes to supersphplot
! * major revamp of rates - split into separate subroutines (faster)
! version = 'NDSPMHD-3D-v5-1'
! * smoothing length iteration on single particles works
! * conservative2primitive and primitive2conservative
! * eos rehashed
! * setup is on primitive variables
! * psep sent to set_uniform_cartesian
! * external forces subroutine and option instead of itoystar
! * main loop in subroutine evolve
! * potential energy from external forces
! * symplectic integrator (not yet for MHD, not as good as leapfrog)
! * ihvar = 3 predictor step only in h update
! * setup_unifsph
! version = 'NDSPMHD-3D-v5-0'
! *** saved as working 3D version ***
! * equations use general alternative formulation
! * compiles in 3D
! * set_ghosts totally rewritten -> works in up to 3D
! * 3D neighbour finding fixed - density calculated OK in 3D
! * fixed particle boundaries in > 1D (itype)
! * handles zero pressure
! * fix for negative thermal energies
! * my symmetrisation of vsig - but not for MHD yet
! * can have different boundary options in different dimensions
! * set_uniform cartesian can be called multiple times
! * shock setups in 3D
! version = 'NDSPMHD-v4-0' ! make sure there are no .'s in version name
! *** versioning now done with CVS ***
! use 'make tag' to tag a working copy of the code in CVS
! version = 'SPMHD-ND-v3.6_22_10_2003'
! * bug fix in thermal conductivity in thermal energy equation
! * bug fix in source term for artificial stress in total energy equation
! * something goes wrong with mshk2 when hfact = 1.3 (art. stress?)
! * use inflow/outflow boundaries in 1D (boundaryND_1D).
! * bug fix in this for outflow boundaries (linkND)
! * another bug fix in inflow/outflow (copies forces etc when reshuffling)
! * setkern not called with ikernel (set inside setkern itself)
! * new header
! * uses dhdt in step/rates
! * uniform setups work in 1D
! * setup_wave2D, setup_wave_x_ND, initialise calls ghosts properly
! *** this is the version used for 1D results in revised papers ***
! version = 'SPMHD-ND-v3.5_10_10_2003'
! * timestep criterion agrees with paper I (and version SPMHD-f90-v2.5)
! * found the bug in iener=3!!! (rates)
! * fixed offset bug in set_uniform_cartesian
! * bug with -ve densities in unused ghosts
! * changes to sametime/multirun
! version = 'SPMHD-ND-v3.4_07_08_2003'
! * interpolate_kernel function - setkern,density,rates,density_partial
! * minor changes to step
! * no explicit references to idim - bugs with this fixed (nlistdim)
! * smoothing in blastwave, random particle setups
! * smooth.f90 -> SPH smoothing of initial conditions
! * eos can do either individual elements or whole arrays
! * experimenting with namelist input
! * gravity by direct summation
! * multirun compiles with SPH source files
! * set_uniform_spherical
! * setkern cleaned up
! * polytrope setup
! * BUG IN LINK FOR NO BOUNDARIES
! * beginnings of maketree
! * trace and idebug not itrace
! version = 'SPMHD-ND-v3.3_22_07_2003'
! * density iteration - only iterates particles which are not converged
! (density_partial.f90, iterate_density.f90) something wrong however
! * toy star: setup using different particle masses - force bug fixed
! * initialise: no longer reads file runname, opens .dat,.ev earlier
! * hconst removed, different particle masses didn't work
! * toy star solution in supersphplot
! * toy star results obtained using this code
! * grad h terms bug fixed in ND
! * setup for 2D MHD blast wave
! * log plots in supersphplot
! version = 'SPMHD-ND-v3.2-25-06-03'
! * rendering for supersphplot
! * rates cleaned up and tried to split into separate subroutines
! * bug fix in MHD dissipation terms, also in boundaryND for periodic bc's
! * initialise can read rootname off command line or from file
! * optimization of density/rates - divisions done only once
! * catches vsig error
! * write header in two parts
! * bug in density iteration fixed (+ Newton Raphson)
! * toy star setup 1-x^2
! * new Makefiles - make .tar.gz's, saves versions, compiles in 1 or 2D
! * processplot.pl to script movie generation
! * supersphplot split into subroutines, new Makefile
! * NB: there is *something* wrong with the iner=3 results for 1D shocks
! version = 'SPMHD-ND-v3.1-16-05-03'
! * uniform density setups in ND (setup_unifdis)
! * 2D shock setup
! * 1D shock setup for ND code - bug fix in get_neighbour_list for 1D
! * ghost particle setup rehashed in ND + bug fixed (dx > 0 only now)
! * supersphplot changes
! * some changes to linkND
! * rates doesn't calculate interactions between ghost particles
! * orszag-tang vortex seems to work OK
! * bug fix in Joe's anticlumping term (2D) - uses 1/hfact, define hconst
! version = 'SPMHD-ND-v3.0-29-04-03'
! * ND (not fully tested)
! * setup for Orszag-Tang vortex in 2D
! * ghost particles setup is simpler - must call set_ghost_particles
! before the call to link (link finds ghosts rather than ghosts creating
! new link list cells)
! * allocate reallocates arrays if particle # > array size
! * ghost particles work in 2D, need to sort out corners in 3D
! * bug in prvaniso with anticlumping term
! * need to fix periodic BC's, 1D problems
! version = 'SPMHD-f90-v2.5-lots-of-options'
! * isothermal MHD
! version = 'SPMHD-f90-v2.4-lots-of-options'
! * allocation of arrays is done in allocate.f90
! * link lists can handle empty cells
! * bug in link fixed (ibound=1)
! * bug in density/rates when only one particle in a link list cell fixed
! * inflow/outflow boundaries bug fixed (Srojeen test + MHD shocks)
! * pressure on fixed particles near edges fixed
! ** this version used to obtain 1D results with inflow/outflow
! version = 'SPMHD-f90-v2.3-lots-of-options'
! * inflow/outflow boundaries work in 1D
! * wrinsph writes an infile.
! * batch running - domulti.pl, multirun.f90, sametime.f90
! ** this is the version used to obtain all the
! results in the 1D case
! version = 'SPMHD-f90-v2.2-lots-of-options'
! * time step condition uses correct signal velocity (previously too slow)
! * energy equation includes grad h terms and Joe's anti-clumping term
! * initialise: ghost particles set to agree initially (t=0 graphs look better)
! * equation_of_state no longer calculates alfven speed (unnecessary)
! * limits/plots improved in supersphplot
! * bug in periodic boundaries fixed (causing spurious currents)
! * fast/slow wave setups properly done
! * inflow/outflow boundary conditions (boundary1D) - not tested
! version = 'SPMHD-f90-v2.1-lots-of-options'
! * MHD works!! Uses AV in direction of vunit for 1.5D problems
! * Joes anticlumping term fixed - now works as it should
! * all setups in f90
! version = 'SPMHD-f90-v2.0-lots-of-options'
! * fortran 90 (uses modules, allocatable arrays)
! version = 'SPMHD-v1.6-lots-of-options'
! * follows Joe's notes - viscosity term in induction equation
! * bug fix in ghost1D - reflects all velocities, not just vx.
! version = 'SPMHD-v1.5-lots-of-options'
! * setup for advection test(ok) and ryu/jones shock tubes (still crap)
! * spits out current density
! * evolves h in evolution equations (if ihvar = 2)
! * Joes clumping correction only in x-direction (line between particles)
! * fixed bug in periodic boundary conditions (link1D/step1D)
! version = 'SPMHD-v1.4-lots-of-options'
! * reverted to old h update (blast wave didn't work)
! * can evolve using either B/rho or B
! (setup is now always just the basic field B)
! * bug check in read_infile
! * neighbour array size checked in get_neigh...(set in COMMONS/linklist)
! * now uses correct signal speed!
! * anisotropic forces with vectors
! version = 'SPMHD-v1.3-lots-of-options'
! * ensures h proportional to 1/rho
! * bug in iav = 2 fixed
! * uses COMMONS/dimen_mhd separate to SUPERSPH1D
! * alternative induction equations (incl. v*div B term)
! * more mag force options (6=Morris' hybrid force)
! * alternative kernel is used only when idivBzero=2
! * iener=3 seems to work with MHD
! version = 'SPMHD-v1.2-lots-of-options'
! * kernel table changed to try the thomas/couchman anti-clumping kernel
! (stores grkern instead of grkern/(r/h) ).
! * changed magnetic field variable to be of any dimension e.g. Brho(3,idim)
! * fixed bug in adjust_bound and in sum_density.
! * Velocity variable also changed to be of arbitrary dimension.
! * fixed bug in read_infile (didn't read hfact)
! * however there is still a bug in iav = 2
! version = 'SPMHD-v1.1-lots-of-options'
! * polytropic equation of state
! * prints out loads of MHD parameters (use evsupersph.f)
! * setup for Brio/Wu shock tubes
! * MHD for iener = 3 in step
! * Joe's clumping correction term (doesn't fix things though)
! version = 'SPMHD-v1.0-lots-of-options'
! * MHD
! version = 'SPH-v1.1-lots-of-options'
! * individual common blocks
trace = .false. ! set tracing flow (prints entry into subroutine)
! trace = .true.
idebug = 'none'
! idebug = 'fixed'
! idebug = 'density'
! idebug = 'neighb'
! idebug = 'link'
! idebug = 'divB'
! itemp = 40000 ! debug one particular particle
call logun ! set logical unit numbers to use for input/output
!
!--get runname(s) off command line
!
iprev = 1
i = 1
do while (runname(iprev)(1:1).ne.' ' .and. i.le.maxruns)
call get_command_argument(i,runname(i))
!!print*,i,runname(i)
iprev = i
i = i + 1
if (i.gt.maxruns .and. runname(iprev)(1:1).ne.' ') then
print*,'WARNING: number of runs >= array size: setting nruns = ',maxruns
endif
enddo
if (i.gt.maxruns .and. runname(maxruns)(1:1).ne.' ') then
nruns = maxruns
else
nruns = iprev - 1
endif
print*,'number of runs = ',nruns
!
!--If nothing on command exit and print usage
!
if (runname(1)(1:1).eq.' ') then
nruns = 1
10 write(6,*) 'Enter name of run:'
read(*,*,err=10) runname(1)
! stop 'Usage: spmhd runname '
endif
!
!--for each runname, perform a simulation
!
do irun = 1,nruns
rootname = runname(irun)
print*,'run = ',irun,' runname = ',rootname
call initialise ! read files and parameters and setup particles
!
!--now call the main timestepping loop
!
call evolve
!
!--close all open files and exit
!
if (iprint.ne.6) then
close(unit=iprint)
endif
close(unit=ievfile)
close(unit=idatfile)
enddo
stop
end program ndspmhd
!!---------------------------------------------------------------------
!! This subroutine performs a graceful exit if the program has crashed
!!---------------------------------------------------------------------
subroutine quit
use loguns
use timestep
implicit none
write(iprint,*) 'performing graceful exit...'
call output(time,-nsteps) ! dump particles before crashing
!
!--close all open files and exit
!
if (iprint.ne.6) then
close(unit=iprint)
endif
close(unit=ievfile)
close(unit=idatfile)
stop
end subroutine quit