diff --git a/fvmbook/chap22/corner/setplot.py b/fvmbook/chap22/corner/setplot.py new file mode 100644 index 00000000..05870e1d --- /dev/null +++ b/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 + + diff --git a/fvmbook/chap22/corner/setrun.py b/fvmbook/chap22/corner/setrun.py new file mode 100644 index 00000000..a4998ef4 --- /dev/null +++ b/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() +