#!/usr/bin/env python """ Use Maccormack's technique to solve the one-dimensional subsonic-supersonic isentropic nozzle flow. Provision of nozzle profile with external file is possible. Calculate analytic solution for ideal-gas and non-ideal-gas gamma. Author: Bob Wimmer Date: January 17, 2009 """ from numpy import * import Gnuplot, Gnuplot.funcutils import time #global definitions gamma = 1.4 #ratio of specific heats, gamma=5/3 for ideal mono-atomic gas C = 0.999 #Courant number (0.5 is a good choice) nx = 60 #number of grid intervals L = 3. #Length of nozzle in meters x = linspace(0.,L,nx+1) #Array along the nozzle dx = L/float(nx) #x-interval along nozzle A = zeros(nx+1) #Profile of the nozzle's cross section R = 8.314 #Universal gas constant in J/(mol K) T_env = 293.15 #degrees Kelvin (room temperature) t_scale = L/sqrt(gamma*R*T_env) #time scale res_acc = 1.e-10 #Accuracy in reduals to which calculations are carried out #Provide cross section profile via externam file infilename = 'A_prof_4.dat' infile = open(infilename,'r') i = 0 for line in infile: #print line bla = line.split(' ') A[i] = float(bla[1]) i = i+1 infile.close() #Profide cross section profile via analytic expression for i in xrange(0,nx+1,1): A[i] = 1. + 2.2*(x[i] - 1.5)**2 #print x[i], A[i] #A[i] = 1. + x[i]**1.5+2.2*(x[i] - 1.5)**2 - 0.5*(x[i]-1.5)**4 A_min = A.min() #This is the minimum cross section area, the throat of the nozzle for i in xrange(0,nx+1,1): if (A[i] == A_min): j_min = i #print "j_mi = " + str(j_min) class A_of_Mach(object): """Calculate cross section as a function of Mach number, M. Calculate derivative with respect to Mach number, M. This function is defined as a class because we need to pass cross sectional area as a parameter. """ def __init__(self,A=1.): self.A = A #This is how I can change the parameter def __call__(self,M=1.): M_sq = M**2 gexp = (gamma + 1.)/(gamma - 1.) paren = 1. + (gamma - 1.)/2.*M_sq f_of_M = (paren**gexp)/M_sq - ((gamma + 1.)/2.)**gexp*self.A**2 df = -2.*paren**gexp/M_sq/M + (gamma +1.)/M*paren**(2./(gamma -1.)) #print f_of_M, df return f_of_M, df def rtsafe(func,x1,x2,xacc): """Safe version of a Newton-Raphson algorithm that combines this with bisection. Call with function func to be solved (we want to solve f(x) = 0 for x)) x1 and x2 are values that bracket the solution and must be passed by user to initialize xacc gives the accuracy that should be attempted. Don't be too greedy! """ maxit = 100 fl, df = func(x1) fh, df = func(x2) #print 333, x1, fl, x2, fh, fl*fh if (fl*fh > 0): print "error with guess_lo = ", x1, ", f= ", fl, " guess_hi = ", x2, ", f= ",fh return "This will not work as initial guesses do not bracket solution" if (fl < 0): xl = x1 xh = x2 else: xh = x1 xl = x2 rtsafe = 0.5*(xl+xh) dxold = fabs(xh-xl) dx = dxold f, df = func(rtsafe) for i in xrange(0,maxit,1): if ((((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) > 0.0) or (fabs(2.0*f) > fabs(dxold*df))): """use bisection""" dxold = dx dx = 0.5*(xh - xl) rtsafe = xl+dx if (xl == rtsafe): return rtsafe else: """use Newton-Raphson step""" dxold=dx dx=f/df temp=rtsafe rtsafe -= dx if (temp == rtsafe): return rtsafe #no appreciable change ===> convergence if (fabs(dx) < xacc): return rtsafe #convergence criterion met f, df = func(rtsafe) if (f < 0.): xl = rtsafe else: xh = rtsafe print "maximum iterations reached in rtsafe! value for x will not be good" return rtsafe #Calculate analytic solution #declare arrays M_an = zeros(nx+1) M_an_sq = zeros(nx+1) T_an = M_an.copy() v_an = M_an.copy() p_an = M_an.copy() rho_an = M_an.copy() j = 0 #determine analytic solutions for gamma = 5/3 if (gamma == 5./3.): p = [1., 12., 54., 0, 81.] #polynomial of order 4 has four (complex) solutions for i in xrange(0,nx+1,1): p[3]=(108. - 256.*(A[i]/A_min)**2) sol= roots(p) sol= real_if_close(sol)#this ensures that we dont have numerically insignificant imginary parts in the two real solutions #Make sure we use the right solution. Need the larger one after the throat if A[i] == A_min: j = 1 M_an_sq[i] = real(sol[3-j]) M_an[i] = sqrt(M_an_sq[i]) check = 1./M_an_sq[i]*(0.75*(1.+1./3.*M_an_sq[i]))**4 #print i + 1, j, A[i], sqrt(check), M_an[i] rho_an[i] = (1.+1./3.*M_an_sq[i])**(-3./2.) T_an[i] = 1./(1.+1./3.*M_an_sq[i]) p_an[i] = (1.+1./3.*M_an_sq[i])**(-5./2.) else: AM = A_of_Mach() guess_lo = 0.01#0.01*rtsafe(AM,1.e-3,1.,1.e-4) guess_hi = 0.99#200.*guess_lo #print guess_lo, guess_hi for i in xrange(0,nx+1,1): AM.A = A[i] agl = AM(guess_lo)[0] agh = AM(guess_hi)[0] k = 0 while ((agl*agh > 0) and (k < 10)): guess_hi = 1.5*guess_hi guess_lo = 0.95*guess_lo agl = AM(guess_lo)[0] agh = AM(guess_hi)[0] k = k+1 #print 11111, guess_lo, agl, guess_hi, agh #print i, A[i], guess_lo, agl, guess_hi, agh #print i, A[i], rtsafe(AM,guess_lo,guess_hi,1.e-4) if (A[i] == A_min): M_an[i] = 1. else: M_an[i] = rtsafe(AM,guess_lo,guess_hi,1.e-4) guess_lo = M_an[i] guess_hi = 1.5*guess_lo #print M_an[i] if (i> 2): guess_hi = M_an[i] + 1.5*(M_an[i]-M_an[i-1]) agl = AM(guess_lo)[0] agh = AM(guess_hi)[0] #print guess_lo, agl, guess_hi, agh M_an_sq[i] = M_an[i]**2 p_an[i] = (1.+(gamma-1.)/2.*M_an_sq[i])**(-gamma/(gamma - 1.)) rho_an[i] = (1.+(gamma-1.)/2.*M_an_sq[i])**(-1./(gamma - 1.)) T_an[i] = (1.+(gamma-1.)/2.*M_an_sq[i])**-1. """mass flux must be constant in stationary solution. Therefore, can use v_an[0] = M_an[0]*sqrt(gamma*T_an[0]) as the value for the analytic solution of v(x)""" v_an[0:nx+1] = M_an[0:nx+1]*sqrt(T_an[0:nx+1]) mf_an = rho_an.copy() mf_an[0:nx+1] = rho_an[0:nx+1]*v_an[0:nx+1]*A[0:nx+1] a = Gnuplot.Gnuplot() #depending on your system, you may need to insert persist=1 in the parenthesis data_rho_an = Gnuplot.Data(x,rho_an,using=(1,2)) data_p_an = Gnuplot.Data(x,p_an,using=(1,2)) data_T_an = Gnuplot.Data(x,T_an,using=(1,2)) data_M_an = Gnuplot.Data(x,M_an,using=(1,2)) data_MF_an = Gnuplot.Data(x,mf_an,using=(1,2)) a.plot(data_M_an,data_rho_an,data_p_an,data_T_an, data_MF_an) def solver(user_action=None): t0 = time.clock() #start timer to find out how much CPU time is used in solver """Define arrays etc.""" rho = zeros(nx+1) v = rho.copy() T = v.copy() dta = v.copy() cs = v.copy() drho_dt = zeros(nx+1) dv_dt = drho_dt.copy() dT_dt = dv_dt.copy() rho_p = zeros(nx+1) v_p = rho_p.copy() T_p = v_p.copy() drho_p_dt = zeros(nx+1) dv_p_dt = drho_p_dt.copy() dT_p_dt = dv_dt.copy() drho_dt_ave = zeros(nx+1) dv_dt_ave = drho_p_dt.copy() dT_dt_ave = dv_dt.copy() rho_pp = zeros(nx+1) v_pp = rho_pp.copy() T_pp = v_pp.copy() mass_flow = v.copy() """Set up initial conditions""" # for i in xrange(0,nx+1,1): # rho[i] = 1. - 0.314*x[i] # T[i] = 1. - 0.231*x[i] # cs[i] = sqrt(T[i]) # v[i] = (0.1 + 1.09*x[i])*cs[i] #vectorized version for initial conditions: rho[0:nx+1] = 1. - 0.3146*x[0:nx+1] T[0:nx+1] = 1. - 0.2314*x[0:nx+1] cs[0:nx+1] = sqrt(T[0:nx+1]) v[0:nx+1] = (0.1 + 1.09*x[0:nx+1])*cs[0:nx+1] #Could also use quite different initial conditions... This does not always work... #rho[0:nx+1] = 0. + exp(-x[0:nx+1]/5.5) #T[0:nx+1] = 0. + exp(-x[0:nx+1]/5.5) #cs[0:nx+1] = sqrt(T[0:nx+1]) #v[0:nx+1] = 0. + exp(-x[0:nx+1]/5.5) #print v """Set up boundary conditions""" rho[0] = 1. T[0] = 1. v[0] = 2.*v[1] - v[2] #one floating boundary condition at inflow rho[nx] = 2.*rho[nx-1] - rho[nx-2] v[nx] = 2.*v[nx-1] - v[nx-2] T[nx] = 2.*T[nx-1] - T[nx-2] #print 'initial values at i = 16: ', x[15], A[15], rho[15], v[15], T[15] #print 'initial values at i = 17: ', x[16], A[16], rho[16], v[16], T[16] t = 0. dres = 1. oldres = 1. while dres > res_acc: """compute time step""" for i in xrange(0,nx+1,1): dta[i] = fabs(C*dx/(cs[i]+v[i])) #print '*** ', t, dta[i], dx, (cs[i] + v[i]), dx/(cs[i]+v[i]) dt = dta.min() #print 'final dt is: ', dta, '\n******', dt t_old = t; t += dt * t_scale """predictor step""" #prepare derivatives using forward differences for i in xrange(0,nx,1): drho_dt[i] = -rho[i]*(v[i+1] - v[i])/dx - rho[i]*v[i]*(log(A[i+1])-log(A[i]))/dx - v[i]*(rho[i+1]-rho[i])/dx dv_dt[i] = -v[i]*(v[i+1]-v[i])/dx - ((T[i+1]-T[i])/dx + T[i]/rho[i]*(rho[i+1]-rho[i])/dx)/gamma dT_dt[i] = -v[i]*(T[i+1]-T[i])/dx - (gamma -1)*T[i]*((v[i+1]-v[i])/dx + v[i]*(log(A[i+1]) - log(A[i]))/dx) #Now calculate predictor values rho_p[i] = rho[i] + drho_dt[i]*dt v_p[i] = v[i] + dv_dt[i]*dt T_p[i] = T[i] + dT_dt[i]*dt #print 'x values at i = 16: ', rho_p[15], v_p[15], T_p[15] #print 'values at i = 17: ', rho[16], v[16], T[16] #print 'deriv values at i = 16: ', drho_dt[15], dv_dt[15], dT_dt[15] #print 'deriv values at i = 17: ', drho_dt[16], dv_dt[16], dT_dt[16] v_p[0] = 2.*v_p[1] - v_p[2] rho_p[0] = 1. T_p[0] = 1. rho_p[nx] = 2.*rho_p[nx-1] - rho_p[nx-2] T_p[nx] = 2.*T_p[nx-1] - T_p[nx-2] v_p[nx] = 2.*v_p[nx-1] - v_p[nx-2] #print rho_p #print v_p #print T_p """corrector step""" #prepare derivatives using rearward differences for i in xrange(1,nx+1,1): drho_p_dt[i] = -rho_p[i]*(v_p[i] - v_p[i-1])/dx - rho_p[i]*v_p[i]*(log(A[i])-log(A[i-1]))/dx - v_p[i]*(rho_p[i]-rho_p[i-1])/dx dv_p_dt[i] = -v_p[i]*(v_p[i]-v_p[i-1])/dx - ((T_p[i]-T_p[i-1])/dx + T_p[i]/rho_p[i]*(rho_p[i]-rho_p[i-1])/dx)/gamma dT_p_dt[i] = -v_p[i]*(T_p[i]-T_p[i-1])/dx - (gamma -1)*T_p[i]*((v_p[i]-v_p[i-1])/dx + v_p[i]*(log(A[i]) - log(A[i-1]))/dx) #print 'dx_p_dt values at i = 16: ', drho_p_dt[15], dv_p_dt[15], dT_p_dt[15] """time propagator""" for i in xrange(0,nx+1,1): drho_dt_ave[i] = 0.5*(drho_dt[i] + drho_p_dt[i]) dv_dt_ave[i] = 0.5*(dv_dt[i] + dv_p_dt[i]) dT_dt_ave[i] = 0.5*(dT_dt[i] + dT_p_dt[i]) rho_pp[i] = rho[i] + drho_dt_ave[i]*dt v_pp[i] = v[i] + dv_dt_ave[i]*dt T_pp[i] = T[i] + dT_dt_ave[i]*dt #print 'dx_dt_ave values at i = 16: ', drho_dt_ave[15], dv_dt_ave[15], dT_dt_ave[15] #print 'final propagated values at i = 16 are: ', rho_pp[15], v_pp[15], T_pp[15] """insert boundary conditions""" rho[0] = 1. T[0] = 1. cs[0] = 1. #print 'at lower boundary: ', v_pp[1], v_pp[2] #print 'at upper boundary: ', rho_pp[nx-1], rho_pp[nx-2] #print 'at upper boundary: ', v_pp[nx-1], v_pp[nx-2] #print 'at upper boundary: ', T_pp[nx-1], T_pp[nx-2] v_pp[0] = 2.*v_pp[1] - v_pp[2] #one floating boundary condition at inflow rho_pp[nx] = 2.*rho_pp[nx-1] - rho_pp[nx-2] v_pp[nx] = 2.*v_pp[nx-1] - v_pp[nx-2] T_pp[nx] = 2.*T_pp[nx-1] - T_pp[nx-2] cs[nx] = sqrt(T[nx]) #print 'boundary conditions are: at 0: ',v[0], ' at 30: ', rho[30], v[30], T[30] """The above values are the new ones (at the next time step)""" rho, v, T = rho_pp, v_pp, T_pp #update for i in xrange(0,nx+1,1): cs[i] = sqrt(T[i]) """Determine residuals""" rho_res = 0. v_res = 0. T_res = 0. for i in xrange(0,nx+1,1): rho_res += drho_dt_ave[i]**2 v_res += dv_dt_ave[i]**2 T_res += dT_dt_ave[i]**2 #print 'residuals: ', rho_res, v_res, T_res dres = fabs(oldres - rho_res) oldres = rho_res if user_action is not None: user_action(rho, v, T, rho_res, v_res, T_res, dres, x, t) #I can do whatever I'd like here t1 = time.clock() return t1-t0 , t #solver does NOT return solution u! #But it returns how long it worked (t1-t0) and the current model time # ---------------- end of solver ------------------------------------------- def test_solver(N): """ Very simple test case. Store the solution at every N time levels. Measure how long the solver actually works """ g = Gnuplot.Gnuplot() #depending on your system, you may need to insert persist=1 in the parenthesis rho_solutions = [] v_solutions = [] T_solutions = [] timesteps = [] rho_residuals = [] v_residuals = [] T_residuals = [] # Need time_level_counter as global variable since # it is assigned in the action function (that makes # a variable local to that block otherwise). global time_level_counter time_level_counter = 0 def action(rho, v, T, rho_res, v_res, T_res, dres, x, t): #this is where I can access solution u global time_level_counter if time_level_counter % N == 0: #only do this every N times timesteps.append(t) rho_residuals.append(rho_res) v_residuals.append(v_res) T_residuals.append(T_res) rho_solutions.append(rho.copy()) v_solutions.append(v.copy()) T_solutions.append(T.copy()) p = zeros(nx+1) p[0:nx+1] = rho[0:nx+1]*T[0:nx+1] M = zeros(nx+1) """calculate mass flow""" mass_flow = M.copy() #mass_flow[0:nx+1] = rho[0:nx+1]*v[0:nx+1]*A[0:nx+1] #mass flow itself mass_flow[0:nx+1] = rho[0:nx+1]*v[0:nx+1]*A[0:nx+1]/mf_an[0:nx+1]*2.5 #mass flow / analytic * 2.5 M[0:nx+1] = v[0:nx+1]/sqrt(T[0:nx+1]) g.reset() g('set ylabel "y-axis [arb. units]"') #this is how you access gnupot commands g('set xlabel "x-axis [arb. units]"') string = 'set title "time: ' + str(t) +'"' #g('set title "blabla"') g(string) data_rho = Gnuplot.Data(x,rho,using=(1,2)) data_v = Gnuplot.Data(x,v,using=(1,2)) data_T = Gnuplot.Data(x,T,using=(1,2)) data_p = Gnuplot.Data(x,p,using=(1,2)) data_M = Gnuplot.Data(x,M,using=(1,2)) data_rho_an = Gnuplot.Data(x,rho_an,using=(1,2), with_ = 'line', title = 'rho analytic') data_p_an = Gnuplot.Data(x,p_an,using=(1,2), with_ = 'line', title = 'p analytic') data_T_an = Gnuplot.Data(x,T_an,using=(1,2), with_ = 'line', title = 'T analytic') data_M_an = Gnuplot.Data(x,M_an,using=(1,2), with_ = 'line', title = 'M analytic') data_MF = Gnuplot.Data(x,mass_flow,using=(1,2)) data_MF_an= Gnuplot.Data(x,2.5*ones(nx+1),using=(1,2), with_ = 'line') data_v_an= Gnuplot.Data(x,v_an,using=(1,2), with_ = 'line', title = 'v analytic') #data_MF_an= Gnuplot.Data(x,2.5*mf_an,using=(1,2), with = 'line') #data_MF_an= Gnuplot.Data(x,mf_an,using=(1,2), with = 'line') g.plot(data_rho,data_v,data_p,data_T,data_M,data_rho_an,data_p_an,data_T_an, data_M_an, data_MF, data_MF_an, data_v_an) #g.plot(data_rho,data_v,data_T) #z = time.sleep(0.1) #wait a little so you can look at it #print rho_res, dres time_level_counter += 1 #call the solver here, note that the solution, u, is not returned #but cpu time is cpu, t = solver(user_action=action) print 'CPU time:', cpu, ' sec' print 'time: ', t return timesteps, rho_solutions, v_solutions, T_solutions, rho_residuals, v_residuals, T_residuals #test_solver has access to solutions through 'action' #and returns solution array to main # --------------------------------- end test_solver ----------------------------- def main(): # open output file #output = open("1dwave-sol.dat","w") timesteps, rho_solutions, v_solutions, T_solutions, rho_residuals, v_residuals, T_residuals = test_solver(2) #this is where I call the solver #print solutions #output.write(str(solutions)) #note: the output is an array #output.close() raw_input('Please press the return key to continue...\n') main()