Source code for pyro.compressible_rk.problems.convection
"""A heat source in a layer at some height above the bottom will driveconvection in an adiabatically stratified atmosphere."""importnumpyasnpfrompyro.utilimportmsgDEFAULT_INPUTS="inputs.convection"PROBLEM_PARAMS={"convection.dens_base":10.0,# density at the base of the atmosphere"convection.scale_height":4.0,# scale height of the isothermal atmosphere"convection.y_height":2.0,"convection.thickness":0.25,"convection.e_rate":0.1,"convection.dens_cutoff":0.01}
[docs]definit_data(my_data,rp):""" initialize the bubble problem """ifrp.get_param("driver.verbose"):msg.bold("initializing the bubble problem...")# get the density, momenta, and energy as separate variablesdens=my_data.get_var("density")xmom=my_data.get_var("x-momentum")ymom=my_data.get_var("y-momentum")ener=my_data.get_var("energy")gamma=rp.get_param("eos.gamma")grav=rp.get_param("compressible.grav")scale_height=rp.get_param("convection.scale_height")dens_base=rp.get_param("convection.dens_base")dens_cutoff=rp.get_param("convection.dens_cutoff")# initialize the components, remember, that ener here is# rho*eint + 0.5*rho*v**2, where eint is the specific# internal energy (erg/g)xmom[:,:]=0.0ymom[:,:]=0.0dens[:,:]=dens_cutoff# create a seeded random number generatorrng=np.random.default_rng(12345)# set the density to be stratified in the y-directionmyg=my_data.gridp=myg.scratch_array()pres_base=scale_height*dens_base*abs(grav)forjinrange(myg.jlo,myg.jhi+1):profile=1.0-(gamma-1.0)/gamma*myg.y[j]/scale_heightifprofile>0.0:dens[:,j]=max(dens_base*(profile)**(1.0/(gamma-1.0)),dens_cutoff)else:dens[:,j]=dens_cutoffifj==myg.jlo:p[:,j]=pres_baseelifdens[0,j]<=dens_cutoff+1.e-30:p[:,j]=p[:,j-1]else:#p[:, j] = p[:, j-1] + 0.5*myg.dy*(dens[:, j] + dens[:, j-1])*gravp[:,j]=pres_base*(dens[:,j]/dens_base)**gamma# set the ambient conditionsmy_data.set_aux("ambient_rho",dens_cutoff)my_data.set_aux("ambient_u",0.0)my_data.set_aux("ambient_v",0.0)my_data.set_aux("ambient_p",p.v().min())# set the energy (P = cs2*dens) -- assuming zero velocityener[:,:]=p[:,:]/(gamma-1.0)# pairs of random numbers between [-1, 1]vel_pert=2.0*rng.random(size=(myg.qx,myg.qy,2))-1cs=np.sqrt(gamma*p/dens)# make vel_pert have M < 0.05vel_pert[:,:,0]*=0.05*csvel_pert[:,:,1]*=0.05*csidx=dens>2*dens_cutoffxmom[idx]=dens[idx]*vel_pert[idx,0]ymom[idx]=dens[idx]*vel_pert[idx,1]ener[:,:]+=0.5*(xmom[:,:]**2+ymom[:,:]**2)/dens[:,:]
[docs]defsource_terms(myg,U,ivars,rp):"""source terms to be added to the evolution"""S=myg.scratch_array(nvar=ivars.nvar)y_height=rp.get_param("convection.y_height")dist=np.abs(myg.y2d-y_height)e_rate=rp.get_param("convection.e_rate")thick=rp.get_param("convection.thickness")S[:,:,ivars.iener]=U[:,:,ivars.idens]*e_rate*np.exp(-(dist/thick)**2)returnS
[docs]deffinalize():""" print out any information to the user at the end of the run """