missing setrun and setplot for corner problem

fvmbook/chap22/corner/setplot.py

@@ -0,0 +1,125 @@
+
+""" 
+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.
+    
+""" 
+
+#--------------------------
+def setplot(plotdata):
+#--------------------------
+
+    """ 
+    Specify what is to be plotted at each frame.
+    Input:  plotdata, an instance of clawpack.visclaw.data.ClawPlotData.
+    Output: a modified version of plotdata.
+    
+    """ 
+
+
+    from clawpack.visclaw import colormaps
+    from numpy import linspace
+
+    plotdata.clearfigures()  # clear any old figures,axes,items data
+
+    def plot_corner(current_data):
+        from pylab import plot
+        plot([.0,.0],[-1,0],'r',linewidth=2)
+        plot([.0,1],[0,.55],'r',linewidth=2)
+
+    def sigmatr(current_data):
+        # trace of sigma
+        q = current_data.q
+        return q[0,:,:] + q[1,:,:]
+
+
+
+    # Figure for trace of sigma
+    plotfigure = plotdata.new_plotfigure(name='trace(sigma)', figno=0)
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = 'auto'
+    plotaxes.ylimits = 'auto'
+    plotaxes.title = 'trace(sigma)'
+    plotaxes.scaled = True
+    plotaxes.afteraxes = plot_corner
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
+    plotitem.plot_var = sigmatr
+    plotitem.pcolor_cmap = colormaps.red_yellow_blue
+    plotitem.pcolor_cmin = -1.
+    plotitem.pcolor_cmax = 1.
+    plotitem.add_colorbar = True
+    plotitem.show = True       # show on plot?
+
+
+
+    # Figure for shear stress
+    plotfigure = plotdata.new_plotfigure(name='shear', figno=1)
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = 'auto'
+    plotaxes.ylimits = 'auto'
+    plotaxes.title = 'shear stress'
+    plotaxes.scaled = True
+    plotaxes.afteraxes = plot_corner
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
+    plotitem.plot_var = 2  # sigma_12
+    plotitem.pcolor_cmap = colormaps.red_yellow_blue
+    plotitem.pcolor_cmin = -0.2
+    plotitem.pcolor_cmax = 0.2
+    plotitem.add_colorbar = True
+    plotitem.show = True       # show on plot?
+
+
+    # Figure for contours 
+    plotfigure = plotdata.new_plotfigure(name='contours', figno=2)
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = 'auto'
+    plotaxes.ylimits = 'auto'
+    plotaxes.title = 'pressure(black) and shear(green)'
+    plotaxes.scaled = True
+    plotaxes.afteraxes = plot_corner
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_contour')
+    plotitem.plot_var = sigmatr
+    plotitem.contour_levels = linspace(-2,8,50)
+    plotitem.contour_colors = 'k'
+    plotitem.show = True       # show on plot?
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_contour')
+    plotitem.plot_var = 2  # sigma_12
+    plotitem.contour_levels = linspace(-0.4,0.4,30)
+    plotitem.contour_colors = 'g'
+    plotitem.show = True       # show on plot?
+
+
+
+    # Parameters used only when creating html and/or latex hardcopy
+    # e.g., via clawpack.visclaw.frametools.printframes:
+
+    plotdata.printfigs = True                # print figures
+    plotdata.print_format = 'png'            # file format
+    plotdata.print_framenos = 'all'          # list of frames to print
+    plotdata.print_fignos = 'all'            # list of figures to print
+    plotdata.html = True                     # create html files of plots?
+    plotdata.html_homelink = '../README.html'   # pointer for top of index
+    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
+
+

fvmbook/chap22/corner/setrun.py

