Blog: meeting update: problem.f90

File problem.f90, 6.0 KB (added by Erica Kaminski, 7 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 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
28MODULE 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
44CONTAINS
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
164END MODULE Problem
165