setplot2.py.html | |
Source file: setplot2.py | |
Directory: /var/www/html/clawpack/links/shockvacuum10/shockvacuum | |
Converted: Mon Aug 16 2010 at 14:13:32 using clawcode2html | |
This documentation file will not reflect any later changes in the source file. |
""" Set up the plot figures, axes, and items to be done for each frame. This module is imported by the plotting routines and then the function setplot is called to set the plot parameters. This version requires Clawpack 4.4. $ svn co http://kingkong.amath.washington.edu/svn/claw4/trunk claw4 Or download the tarfile (see documentation). For documentation, see ... http://kingkong.amath.washington.edu/claw4/www/users/ """ import os import numpy as np import pylab from pyclaw.data import Data import similarity_solution as ss reload(ss) # in case it has been changed by user # Global parameters used in functions below: #------------------------------------------- # Specify where the data comes from for plotting. # outdir should be set already! #outdir = '_output' # Read the values from setprob.data, which should reflect the values set # at the top of setrun.py. setprob_file = os.path.join(outdir, 'setprob.data') params = Data(setprob_file) ieos = params.ieos # = 1 for Euler, 2 for polytropic n = params.npoly nstar = params.nstar # = alpha K = params.sK grav = params.grav gamma = 1. + 1./n gamma_star = 1. + 1./nstar tfontsize = 20 # fontsize for some titles plot_simsol = True # Set to False to suppress similarity solution plots. if plot_simsol & (ieos==1): if (n==1) & (nstar==1.): fname_simsol = 'simsol_euler_n1.dat' elif (n==1.1) & (nstar==1.1): fname_simsol = 'simsol_euler_n1.1.dat' elif (n==1.5) & (nstar==1.5): fname_simsol = 'simsol_euler_n3o2.dat' #fname_simsol = 'simsol_euler_n1.5_alpha1.5.dat' elif (n==1.) & (nstar==1.5): fname_simsol = 'simsol_euler_n1_alpha1.5.dat' elif (n==1.5) & (nstar==3.): fname_simsol = 'simsol_euler_n3o2_alpha3.dat' elif n==3: fname_simsol = 'simsol_euler_n3.dat' else: print '*** Cannot find a data file for similarity solution' plot_simsol = False print "Assuming computation was done with ieos = %i, n = %g" % (ieos,n) if ieos==1: print " and nstar = ",nstar #------------------------------------------ # Read in t and xs from fort.xs try: fname = outdir+'/fort.xs' txs_array = np.loadtxt(fname, dtype=float) except: print "*** Error reading ",fname txs_array = None print "exists? ",os.path.isfile(fname) raise if plot_simsol: if ieos==2: # compute the similarity solution for polytropic: #vstar = 4.8905 # for n=2, from numerics with order=2, mx=2400 vstar = 3.7418 # for n=1, mx=1200 vstar = 4.89 # for n=1. mx=4800 vstar = 4.5 # for n=3. mx=4800 vstar = 4.89 # for n=1. mx=2400 tfinal = 2.59125913623e-07 # for polytropic npoly=1, mx=2400 tfinal = 2.4477500128e-07 # for polytropic npoly=1, mx=4800 tfinal = 2.641111388e-07 # for polytropic npoly=3, mx=4800 tfinal = 2.4967677e-07 # for polytropic npoly=1, mx=2400 tfinal = 2.447750014889652e-07 # for polytropic npoly=1, mx=4800 # xreset = -1e-4: tfinal = 4.0066713534085611e-06 # for polytropic npoly=1, mx=4800, tfinal, vstar = ss.estimate_params_polytropic(txs_array) print 'tfinal = ',tfinal #vstar = 4.54 # for n=1.5 mx=2400 #vstar = 4.16 # for n=3.0 mx=1200 #vstar = 4.89 # for n=1.0 mx=1200 #vstar = 4.92 # for n=1.0 mx=1200 # Reset vstar based on Riemann invariant to better match... vstar = 4.893 print 'Reset vstar = ',vstar pylab.figure(23) pylab.clf() xs = txs_array[1:,1] t_xs = txs_array[1:,0] - tfinal pylab.loglog(-t_xs, -xs/vstar, 'b', linewidth=3) pylab.title('Computed shock location') pylab.xlabel(r'$t_f-t_c$', fontsize=tfontsize) pylab.ylabel(r'$-x_s/v_*$', fontsize=tfontsize) pylab.axis('scaled') pylab.xticks(np.logspace(-15,-5,6), fontsize=tfontsize) pylab.yticks(np.logspace(-15,-5,6), fontsize=tfontsize) nypts = 1000 # number of points to evaluate sim. soln. ymax = 1000. print "Computing similarity solution using tfinal = ",tfinal simsol = ss.compute_simsol_polytropic(ymax, nypts, grav, n, vstar, tfinal) else: # Read in the tabulated similarity solution for Euler: print "Reading similarity solution from file ",fname_simsol #print " using tfinal = ",tfinal simsol = ss.read_tabular_simsol(fname_simsol) simsol.grav = grav tfinal, C_y = ss.estimate_params(simsol.beta, simsol.ys, txs_array) simsol.C_y = C_y simsol.tfinal = tfinal print "Estimated tfinal = %22.15e" % tfinal print "Estimated C_y = %22.15e" % C_y simsol.K = K simsol.gamma = gamma C_rho = (grav / ((nstar+1)*K))**nstar simsol.C_rho = C_rho print "C_rho = ",C_rho ftreset = os.path.join(outdir, 'fort.treset') if os.path.isfile(ftreset): treset = float(open(ftreset,'r').readline().split()[0]) else: treset = 0. #--------------------- def setplot(plotdata): #--------------------- """ Specify what is to be plotted at each frame. Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. Output: a modified version of plotdata. """ from pyclaw.plotters.data import ClawPlotData from pyclaw.plotters.frametools import var_limits if not isinstance(plotdata,ClawPlotData): print '*** Error in setplot: plotdata must be a ClawPlotData object' return plotdata.clearfigures() # clear any old figures,axes,items data plotdata.outdir = outdir # where to find solution plotdata.beforeframe = beforeframe_simsol #plotdata.afterframe = None plotdata.afterframe = afterframe_savefigs if plotdata.mode() == 'printframes': # call function defined below to set html and latex parameters: plotdata = setplot_print(plotdata) # limits of x-axis, set to 'auto' or set to an interval [a,b] xlim = 'auto' # determine min and max of each component of q over all frames # and use this to set limits on y-axis for plots so they # don't change from one frame to the next: # set of variables that will be plotted: vars = [2, velocity, soundspeed, vp2nc, vm2nc, h25, h30] use_global_limits = False if use_global_limits: # scan all frames to determine limits: varmin,varmax,varlim = var_limits(plotdata,vars,padding=0.1) # v+2nc and v-2nc appear on same plot, and we always have v-2nc < v+2nc: vpm2nc_lim = [varmin[vm2nc], varmax[vp2nc]] # limits for powers of h: hpowers_lim = 'auto' else: # choose all axes limits automatically (may change each frame): varlim = {} for var in vars: varlim[var] = 'auto' vpm2nc_lim = 'auto' hpowers_lim = 'auto' plotstyle = 'o' color = 'b' def add_xshock(current_data): import pylab xshock = current_data.user.xshock if xshock is not None: pylab.xlabel("shock location = %8.2e" % xshock, fontsize=15) #=== Figure 1 ==== title = 'rho' plotfigure = plotdata.new_plotfigure(name=title, figno=1) #plotfigure.kwargs = {'figsize':(12,6)} plotfigure.kwargs = {'figsize':(8,4)} plotfigure.show = True plotaxes = plotfigure.new_plotaxes(name=title) plotaxes.title = title plotaxes.xlimits = xlim plotaxes.ylimits = 'auto' if plot_simsol: plotaxes.afteraxes = plot_rho_simsol # add plot of similarity soln plotitem = plotaxes.new_plotitem('rho','1d_plot') plotitem.plot_var = 0 plotitem.plotstyle = 'b-' plotitem.kwargs = {'linewidth': 3} #=== Figure 2 ==== title = 'velocity' plotfigure = plotdata.new_plotfigure(name=title, figno=2) plotfigure.kwargs = {'figsize':(8,4)} plotfigure.show = True plotaxes = plotfigure.new_plotaxes(name=title) plotaxes.title = title plotaxes.xlimits = xlim plotaxes.ylimits = 'auto' plotitem = plotaxes.new_plotitem('velocity','1d_plot') plotitem.plot_var = velocity plotitem.plotstyle = 'b-' plotitem.kwargs = {'linewidth': 3} if plot_simsol: plotaxes.afteraxes = plot_v_simsol # add plot of similarity soln #plotitem = plotaxes.new_plotitem('polytropic','1d_plot') #plotitem.outdir = '_output_npoly1' #plotitem.plot_var = velocity #plotitem.plotstyle = 'g-' #=== Figure 22 ==== title = 'v and c' plotfigure = plotdata.new_plotfigure(name=title, figno=22) plotfigure.show = False plotfigure.kwargs = {'figsize':(8,4)} plotaxes = plotfigure.new_plotaxes(name=title) plotaxes.title = title plotaxes.xlimits = xlim plotaxes.ylimits = [-1., 5.] def add_legend(current_data): import pylab tsim = current_data.user.tsim pylab.legend(('sound speed','velocity')) pylab.title("v and c at time t = %14.8e" % tsim, fontsize=tfontsize) plotaxes.afteraxes = add_legend plotitem = plotaxes.new_plotitem('soundspeed','1d_plot') plotitem.plot_var = soundspeed plotitem.plotstyle = 'b-' plotitem.kwargs = {'linewidth': 3} plotitem = plotaxes.new_plotitem('velocity','1d_plot') plotitem.plot_var = velocity plotitem.plotstyle = 'k-' plotitem.kwargs = {'linewidth': 2} #=== Figure 3 ==== title = 'Sound speed' plotfigure = plotdata.new_plotfigure(name=title, figno=3) plotfigure.kwargs = {'figsize':(8,4)} plotfigure.show = True plotaxes = plotfigure.new_plotaxes(name=title) plotaxes.title = title plotaxes.xlimits = xlim plotaxes.ylimits = 'auto' plotitem = plotaxes.new_plotitem('soundspeed','1d_plot') plotitem.plot_var = soundspeed plotitem.plotstyle = 'b-' plotitem.kwargs = {'linewidth': 3} if plot_simsol: plotaxes.afteraxes = plot_c_simsol # add plot of similarity soln #=== Figure 4 ==== title = 'internal energy per unit mass' plotfigure = plotdata.new_plotfigure(name=title, figno=4) plotfigure.kwargs = {'figsize':(8,4)} plotfigure.show = True plotaxes = plotfigure.new_plotaxes(name=title) plotaxes.title = title plotaxes.xlimits = xlim plotaxes.ylimits = 'auto' plotitem = plotaxes.new_plotitem('energy','1d_plot') plotitem.plot_var = energy plotitem.plotstyle = 'bo-' #=== Figure 5 ==== title = 'Riemann invariants' plotfigure = plotdata.new_plotfigure(name=title, figno=5) plotfigure.kwargs = {'figsize':(8,4)} plotfigure.show = True plotaxes = plotfigure.new_plotaxes(name=title) plotaxes.title = title plotaxes.xlimits = xlim plotaxes.ylimits = vpm2nc_lim if plot_simsol: plotaxes.afteraxes = plot_vpm2nc plotitem = plotaxes.new_plotitem('vp2nc','1d_plot') plotitem.plot_var = vp2nc # function defined below plotitem.plotstyle = 'o' plotitem.color = 'g' # green for all grids plotitem = plotaxes.new_plotitem('vm2nc','1d_plot') # needs distinct name plotitem.plot_var = vm2nc # function defined below plotitem.plotstyle = 'o' plotitem.color = 'r' # red for all grids #=== Figure 8 ==== title = "v vs. y" plotfigure = plotdata.new_plotfigure(name='v vs. y', figno=8) plotfigure.kwargs = {'figsize':(8,4)} plotfigure.show = False plotaxes = plotfigure.new_plotaxes(name=title, type='empty') plotaxes.title = title # plot is done in afteraxes, plotting against y instead of x: if plot_simsol: plotaxes.afteraxes = plot_v_vs_y #=== Figure 9 ==== title = "c vs. y" plotfigure = plotdata.new_plotfigure(name='c vs. y', figno=9) plotfigure.kwargs = {'figsize':(8,4)} plotfigure.show = False plotaxes = plotfigure.new_plotaxes(name=title, type='empty') plotaxes.title = title # plot is done in afteraxes, plotting against y instead of x: if plot_simsol: plotaxes.afteraxes = plot_c_vs_y return plotdata #------------------------------------------------------ def afterframe_savefigs(current_data): """ Save the figures needed for the paper... may need to change the filenames below to be appropriate for this run. """ from pylab import figure, savefig n = current_data.frameno if ieos==1: nlist = [30, 40, 50] # for Euler fname = 'euler_n%3.1f_nstar%3.1f_' % (simsol.npoly,simsol.nstar) else: nlist = [0, 1, 2, 25, 30, 35, 40, 45, 50] # for polytropic fname = 'poly_n%3.1f_' % simsol.npoly if n in nlist: figure(1) # rho savefig(fname + 'rho%s.eps' % n) figure(2) # velocity savefig(fname + 'v%s.eps' % n) figure(5) # sound speed savefig(fname + 'c%s.eps' % n) if n < 10: figure(22) # v and c savefig(fname + 'vc%s.eps' % n) print 'Saved figures for frame ',n return #------------------------------------------------------ # Define the functions to be plotted. # q is the full q array (all columns, i.e. what was called data in Matlab) # on a single grid. This is called from the plotting routine for # each grid in an AMR calculation. # Note that indexing starts at 0 in Python so e.g. q[:,0]= density #------------------------------------------------------ def surface(cd): grid = cd.grid x = cd.x B = x # need to fix for general polytropic! surface = cd.q[:,0] + B return surface #------------------------------------------------------ def energy(cd): from numpy import where h = cd.q[:,0] hu = cd.q[:,1] E = cd.q[:,2] energy = where(h>0, (E-0.5*hu**2/h)/h, 0) return energy #------------------------------------------------------ def velocity(cd): from numpy import where h = cd.q[:,0] hu = cd.q[:,1] velocity = where(h>0, hu/h, 0) # h = 0 ==> v = 0 return velocity #------------------------------------------------------ def soundspeed(cd): from numpy import sqrt, where if ieos==2: rho = cd.q[:,0] #c = sqrt(grav*rho) c = sqrt(K*gamma*rho**(gamma-1)) else: e = energy(cd) c = sqrt(gamma*(gamma-1)*e) return c #------------------------------------------------------ def vp2nc(cd): c = soundspeed(cd) v = velocity(cd) vp2nc = v + 2.*n*c return vp2nc #------------------------------------------------------ def vm2nc(cd): c = soundspeed(cd) v = velocity(cd) vm2nc = v - 2.*n*c return vm2nc #------------------------------------------------------ def vmc(cd): c = soundspeed(cd) v = velocity(cd) vmc = v - c return vmc #------------------------------------------------------ def h25(cd): h25 = cd.q[:,0] ** 2.5 return h25 #------------------------------------------------------ def h30(cd): h30 = cd.q[:,0] ** 3.0 return h30 #------------------------------------------------------ def beforeframe_simsol(current_data): # Compute the x,c,u from the similarity solution at the current time. # for adding to plots of computed soln vs. x. # e.g. in plot_v_simsol. from numpy import hstack t = current_data.t # time from this frame if (t==0.) or (t>1e-5): t = t - treset if ieos==2: (x,c,v,rho) = ss.compute_xcv_polytropic(t,simsol) else: (x,c,v,rho) = ss.compute_xcv_euler(t,simsol) current_data.user.xshocksim = x[0] print "from simsol, xshock, cshock = ",x[0],c[0] if hasattr(simsol,'C_rho'): rhos = simsol.C_rho * (-x[0])**nstar else: print "*** simsol lacks attribute C_rho" rhos = 0. # need to fix print "+++ C_rho, x[0] = ",simsol.C_rho, x[0] print "+++ rhos = ",rhos #drho = (simsol.gamma - 1.) / (simsol.gamma + 1.) #drho2 = y[0]**(-nstar) * rhobar[0] #print "Delta rho = ", drho, drho2 x = hstack((0., x[0], x)) c = hstack((0., 0., c)) # not correct since c !=0 to right of shock, # but ok for plotting. v = hstack((0., 0., v)) rho = hstack((0., rhos, rho)) current_data.user.xsim = x current_data.user.csim = c current_data.user.vsim = v current_data.user.rhosim = rho tsim = t - tfinal # shift by treset has already been done current_data.user.tsim = tsim print "+++ t, treset, tfinal, tsim: ",t,treset,tfinal,tsim # compute xshock, current shock location: if (t>0.) and (t< 1e-5): xshock = np.interp(t, txs_array[:,0], txs_array[:,1]) current_data.user.xshock = xshock else: current_data.user.xshock = None print "xshock = ",current_data.user.xshock return current_data #------------------------------------------------------ def plot_v_vs_y(cd): # plot velocity vs. y instead of x # Similarity solution assumed to be precomputed in top of this module. from pylab import figure, plot,axes, hold, array, clf, isnan, title,xlim t = cd.t if ieos==2: ycell = ss.mapx2y_polytropic(cd.x, t, simsol) else: ycell = ss.mapx2y_euler(cd.x, t, simsol) v = velocity(cd) # from computed results clf() plot(ycell,v,'bo') title("velocity vs. y at t = %14.8f " % cd.t) # add plot of similarity solution: hold(True) ysim = simsol.y vsim = cd.user.vsim plot(ysim, vsim, 'r-', linewidth=2) xlim([ycell.min(),ycell.max()]) # -------------------------------------------------- def plot_c_vs_y(cd): # plot soundspeed vs. y instead of x # Similarity solution assumed to be precomputed in top of this module. from pylab import figure, plot,axes, hold, array, clf, isnan, title,xlim t = cd.t if ieos==2: ycell = ss.mapx2y_polytropic(cd.x, t, simsol) else: ycell = ss.mapx2y_euler(cd.x, t, simsol) c = soundspeed(cd) # from computed results clf() plot(ycell,c,'bo') title("soundspeed vs. y at t = %14.8f " % cd.t) # add plot of similarity solution: hold(True) ysim = simsol.y csim = cd.user.csim[2:] plot(ysim, csim, 'r-', linewidth=2) xlim([ycell.min(),ycell.max()]) # -------------------------------------------------- def plot_rho_simsol(current_data): # Plot rho(x,t) on top of density plot. # Similarity solution assumed to be precomputed in top of this module. # Also assumes x,u have been computed from simsol in beforeframe. from pylab import plot,xlim, legend,title # Only plot for times after t was reset: t = current_data.t tsim = current_data.user.tsim x = current_data.user.xsim rho = current_data.user.rhosim plot(x,rho,'r--', linewidth=2) xlim([current_data.xlower, 0.]) #rhoss = simsol.C_rho * (-x)**nstar #plot(x,0.5*rhoss,'g--', linewidth=2) if current_data.user.xshock is not None: xs = current_data.user.xshock xssim = current_data.user.xshocksim legend(('computed soln, xshock = %9.3e' % xs, \ 'similarity soln, xshock = %9.3e' % xssim),loc='upper right') title("density at time t = %14.8e" % tsim, fontsize=tfontsize) # -------------------------------------------------- def plot_v_simsol(current_data): # Plot v(x,t) on top of velocity plot. # Similarity solution assumed to be precomputed in top of this module. # Also assumes x,u have been computed from simsol in beforeframe. from pylab import plot,xlim, legend,title # Only plot for times after t was reset: t = current_data.t tsim = current_data.user.tsim x = current_data.user.xsim v = current_data.user.vsim plot(x,v,'r--', linewidth=2) xlim([current_data.xlower, 0.]) if current_data.user.xshock is not None: xs = current_data.user.xshock xssim = current_data.user.xshocksim legend(('computed soln, xshock = %9.3e' % xs, \ 'similarity soln, xshock = %9.3e' % xssim),loc='lower left') title("velocity at time t = %14.8e" % tsim, fontsize=tfontsize) #------------------------------------------------------ def plot_c_simsol(current_data): # Plot c(x,t) on top of velocity plot. # Similarity solution assumed to be precomputed in top of this module. # Also assumes x,c have been computed from simsol in beforeframe. from pylab import plot,xlim, legend,title, axes,ylim,xticks # Only plot for times after t was reset: t = current_data.t tsim = current_data.user.tsim x = current_data.user.xsim c = current_data.user.csim plot(x,c,'r--', linewidth=2) xlim([current_data.xlower, 0.]) if current_data.user.xshock is not None: xs = current_data.user.xshock xssim = current_data.user.xshocksim legend(('computed soln, xshock = %9.3e' % xs, \ 'similarity soln, xshock = %9.3e' % xssim),loc='lower left') title("sound speed at time t = %14.8e" % tsim, fontsize=tfontsize) if 0: ylims = ylim() # sub-plot with zoomed view: # set up for polytropic with n=1 a = axes([.2,.2,.3,.3]) c1 = current_data.var x1 = current_data.x plot(x1,c1,'b-', linewidth=2) plot(x,c,'r--', linewidth=2) xlim([1.05*xs, 0.95*xs]) ylim([0., 0.5*ylims[1]]) title("Zoom near shock") xticks([1.05*xs, xs, 0.95*xs]) #------------------------------------------------------ def plot_vpm2nc(current_data): # Plot v+2*n*c, v-2*n*c from similarity solution # Similarity solution assumed to be precomputed in top of this module. # Also assumes x,c,u have been computed from simsol in beforeframe. from pylab import plot,xlim n = simsol.npoly x = current_data.user.xsim c = current_data.user.csim v = current_data.user.vsim plot(x,v+2*n*c,'g--', linewidth=2) plot(x,v-2*n*c,'r--', linewidth=2) xlim([current_data.xlower, 0.]) return #------------------- def setplot_print(plotdata): #------------------- """ Parameters used only when creating html and/or latex hardcopy e.g., via pyclaw.plotters.frametools.printframes """ plotdata.printfigs = True # print figures plotdata.print_format = 'png' # file format plotdata.print_framenos = 'all' plotdata.print_fignos = 'all' plotdata.html = True # create html files of plots? plotdata.latex = True # create latex file of plots? plotdata.latex_figsperline = 2 # layout of plots plotdata.latex_framesperline = 1 # layout of plots plotdata.latex_makepdf = False # also run pdflatex? return plotdata