-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsetup_toystar1D_pmass.f90
189 lines (177 loc) · 5.47 KB
/
setup_toystar1D_pmass.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
!------------------------------------------------------------------------------!
! 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: !
!------------------------------------------------------------------------------!
!!----------------------------------------------------------------
!! Set up a toy star (with oscillations) in 1D
!! Starts with particles on a uniform density grid
!!----------------------------------------------------------------
SUBROUTINE setup
!
!--include relevant global variables
!
USE dimen_mhd
USE debug
USE loguns
USE bound
USE eos
USE options
USE part
USE setup_params
!
!--define local variables
!
IMPLICIT NONE
INTEGER :: i,j
INTEGER, PARAMETER :: itsmax = 100
REAL, PARAMETER :: tol = 1.e-7
INTEGER :: imax,its,ipart,norder
REAL :: massp,totmass,radstar,sigma,denszero
REAL :: ampl,wk,xlambda,dxmax,denom,gam1,dgam1
REAL :: H,C,A ! toy star parameters
REAL :: Gn
LOGICAL :: oscills
CHARACTER(LEN=40) :: tstarfile
!
!--allow for tracing flow
!
IF (trace) WRITE(iprint,*) ' Entering subroutine setup (toystar)'
IF (ndim.GT.1) STOP 'Toy star not implemented for ndim > 1'
!
!--initially set up a uniform density grid
!
ibound = 0 ! no boundaries
nbpts = 0 ! no fixed particles
gam1 = gamma - 1.
dgam1 = 1./gam1
!
!--toy star parameters
!
H = 1.0 ! choose H=1 so central density is unity
C = 1.0 ! C > 0
A = 0.05/SQRT(2.) !
sigma = 1./SQRT(2.)
oscills = .true. ! with oscillations ?
norder = 3 ! mode of oscillation
!
!--read parameters from file
!
tstarfile = rootname(1:LEN_TRIM(rootname))//'.tstar'
OPEN(UNIT=ireadf,ERR=11,FILE=tstarfile,STATUS='old',FORM='formatted')
READ(ireadf,*,ERR=12) H,C,A
READ(ireadf,*,ERR=12) sigma
READ(ireadf,*,ERR=12) norder
CLOSE(unit=ireadf)
oscills = .true.
GOTO 13
11 CONTINUE
WRITE(iprint,*) tstarfile,' not found, using default options '
GOTO 13
12 CONTINUE
WRITE(iprint,*) ' error reading ',tstarfile
13 CONTINUE
IF (norder.LT.0) oscills = .false.
!
!--set dependent parameters
!
totmass = 4./3.*SQRT(H**3/C)
radstar = SQRT(H/C)
xmin(1) = -radstar ! set position of boundaries
xmax(1) = radstar
imax = INT((xmax(1)-xmin(1))/psep)
massp = psep
gam1 = gamma - 1.
iexternal_force = 1
!
!--allocate memory here
!
CALL alloc(imax)
WRITE(iprint,*) ' Toy star : total mass = ',totmass
WRITE(iprint,*) ' toy star radius = ',radstar
WRITE(iprint,*) ' particle mass = ',massp
WRITE(iprint,*) ' A = ',A,' C = ',C,' H = ',H
IF (oscills) WRITE(iprint,*) ' -> with oscills, order = ',norder
WRITE(iprint,*)
!
!--setup uniform density distribution
!
DO i=1,imax
x(1,i) = xmin(1) + (i-1)*psep + 0.5*psep
dens(i) = (H - C*x(1,i)**2)**dgam1
IF (oscills) THEN
vel(1,i) = A*Gn(x(1,i),norder)
ELSE
vel(1,i) = A*x(1,i)
ENDIF
IF (ndimV.GT.1) vel(2:3,i) = 0.
pmass(i) = psep*dens(i)
uu(i) = polyk*dens(i)**(gam1)/gam1
IF (imhd.GE.1) THEN
Bfield(1,i) = 0.0 !SQRT(1.5)
Bfield(2,i) = sigma*dens(i)
Bfield(3,i) = 0.0
ELSE
Bfield(:,i) = 0.
ENDIF
! PRINT*,x(1,i),dens(i),pmass(i),hh(i)
ENDDO
npart = imax
WRITE(iprint,*) ' finished toy star setup, npart = ',npart
ntotal = npart
RETURN
END
!
!--function to evaluate the Gegenbauer polynomial of index n given x
! (the parameter lambda = 3/2 in this case)
!
FUNCTION Gn(x,n)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: x
INTEGER :: i
REAL :: Gn,Gnminus1,Gnminus2
REAL :: fnorm
fnorm = 2.*(n+1)*(n+2)/REAL(2.*n + 3.)
! PRINT*,' fnorm = ',fnorm
!
!--specify first two Gegenbauer polynomials
!
Gnminus2 = 1.
Gnminus1 = 3.*x
SELECT CASE (n)
CASE (0)
Gn = Gnminus2
CASE (1)
Gn = Gnminus1
CASE (2:)
DO i=2,n
Gn = ((2*i+1)*x*Gnminus1 - (i+1)*Gnminus2)/REAL(i)
Gnminus2 = Gnminus1
Gnminus1 = Gn
ENDDO
END SELECT
Gn = Gn/fnorm
END FUNCTION Gn
!
! use this routine to modify the dump upon code restart
!
subroutine modify_dump()
implicit none
end subroutine modify_dump