GNU-Darwin Web
i/demo1
 Home Manual Packages Global Index Keywords Quick Reference ``` /* DEMO1.I 1-D hydro code written in Yorick language \$Id: demo1.i,v 1.1 1993/08/27 18:50:06 munro Exp \$ */ /* Copyright (c) 1994. The Regents of the University of California. All rights reserved. */ /* ------------- Set a few parameters --------------------------- */ /* Set initial density and temperature. */ rho0= 1.1845e-3; /* g/cc dry air at NTP */ RT0= 298.16; /* Kelvin at NTP */ /* Set number of zones for hydro calculation, and total length of column. */ n_zones= 200; column_length= 100.0/*cm*/; /* A conversion factor. */ R_gas= 8.314e7; /* erg/K/mole ideal gas constant */ /* Gamma-law gas equation of state parameters. */ gamma= 1.4; /* Cp/Cv for air is 7/5 */ gammaM1= gamma - 1; A_bar= 1.1845e-3/*g/cc dry air*/ * 24.5e3/*cc/mole*/; K_to_v2= 1.e-6*R_gas/A_bar; /* (cm/ms)^2/Kelvin for dry air */ /* Local temperature units are (cm/ms)^2. -- Note dependence of temperature units on A_bar. */ RT0*= K_to_v2; /* ------------- Set initial conditions --------------------------- */ /* Set initial conditions for hydro calculation. */ func reset { extern RT, z, v, M, p; RT= array(RT0, n_zones); /* temperature */ z= span(0, column_length, n_zones+1); /* zone boundary coordinates */ v= array(0.0, dimsof(z)); /* zone bndy velocities */ /* The column consists initially of n_zones zones of equal size containing equal amounts of gas. */ M= rho0*abs(z(2)-z(1)); /* mass/area/zone */ /* The pressure array includes a pressure before and after the column, p(0) and p(n_zones+1), which will be used to set pressure boundary conditions. (Velocity boundary conditions will be applied by setting v(0) and v(n_zones).) Note that the units of this pressure are g*(cm/ms)^2/cc. */ p= array(rho0*RT0, n_zones+2); } /* ------------- Set boundary conditions --------------------------- */ func sound /* DOCUMENT sound Set up the initial conditions for evolve to launch a weak sound wave. SEE ALSO: shock, evolve */ { extern bc0_v, bc0_time, bc0_p, bc0_Z; bc0_v= 0.1*sin(span(0, 2*pi, 100)); bc0_time= span(0, 1.0, 100); bc0_p= bc0_Z= []; reset; } func shock /* DOCUMENT sound Set up the initial conditions for evolve to launch a strong wave, which steepens into a shock as it propagates. SEE ALSO: sound, evolve */ { extern bc0_v, bc0_time, bc0_p, bc0_Z; bc0_v= 10.0*sin(span(0, 2*pi, 100)); bc0_time= span(0, 1.0, 100); bc0_p= bc0_Z= []; reset; } sound; func nobc { extern bcN_v, bcN_p, bcN_Z; bcN_v= bcN_p= bcN_Z= []; } func hardbc { extern bcN_v, bcN_p, bcN_Z; bcN_v= 0; bcN_p= bcN_Z= []; } func softbc { extern bcN_v, bcN_p, bcN_Z; bcN_p= rho0*RT0; bcN_v= bcN_Z= []; } func matchbc { extern bcN_v, bcN_p, bcN_Z; softbc; bcN_Z= rho0*sqrt((gammaM1+1)*RT0); } /* ------------- Define the main function --------------------------- */ /* The DOCUMENT comment will be printed in response to: help, evolve */ func evolve (time1, time0) /* DOCUMENT evolve, time1 or evolve, time1, time0 Step the hydro calculation forward to TIME1, starting with the initial conditions in the RT, z, and v arrays at time TIME0 (default 0.0 if omitted). The calculation also depends on the constants M (mass/area/zone) and gammaM1 (gamma-1 for the gamma-law equation of state). The pressure array p is updated in addition to the primary state arrays RT, z, and v. Boundary conditions are specified by setting either a boundary pressure or a boundary velocity at each end of the fluid column. bc0_v - Boundary velocity at z(0), or [] if z(0) has pressure BC. bc0_p - Boundary pressure beyond z(0). bc0_time - If bc0_v or bc0_p is an array, bc0_time is an array of the same length specifying the corresponding times for time dependent BCs. bc0_Z - Acoustic impedance at z(0) if bc0_v is nil (default is 0). bcN_v, bcN_p, bcN_time, and bcN_Z have the same meanings for the z(n_zones) boundary. The worker routines OutputResults and TakeStep must be supplied. */ { if (is_void(time0)) time0= 0.0; for (time=time0 ; ; time+=dt) { dt= GetTimeStep(); SetBoundaryValues, time, dt; OutputResults, time, dt; if (time >= time1) break; TakeStep, dt; } } /* ----------------- compute time step --------------------- */ func GetTimeStep (dummy) { dz= abs(z(dif)); dv= abs(v(dif)); cs= sqrt((gammaM1+1)*RT); return min( dz / max(courant*cs, accuracy*dv) ) * dt_multiplier; } /* Set reasonable default values for courant and accuracy parameters. */ courant= 2.0; /* number of cycles for sound signal to cross one zone -- must be >=2.0 for numerical stability */ accuracy= 3.0; /* number of cycles for zone volume to change by a factor of ~2 -- must be >1.0 to avoid possible collapse to zero volume */ dt_multiplier= 1.0; /* ----------------- set boundary conditions ------------------ */ func SetBoundaryValues (time, dt) { vtime= time + 0.5*dt; /* velocity is 1/2 step ahead of pressure, z */ /* boundary at z(0) */ if (!is_void(bc0_v)) { /* velocity BC */ v(1)= BCinterp(bc0_v, bc0_time, vtime); } else if (!is_void(bc0_p)) { /* pressure BC */ p(1)= BCinterp(bc0_p, bc0_time, time); /* acoustic impedance BC */ if (!is_void(bc0_Z)) p(1)-= sign(z(2)-z(1))*bc0_Z*v(1); } /* boundary at z(n_zones) (written here as z(0)) */ if (!is_void(bcN_v)) { v(0)= BCinterp(bcN_v, bcN_time, vtime); } else if (!is_void(bcN_p)) { p(0)= BCinterp(bcN_p, bcN_time, time); if (!is_void(bcN_Z)) p(0)+= sign(z(0)-z(-1))*bcN_Z*v(0); } } func BCinterp (values, times, time) { if (numberof(times)<2) return values(1); else return interp(values, times, time); } /* ----------------- produce output ----------------------- */ func OutputVPlot (time, dt) { extern cycle_number; if (time==time0) cycle_number= 0; else cycle_number++; if (!(cycle_number%output_period)) { fma; plg, v, z; zx= sqrt((gammaM1+1)*RT0)*time; if (zx>max(z)) { /* make a dotted "ruler" at the location a sound wave would have reached */ zx= 2*max(z)-zx; if (zx