u/erica/2DShockedClumpsSNR: problem.f90

File problem.f90, 10.2 KB (added by Erica Kaminski, 8 years ago)
Line 
1!#########################################################################
2!
3! Copyright (C) 2003-2012 Department of Physics and Astronomy,
4! University of Rochester,
5! Rochester, NY
6!
7! problem.f90 of module SingleClump is part of AstroBEAR.
8!
9! AstroBEAR is free software: you can redistribute it and/or modify
10! it under the terms of the GNU General Public License as published by
11! the Free Software Foundation, either version 3 of the License, or
12! (at your option) any later version.
13!
14! AstroBEAR is distributed in the hope that it will be useful,
15! but WITHOUT ANY WARRANTY; without even the implied warranty of
16! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17! GNU General Public License for more details.
18!
19! You should have received a copy of the GNU General Public License
20! along with AstroBEAR. If not, see <http://www.gnu.org/licenses/>.
21!
22!#########################################################################
23!> @dir MultiClumps
24!! @brief Contains files necessary for the Shape Tester problem
25
26!> @file problem.f90
27!! @brief Main file for module Problem
28
29!> @defgroup MultiClumps Shape Tester Module
30!! @brief Module for setting up orbiting particles
31!! @ingroup Modules
32
33!> MultiClump Module
34!! @ingroup MultiClumps
35MODULE Problem
36 USE DataDeclarations
37 USE ParticleDeclarations
38 USE Clumps
39 USE Winds
40 USE Ambients
41! USE VectorPerturbation
42! USE Histograms
43! USE Totals
44! USE Fields
45 IMPLICIT NONE
46 SAVE
47
48 PUBLIC ProblemModuleInit, ProblemGridInit, &
49 ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
50 REAL(KIND=qPREC) :: wind_thickness, wind_vel, wind_refinement_time
51 TYPE(WindDef), POINTER :: Wind
52CONTAINS
53
54 !> Initializes module variables
55 SUBROUTINE ProblemModuleInit()
56 INTEGER :: ClumpTracer, AltClumpTracer,i,j,dir,edge
57 Logical :: ClumpToTrace
58 TYPE(ClumpDef), POINTER :: Clump
59 !TYPE(HistogramDef), POINTER :: Histogram
60 !TYPE(TotalDef), POINTER :: KE,TE,BE,ENSTROPHY
61 !LOGICAL :: lCooling!, lMagClump, lMHDPerturbation
62 REAL(KIND=qPREC) :: density, temperature, velocity, pressure, smooth_distance, position(3), radius, B(3)
63 REAL(KIND=qPREC) :: temp_cgs, mach, cs_cgs
64 REAL(KIND=qPREC) :: theta=0
65 REAL(KIND=qPREC) :: phi=0
66 REAL(KIND=qPREC) :: B_tor=0
67 REAL(KIND=qPREC) :: B_pol=0
68 REAL(KIND=qPREC) :: B_phi=0
69 REAL(KIND=qPREC) :: B_theta=0
70 TYPE(AmbientDef), POINTER :: Ambient
71 REAL(KIND=qPREC) :: rhoOut, tempOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut, Pout
72 !INTEGER :: kmin, kmax
73 !REAL(KIND=qPREC) :: beta, field_amp, kscale
74 !INTEGER :: nwaves=0, nMHDwaves=0
75 NAMELIST /AmbientData/ rhoOut, tempOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut
76 NAMELIST /ProblemData/ ClumpToTrace!, lMagClump, lMHDPerturbation
77 NAMELIST /WindData/ density, temperature, velocity, B, wind_thickness, wind_refinement_time, dir, edge
78 NAMELIST /ClumpData/ density, temperature, radius, smooth_distance, theta, phi, B_tor, B_pol, B_theta, B_phi, position
79 !NAMELIST /PerturbData/ nwaves, nMHDwaves
80 !NAMELIST /MHDPerturbationData/ kmin, kmax, beta, field_amp, kscale
81
82 OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
83 READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
84 dir=1
85 edge=1
86 CALL AddTracer(ClumpTracer, 'ClumpTracer')
87 CALL AddTracer(AltClumpTracer, 'AltClumpTracer')
88 CALL CreateAmbient(Ambient)
89 READ(PROBLEM_DATA_HANDLE,NML=AmbientData)
90 cs_cgs = sqrt((gamma*boltzmann*tempout)/(xmu*amu))
91! print *, 'cs_cgs=', cs_cgs
92 Ambient%density=rhoOut/nscale
93 pOut=(rhoOut*Boltzmann*tempOut)
94 Ambient%pressure=pOut/pScale
95 Ambient%B(:)=(/BxOut, ByOut, BzOut/)
96 !Ambient%v(:)=(/vxOut, vyOut, vzOut/)
97 Ambient%velocity(:)=(/vxOut, vyOut, vzOut/)
98
99 READ(PROBLEM_DATA_HANDLE,NML=ClumpData)
100 CALL CreateClump(Clump)
101 Clump%density=density/nscale
102 Clump%temperature=temperature
103 Clump%position=position
104 IF (ClumpToTrace) THEN
105 Clump%iTracer=AltClumpTracer
106 ELSE
107 Clump%iTracer=ClumpTracer
108 END IF
109 Clump%thickness=smooth_distance
110 Clump%radius=radius
111 !Clump%Magnetized=lMagClump
112 Clump%B_toroidal=B_tor
113 Clump%B_poloidal=B_pol
114 Clump%theta=theta*pi/180d0
115 Clump%phi=phi*pi/180d0
116 Clump%B_theta=B_theta*pi/180d0
117 Clump%B_phi=B_phi*pi/180d0
118 CALL UpdateClump(Clump)
119
120 !IF (lMHDPerturbation) THEN
121 ! READ(PROBLEM_DATA_HANDLE, NML=MHDPerturbationData)
122 ! CALL CreateVectorPerturbationObj(Clump%MagneticPerturbation)
123 ! CALL CreateSolenoidalSpectra(Clump%MagneticPerturbation, kmin, kmax, beta, field_amp, kscale)
124 !END IF
125
126
127 READ(PROBLEM_DATA_HANDLE,NML=WindData)
128 CALL CreateWind(Wind)
129 !!!!! Calculate wind params based on R-H jump conditions
130 mach = velocity/cs_cgs
131 density=(4.0*rhoOut)/(1.0+(3.0/Mach**2))
132 pressure=((5.0*(Mach**2)-1.0)*pOut)/4
133 temperature=pressure/(density*boltzmann)
134 velocity=((3.0*((Mach**2)-1.0))/(4.0*Mach))*cs_cgs
135 IF(mpi_id==0)THEN
136 print *, 'Supersonic inflow boundary conditions:'
137 print *, 'mach = ', mach
138 print *, 'density(cm^-3) =', density
139 print *, 'pressure (Ba) =', pressure
140 print *, 'temperature (K) =', temperature
141 print *, 'velocity (cm/s) =', velocity
142 END IF
143 Wind%density=density/nscale
144 Wind%temperature=temperature
145 Wind%velocity=velocity/VelScale
146 Wind%B=B
147 wind_vel=velocity/VelScale
148 CALL AddTracer(Wind%iTracer, 'Wind_Tracer')
149 Wind%dir=dir
150 Wind%edge=edge
151
152! realspeed=(wind%velocity*VelScale)
153! IF(mpi_id==0)THEN
154! PRINT*, "realdensity = ", realden*1d3, "mg/cc"
155! PRINT*, "realtemp = ", Boltzmann*realtemp/1.602d-12, "ev"
156! PRINT*, "realradius = ", realradius, "mm"
157! PRINT*, "realbeta = ", realbeta
158! PRINT*, "realspeed = ", realspeed/1d5, "km/s"
159! PRINT*, "timescale = ", timescale
160 !PRINT*, "Reynolds# = ", wind%velocity*clump%radius/resistivity
161! END IF
162 !STOP
163
164 !CALL CreateTotal(KE)
165 !KE%Field%id=PKE_Field
166 !KE%Field%component=GASCOMP
167 !KE%Field%name='PKE'
168
169 !CALL CreateTotal(TE)
170 !TE%Field%id=PiE_Field
171 !TE%Field%component=GASCOMP
172 !TE%Field%name='PTE'
173
174 !CALL CreateTotal(BE)
175 !BE%Field%id=BE_Field
176 !BE%Field%component=GASCOMP
177 !BE%Field%name='BE'
178
179! CALL CreateHistogram(Histogram)
180! Histogram%Field%id=MixingRatio12_Field
181! Histogram%Field%component=GASCOMP
182! Histogram%minvalue=0d0
183! Histogram%maxvalue=1d0
184! Histogram%nbins=100
185! Histogram%scale=LINEARSCALE
186! Histogram%WeightField=BINBYMASS
187
188 !CALL CreateTotal(ENSTROPHY)
189 !ENSTROPHY%Field%iD=Enstrophy_Field
190 !ENSTROPHY%Field%Component=GASCOMP
191 !ENSTROPHY%Field%name='ENSTROPHY'
192 !Processing_mbc=max(Processing_mbc,1)
193 CLOSE(PROBLEM_DATA_HANDLE)
194
195 END SUBROUTINE ProblemModuleInit
196
197 !> Applies initial conditions
198 !! @param Info Info object
199 SUBROUTINE ProblemGridInit(Info)
200 TYPE(InfoDef) :: Info
201 END SUBROUTINE ProblemGridInit
202
203 !> Applies Boundary conditions
204 !! @param Info Info object
205 SUBROUTINE ProblemBeforeStep(Info)
206 TYPE(InfoDef) :: Info
207 INTEGER :: dir,edge, ierr
208 REAL(KIND=qPREC) :: cs_cgs, pOut, mach, pressure
209 REAL(KIND=qPREC) :: density, temperature, velocity, B(3), wind_thickness, wind_refinement_time
210 REAL(KIND=qPREC) :: rhoOut, tempOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut
211 NAMELIST /AmbientData/ rhoOut, tempOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut
212 NAMELIST /WindData/ density, temperature, velocity, B, wind_thickness, wind_refinement_time, dir, edge
213 !NAMELIST /WindData/ density, temperature, velocity, dir, edge
214
215 OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', IOSTAT=iErr, STATUS="OLD", POSITION="REWIND")
216 READ(PROBLEM_DATA_HANDLE,NML=AmbientData)
217 READ(PROBLEM_DATA_HANDLE,NML=WindData)
218! PRINT *, 'iErr=', iErr
219! STOP
220
221
222 !!!!! Calculate wind params based on R-H jump conditions
223 cs_cgs = sqrt((gamma*boltzmann*tempout)/(xmu*amu))
224 mach = velocity/cs_cgs
225 pOut=(rhoOut*Boltzmann*tempOut)
226 density=(4.0*rhoOut)/(1.0+(3.0/Mach**2))
227 pressure=((5.0*Mach**2-1.0)*pOut)/4
228 temperature=pressure/(density*boltzmann)
229 velocity=((3.0*((Mach**2)-1.0))/(4.0*Mach))*cs_cgs
230 IF(mpi_id==0)THEN
231 print *, 'Supersonic inflow boundary conditions:'
232 print *, 'mach = ', mach
233 print *, 'density(cm^-3) =', density
234 print *, 'pressure (Ba) =', pressure
235 print *, 'temperature (K) =', temperature
236 print *, 'velocity (cm/s) =', velocity
237 END IF
238
239
240 CALL CreateWind(Wind)
241 Wind%density=density/nscale
242 Wind%temperature=temperature
243 Wind%velocity=velocity/VelScale
244! Wind%B=B
245 wind_vel=velocity/VelScale
246 CALL AddTracer(Wind%iTracer, 'Wind_Tracer')
247 Wind%dir=dir
248 Wind%edge=edge
249 CLOSE(PROBLEM_DATA_HANDLE)
250
251 END SUBROUTINE ProblemBeforeStep
252
253 !> Could be used to update grids pre-output
254 !! @param Info Info Object
255 SUBROUTINE ProblemAfterStep(Info)
256 TYPE(InfoDef) :: Info
257 END SUBROUTINE ProblemAfterStep
258
259 !> Could be used to set force refinement
260 !! @param Info Info object
261 SUBROUTINE ProblemSetErrFlag(Info)
262 TYPE(InfoDef) :: Info
263 INTEGER :: ip(3,2)
264 REAL(KIND=qPREC) :: wind_pos, dx
265
266 IF (levels(Info%level)%tnow <= wind_refinement_time) THEN
267 dx=levels(Info%level)%dx
268 wind_pos=wind_vel/2d0*levels(Info%level)%tnow
269 ip(:,1)=1
270 ip(:,2)=Info%mX
271 ip(wind%dir,1)=ceiling((wind_pos-wind_thickness-Info%xBounds(wind%dir,1))/dx)
272 ip(wind%dir,2)=ceiling((wind_pos+wind_thickness-Info%xBounds(wind%dir,1))/dx)
273 ip(wind%dir,1)=max(ip(wind%dir,1), 1)
274 ip(wind%dir,2)=min(ip(wind%dir,2),Info%mX(wind%dir))
275 IF (ip(wind%dir,2) >= ip(wind%dir,1)) Info%ErrFlag(ip(1,1):ip(1,2), ip(2,1):ip(2,2), ip(3,1):ip(3,2))=1
276 END IF
277 END SUBROUTINE ProblemSetErrFlag
278
279 SUBROUTINE ProblemBeforeGlobalStep(n)
280 INTEGER :: n
281 END SUBROUTINE ProblemBeforeGlobalStep
282
283
284END MODULE Problem
285