-
Notifications
You must be signed in to change notification settings - Fork 0
/
hyloop.pro
executable file
·240 lines (217 loc) · 9.51 KB
/
hyloop.pro
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
;+
;NAME:
; HyLoop <Based on msu_loopmodel3
;
;PATC taken out
;Viscosity
;Saturated Cond. flux
;
;PUPROSE:
; Front end for numerical loop model (SHrEC). Takes
; an initial state and propagates forward in time.
; Results are saved to disk at 1s time intervals.
;CALLING SEQUENCE:
; loopmodel, g, A, x, heat, interval, [, state=state, lhist=lhist, $
; T0=T0, outfile=outfile, /src, /fal, /uri, depth=depth]
;INPUT PARAMETERS:
; interval - amount of time (s) to run the model forward.
; This gets rounded to the nearest second.
;OPTIONAL KEYWORD INPUTS:
; lhist - An array of state vectors giving the previous
; history of the loop. If lhist is specified, then
; the initial conditions are read from the last
; entry of lhist and state (if given) is ignored.
; On output, lhist contains a history of states.
; state - Initial state structure. Must be supplied if
; lhist is not. On output, contains the final state.
; state.e - internal energy (N volume elements)
; state.n_e - electron density (N volume elements)
; state.v - velocity in x direction at the
; interstices between the volume
; elements (N-1 elements)
; state.time - time of specified state.
; T0 - Temperature at the loop boundaries (chromospheric
; footpoints). If not specified, it is read
; from the zeroth volume element of the state vector.
; outfile - File to store results (lhist, g, A, x, heat). If
; the display device is set to the z-buffer, then
; real time graphical output is suppressed, and the
; present state of the system (state, g, A, x, heat)
; is output as a save file called outfile+'.snap'.
; src - optional switch to cut off radiative losses below
; T=T0, thus holding the chromospheric temperature
; above the boundary condition temperature.
; depth - optional desired chromospheric depth. If this keyword
; is nonzero, then corrections are applied to try to
; maintain a chromosphere of a desired depth (add or
; remove mass at the 1.1*T0 isotherm, which slides
; the isotherm toward desired depth). If
; depth <= 1, then the default depth of 1e6 cm is used.
; If depth is greater than 1, then the value of the
; keyword is the desired chromospheric depth in cm.
; safety - timestepping safety factor; see loop6.pro.
;OUTPUTS:
; lhist - array of state structures, giving the loop history
; as a function of time.
; state - final state of the loop.
;HISTORY:
; 1999-Feb-4 written by Charles C. Kankelborg
; 2009:SHreC
pro hyloop, loop, interval, $
T0=T0, FILE_PREFIX=FILE_PREFIX, $
FILE_EXT=FILE_EXT, src=src, uri=uri, fal=fal, $
depth=depth, safety=safety,QUIET=QUIET, $ ;Begin regrid keywords
REGRID=REGRID, GRID_SAFETY=GRID_SAFETY, $
PERCENT_DIFFERENCE=PERCENT_DIFFERENCE, $
MAX_STEP=MAX_STEP, $
QUADRATIC=QUADRATIC, LSQUADRATIC=LSQUADRATIC, $
SPLINE=SPLINE, SHOWME=SHOWME,$
DELTA_T=DELTA_T, DEBUG=DEBUG, $
HEAT_FUNCTION=HEAT_FUNCTION, $
WINDOW_REGRID=WINDOW_REGRID, $
WINDOW_STATE=WINDOW_STATE, $
SO=SO, E_H=E_H, $
NOVISC=NOVISC, $
NO_SAT=NO_SAT,$
NT_BEAM=NT_BEAM,$
NT_PART_FUNC=NT_PART_FUNC,$ ; PATC section
NT_BREMS=NT_BREMS,$
NT_DELTA_E=NT_DELTA_E,$
NT_DELTA_MOMENTUM=NT_DELTA_MOMENTUM ,$
PATC_heating_rate=PATC_heating_rate,$
extra_flux=extra_flux,$
NO_V_HEAT=NO_V_HEAT , DT_COND_MIN=DT_COND_MIN,$
CONSTANT_CHROMO=CONSTANT_CHROMO, $
SLIDE_CHROMO=SLIDE_CHROMO
;if size(nt_beam, /TYPE) eq 0 then nt_beam=0d0
;if size(nt_brems, /TYPE) eq 0 then nt_brems=0d0
;if size(PATC_heating_rate, /TYPE) eq 0 then PATC_heating_rate=0d0
;if size(extra_flux, /TYPE) eq 0 then extra_flux=0d0
;if size(DELTA_MOMENTUM, /TYPE) eq 0 then DELTA_MOMENTUM=0d0
COMPILE_OPT STRICTARR
;initialize computer time
t_start = systime(1)
;Initialize System variables if necessary
defsysv, !shrec_set, exists=test_set
if test_set le 0 then set_shrec_sys_variables
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Check and set keywords
if size(loop, /TYPE) eq 0 then begin
print,'Hey! What are the initial conditions?'
message,'No loop input!'
goto, end_jump
endif
if not keyword_set(HEAT_FUNCTION)then HEAT_FUNCTION='get_constant_heat'
;file in which to save results
if not keyword_set(FILE_PREFIX) then FILE_PREFIX =''
if not keyword_set(FILE_EXT) then FILE_EXT='loop'
;timing information
;time interval (s) between saved data
if not keyword_set(DELTA_T) then delta_t = 1.0
if not keyword_set(T0) then T0 = loop.state.e[0]/(3*loop.state.n_e[0]*!shrec_kB)
if not keyword_set(REGRID) then REGRID=0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
old_loop=loop
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
num_iterations = long(interval/DELTA_T+0.5) ;round to nearest second
i1 = long(n_elements(old_loop)) ;where are we starting?
loop=old_loop[i1-1]
for i=0ul, num_iterations-1l do begin
case REGRID of
1: regrid_step,loop, SHOWME=SHOWME, GRID_SAFETY=GRID_SAFETY, $
MAX_STEP=MAX_STEP,PERCENT_DIFFERENCE=PERCENT_DIFFERENCE, $
ENERGY_PD=ENERGY_PD, PARTICLE_PD=PARTICLE_PD, $
MOMENTUM_PD=MOMENTUM_PD, VOLUME_PD=VOLUME_PD, $
ENERGY_CHANGE=ENERGY_CHANGE, $
MOMENTUM_CHANGE=MOMENTUM_CHANGE, $
VOLUME_CHANGE=VOLUME_CHANGE ,$
WINDOW=WINDOW_REGRID
2: regrid_step2,loop, SHOWME=SHOWME, GRID_SAFETY=GRID_SAFETY, $
MAX_STEP=MAX_STEP,PERCENT_DIFFERENCE=PERCENT_DIFFERENCE, $
ENERGY_PD=ENERGY_PD, PARTICLE_PD=PARTICLE_PD, $
MOMENTUM_PD=MOMENTUM_PD, VOLUME_PD=VOLUME_PD, $
ENERGY_CHANGE=ENERGY_CHANGE, $
MOMENTUM_CHANGE=MOMENTUM_CHANGE, $
VOLUME_CHANGE=VOLUME_CHANGE ,$
WINDOW=WINDOW_REGRID
else: begin
MAX_STEP=0
ENERGY_PD=0
PARTICLE_PD=0
MOMENTUM_PD=0
VOLUME_PD=0
ENERGY_CHANGE=0
MOMENTUM_CHANGE=0
VOLUME_CHANGE=0
end
endcase
;Get the grid spacing
N = n_elements(loop.state.e) ;figure out grid size
ds = loop.s[1:N-2] - loop.s[0:N-3]
if not keyword_set(QUIET) then Print, 'Number of volume grids: ',N
if keyword_set(CONSTANT_CHROMO) then src=1
;Evolve the loop state in time
shrec, loop, delta_t, T0=T0, $
src=src, fal=fal, uri=uri,$
depth=depth, safety=safety, $
DEBUG=DEBUG, HEAT_FUNCTION=HEAT_FUNCTION, $
SO=SO, E_H=E_H,$
NOVISC=1,$ ;NOVISC, $
NO_SAT=NO_SAT,$
NT_BEAM=NT_BEAM,NT_PART_FUNC=NT_PART_FUNC,$
NT_BREMS=NT_BREMS,$
NT_DELTA_E=NT_DELTA_E,$
NT_DELTA_MOMENTUM=NT_DELTA_MOMENTUM,$
PATC_heating_rate=PATC_heating_rate,$
extra_flux=extra_flux , DT_COND_MIN=DT_COND_MIN, $
CONSTANT_CHROMO=CONSTANT_CHROMO, $
SLIDE_CHROMO=SLIDE_CHROMO
;NO_V_HEAT=NO_V_HEAT
;Save a snapshot of the present state.
flux=!constant_heat_flux
save, loop, ENERGY_PD, PARTICLE_PD, $
MOMENTUM_PD, VOLUME_PD, $
ENERGY_CHANGE, $
MOMENTUM_CHANGE, $
VOLUME_CHANGE,$
E_H,$
nt_beam, nt_brems, PATC_heating_rate, extra_flux, $
DELTA_MOMENTUM,flux ,$
filename=strcompress(FILE_PREFIX $
+string(loop.state.time,FORMAT='(I06)') $
+'.'+FILE_EXT,/REMOVE_ALL), $
/COMPRESS
if keyword_set(SHOWME) then $
stateplot2, loop.s, loop.state, /screen, WINDOW=WINDOW_STATE
;stop
if not keyword_set(QUIET) then begin
print,'Computer: ', strupcase(!computer)
print,'Computer Time: ', systime()
print, 'Heat Function: ',HEAT_FUNCTION
print,'Model time: '+strcompress(string(loop.state.time),/REMOVE_ALL)+ $
's Real time: '+strcompress(string(systime(1)-t_start),/REMOVE_ALL),'s'
print,'Min/Max V: '+strcompress(string(min(loop.state.v)/1d5),/REMOVE_ALL)+ $
'/'+strcompress(string(max(loop.state.v)/1d5),/REMOVE_ALL)+'[km/s]'
print, 'T_max: '+strcompress(string(loop.t_max/1d6),/REMOVE_ALL)+'MK'
print,'grid spacing runs from ',min(ds)/100.0,' m to ',max(ds)/100000.0,' km.'
if keyword_set(CONSTANT_CHROMO) then print, 'HyLoop: Constant Chromosphere Set'
if keyword_set(SLIDE_CHROMO) then print, 'HyLoop: Sliding Chromosphere Set'
help, nt_beam
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
n_cells=n_elements(loop.state.n_e)
print, 'Flux: '+strcompress($
string(int_tabulated(loop.s_alt[1:n_cells-2],loop.e_h)),$
/REMOVE_ALL)+' erg cm^-2 sec^-1'
if not keyword_set(depth) then $
print,'keyword DEPTH = OFF' Else $
print,'keyword DEPTH =',depth
if keyword_set(NOVISC) then $
print,'Viscosity model = OFF' Else $
print,'Viscosity model = ON'
print,'########################################################################'
endif
endfor
end_jump:
t_end= systime(1) ;get computer time
run_time=t_end-t_start
end