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 HydroStaticStar 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 HydroStaticStar
|
---|
24 | !! @brief Contains files necessary for the Hydrostatic Equilibrium calculations
|
---|
25 |
|
---|
26 | !> @file problem.f90
|
---|
27 | !! @brief Main file for module Problem
|
---|
28 | MODULE Problem
|
---|
29 | USE DataDeclarations
|
---|
30 | USE Ambients
|
---|
31 | !USE PointGravitySrc
|
---|
32 | !USE ParticleDeclarations !n
|
---|
33 | USE Clumps
|
---|
34 | USE Profiles
|
---|
35 | USE ProfileFunctions
|
---|
36 | USE Winds
|
---|
37 | USE Refinements
|
---|
38 | IMPLICIT NONE
|
---|
39 | SAVE
|
---|
40 |
|
---|
41 | PUBLIC ProblemModuleInit, ProblemGridInit, &
|
---|
42 | ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
|
---|
43 |
|
---|
44 | CONTAINS
|
---|
45 |
|
---|
46 | !> Initializes module variables
|
---|
47 | SUBROUTINE ProblemModuleInit()
|
---|
48 |
|
---|
49 | TYPE(AmbientDef), POINTER :: Ambient
|
---|
50 | TYPE(WindDef), POINTER :: Wind
|
---|
51 | TYPE(ClumpDef), POINTER :: Clump
|
---|
52 | INTEGER :: edge
|
---|
53 | INTEGER :: i, nEntries !soft_function
|
---|
54 | REAL(KIND=qPREC), DIMENSION(3) :: xloc=(/0,0,0/), velocity= (/0,0,0/)
|
---|
55 | REAL(KIND=qPREC) :: rhoOut, pOut, tempOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut, r_cutoff
|
---|
56 |
|
---|
57 |
|
---|
58 | NAMELIST /ClumpData/ xloc, velocity!, soft_radius, soft_function!, mass
|
---|
59 | NAMELIST /AmbientData/ rhoOut, pOut, tempOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut
|
---|
60 |
|
---|
61 | !CALL ClearAllRefinements()
|
---|
62 |
|
---|
63 | OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
|
---|
64 | READ(PROBLEM_DATA_HANDLE,NML=AmbientData)
|
---|
65 | CALL CreateAmbient(Ambient)
|
---|
66 | Ambient%density=rhoOut/rscale
|
---|
67 | !Ambient%density=.1
|
---|
68 | !Ambient%pressure=10
|
---|
69 | Ambient%pressure=pOut/pscale
|
---|
70 | ! Ambient%PersistInBoundaries=.true.
|
---|
71 | CALL UpdateAmbient(Ambient)
|
---|
72 |
|
---|
73 | READ(PROBLEM_DATA_HANDLE,NML=ClumpData)
|
---|
74 | CALL CreateClump(Clump)
|
---|
75 | Clump%position=xloc
|
---|
76 | Clump%velocity=velocity
|
---|
77 | !Clump%density=1 these things are set using the profile object
|
---|
78 | !Clump%temperature=1 ""
|
---|
79 | Clump%thickness=0d0
|
---|
80 | Clump%subsample=1
|
---|
81 |
|
---|
82 | CLOSE(PROBLEM_DATA_HANDLE)
|
---|
83 |
|
---|
84 | OPEN (UNIT=PROBLEM_DATA_HANDLE, FILE='progenitor_profiles.dat', STATUS="OLD", POSITION='REWIND')
|
---|
85 | ! PRINT *, 'OPENED FILE OKAY'
|
---|
86 | READ (PROBLEM_DATA_HANDLE, *) nEntries
|
---|
87 | ! PRINT *, 'nEntries=', nEntries
|
---|
88 | READ (PROBLEM_DATA_HANDLE, *) !Skip the header column
|
---|
89 | CALL CreateProfile(Clump%Profile, nEntries, (/Mass_Field, P_Field/), RADIAL) !Mass field - is gas density
|
---|
90 | DO i=1,nEntries
|
---|
91 | READ(PROBLEM_DATA_HANDLE, *) Clump%profile%data(i,:) ! position, density, pressure in cgs
|
---|
92 | ! PRINT *, 'CLUMP%PROFILE%DATA(I,1)=radius (cm)', 'CLUMP%PROFILE%DATA(I,2)=density (g/cc)', 'CLUMP%PROFILE%DATA(I,3)=press(dyn/cc)'
|
---|
93 | ! PRINT *, CLUMP%PROFILE%DATA(I,1), CLUMP%PROFILE%DATA(I,2), CLUMP%PROFILE%DATA(I,3)
|
---|
94 | Clump%profile%data(i,:)=Clump%profile%data(i,:) / (/lScale, rScale, pScale/) !rescale to cu
|
---|
95 | ! PRINT *, 'CLUMP%PROFILE%DATA(I,1)=radius (cu)', 'CLUMP%PROFILE%DATA(I,2)=density (cu)', 'CLUMP%PROFILE%DATA(I,3)=press(cu)'
|
---|
96 | ! PRINT *, CLUMP%PROFILE%DATA(I,1), CLUMP%PROFILE%DATA(I,2), CLUMP%PROFILE%DATA(I,3)
|
---|
97 | END DO
|
---|
98 | CLOSE(PROBLEM_DATA_HANDLE)
|
---|
99 |
|
---|
100 |
|
---|
101 | OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='profile_init.out')
|
---|
102 | WRITE(PROBLEM_DATA_HANDLE,'(3E25.15)') (/(clump%profile%data(i,:)*(/lscale, rScale, pScale/), i=1,nEntries)/)
|
---|
103 | CLOSE(PROBLEM_DATA_HANDLE)
|
---|
104 |
|
---|
105 | pOut=pOut/pscale
|
---|
106 |
|
---|
107 | CALL Profile_SelfGravityHSE_Outside(Clump%Profile, pOut) !Use scaled ambient pressure from profile data file to recalculate pressure profile in HSE
|
---|
108 | CALL UpdateProfile(Clump%Profile) !don't think this is actually needed -- subroutine is a stub in the code and data values are updated in prev func
|
---|
109 |
|
---|
110 | OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='profile.out')
|
---|
111 | WRITE(PROBLEM_DATA_HANDLE,'(3E25.15)') (/(clump%profile%data(i,:)*(/lscale, rScale, pScale/), i=1,nEntries)/)
|
---|
112 | CLOSE(PROBLEM_DATA_HANDLE)
|
---|
113 |
|
---|
114 | r_cutoff=Clump%profile%data(nEntries,1)
|
---|
115 | Clump%radius=r_cutoff
|
---|
116 | print *, 'clump%radius=', clump%radius, 'clump%radius (cm)=', clump%radius*lscale
|
---|
117 | CALL UpdateClump(Clump)
|
---|
118 |
|
---|
119 | DO i=1,nDim
|
---|
120 | DO edge=1,2
|
---|
121 | IF (Gmthbc(i,edge) == 1) THEN
|
---|
122 | CALL CreateWind(Wind)
|
---|
123 | Wind%density=rhoOut/rScale
|
---|
124 | wind%temperature=tempOut/tempscale
|
---|
125 | wind%type=OUTFLOW_ONLY
|
---|
126 | Wind%dir=i
|
---|
127 | Wind%edge=edge
|
---|
128 | END IF
|
---|
129 | END DO
|
---|
130 | END DO
|
---|
131 |
|
---|
132 | END SUBROUTINE ProblemModuleInit
|
---|
133 |
|
---|
134 | !> Applies initial conditions
|
---|
135 | !! @param Info Info object
|
---|
136 | SUBROUTINE ProblemGridInit(Info)
|
---|
137 | TYPE(InfoDef) :: Info
|
---|
138 | END SUBROUTINE ProblemGridInit
|
---|
139 |
|
---|
140 | !> Applies Boundary conditions
|
---|
141 | !! @param Info Info object
|
---|
142 | SUBROUTINE ProblemBeforeStep(Info)
|
---|
143 | TYPE(InfoDef) :: Info
|
---|
144 | END SUBROUTINE ProblemBeforeStep
|
---|
145 |
|
---|
146 | !> Could be used to update grids pre-output
|
---|
147 | !! @param Info Info Object
|
---|
148 | SUBROUTINE ProblemAfterStep(Info)
|
---|
149 | TYPE(InfoDef) :: Info
|
---|
150 | END SUBROUTINE ProblemAfterStep
|
---|
151 |
|
---|
152 | !> Could be used to set force refinement
|
---|
153 | !! @param Info Info object
|
---|
154 | SUBROUTINE ProblemSetErrFlag(Info)
|
---|
155 | TYPE(InfoDef) :: Info
|
---|
156 | INTEGER :: i,j,k
|
---|
157 | ! Info%ErrFlag(:,:,:)=1 ! refines over entire grid
|
---|
158 | END SUBROUTINE ProblemSetErrFlag
|
---|
159 |
|
---|
160 | SUBROUTINE ProblemBeforeGlobalStep(n)
|
---|
161 | INTEGER :: n
|
---|
162 | END SUBROUTINE ProblemBeforeGlobalStep
|
---|
163 |
|
---|
164 | END MODULE Problem
|
---|
165 |
|
---|