@@ -0,0 +1,243 @@
+""" 
+Module to set up run time parameters for Clawpack -- classic code.
+
+The values set in the function setrun are then written out to data files
+that will be read in by the Fortran code.
+    
+""" 
+
+import os
+import numpy as np
+
+#------------------------------
+def setrun(claw_pkg='classic'):
+#------------------------------
+
+    """ 
+    Define the parameters used for running Clawpack.
+
+    INPUT:
+        claw_pkg expected to be "classic" for this setrun.
+
+    OUTPUT:
+        rundata - object of class ClawRunData 
+    
+    """ 
+
+    from clawpack.clawutil import data 
+
+
+    assert claw_pkg.lower() == 'classic',  "Expected claw_pkg = 'classic'"
+
+    num_dim = 2
+    rundata = data.ClawRunData(claw_pkg, num_dim)
+
+    #------------------------------------------------------------------
+    # Problem-specific parameters to be written to setprob.data:
+    #------------------------------------------------------------------
+    # Sample setup to write one line to setprob.data ...
+    probdata = rundata.new_UserData(name='probdata',fname='setprob.data')
+    probdata.add_param('rho1',     1.,  'rho1')
+    probdata.add_param('lam1',     4.,  'lam1')
+    probdata.add_param('mu1',     0.5,  'mu1')
+    probdata.add_param('rho2',     1.,  'rho2')
+    probdata.add_param('lam2',     2.,  'lam2')
+    probdata.add_param('mu2',     1.0,  'mu2')
+
+    #------------------------------------------------------------------
+    # Standard Clawpack parameters to be written to claw.data:
+    #------------------------------------------------------------------
+
+    clawdata = rundata.clawdata  # initialized when rundata instantiated
+
+
+    # ---------------
+    # Spatial domain:
+    # ---------------
+
+    # Number of space dimensions:
+    clawdata.num_dim = num_dim
+
+    # Lower and upper edge of computational domain:
+    clawdata.lower[0] = -1.000000e+00          # xlower
+    clawdata.upper[0] = 1.000000e+00          # xupper
+    clawdata.lower[1] = -1.000000e+00          # ylower
+    clawdata.upper[1] = 1.000000e+00          # yupper
+
+    # Number of grid cells:
+    clawdata.num_cells[0] = 200      # mx
+    clawdata.num_cells[1] = 200      # my
+
+
+    # ---------------
+    # Size of system:
+    # ---------------
+
+    # Number of equations in the system:
+    clawdata.num_eqn = 5
+
+    # Number of auxiliary variables in the aux array (initialized in setaux)
+    clawdata.num_aux = 5
+
+    # Index of aux array corresponding to capacity function, if there is one:
+    clawdata.capa_index = 0
+
+
+    # -------------
+    # Initial time:
+    # -------------
+
+    clawdata.t0 = 0.000000
+
+
+    # Restart from checkpoint file of a previous run?
+    # Note: If restarting, you must also change the Makefile to set:
+    #    RESTART = True
+    # If restarting, t0 above should be from original run, and the
+    # restart_file 'fort.qNNNN' specified below should be in 
+    # the OUTDIR indicated in Makefile.
+
+    clawdata.restart = False               # True to restart from prior results
+    clawdata.restart_file = 'fort.q0006'   # File to use for restart data
+
+
+    # -------------
+    # Output times:
+    #--------------
+
+    # Specify at what times the results should be written to fort.q files.
+    # Note that the time integration stops after the final output time.
+
+    clawdata.output_style = 1
+
+    if clawdata.output_style==1:
+        # Output ntimes frames at equally spaced times up to tfinal:
+        # Can specify num_output_times = 0 for no output
+        clawdata.num_output_times = 10
+        clawdata.tfinal = 0.500000
+        clawdata.output_t0 = True  # output at initial (or restart) time?
+
+    elif clawdata.output_style == 2:
+        # Specify a list or numpy array of output times:
+        # Include t0 if you want output at the initial time.
+        clawdata.output_times =  [0., 0.1]
+
+    elif clawdata.output_style == 3:
+        # Output every step_interval timesteps over total_steps timesteps:
+        clawdata.output_step_interval = 2
+        clawdata.total_steps = 4
+        clawdata.output_t0 = True  # output at initial (or restart) time?
+
+
+    clawdata.output_format = 'ascii'      # 'ascii', 'binary', 'netcdf'
+
+    clawdata.output_q_components = 'all'   # could be list such as [True,True]
+    clawdata.output_aux_components = 'none'  # could be list
+    clawdata.output_aux_onlyonce = True    # output aux arrays only at t0
+
+
+    # ---------------------------------------------------
+    # Verbosity of messages to screen during integration:  
+    # ---------------------------------------------------
+
+    # The current t, dt, and cfl will be printed every time step
+    # at AMR levels <= verbosity.  Set verbosity = 0 for no printing.
+    #   (E.g. verbosity == 2 means print only on levels 1 and 2.)
+    clawdata.verbosity = 1
+
+
+
+    # --------------
+    # Time stepping:
+    # --------------
+
+    # if dt_variable==True:  variable time steps used based on cfl_desired,
+    # if dt_variable==False: fixed time steps dt = dt_initial always used.
+    clawdata.dt_variable = True
+
+    # Initial time step for variable dt.  
+    # (If dt_variable==0 then dt=dt_initial for all steps)
+    clawdata.dt_initial = 1.000000e-02
+
+    # Max time step to be allowed if variable dt used:
+    clawdata.dt_max = 1.000000e+99
+
+    # Desired Courant number if variable dt used 
+    clawdata.cfl_desired = 0.900000
+    # max Courant number to allow without retaking step with a smaller dt:
+    clawdata.cfl_max = 1.000000
+
+    # Maximum number of time steps to allow between output times:
+    clawdata.steps_max = 500
+
+
+    # ------------------
+    # Method to be used:
+    # ------------------
+
+    # Order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
+    clawdata.order = 2
+
+    # Use dimensional splitting? (not yet available for AMR)
+    clawdata.dimensional_split = 'unsplit'
+
+    # For unsplit method, transverse_waves can be 
+    #  0 or 'none'      ==> donor cell (only normal solver used)
+    #  1 or 'increment' ==> corner transport of waves
+    #  2 or 'all'       ==> corner transport of 2nd order corrections too
+    clawdata.transverse_waves = 2
+
+
+    # Number of waves in the Riemann solution:
+    clawdata.num_waves = 4
+
+    # List of limiters to use for each wave family:  
+    # Required:  len(limiter) == num_waves
+    # Some options:
+    #   0 or 'none'     ==> no limiter (Lax-Wendroff)
+    #   1 or 'minmod'   ==> minmod
+    #   2 or 'superbee' ==> superbee
+    #   3 or 'vanleer'  ==> van Leer
+    #   4 or 'mc'       ==> MC limiter
+    clawdata.limiter = ['mc', 'mc', 'mc', 'mc']
+
+    clawdata.use_fwaves = False    # True ==> use f-wave version of algorithms
+
+    # Source terms splitting:
+    #   src_split == 0 or 'none'    ==> no source term (src routine never called)
+    #   src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, 
+    #   src_split == 2 or 'strang'  ==> Strang (2nd order) splitting used,  not recommended.
+    clawdata.source_split = 0
+
+
+    # --------------------
+    # Boundary conditions:
+    # --------------------
+
+    # Number of ghost cells (usually 2)
+    clawdata.num_ghost = 2
+
+    # Choice of BCs at xlower and xupper:
+    #   0 or 'user'     => user specified (must modify bcNamr.f to use this option)
+    #   1 or 'extrap'   => extrapolation (non-reflecting outflow)
+    #   2 or 'periodic' => periodic (must specify this at both boundaries)
+    #   3 or 'wall'     => solid wall for systems where q(2) is normal velocity
+
+    clawdata.bc_lower[0] = 'extrap'   # at xlower
+    clawdata.bc_upper[0] = 'extrap'   # at xupper
+
+    clawdata.bc_lower[1] = 'extrap'   # at ylower
+    clawdata.bc_upper[1] = 'extrap'   # at yupper
+
+    return rundata
+
+    # end of function setrun
+    # ----------------------
+
+
+if __name__ == '__main__':
+    # Set up run-time parameters and write all data files.
+    import sys
+    rundata = setrun(*sys.argv[1:])
+    rundata.write()
+