1 /*
2 * $Id: demo1.i,v 1.1 2005-09-18 22:05:55 dhmunro Exp $
3 * 1-D hydro code written in Yorick language
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 /* ------------- Set a few parameters --------------------------- */
12
13 /* Set initial density and temperature. */
14 rho0= 1.1845e-3; /* g/cc dry air at NTP */
15 RT0= 298.16; /* Kelvin at NTP */
16
17 /* Set number of zones for hydro calculation, and
18 total length of column. */
19 n_zones= 200;
20 column_length= 100.0/*cm*/;
21
22 /* A conversion factor. */
23 R_gas= 8.314e7; /* erg/K/mole ideal gas constant */
24
25 /* Gamma-law gas equation of state parameters. */
26 gamma= 1.4; /* Cp/Cv for air is 7/5 */
27 gammaM1= gamma - 1;
28 A_bar= 1.1845e-3/*g/cc dry air*/ * 24.5e3/*cc/mole*/;
29 K_to_v2= 1.e-6*R_gas/A_bar; /* (cm/ms)^2/Kelvin
30 for dry air */
31 /* Local temperature units are (cm/ms)^2. --
32 Note dependence of temperature units on A_bar. */
33 RT0*= K_to_v2;
34
35 /* ------------- Set initial conditions --------------------------- */
36
37 /* Set initial conditions for hydro calculation. */
38 func reset
39 {
40 extern RT, z, v, M, p;
41
42 RT= array(RT0, n_zones); /* temperature */
43 z= span(0, column_length, n_zones+1); /* zone
44 boundary coordinates */
45 v= array(0.0, dimsof(z)); /* zone bndy velocities */
46
47 /* The column consists initially of n_zones zones
48 of equal size containing equal amounts of gas. */
49 M= rho0*abs(z(2)-z(1)); /* mass/area/zone */
50
51 /* The pressure array includes a pressure before and
52 after the column, p(0) and p(n_zones+1), which
53 will be used to set pressure boundary conditions.
54 (Velocity boundary conditions will be applied by
55 setting v(0) and v(n_zones).)
56 Note that the units of this pressure are
57 g*(cm/ms)^2/cc. */
58 p= array(rho0*RT0, n_zones+2);
59 }
60
61 /* ------------- Set boundary conditions --------------------------- */
62
63 func sound
64 /* DOCUMENT sound
65 Set up the initial conditions for evolve to launch a weak sound wave.
66
67 SEE ALSO: shock, evolve
68 */
69 {
70 extern bc0_v, bc0_time, bc0_p, bc0_Z;
71 bc0_v= 0.1*sin(span(0, 2*pi, 100));
72 bc0_time= span(0, 1.0, 100);
73 bc0_p= bc0_Z= [];
74 reset;
75 }
76
77 func shock
78 /* DOCUMENT sound
79 Set up the initial conditions for evolve to launch a strong wave, which
80 steepens into a shock as it propagates.
81
82 SEE ALSO: sound, evolve
83 */
84 {
85 extern bc0_v, bc0_time, bc0_p, bc0_Z;
86 bc0_v= 10.0*sin(span(0, 2*pi, 100));
87 bc0_time= span(0, 1.0, 100);
88 bc0_p= bc0_Z= [];
89 reset;
90 }
91
92 sound;
93
94 func nobc {
95 extern bcN_v, bcN_p, bcN_Z;
96 bcN_v= bcN_p= bcN_Z= [];
97 }
98
99 func hardbc {
100 extern bcN_v, bcN_p, bcN_Z;
101 bcN_v= 0;
102 bcN_p= bcN_Z= [];
103 }
104
105 func softbc {
106 extern bcN_v, bcN_p, bcN_Z;
107 bcN_p= rho0*RT0;
108 bcN_v= bcN_Z= [];
109 }
110
111 func matchbc {
112 extern bcN_v, bcN_p, bcN_Z;
113 softbc;
114 bcN_Z= rho0*sqrt((gammaM1+1)*RT0);
115 }
116
117 /* ------------- Define the main function --------------------------- */
118 /* The DOCUMENT comment will be printed in response to: help, evolve */
119
evolve(time1,time0)120 func evolve(time1, time0)
121 /* DOCUMENT evolve, time1
122 or evolve, time1, time0
123 Step the hydro calculation forward to TIME1,
124 starting with the initial conditions in the
125 RT, z, and v arrays at time TIME0 (default 0.0
126 if omitted). The calculation also depends on
127 the constants M (mass/area/zone) and gammaM1
128 (gamma-1 for the gamma-law equation of state).
129 The pressure array p is updated in addition to
130 the primary state arrays RT, z, and v.
131
132 Boundary conditions are specified by setting
133 either a boundary pressure or a boundary
134 velocity at each end of the fluid column.
135 bc0_v - Boundary velocity at z(0), or []
136 if z(0) has pressure BC.
137 bc0_p - Boundary pressure beyond z(0).
138 bc0_time - If bc0_v or bc0_p is an array,
139 bc0_time is an array of the same
140 length specifying the corresponding
141 times for time dependent BCs.
142 bc0_Z - Acoustic impedance at z(0) if bc0_v
143 is nil (default is 0).
144 bcN_v, bcN_p, bcN_time, and bcN_Z have the same
145 meanings for the z(n_zones) boundary.
146
147 The worker routines OutputResults and
148 TakeStep must be supplied.
149 */
150 {
151 if (is_void(time0)) time0= 0.0;
152
153 for (time=time0 ; ; time+=dt) {
154 dt= GetTimeStep();
155 SetBoundaryValues, time, dt;
156 OutputResults, time, dt;
157 if (time >= time1) break;
158 TakeStep, dt;
159 }
160 }
161
162 /* ----------------- compute time step --------------------- */
163
GetTimeStep(dummy)164 func GetTimeStep(dummy)
165 {
166 dz= abs(z(dif));
167 dv= abs(v(dif));
168 cs= sqrt((gammaM1+1)*RT);
169 return min( dz / max(courant*cs, accuracy*dv) ) * dt_multiplier;
170 }
171 /* Set reasonable default values for courant and
172 accuracy parameters. */
173 courant= 2.0; /* number of cycles for sound signal
174 to cross one zone -- must be >=2.0
175 for numerical stability */
176 accuracy= 3.0; /* number of cycles for zone volume to
177 change by a factor of ~2 -- must be
178 >1.0 to avoid possible collapse to
179 zero volume */
180 dt_multiplier= 1.0;
181
182 /* ----------------- set boundary conditions ------------------ */
183
SetBoundaryValues(time,dt)184 func SetBoundaryValues(time, dt)
185 {
186 vtime= time + 0.5*dt; /* velocity is 1/2 step
187 ahead of pressure, z */
188
189 /* boundary at z(0) */
190 if (!is_void(bc0_v)) {
191 /* velocity BC */
192 v(1)= BCinterp(bc0_v, bc0_time, vtime);
193
194 } else if (!is_void(bc0_p)) {
195 /* pressure BC */
196 p(1)= BCinterp(bc0_p, bc0_time, time);
197 /* acoustic impedance BC */
198 if (!is_void(bc0_Z))
199 p(1)-= sign(z(2)-z(1))*bc0_Z*v(1);
200 }
201
202 /* boundary at z(n_zones) (written here as z(0)) */
203 if (!is_void(bcN_v)) {
204 v(0)= BCinterp(bcN_v, bcN_time, vtime);
205
206 } else if (!is_void(bcN_p)) {
207 p(0)= BCinterp(bcN_p, bcN_time, time);
208 if (!is_void(bcN_Z))
209 p(0)+= sign(z(0)-z(-1))*bcN_Z*v(0);
210 }
211 }
212
BCinterp(values,times,time)213 func BCinterp(values, times, time)
214 {
215 if (numberof(times)<2) return values(1);
216 else return interp(values, times, time);
217 }
218
219 /* ----------------- produce output ----------------------- */
220
OutputVPlot(time,dt)221 func OutputVPlot(time, dt)
222 {
223 extern cycle_number;
224 if (time==time0) cycle_number= 0;
225 else cycle_number++;
226 if (!(cycle_number%output_period)) {
227 fma;
228 plg, v, z;
229 zx= sqrt((gammaM1+1)*RT0)*time;
230 if (zx>max(z)) {
231 /* make a dotted "ruler" at the location a sound wave
232 would have reached */
233 zx= 2*max(z)-zx;
234 if (zx<min(z)) zx= 2*min(zx)-zx;
235 }
236 pldj, zx, min(bc0_v), zx, max(bc0_v), type="dot";
237 }
238 }
239 output_period= 8; /* about 50 frames for wave to
240 transit 200 zones */
241
242 /* OutputResults can be switched among several possibilities */
243 OutputResults= OutputVPlot;
244
245 /* ----------------- 1-D hydro worker ----------------------- */
246
TakeStep1(dt)247 func TakeStep1(dt)
248 {
249 /* velocities 1/2 step ahead of coordinates */
250 z+= dt * v;
251 dv= v(dif);
252 dz= z(dif);
253
254 /* Compute artificial viscosity. */
255 q= array(0.0, n_zones+2);
256 q(2:-1)= q_multiplier * M * max(-dv/dz, 0.0)*abs(dv);
257
258 /* Apply 1st law of thermodynamics to compute
259 temperature change from work p*v(dif) done
260 on zone. Note that p is not time-centered
261 properly here. */
262 RT-= (dt * gammaM1/M) * (p+q)(2:-1) * dv;
263
264 /* Update the pressures from the updated densities
265 M/z(dif) and temperatures RT. Note that p(0)
266 and p(-1) updated by SetBoundaryValues. */
267 p(2:-1)= M * RT/dz;
268
269 /* Apply Newton's 2nd law to update velocities. */
270 v-= (dt/M) * (p+q)(dif);
271 }
272 q_multiplier= 1.0;
273
274 TakeStep= TakeStep1;
275
276 /* ----------------- simple interface ----------------------- */
277
278 func demo1
279 /* DOCUMENT demo1
280 run the 1-D hydrocode demonstration. An X window should pop up
281 in which a movie of a wave propagating down a shock tube appears.
282 Use the 'sound' and 'shock' commands to set two different interesting
283 initial conditions; the default is 'sound'.
284
285 SEE ALSO: sound, shock, evolve
286 */
287 {
288 window, 0, wait=1;
289 reset;
290 evolve, 5;
291 }
292
293 /* ----------------- end of demo1.i ----------------------- */
294