#!/net/python/bin/python from math import log #These abbreviations/variables are used throughout # P - Pressure # p - Density # v - Velocity # pv - Momentum # E - Energy # s - Specific Entropy #here are the analytic values that were calculated P_A=[1,0.294,0.294,0.1] p_A=[1,0.48,0.23,0.125] v_A=[0,0.84,0.84,0] pv_A=[p_A[i]*v_A[i] for i in range(4)] E_A=[p_A[i]*pv_A[i]/2+3*(1/2.0)*P_A[i] for i in range(4)] s_A=[log(P_A[i]/(p_A[i]**(5.0/3))) for i in range(4)] ######################################################################## # Initialization of class Shock # ######################################################################## class Shock: def __init__(self,size=100,xstep=0.01,tstep=0.001): self.data=size*[[0.0,0.0,0.0]] #each point has three values: (Density,Momemtum,Energy) self.totmass=0 #total mass of the system self.totenergy=0 #total energy of the system self.totmomentum=0 self.Pdifference=0 self.size=size #number of bins size dx self.xstep=xstep self.tstep=tstep self.range=size*xstep self.steps=0 # steps will increase after each # tstep so we know how many steps have been taken def reset(self,Dens_L=1,Dens_R=0.125,Pres_L=1,Pres_R=0.1): """resets state to initial conditions, i.e. a loaded shock tube""" self.steps=0 for i in range(self.size/2): self.data[i]=[Dens_L,0,3.0*Pres_L/2] #initialize [Density,Momentum,Energy] for left side. Notice v=0 => pv=0 & E=(3/2)P for i in range(self.size/2,self.size): self.data[i]=[Dens_R,0,3.0*Pres_R/2] #set initial values for energy momentum, and mass: self.totmass=(1.0/2)*(self.size*Dens_L+self.size*Dens_R) En=0 Mom=0 for i in range(self.size): En+=self.data[i][2] Mom+=self.data[i][1] self.totenergy=En self.totmomentum=Mom self.Pdifference=Pres_L-Pres_R ######################################################################## # Analytic Solutions # ######################################################################## def analytic(self,P0,P1,P2,P3,vr=-1.291,vcd=0.841,vs=1.844): """An algorithm for creating a list containing my analytic solution of the SOD conditions""" dx=self.xstep dt=self.tstep x0=int(self.size/2) xr=int(x0+vr*dt*self.steps/dx) xcd=int(x0+vcd*dt*self.steps/dx) xs=int(x0+vs*dt*self.steps/dx) if xr<0: xr=0 if xcd>self.size: xcd=self.size if xs>self.size: xs=self.size D=[] for x in range(0,xr): D.append(P0) for x in range(xr,x0): y=(P0-P1)*(x-xr)**2/(x0-xr)**2-2*(P0-P1)*(x-xr)/(x0-xr)+P0 #quadratic fit for the rarefaction region D.append(y) for x in range(x0,xcd): D.append(P1) for x in range(xcd,xs): D.append(P2) for x in range(xs,self.size): D.append(P3) return D def analyticpressure(self): """Returns a list containing the analytic solution for Sod's initial conditions for pressure""" return self.analytic(P_A[0],P_A[1],P_A[2],P_A[3]) def analyticdensity(self): """Returns a list containing the analytic solution for Sod's initial conditions for density""" return self.analytic(p_A[0],p_A[1],p_A[2],p_A[3]) def analyticvelocity(self): return self.analytic(v_A[0],v_A[1],v_A[2],v_A[3]) def analyticmomentum(self): return self.analytic(pv_A[0],pv_A[1],pv_A[2],pv_A[3]) def analyticenergy(self): return self.analytic(E_A[0],E_A[1],E_A[2],E_A[3]) def analyticentropy(self): return self.analytic(s_A[0],s_A[1],s_A[2],s_A[3]) ######################################################################## # Functions for grabbing data # ######################################################################## def time(self): return self.tstep*self.steps def xval(self): """returns a list of x values for our data""" X=[] for i in range(self.size): X.append(i*self.xstep) return X def density(self,n=-1): """returns p[n] or a list of density values""" if n in range(self.size): return self.data[n][0] else: return [D[0] for D in self.data] def energy(self,n=-1): """returns E[n] or a list of energy values""" if n in range(self.size): return self.data[n][2] else: return [E[2] for E in self.data] def momentum(self,n=-1): """returns pv[n] or a list of momentum values""" if n in range(self.size): return self.data[n][1] else: return [P[1] for P in self.data] def pressure(self,n=-1): """returns P[n] or a list of pressure values""" if n in range(self.size): return (2.0/3)*self.data[n][2]-(1.0/3)*self.data[n][1]*self.data[n][1]/self.data[n][0] else: return [(2.0/3)*D[2]-(1.0/3)*D[1]*D[1]/D[0] for D in self.data] #P=(2/3)E-(1/3)p(v^2) def velocity(self,n=-1): """returns v[n] or a list of velocity values""" if n in range(self.size): return self.data[n][1]/self.data[n][0] else: return [D[1]/D[0] for D in self.data] #v=pv/p def entropy(self,n=-1): if n in range(self.size): return log(self.pressure(n)/self.density(n)**(5.0/3)) else: return [log(self.pressure(i)/self.density(i)**(5.0/3)) for i in range(self.size)] ######################################################################## # Conservation Checks # ######################################################################## def masscons(self): """Returns fraction of mass lost in code""" return (self.totmass-sum(self.density()))/self.totmass def energycons(self): """Returns fraction of energy lost in code""" return (self.totenergy-sum(self.energy()))/self.totenergy def momentumcons(self): """Returns fraction of momentum lost in code""" mom=100*(self.totmomentum+self.Pdifference*self.time()) return (mom-sum(self.momentum()))/mom ######################################################################## # Time Step Functions # ######################################################################## def order1lax(self,n=1): """Takes our present state and performs n steps using lax method""" #Note: does not work very well D=[A[:] for A in self.data] #A new list to store time stepped values dt=self.tstep #length of time step dx=self.xstep #length of distance step for num in range(n): #now perform algorithm on all but the end points # so boundary conditions can be ignored (until shock hits boundaries) P=self.pressure() v=self.velocity() p=self.density() pv=self.momentum() E=self.energy() self.steps+=1 for k in range(1,self.size-1): ##Compute using the following algorithm: ## ## dq/dt=(q[t+1,x]-(q[t,x+1]+q[t,x-1])/2)/dt ## ## dq/dx=(q[t,x+1]-q[t,x-1])/(2dx) ## ##Note: [0] refers to density ## [1] refers to momentum ## [2] refers to energy ##now perform the time step # p=p0-dt(d/dx)(pv) D[k][0]=(p[k+1]+p[k-1])/2-dt*(pv[k+1]-pv[k-1])/(2*dx) # pv=pv0-dt(d/dx)(pv^2+P) D[k][1]=(pv[k+1]+pv[k-1])/2-dt*((P[k+1]+pv[k+1]*v[k+1])-(P[k-1]+pv[k-1]*v[k-1]))/(2*dx) # E=E0-dt(d/dx)[v(E+P)] D[k][2]=(E[k+1]+E[k-1])/2-dt*(v[k+1]*(E[k+1]+P[k+1])-v[k-1]*(E[k-1]+P[k-1]))/(2*dx) self.data=[A[:] for A in D] def order2lax(self,n=1,art_viscosity=True): """Takes our present state and performs n steps using 2nd order lax method""" D=[A[:] for A in self.data] #A new list to store time stepped values dt=self.tstep #length of time step dx=self.xstep #length of distance step for num in range(n): #now perform algorithm on all but the end points # so boundary conditions can be ignored (until shock hits boundaries) #first we update the value of mid: P=self.pressure() v=self.velocity() p=self.density() pv=self.momentum() E=self.energy() self.steps+=1 #define holders for the middle values pM=[0,0] pvM=[0,0] EM=[0,0] PM=[0,0] vM=[0,0] for k in range(1,self.size-1): ##Compute using the following algorithm: ## ##First use a lax step to get a t+1/2 value ## then another step on these values to achieve second order #define artificial viscosity # for [x-1,x,x+1] Q=[0,0,0] if art_viscosity: q=1 if v[k-1]>v[k]: Q[0]=p[k-1]*(q*(v[k]-v[k-1]))**2 else: Q[0]=0 if v[k-1]>v[k+1]: Q[1]=p[k]*(q*(v[k+1]-v[k-1])/2)**2 else: Q[1]=0 if v[k]>v[k+1]: Q[2]=p[k+1]*(q*(v[k+1]-v[k]))**2 else: Q[2]=0 else: pass ##now perform the time step # p[k-1/2] pM[0]=(p[k]+p[k-1])/2-dt*(pv[k]-pv[k-1])/(2*dx) # pv[k-1/2] pvM[0]=(pv[k]+pv[k-1])/2-dt*((P[k]+Q[1]+pv[k]*v[k])-(P[k-1]+Q[0]+pv[k-1]*v[k-1]))/(2*dx) # E[k-1/2] EM[0]=(E[k]+E[k-1])/2-dt*(v[k]*(E[k]+P[k]+Q[1])-v[k-1]*(E[k-1]+P[k-1]+Q[0]))/(2*dx) # p[k+1/2] pM[1]=(p[k+1]+p[k])/2-dt*(pv[k+1]-pv[k])/(2*dx) # pv[k+1/2] pvM[1]=(pv[k+1]+pv[k])/2-dt*((P[k+1]+Q[2]+pv[k+1]*v[k+1])-(P[k]+Q[1]+pv[k]*v[k]))/(2*dx) # E[k+1/2] EM[1]=(E[k+1]+E[k])/2-dt*(v[k+1]*(E[k+1]+P[k+1]+Q[2])-v[k]*(E[k]+P[k]+Q[1]))/(2*dx) # calculate values for v[k+-1/2] and P[k+-1/2] for t+1/2 data vM=[pvM[i]/pM[i] for i in (0,1)] PM=[(2.0/3)*EM[i]-(1.0/3)*pvM[i]*pvM[i]/pM[i] for i in (0,1)] # define artificial viscosity for the t+1/2 data if art_viscosity and vM[0]>vM[1]: QM=[rho*(q*(vM[1]-vM[0]))**2 for rho in pM] else: QM=[0,0] # p[t+1] D[k][0]=p[k]-dt*(pvM[1]-pvM[0])/dx # pv[t+1] D[k][1]=pv[k]-dt*((PM[1]+QM[1]+pvM[1]*vM[1])-(PM[0]+QM[0]+pvM[0]*vM[0]))/dx # E=[t+1] D[k][2]=E[k]-dt*(vM[1]*(EM[1]+PM[1]+QM[1])-vM[0]*(EM[0]+PM[0]+QM[0]))/dx self.data=[A[:] for A in D] def upwind(self,n=1): """Takes our present state and performs n steps using lax method""" D=[A[:] for A in self.data] #A new list to store time stepped values dt=self.tstep #length of time step dx=self.xstep #length of distance step for num in range(n): #now perform algorithm on all but the end points # so boundary conditions can be ignored (thus we cannot let waves hit the boundaries) self.steps+=1 P=self.pressure() v=self.velocity() p=self.density() pv=self.momentum() E=self.energy() for k in range(1,self.size-1): #set s equal to the sign of v if v[k]<0: s=-1 else: s=1 # s=1 for v=0 # this last condition turns out to be very important. # for the v=0 preshock region, we want it to be affected # by the shock to the left, so we treat it as if the # velocity is positive. ##now perform the upwind time step # if s=-1 then second term is -dt*(Q[k+1]-Q[k])/dx # if s=+1 then second term is -dt*(Q[k]-Q[k-1])/dx # p=p0-dt(d/dx)(pv) D[k][0]=p[k]-s*(dt/dx)*(pv[k]-pv[k-s]) # pv=pv0+dt(d/dx)(pv*v+P) D[k][1]=pv[k]-s*(dt/dx)*((pv[k]*v[k])-(pv[k-s]*v[k-s]))-(dt/(2*dx))*(P[k+1]-P[k-1]) # E=E0+dt(d/dx)(v(E+P)) D[k][2]=E[k]-s*(dt/dx)*(v[k]*(E[k]+P[k])-v[k-s]*(E[k-s]+P[k-s])) self.data=[A[:] for A in D] ######################################################################## # Let the fun begin... # ######################################################################## if __name__=='__main__': # Import graphing utilities import pylab x=Shock(120,0.01,0.001) x.reset() x.upwind(250) print 'Upwind method:\nFraction of mass lost: '+ str(x.masscons()) print 'Franction of energy lost: '+str(x.energycons()) print 'Fraction of momentum lost: '+str(x.momentumcons()) pylab.figure(1) pylab.subplot(311) pylab.title('Momentum') pylab.plot(x.xval(),x.momentum()) pylab.plot(x.xval(),x.analyticmomentum()) pylab.subplot(312) pylab.title('Velocity') pylab.plot(x.xval(),x.velocity()) pylab.plot(x.xval(),x.analyticvelocity()) pylab.subplot(313) pylab.title('Pressure') pylab.plot(x.xval(),x.pressure()) pylab.plot(x.xval(),x.analyticpressure()) pylab.axis([0,1.2,0,1.1]) pylab.show() pylab.figure(1) pylab.subplot(311) pylab.title('Density') pylab.plot(x.xval(),x.density()) pylab.plot(x.xval(),x.analyticdensity()) pylab.axis([0,1.2,0,1.1]) pylab.subplot(312) pylab.title('Entropy') pylab.plot(x.xval(),x.entropy()) pylab.plot(x.xval(),x.analyticentropy()) pylab.axis([0,1.2,-0.2,1.4]) pylab.subplot(313) pylab.title('Energy') pylab.plot(x.xval(),x.energy()) pylab.plot(x.xval(),x.analyticenergy()) pylab.show()