u/erica/memoryBaseline: problem.f90

File problem.f90, 14.3 KB (added by Erica Kaminski, 11 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 MolecularCloudFormation 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 MolecularCloudFormation
24!! @brief Contains files necessary for the Molecular Cloud Formation problem
25
26!> @file problem.f90
27!! @brief Main file for module Problem
28
29!> @defgroup MolecularCloudFormation Molecular Cloud Formation Module
30!! @brief Module for setting up orbiting particles
31!! @ingroup Modules
32
33!> Molecular Cloud Formation Module
34!! @ingroup MolecularCloudFormation
35MODULE Problem
36 USE DataDeclarations
37 USE ParticleDeclarations
38 USE ProcessingDeclarations
39 USE CollidingFlows
40 USE Winds
41 USE Ambients
42 USE Histograms
43 USE Fields
44 USE Totals
45 USE PDFs
46 USE Projections
47 USE Clumps
48 USE Shapes
49 USE Refinements
50 USE IICooling
51 IMPLICIT NONE
52 SAVE
53
54 REAL(KIND=qPREC), DIMENSION(0:MaxDepth) :: cells_per_cooling_length=0
55 REAL(KIND=qPREC), DIMENSION(0:MaxDepth) :: InterfaceWidth=0
56 TYPE(CollidingFlowDef), POINTER :: CollidingFlow
57 TYPE(ShapeDef), POINTER :: DerefineShape=>Null()
58 PUBLIC ProblemModuleInit, ProblemGridInit, &
59 ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
60 REAL(KIND=qPREC) :: IIScaleCool, InterfaceTime, IIScaleHeat, derefine_radius=1e20
61 INTEGER :: derefine_dist=0
62 TYPE(AmbientDef), POINTER :: Ambient
63 LOGICAL :: DeRefineOutSide=.false.
64 LOGICAL :: lClumps=.false.
65 LOGICAL :: lPlaceClumps=.true.
66 REAL(KIND=qPREC) :: ClumpJeansFact = .5d0 !Make them 'X=.5' jeans lengths in diameter.
67 REAL(KIND=qPREC) :: MeanDensitywClumps=3d0
68 REAL(KIND=qPREC) :: ClumpChi=10
69 REAL(KIND=qPREC) :: separation_param=.3d0
70 REAL(KIND=qPREC) :: TShutOff=1d30
71 REAL(KIND=qPREC) :: RampTime=1d0
72
73CONTAINS
74
75 !> Initializes module variables
76 SUBROUTINE ProblemModuleInit()
77 INTEGER :: nRegions, nWaves
78 REAL(KIND=qPREC) :: density, temperature, smooth_distance, velocity, position(3), size_param(3), psi, theta, phi, wavevector(2), amplitude, phase, vel(3), a, alpha, dk, kmag, beta, mach, interface_dist, clumprho, clumptemp, clumpff, clumpradius, clumpnumbervolume, region(3,2), min_separation, ram_dens, ram_press, ram_temp, phi2
79 REAL(KIND=qPREC), DIMENSION(:,:), POINTER :: positions
80 INTEGER :: i, j, type, subsample=1, smooth_function, edge, kx, ky, interface_func
81 TYPE(WindDef), POINTER :: Wind
82 TYPE(HistogramDef), POINTER :: HISTOGRAM
83 TYPE(PDFDef), POINTER :: PDF
84 TYPE(ProjectionDef), POINTER :: Projection
85 INTEGER :: nWinds, nbins=100, clumprefinelevel
86 INTEGER, DIMENSION(10) :: myseed
87 INTEGER :: kmax,nclumps,nParticleCells=0
88 TYPE(ClumpDef), POINTER :: Clump
89 TYPE(ShapeDef), POINTER :: LeftClumpRegion, RightClumpRegion, TotalShape
90 LOGICAL :: lPersist=.true.
91 REAL(KIND=qpREC) :: TotalWidth, FlowDiameter(2), ShearAngle, MaxDiameter
92 TYPE(RefinementDef), POINTER :: Refinement
93 INTEGER :: MustRefineDepth=MAXDEPTH
94 INTEGER :: CanRefineDepth=MAXDEPTH
95 REAL(KIND=qPREC) :: MustRefineWidth=0d0, CanRefineWidth=1d30
96 NAMELIST /ProblemData/ density, velocity, mach, temperature, smooth_distance, smooth_function, subsample, &
97 myseed, A, alpha, kmax, beta, &
98 interface_func, interface_dist, nbins, lClumps, ClumpJeansFact, MeanDensitywClumps, &
99 ClumpChi, separation_param, nParticleCells,CellsPerJeansLength, clumprefinelevel, IICoolPar, TShutOff, RampTime, lPlaceClumps, &
100 TotalWidth, FlowDiameter, ShearAngle, MustRefineWidth, MustRefineDepth, CanRefineWidth, CanRefineDepth
101
102 OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
103 READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
104 CLOSE(PROBLEM_DATA_HANDLE)
105
106 write(*,*) 'Shear angle=', ShearAngle
107 ! Rescale various physical parameters to computational units
108 FlowDiameter=FlowDiameter*pc/lScale
109 density=density/nScale
110 temperature=temperature/TempScale
111 velocity=velocity*1e5/VelScale
112 TotalWidth=TotalWidth*pc/lScale
113 RampTime=RampTime*1e6*yr/TimeScale
114 TShutOff=TShutOff*1e6*yr/TimeScale
115 InterfaceTime=InterfaceTime*1e6*yr/TimeScale
116
117 IF (nParticleCells > 0) DefaultParticleRadius=nParticleCells
118
119
120 IF (temperature == 0d0) THEN
121 temperature=GetIICoolEqTemp(density*nScale)/TempScale
122 IF (MPI_ID == 0) write(*,*) 'Background Density = ', density*nScale, ' particles/cc'
123 IF (MPI_ID == 0) write(*,*) 'Equilibrium temp = ', temperature*TempScale, ' K'
124 IF (MPI_ID == 0) write(*,*) 'sound speed = ', sqrt(gamma*temperature), sqrt(gamma*temperature)*velscale/1e5, 'km/s'
125 END IF
126 IF (mach /= 0d0) velocity=mach*sqrt(gamma*temperature)
127
128 IF (MPI_ID == 0) write(*,*) 'Background Jeans Length =', JeansLengthFromRhoTemp(density, temperature*TempScale)*lScale/pc, ' pc'
129 IF (MPI_ID == 0) write(*,*) 'Background free fall time =', sqrt(3d0*pi/32d0/ScaleGrav/density)*TimeScale/yr/1e6, 'myr'
130 IF (MPI_ID == 0) write(*,*) 'flow velocity = ', velocity*velscale/1e5, 'km/s'
131
132 CALL CreateAmbient(Ambient)
133 Ambient%density=density
134 Ambient%pressure=density*temperature
135 IF (beta /= 0d0) THEN
136 Ambient%B=sqrt(2d0*temperature*density/beta)*(/1d0,0d0,0d0/)
137 IF (MPI_ID == 0) write(*,'(A,3E25.15)') 'Ambient%B=',Ambient%B
138 END IF
139
140 DO i=1,nDim
141 DO edge=1,2
142 IF (Gmthbc(i,edge) == 1) THEN
143 CALL CreateWind(Wind)
144 Wind%dir=i
145 Wind%edge=edge
146 Wind%type=OUTFLOW_ONLY
147 Wind%density=Ambient%density
148 Wind%temperature=Ambient%Pressure/Ambient%Density
149 Wind%B=Ambient%B
150 END IF
151 END DO
152 END DO
153
154
155
156 CALL CreateCollidingFlow(CollidingFlow)
157
158 MaxDiameter=maxval(FlowDiameter)
159 CollidingFlow%SubSample = subsample
160 CollidingFlow%density = density
161 CollidingFlow%velocity = velocity
162 CollidingFlow%Temperature = Temperature
163 CollidingFlow%PersistInBoundaries = lPersist
164 CollidingFlow%smooth_function = smooth_function
165 CollidingFlow%smooth_distance = smooth_distance*MaxDiameter
166 CollidingFlow%interface_func = interface_func
167 CollidingFlow%interface_dist = interface_dist*MaxDiameter
168 CollidingFlow%tShutOff = TShutOff
169 CollidingFlow%RampTime = RampTime
170 CALL AddTracer(CollidingFlow%iTracer(1),'FlowTracer1')
171 CALL AddTracer(CollidingFlow%iTracer(2),'FlowTracer2')
172
173 CALL SetShapeType(CollidingFlow%Shape, ELLIPTICAL_PRISM, half*(/FlowDiameter(:),2d0*(GxBounds(1,2)-GxBounds(1,1))/))
174 CALL SetShapeOrientation(CollidingFlow%Shape, 0d0, Pi/2d0, 0d0)
175 CollidingFlow%Shape%Position=half*SUM(GxBounds(:,:), 2)
176 CALL SetShapeBounds(CollidingFlow%Shape)
177
178 IF (derefine_radius > 0d0) THEN
179 CALL CreateShape(DerefineShape)
180 CALL SetShapeType(DerefineShape, ELLIPTICAL_PRISM, (/derefine_radius*half*FlowDiameter(:), 2d0*(GxBounds(1,2)-GxBounds(1,1))/))
181 CALL SetShapeOrientation(DerefineShape, 0d0, Pi/2d0, 0d0)
182 DerefineShape%position=half*SUM(GxBounds(:,:),2)
183 CALL SetShapeBounds(DerefineShape)
184 END IF
185 IF (MPI_ID == 0) THEN
186 write(*,*) 'mass flux = ', 2d0*velocity*velScale * (SUM(FlowDiameter(1:nDim-1)**2))*Pi/4d0*density*rScale*lScale**2/(mSolar/(1e6*yr)), 'm_sun/Myr'
187 ram_press=velocity**2*TempScale*density*nScale
188 write(*,*) 'ram pressure = ', ram_press, 'particles K/cc'
189 ram_dens=GetIICoolEqDensity(ram_press)
190 write(*,*) 'ram density = ', ram_dens, 'particles per cc'
191 ram_temp=ram_press/ram_dens
192 write(*,*) 'ram temp = ', ram_temp, 'K'
193 write(*,*) 'ram density jeans length =', JeansLengthFromRhoTemp(ram_dens/nScale, ram_temp)*lScale/pc, ' pc'
194 END IF
195 ! Then set up interface
196 CollidingFlow%InterfaceObj%Position=CollidingFlow%Shape%Position
197
198 IF (nDim == 2) THEN
199 CALL SetInterfaceOrientation(CollidingFlow%InterfaceObj, Pi/2d0,ShearAngle/180d0*Pi)
200 ELSE
201 CALL SetInterfaceOrientation(CollidingFlow%InterfaceObj, Pi/2d0-ShearAngle/180d0*Pi,0d0)
202 END IF
203 CALL RANDOM_SEED(size=j)
204 CALL RANDOM_SEED(PUT=myseed(1:j))
205 nWaves=0
206 DO kx=1, kMax
207 DO ky=1, kMax
208 kmag=sqrt(real(kx**2+ky**2))
209 if (kmag < kMax) THEN
210 nwaves=nwaves+1
211 END if
212 END DO
213 END DO
214 dk=2d0*Pi/MaxDiameter
215 CALL InitInterfaceWaves(CollidingFlow%InterfaceObj, nWaves)
216 DO kx=1, kMax
217 DO ky=1, kMax
218 kmag=sqrt(real(kx**2+ky**2))
219 if (kmag < kMax) THEN
220 CALL RANDOM_NUMBER(phase)
221 phase=phase*2d0*Pi
222 wavevector=(/kx*dk,ky*dk/)
223 amplitude=A*kmag**alpha*MaxDiameter
224 CALL AddInterfaceWave(CollidingFlow%InterfaceObj, wavevector, phase, amplitude)
225 ! write(*,*) wavevector, amplitude, phase
226 END if
227 END DO
228 END DO
229
230 CALL ClearAllRefinements()
231
232
233 !First create refinement object to refine around interface
234 CALL CreateRefinement(Refinement)
235 CALL CreateShape(Refinement%Shape)
236 CALL SetShapeType(Refinement%Shape, ELLIPTICAL_PRISM, half*(/FlowDiameter(1)/cos(ShearAngle/180d0*Pi), FlowDiameter(2), MustRefineWidth/))
237 IF (nDim == 2) THEN
238 CALL SetShapeOrientation(Refinement%Shape, 0d0, Pi/2d0, ShearAngle/180d0*Pi)
239 ELSE
240 CALL SetShapeOrientation(Refinement%Shape, 0d0, Pi/2d0-ShearAngle/180d0*Pi, 0d0)
241 END IF
242 Refinement%Shape%Position=CollidingFlow%Shape%Position
243 CALL UpdateShape(Refinement%Shape)
244 Refinement%type=REFINE_INSIDE
245! Refinement%field=Mass_Field
246! Refinement%tolerance=1
247! Refinement%scale=LOGSCALE
248 Refinement%BufferCells=1
249 Refinement%MaxLevel=MustRefineDepth
250
251 CALL CreateRefinement(Refinement)
252 CALL CreateShape(Refinement%Shape)
253 CALL SetShapeType(Refinement%Shape, ELLIPTICAL_PRISM, half*(/FlowDiameter(1)/cos(ShearAngle/180d0*Pi), FlowDiameter(2), CanRefineWidth/))
254 IF (nDim == 2) THEN
255 CALL SetShapeOrientation(Refinement%Shape, 0d0, Pi/2d0, ShearAngle/180d0*Pi)
256 ELSE
257 CALL SetShapeOrientation(Refinement%Shape, 0d0, Pi/2d0-ShearAngle/180d0*Pi, 0d0)
258 END IF
259 Refinement%Shape%Position=CollidingFlow%Shape%Position
260 CALL SetShapeBounds(Refinement%Shape)
261 Refinement%field=Mass_Field
262 Refinement%tolerance=1
263 Refinement%scale=LOGSCALE
264 Refinement%BufferCells=2
265 Refinement%MaxLevel=CanRefineDepth
266
267
268 CALL AddRefinementThreshold(JeansLength_Field, LESSTHAN, (/(CellsPerJeansLength*levels(i)%dx,i=0,MaxLevel)/))
269
270
271! CALL AddDiagnosticVar(ErrFlag_Field)
272
273
274! Uncomment this to turn off processing stuff
275 RETURN
276
277 CALL CreateShape(TotalShape)
278 CALL SetShapeType(TotalShape, ELLIPTICAL_PRISM, half*(/FlowDiameter(:),TotalWidth/))
279 CALL SetShapeOrientation(TotalShape, 0d0, Pi/2-ShearAngle/180d0*Pi,0d0)
280 TotalShape%position(:) = CollidingFlow%Shape%position
281 CALL SetShapeBounds(TotalShape)
282
283 CALL AddAllTotals(GASCOMP, TotalShape)
284 CALL AddAllTotals(PARTICLECOMP,TotalShape)
285 CALL AddAllTotals(BOTHCOMP,TotalShape)
286
287 CALL CreatePDF(PDF)
288 PDF%Field(1)%iD=Mass_Field
289 PDF%Field(1)%name='density'
290 PDF%Field(1)%component=GASCOMP
291 PDF%Field(2)%iD=VMag_Field
292 PDF%Field(2)%name='velocity'
293 PDF%Field(2)%component=GASCOMP
294 PDF%minvalue=(/.01,.001/)
295 PDF%maxvalue=(/1e7,1e4/)
296 PDF%nbins=(/400,400/)
297 PDF%Scale=(/LOGSCALE,LOGSCALE/)
298 PDF%WeightField=BINBYVOLUME
299 PDF%Shape=>TotalShape
300
301
302 CALL CreateHistogram(Histogram)
303 Histogram%Field%iD=1
304 Histogram%Field%name='density'
305 Histogram%Field%component=GASCOMP
306 Histogram%minvalue=.1d0
307 Histogram%maxvalue=1d8
308 Histogram%nbins=nbins
309 Histogram%scale=LOGSCALE
310 Histogram%shape=>TotalShape
311
312
313 CALL CreateHistogram(Histogram)
314 Histogram%Field%iD=MixingRatio12_Field
315 Histogram%Field%name='Mixing_Ratio'
316 Histogram%Field%component=GASCOMP
317 Histogram%minvalue=0d0
318 Histogram%maxvalue=1d0
319 Histogram%nbins=nbins
320 Histogram%scale=LINEARSCALE
321 Histogram%WeightField=BINBYVOLUME
322 Histogram%shape=>TotalShape
323
324 CALL CreatePDF(PDF)
325 PDF%Field(:)%iD=(/Mass_Field, P_Field/)
326 PDF%Field(1)%name='density'
327 PDF%Field(2)%name='pressure'
328 PDF%Field(:)%component=GASCOMP
329 PDF%minvalue=(/.01,100.0/)
330 PDF%maxvalue=(/1e7,1e6/)
331 PDF%nbins=(/400,400/)
332 PDF%Scale=(/LOGSCALE,LOGSCALE/)
333 PDF%WeightField=BINBYMASS
334 PDF%Shape=>TotalShape
335
336
337 CALL CreateProjection(projection)
338 Projection%Field(1)%iD=Mass_Field
339 Projection%Field(1)%component=BOTHCOMP
340 Projection%dim=1
341! Projection%Shape=>TotalShape
342
343 CALL CreateProjection(projection)
344 Projection%Field(1)%iD=Mass_Field
345 Projection%Field(1)%component=BOTHCOMP
346 Projection%dim=2
347
348 CALL CreateProjection(projection)
349 Projection%Field(1)%iD=Mass_Field
350 Projection%Field(1)%component=BOTHCOMP
351 Projection%dim=3
352
353 END SUBROUTINE ProblemModuleInit
354
355
356 !> Applies initial conditions
357 !! @param Info Info object
358 SUBROUTINE ProblemGridInit(Info)
359 TYPE(InfoDef) :: Info
360 END SUBROUTINE ProblemGridInit
361
362 !> Could be used to update grids pre-step
363 !! @param Info Info Object
364 SUBROUTINE ProblemBeforeStep(Info)
365 TYPE(InfoDef) :: Info
366 END SUBROUTINE ProblemBeforeStep
367
368 !> Could be used to update grids pre-output
369 !! @param Info Info Object
370 SUBROUTINE ProblemAfterStep(Info)
371 TYPE(InfoDef) :: Info
372
373 END SUBROUTINE ProblemAfterStep
374
375
376 !> Could be used to set force refinement
377 !! @param Info Info object
378 SUBROUTINE ProblemSetErrFlag(Info)
379 TYPE(InfoDef) :: Info
380 END SUBROUTINE ProblemSetErrFlag
381
382
383 SUBROUTINE ProblemBeforeGlobalStep(n)
384 INTEGER :: n
385 END SUBROUTINE ProblemBeforeGlobalStep
386
387
388 SUBROUTINE SetupTotals()
389
390 END SUBROUTINE SetupTotals
391
392END MODULE Problem
393