-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsetup_rayleightaylor2D_abel.f90
137 lines (131 loc) · 3.43 KB
/
setup_rayleightaylor2D_abel.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
!----------------------------------------------------------------
! Set up a Rayleigh-Taylor instability test as
!----------------------------------------------------------------
subroutine setup
!
!--include relevant global variables
!
use dimen_mhd
use debug
use loguns
use bound
use options
use part
use setup_params, only:psep,pi
use eos, only:gamma
use mem_allocation, only:alloc
use uniform_distributions
use cons2prim, only:primitive2conservative
!
!--define local variables
!
implicit none
integer :: i,ipart
real :: massp,volume,totmass,massmedium
real :: denszero,densmedium,przero,psepmedium,pri
logical, parameter :: equalmass = .false.
!
!--allow for tracing flow
!
if (trace) write(iprint,*) ' entering subroutine setup(RT-Abel)'
if (ndim.ne.2) stop 'error: need ndim=2 for this problem'
!
!--set boundaries
!
ibound(1) = 2 ! periodic in x
ibound(2) = 1 ! fixed in y
nbpts = 0 ! use ghosts not fixed
xmin(1) = 0. ! set position of boundaries
xmax(1) = 0.5
xmin(2) = 0. !- 6.*psep ! set position of boundaries
xmax(2) = 1. !+ 6.*psep
if (iexternal_force.ne.9) stop 'need iexternal force = 9 for this problem'
!
!--set up the uniform density grid
!
denszero = 1.0
densmedium = 2.0
przero = denszero/gamma
psepmedium = psep*(denszero/densmedium)**(1./ndim)
!massmedium = massp*(densmedium/denszero)
volume = product(xmax(:)-xmin(:))
totmass = denszero*volume
!
!--unequal masses setup whole grid
!
call set_uniform_cartesian(1,psep,xmin,xmax,fill=.true.)
!call set_uniform_cartesian(2,psep,xmin,xmax,fill=.true.)
massp = totmass/float(npart)
massmedium = massp*(densmedium/denszero)
npart = ntotal
print*,'npart =',npart
!
!--now assign particle properties
!
do i=1,ntotal
vel(:,i) = 0.
if (x(2,i).gt.0.3 .and. x(2,i).lt.0.7) then
vel(2,i) = 0. !0.01*(1. + cos(8.*pi*(x(1,i)+0.25)))*&
! (1. + cos(2./0.4*pi*(x(2,i)-0.5)))/4.
else
vel(2,i) = 0.
endif
dens(i) = densmedium + (denszero-densmedium)/ &
(1. + exp(-2.*(x(2,i)-0.5)/0.05))
pmass(i) = massmedium + (massp-massmedium)/ &
(1. + exp(-2.*(x(2,i)-0.5)/0.05))
! if (x(2,i).gt.0.5) then
! dens(i) = densmedium
! pmass(i) = massmedium
! else
! dens(i) = denszero
! pmass(i) = massp
! endif
pri = przero - 0.5*dens(i)*(x(2,i) - 0.5)
uu(i) = pri/((gamma-1.)*dens(i))
!--set fixed particles in y
if (x(2,i).gt.0.9 .or. x(2,i).lt.0.1) then
nbpts = nbpts + 1
itype(i) = 1
endif
! Bfield(:,i) = 0.
! Bfield(1,i) = 0.5
enddo
!
!--get rho from a sum and then set u to give a
! smooth pressure
!
write(iprint,*) 'calling density to make smooth pressure jump...'
call primitive2conservative
do i=1,npart
pri = przero - 0.5*rho(i)*(x(2,i) - 0.5)
uu(i) = pri/((gamma-1.)*dens(i))
enddo
!
!--allow for tracing flow
!
if (trace) write(iprint,*) ' exiting subroutine setup'
return
end
subroutine modify_dump
use loguns, only:iprint
use part
use options, only:imhd
use timestep, only:time
use setup_params, only:pi
implicit none
integer :: i
!
!--now assign particle properties
!
write(iprint,*) 'modifying dump with velocities for Rayleigh-Taylor run'
do i=1,ntotal
vel(1,i) = 0.
if (abs(x(2,i)).lt.1./3.) then
vel(2,i) = 0.5*(1. + cos(4.*pi*x(1,i)))*(1. + cos(3.*pi*x(2,i)))/4.
else
vel(2,i) = 0.
endif
enddo
time = 0.
end subroutine modify_dump