Commits

Aron Ahmadia committed e6cea32

figures 7.4 and 7.5

  • Participants
  • Parent commits cd07c1b

Comments (0)

Files changed (21)

 	@echo "Re-running WENO comparison of shock-bubble interaction"
 	(cd source/figure_7.3/experiment && time $(MAKE) data)
 	@cp $^ $(BUILD_DIR)/$@
+
+figure_7.4: source/figure_7.4/plots.py source/figure_7.4/README.txt
+	@mkdir -p $(BUILD_DIR)/$@
+	@cat source/figure_7.4/README.txt
+	@cp source/figure_7.4/README.txt $(BUILD_DIR)/$@
+	@echo "Building figure 7.4 from pyclaw output"
+	(cd source/figure_7.4 && $(MAKE) $@)
+
+figure_7.4/experiment: source/figure_7.4/experiment/psystem.py 
+	@mkdir -p $(BUILD_DIR)/$@	
+	@cat source/figure_7.4/experiment/README.txt
+	@cp source/figure_7.4/experiment/README.txt $(BUILD_DIR)/$@
+	@cp $^ $(BUILD_DIR)/$@
+
+figure_7.5: source/figure_7.5/plots.py source/figure_7.5/README.txt
+	@mkdir -p $(BUILD_DIR)/$@
+	@cat source/figure_7.5/README.txt
+	@cp source/figure_7.5/README.txt $(BUILD_DIR)/$@
+	@echo "Building figure 7.5 from pyclaw output"
+	(cd source/figure_7.5 && $(MAKE) $@)
+
+figure_7.5/experiment: source/figure_7.5/experiment/psystem.py 
+	@mkdir -p $(BUILD_DIR)/$@	
+	@cat source/figure_7.5/experiment/README.txt
+	@cp source/figure_7.5/experiment/README.txt $(BUILD_DIR)/$@
+	@cp $^ $(BUILD_DIR)/$@
+

source/figure_7.4/Makefile

+coarsen_data: experiment/_output/claw_p.pkl0000 experiment/_output/claw_p.ptc0000
+	python coarsen.py 
+
+compress_data: claw_p_coarsened.pkl0000 claw_p_coarsened.ptc0000
+	bzip2 claw_p_coarsened.pkl0000 claw_p_coarsened.ptc0000
+
+figure_7.4: experiment/_output/claw_p_coarsened.pkl0000.bz2 experiment/_output/claw_p_coarsened.ptc0000.bz2 plots.py
+	bunzip2 experiment/_output/claw_p_coarsened.pkl0000.bz2 experiment/_output/claw_p_coarsened.ptc0000.bz2 
+	python plots.py

source/figure_7.4/README.txt

+ABOUT THIS EXPERIMENT
+================
+Figure 7.4 shows the medium and initial condition (both shown zoomed-in) for the cylindrical solitary wave problem.
+
+The experimental output data is archived in compressed form (bzip2).  The uncompressed PETSc data uses 64-bit integer indices, to load it you will need to use a PETSc built with the following configuration option: --with-64-bit-indices=1
+
+RESPONSIBLE AUTHORS
+==============
+Manuel Quezada
+Aron Ahmadia

source/figure_7.4/coarsen.py

+from pyclaw.solution import Solution
+from petclaw.io.petsc import read_petsc
+from petclaw.io.petsc import write_petsc
+from petsc4py import PETSc
+import petclaw
+import numpy as np
+import os
+
+def read_original_solution(frame,file_prefix='claw_p',path='./experiment/_output/_p'):
+    sol = Solution()
+    read_petsc(sol,frame=frame,path=path,file_prefix=file_prefix)
+    return sol
+
+def coarsen_xhalf(sol):
+    # Interpolation
+    sol.state.q_da.setBoundaryType(PETSc.DA.BoundaryType.NONE)
+    sol.state.q_da.setRefinementFactor(refine_x=2,refine_y=2)
+    qda_coarsen = sol.state.q_da.coarsen()
+    gqVec_coarsen = qda_coarsen.createGlobalVec()
+    qda_coarsen.setInterpolationType(interp_type=0)
+    mat,scale = qda_coarsen.getInterpolation(sol.state.q_da)
+    mat.multTranspose(sol.state.gqVec,gqVec_coarsen)
+    gqVec_coarsen = scale*gqVec_coarsen
+    qda_coarsen.setBoundaryType(PETSc.DA.BoundaryType.PERIODIC)
+
+    # get coarsened ng
+    ranges = qda_coarsen.getRanges()
+    ng = np.zeros(2)
+    for i,nrange in enumerate(ranges):
+        ng[i]=nrange[1]-nrange[0]
+
+    # reshape the coarsened vector
+    q_dim = [sol.state.meqn]
+    q_dim.extend(ng)
+    q_coarsen = gqVec_coarsen.getArray().copy().reshape(q_dim,order = 'F')
+    
+    # Create coarsen grid        
+    lower = sol.grid.lower; upper = sol.grid.upper
+    mx = qda_coarsen.getSizes()[0]; my = qda_coarsen.getSizes()[1]
+    x = petclaw.Dimension('x',lower[0],upper[0],mx)
+    y = petclaw.Dimension('y',lower[1],upper[1],my)
+    grid = petclaw.Grid([x,y])
+
+    # Creat coarsened state
+    state = petclaw.State(grid)
+    state.meqn = sol.state.meqn
+    state.t = sol.state.t
+    state.q = q_coarsen
+
+    return petclaw.Solution(state)
+
+if __name__== "__main__":
+    path='./experiment/_output/'
+    first_frame = 0
+    last_frame = 200
+    
+    coarsening_level = 3
+    
+    for frame in range(first_frame,last_frame+1):
+        print 'coarsening frame '+str(frame)
+        if frame<10:
+            filenumber=str(frame)
+            filenumber='000'+filenumber
+        elif frame<100:
+            filenumber=str(frame)
+            filenumber='00'+filenumber
+        elif frame<1000:
+            filenumber=str(frame)
+            filenumber='0'+filenumber
+        elif frame <10000:
+            filenumber=str(frame)
+
+        filename = ["claw_p.pkl%s" %filenumber,"claw_p.ptc%s" %filenumber,"claw_p.ptc%s.info" %filenumber]
+
+        original_sol = Solution()
+        read_petsc(original_sol,frame=frame,path=path,file_prefix='claw_p')
+        coarsened_sol = original_sol
+
+        for i in xrange(coarsening_level):
+            coarsened_sol = coarsen_xhalf(coarsened_sol)
+
+        # Write coarsened solution
+        #if not os.path.exists('./_output/_p/coarsened/'): 
+        #    os.mkdir('./_output/_p/coarsened/')
+        write_petsc(coarsened_sol,frame,'./',file_prefix='claw_p_coarsened')

source/figure_7.4/experiment/README.txt

+ABOUT THIS EXPERIMENT
+================
+Figure 7.4 shows the medium and initial condition (both shown zoomed-in) for the cylindrical solitary wave problem.
+
+The experimental output data is archived in compressed form (bzip2)
+
+Common Environment Software (IBM BlueGene/P) 
+==========================================
+- IBM XL Fortran 11.1 
+- IBM BlueGene/P System Driver V1R4M2_200_2010-100508P
+
+PyClaw Software 
+==========================================
+- Python 2.6
+- PETSc 3.1.8 (compiled with 64-bit integer support)
+- petsc4py 1.2
+- numpy 1.6.1 (patched on Shaheen to use IBM XL compilers)
+- http://github.com/clawpack/pyclaw/commit/XXX 
+   (use "git reset --hard XXX" in a PyClaw repository to switch to this version of
+   the code)
+- Riemann repository 
+- prod.sh script (included in this repository)
+
+RESPONSIBLE AUTHORS
+==============
+Manuel Quezada
+Aron Ahmadia

source/figure_7.4/experiment/_output/claw_p_coarsened.pkl0000

+(dp0
+S'ndim'
+p1
+I2
+sS'maux'
+p2
+I0
+sS'meqn'
+p3
+L1L
+sS'aux_global'
+p4
+(dp5
+sS'write_aux'
+p6
+I00
+sS't'
+p7
+F0.0
+sS'nstates'
+p8
+I1
+s.(dp0
+S'lower'
+p1
+(lp2
+F0.25
+aF0.25
+asS'd'
+p3
+(lp4
+F0.033333333333333333
+aF0.033333333333333333
+asS'level'
+p5
+I1
+sS'stateno'
+p6
+I1
+sS'n'
+p7
+(lp8
+I6000
+aI6000
+asS'names'
+p9
+(lp10
+S'x'
+p11
+aS'y'
+p12
+as.

source/figure_7.4/experiment/_output/claw_p_coarsened.ptc0000

Binary file added.

source/figure_7.4/experiment/prod.sh

+pyclaw_build() {
+    export KSL_PROCESS_MODULES=/project/k47/prod/process/modules
+    export MODULE_HOME=/project/k47/modules
+    export MODULEPATH=$KSL_PROCESS_MODULES/shared/tools:$KSL_PROCESS_MODULES/ppc64/compilers:$KSL_PROCESS_MODULES/ppc64/targets:$MODULE_HOME/ppc64:$MODULE_HOME/ppc450d
+    export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/share/ksl/tcl/8.5.8/ppc64/lib
+    export PYTHONPATH=/project/k47/prod/numpy/dev-aug29/ppc450d/lib/python
+    module load init genie ibm python numpy
+    if [ -z $CLAW ]
+    then 
+        echo "CLAW Environment variable not set!  Please set before compiling." 1>&2
+    fi
+    if [ -z $PYCLAW ]
+    then 
+        echo "PYCLAW Environment variable not set!  Please set before compiling." 1>&2
+    fi
+    if [ -z $RIEMANN ]
+    then 
+        echo "RIEMANN Environment variable not set!  Please set before compiling." 1>&2
+    fi
+    alias make_f2py='make F2PY="f2py -L/opt/ibmcmp/xlsmp/bg/1.7/lib -L/opt/ibmcmp/xlmass/bg/4.4/bglib -L/opt/ibmcmp/xlf/bg/11.1/bglib -lxlf90_r -lxlopt -lxlsmp -lxl -lxlfmath" FC="/bgsys/drivers/ppcfloor/comm/bin/mpixlf77_r -qextname -I/project/k47/prod/python/2.6/ppc64/include -I/project/k47/prod/readline/6.1/ppc64/include -I/project/k47/prod/numpy/dev-aug29/ppc450d/lib/python/numpy/core/include -O3 -qpic=small"'
+    echo "*****************************************************************"
+    echo "make_f2py alias initialized"
+    echo "please use" 
+    echo 'make F2PY="f2py -L/opt/ibmcmp/xlsmp/bg/1.7/lib -L/opt/ibmcmp/xlmass/bg/4.4/bglib -L/opt/ibmcmp/xlf/bg/11.1/bglib -lxlf90_r -lxlopt -lxlsmp -lxl -lxlfmath" FC="/bgsys/drivers/ppcfloor/comm/bin/mpixlf77_r -qextname -I/project/k47/prod/python/2.6/ppc64/include -I/project/k47/prod/readline/6.1/ppc64/include -I/project/k47/prod/numpy/dev-aug29/ppc450d/lib/python/numpy/core/include -O3 -qpic=small"'
+    echo "or" 
+    echo 'make_f2py'
+    echo "*****************************************************************"
+    echo "the following commands are available:"
+    echo "make_f2py    - calls make with the appropriate F2PY variable set"
+    echo "*****************************************************************"
+    echo "Example:"
+    echo 'cd $PYCLAW/apps/advection/1d/constant'
+    echo 'make_f2py'
+    echo "*****************************************************************"
+    echo "numpy/f2py/xlf/python2.6 build environment initialized"
+    echo "*****************************************************************"
+ }
+
+pyclaw_run() {
+    err=0
+    if [ -z $CLAW ]
+    then 
+        echo "CLAW Environment variable not set!  Please set and then re-run this script." 1>&2
+        err=1
+    fi
+    if [ -z $PYCLAW ]
+    then 
+        echo "PYCLAW Environment variable not set!  Please set and then re-run this script." 1>&2
+        err=1
+    fi
+    if [ -z $RIEMANN ]
+    then 
+        echo "RIEMANN Environment variable not set!  Please set and then re-run this script." 1>&2
+        err=1
+    fi
+    if [ err = 1 ]
+    then
+        return 1
+    fi
+    export PATH=$PATH:/opt/share/ksl/python3/3.2/ppc64/bin
+    export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/share/ksl/readline/6.1/ppc64/lib:/opt/share/ksl/python3/3.2/ppc64/lib
+    export KSL_PPC450D_PYTHONPATH=$KSL_PPC450D_PYTHONPATH:/project/k47/prod/petsc4py/dev-jun2011/ppc450d-3.1.8_64/lib/python:/project/k47/prod/numpy/dev-aug29/ppc450d/lib/python:$PYCLAW/src:$RIEMANN
+    export KSL_PPC450D_LD_LIBRARY_PATH=$KSL_PPC450D_LD_LIBRARY_PATH:/project/k47/prod/petsc/3.1.8/ppc450d-bgp_gnu_shared_64/lib:/project/k47/prod/numpy/dev-aug29/ppc450d/lib:/project/k47/prod/petsc4py/dev-jun2011/ppc450d-3.1.8_64/lib/python/petsc4py/lib:/opt/ibmcmp/xlsmp/bg/1.7/lib:/opt/ibmcmp/xlmass/bg/4.4/bglib:/opt/ibmcmp/xlf/bg/11.1/bglib
+    export BGP_PYTHON=/bgsys/drivers/ppcfloor/gnu-linux/bin/python
+    echo "*****************************************************************"
+    echo "the following commands are available:"
+    echo "kslrun    - simplifies launching jobs on Shaheen"
+    echo "*****************************************************************"
+    echo "Example:"
+    echo 'cd $PYCLAW/apps/advection/1d/constant'
+    echo 'kslrun -a k47 -n 4 "$BGP_PYTHON advection.py kernel_language=Fortran use_petsc=True"'
+    echo ''
+    echo 'type kslrun -h for configuration options'
+    echo "*****************************************************************"
+    echo "pyclaw execution environment initialized"
+    echo "*****************************************************************"
+}
+
+echo "pyclaw production execution environment on Shaheen initialized"
+echo "the following commands are available:"
+echo "pyclaw_build - enables compiling/building Fortran kernel modules"
+echo "pyclaw_run   - enables launching jobs using kslrun or loadleveler scripts"

source/figure_7.4/experiment/psystem.py

+#!/usr/bin/env python
+# encoding: utf-8
+
+import numpy as np
+
+# material parameters
+E1=1.;   p1=1.
+E2=5.;   p2=5.
+# interface parameters
+alphax=0.5; deltax=1.
+alphay=0.5; deltay=1.
+# Linearity parameters
+linearity_mat1=2; linearity_mat2=2
+# heterogeneity type
+het_type='checkerboard'
+#het_type='sinusoidal'
+#het_type='smooth_checkerboard'
+sharpness=10
+
+def qinit(state,A,x0,y0,varx,vary):
+    r""" Set initial conditions for q."""
+    x =state.grid.x.center; y =state.grid.y.center
+    # Create meshgrid
+    [yy,xx]=np.meshgrid(y,x)
+    s=A*np.exp(-(xx-x0)**2/(2*varx)-(yy-y0)**2/(2*vary)) #sigma(@t=0)
+    #parameters from aux
+    linearity_mat=state.aux[2,:]
+    E=state.aux[1,:]
+    #initial condition
+    state.q[0,:,:]=np.where(linearity_mat==1,1,0)*s/E+np.where(linearity_mat==2,1,0)*np.log(s+1)/E
+    state.q[1,:,:]=0; state.q[2,:,:]=0
+
+def setaux(x,y):
+    r"""Creates a matrix representing every grid cell in the domain, 
+    whose size is len(x),len(y)
+    Each entry of the matrix contains a vector of size 3 with:
+    The material density p
+    The young modulus E
+    A flag indicating which material the grid is made of.
+    The domain pattern is a checkerboard."""
+
+    aux = np.empty((4,len(x),len(y)), order='F')
+    if het_type == 'checkerboard':
+        # xfrac and yfrac are x and y relative to deltax and deltay resp.
+        xfrac=x-np.floor(x/deltax)*deltax
+        yfrac=y-np.floor(y/deltay)*deltay
+        # create a meshgrid out of xfrac and yfrac
+        [yyfrac,xxfrac]=np.meshgrid(yfrac,xfrac)
+        # density 
+        aux[0,:,:]=p1*(xxfrac<=alphax*deltax)*(yyfrac<=alphay*deltay)\
+            +p1*(xxfrac >alphax*deltax)*(yyfrac >alphay*deltay)\
+            +p2*(xxfrac >alphax*deltax)*(yyfrac<=alphay*deltay)\
+            +p2*(xxfrac<=alphax*deltax)*(yyfrac >alphay*deltay)
+        #Young modulus
+        aux[1,:,:]=E1*(xxfrac<=alphax*deltax)*(yyfrac<=alphay*deltay)\
+            +E1*(xxfrac >alphax*deltax)*(yyfrac >alphay*deltay)\
+            +E2*(xxfrac >alphax*deltax)*(yyfrac<=alphay*deltay)\
+            +E2*(xxfrac<=alphax*deltax)*(yyfrac >alphay*deltay)
+        # linearity of material
+        aux[2,:,:]=linearity_mat1*(xxfrac<=alphax*deltax)*(yyfrac<=alphay*deltay)\
+            +linearity_mat1*(xxfrac >alphax*deltax)*(yyfrac >alphay*deltay)\
+            +linearity_mat2*(xxfrac >alphax*deltax)*(yyfrac<=alphay*deltay)\
+            +linearity_mat2*(xxfrac<=alphax*deltax)*(yyfrac >alphay*deltay)
+    elif het_type == 'sinusoidal' or het_type == 'smooth_checkerboard':
+        [yy,xx]=np.meshgrid(y,x)
+        Amp_p=np.abs(p1-p2)/2; offset_p=(p1+p2)/2
+        Amp_E=np.abs(E1-E2)/2; offset_E=(E1+E2)/2
+        if het_type == 'sinusoidal':
+            frec_x=2*np.pi/deltax; frec_y=2*np.pi/deltay
+            fun=np.sin(frec_x*xx)*np.sin(frec_y*yy)
+        else:
+            fun_x=xx*0; fun_y=yy*0
+            for i in xrange(0,1+int(np.ceil((x[-1]-x[0])/(deltax*0.5)))):
+                fun_x=fun_x+(-1)**i*np.tanh(sharpness*(xx-deltax*i*0.5))
+            for i in xrange(0,1+int(np.ceil((y[-1]-y[0])/(deltay*0.5)))):
+                fun_y=fun_y+(-1)**i*np.tanh(sharpness*(yy-deltay*i*0.5))
+            fun=fun_x*fun_y
+        aux[0,:,:]=Amp_p*fun+offset_p
+        aux[1,:,:]=Amp_E*fun+offset_E
+        aux[2,:,:]=linearity_mat1
+    return aux
+
+def b4step(solver,solutions):
+    r"""put in aux[3,:,:] the value of q[0,:,:] (eps). This is required in rptpv.f"""
+    state = solutions['n'].states[0]   
+    state.aux[3,:,:] = state.q[0,:,:]
+
+    # To set to 0 1st 1/2 of the domain. Used in rect domains with PBC in x
+    if state.aux_global['turnZero_half_2D']==1:
+        if state.t>=state.aux_global['t_turnZero'] and state.t<=state.aux_global['t_turnZero']+1:
+            if state.grid.x.nend <= np.floor(state.grid.x.n/2):
+                state.q[:,:,:]=0
+
+def compute_p(state):
+    state.p[0,:,:]=np.exp(state.q[0,:,:]*state.aux[1,:,:])-1
+
+def compute_F(state):
+    rho = state.aux[0,:,:]; E = state.aux[1,:,:]
+    
+    #Compute the entropy
+    u = state.q[1,:,:]/rho
+    v = state.q[2,:,:]/rho
+
+    nrg=rho * (u**2 + v**2)/2.
+
+    eps = state.q[0,:,:]
+    sigma = np.exp(E*eps) - 1.
+    sigint = (sigma-np.log(sigma+1.))/E
+
+    dx=state.grid.d[0]; dy=state.grid.d[1]
+    
+    state.F[0,:,:] = (sigint+nrg)*dx*dy 
+
+def gauge_pfunction(q,aux):
+    p = np.exp(q[0]*aux[1])-1
+    return [p]
+
+def gauges(radii,thetas):
+    gauges_list=[]
+    for radius in radii:
+        for theta in thetas: 
+            gauges_list.append([radius,theta])
+    return gauges_list
+
+def psystem2D(use_petsc=True,solver_type='classic',iplot=False,htmlplot=False):
+    """
+    Solve the p-system in 2D with variable coefficients
+    """
+    if use_petsc:
+        import petclaw as pyclaw
+    else:
+        import pyclaw
+
+    ####################################
+    ######### MAIN PARAMETERS ##########
+    ####################################
+    # Domain
+    x_lower=0.25; x_upper=200.25
+    y_lower=0.25; y_upper=200.25
+    # Grid cells per layer
+    Ng=240
+    mx=(x_upper-x_lower)*Ng; my=(y_upper-y_lower)*Ng
+    # Initial condition parameters
+    A=5.
+    x0=0.25 # Center of initial perturbation
+    y0=0.25 # Center of initial perturbation
+    varx=5.0; vary=5.0 # Width of initial perturbation
+    # Boundary conditions
+    mthbc_x_lower=pyclaw.BC.reflecting; mthbc_x_upper=pyclaw.BC.outflow
+    mthbc_y_lower=pyclaw.BC.reflecting; mthbc_y_upper=pyclaw.BC.outflow
+    # Turning off 1st half of the domain. Useful in rect domains
+    turnZero_half_2D=0 #flag
+    t_turnZero=50
+    # Regarding time
+    tfinal=200.0
+    nout=200
+    t0=0.0
+    # restart options
+    restart_from_frame = None
+    #solver = pyclaw.ClawSolver2D()
+    solver = pyclaw.SharpClawSolver2D()
+    solver.mwaves = 2
+    solver.limiters = pyclaw.limiters.tvd.MC
+
+    solver.mthbc_lower[0]=mthbc_x_lower
+    solver.mthbc_upper[0]=mthbc_x_upper
+    solver.mthbc_lower[1]=mthbc_y_lower
+    solver.mthbc_upper[1]=mthbc_y_upper
+    solver.mthauxbc_lower[0]=mthbc_x_lower
+    solver.mthauxbc_upper[0]=mthbc_x_upper
+    solver.mthauxbc_lower[1]=mthbc_y_lower
+    solver.mthauxbc_upper[1]=mthbc_y_upper
+
+    solver.fwave = True
+    solver.cfl_max = 2.5
+    solver.cfl_desired = 2.45
+    solver.start_step = b4step
+    solver.dim_split=False
+    solver.max_steps = 1000000000
+    
+    #controller
+    claw = pyclaw.Controller()
+    claw.tfinal = tfinal
+    claw.solver = solver
+
+    if restart_from_frame is not None:
+        claw.solution = pyclaw.Solution(restart_from_frame, format='petsc',read_aux=False)
+        claw.solution.state.mp = 1
+        claw.solution.state.mF = 1
+        grid = claw.solution.grid
+        claw.solution.state.aux = setaux(grid.x.center,grid.y.center)
+        claw.nout = nout - restart_from_frame
+        claw.start_frame = restart_from_frame
+    else:
+        ####################################
+        ####################################
+        ####################################
+        #Creation of grid
+        x = pyclaw.Dimension('x',x_lower,x_upper,mx)
+        y = pyclaw.Dimension('y',y_lower,y_upper,my)
+        grid = pyclaw.Grid([x,y])
+        state = pyclaw.State(grid)
+        state.meqn = 3
+        state.mF = 1
+        state.mp = 1
+        state.t=t0
+        #Set global parameters
+        state.aux_global = {}
+        state.aux_global['turnZero_half_2D'] = turnZero_half_2D
+        state.aux_global['t_turnZero'] = t_turnZero
+
+        state.aux = setaux(grid.x.center,grid.y.center)
+        #Initial condition
+        qinit(state,A,x0,y0,varx,vary)
+
+        claw.solution = pyclaw.Solution(state)
+        claw.nout = nout
+
+    claw.compute_p = compute_p
+    claw.compute_F = compute_F
+    gauges_radii=[10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200]
+    gauges_thetas=[0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90]
+    grid.add_gauges(gauges(gauges_radii,gauges_thetas))
+    solver.compute_gauge_values = gauge_pfunction
+    claw.write_aux_init = False
+
+    #Solve
+    status = claw.run()
+    
+    #strain=claw.frames[claw.nout].state.gqVec.getArray().reshape([grid.ng[0],grid.ng[1],state.meqn])[:,:,0]
+    #return strain
+
+    if iplot:    pyclaw.plot.plotInteractive()
+    if htmlplot: pyclaw.plot.plotHTML()
+
+if __name__=="__main__":
+    import sys
+    from pyclaw.util import _info_from_argv
+    args, kwargs = _info_from_argv(sys.argv)
+    psystem2D(*args,**kwargs)

source/figure_7.4/plots.py

+from pyclaw.solution import Solution
+from petclaw.io.petsc import read_petsc
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pylab as pl
+from matplotlib import rc
+rc('text', usetex=True)
+import matplotlib.cm as cm
+import numpy as np
+import os
+
+def plot_checker(Ka=5,Kb=1,length=10.25):
+    x = np.linspace(0.25,length,1000)
+    y = np.linspace(0.25,length,1000)
+    yy,xx = np.meshgrid(y,x)
+    #checker                                                      
+    xfrac=x-np.floor(x);
+    yfrac=y-np.floor(y);
+    [yyfrac,xxfrac]=np.meshgrid(yfrac,xfrac);
+
+    checker = Ka*(xxfrac<=0.5)*(yyfrac<=0.5)+Ka*(xxfrac>0.5)*(yyfrac>0.5)+Kb*(xxfrac>0.5)*(yyfrac<=0.5)+Kb*(xxfrac<=0.5)*(yyfrac>0.5);
+    pl.pcolormesh(xx,yy,checker,cmap=cm.gray)
+    pl.xlabel('x',fontsize=20); pl.ylabel('y',fontsize=20)
+    pl.xticks(size=20); pl.yticks(size=20)
+    pl.axis('scaled')
+    pl.axis([np.min(x),np.max(x),np.min(y),np.max(y)])
+    pl.savefig('checker.png')
+    pl.close()
+
+
+def plot_p(frame,skip=1,file_prefix='claw_p',path='./_output/_p',title='',axis=None,file_name='fig',axis_slices=None):
+    sol=Solution()
+    read_petsc(sol,frame=frame,path=path,file_prefix=file_prefix)
+    x=sol.grid.x.center; y=sol.grid.y.center
+    mx=len(x); my=len(y)
+    subx=range(0,mx,skip)
+    suby=range(0,my,skip)
+
+    mp=sol.state.meqn
+
+    yy,xx = np.meshgrid(y[suby],x[subx])
+
+    for i in xrange(mp):
+        p=sol.state.q[i,:,:]
+        p_subx=p[subx,:]
+        p_subxy=p_subx[:,suby]
+        
+        pl.pcolormesh(xx,yy,p_subxy,cmap=cm.OrRd)
+        #pl.pcolormesh(xx,yy,p_subxy)
+        pl.title(title,fontsize=20)
+        pl.xticks(size=20); pl.yticks(size=20)
+        cb = pl.colorbar();
+        imaxes = pl.gca(); pl.axes(cb.ax)
+        pl.yticks(fontsize=20); pl.axes(imaxes)
+        if axis is not None:
+            pl.axis([axis[0],axis[1],axis[0],axis[1]])
+        else:
+            pl.axis([np.min(x),np.max(x),np.min(y),np.max(y)])
+        pl.savefig(file_name+'.png')
+        pl.close()
+        if axis_slices is not None:
+            p_45=np.zeros(my)
+            for j in xrange(0,my):
+                p_45[j]=p[j,j]
+            x_45=x*np.sqrt(2)
+            pl.plot(x_45,p_45,'-k')
+            pl.plot(x,p[:,2],'--r')
+            pl.xlabel('Radial distance',fontsize=20)
+            pl.ylabel('Stress',fontsize=20)
+            pl.xticks(size=20); pl.yticks(size=20)
+            pl.axis([axis_slices[0],axis_slices[1],np.min(p),np.max(p)])
+            pl.savefig(file_name+'_slices.png')
+            pl.close()
+
+if __name__== "__main__":
+    plot_initial_condition = True
+    plot_medium = True
+    plot_steg = False
+
+    path='./experiment/_output'
+
+    if plot_initial_condition:
+        print('**********************')
+        print('**********************')
+        print('Plotting initial condition ...')
+        plot_p(frame=0,skip=1,path=path,file_prefix='claw_p_coarsened',title='Initial Condition',axis=[0.25,10.25],file_name='initial_condition')
+
+    if plot_medium:
+        print('**********************')
+        print('**********************')
+        print('Plotting medium ...')
+        plot_checker()
+
+    if plot_steg:
+        print('**********************')
+        print('**********************')
+        print('Plotting stegoton ...')
+        plot_p(frame=200,skip=1,path=path,file_prefix='claw_p_coarsened',title='t=90.0',file_name='stress_checker',axis_slices=[160,190])

source/figure_7.5/Makefile

+coarsen_data: experiment/_output/claw_p.pkl0000 experiment/_output/claw_p.ptc0000
+	python coarsen.py 
+
+compress_data: claw_p_coarsened.pkl0000 claw_p_coarsened.ptc0000
+	bzip2 claw_p_coarsened.pkl0000 claw_p_coarsened.ptc0000
+
+figure_7.4: experiment/_output/claw_p_coarsened.pkl0000.bz2 experiment/_output/claw_p_coarsened.ptc0000.bz2 plots.py
+	bunzip2 experiment/_output/claw_p_coarsened.pkl0000.bz2 experiment/_output/claw_p_coarsened.ptc0000.bz2 
+	python plots.py

source/figure_7.5/README.txt

+ABOUT THIS EXPERIMENT
+================
+Figure 7.5 shows the solution of the p-system in a checkerboard medium, demonstrating 2D cylindrical solitary wave formation
+
+The experimental output data is archived in compressed form (bzip2).  The uncompressed PETSc data uses 64-bit integer indices, to load it you will need to use a PETSc built with the following configuration option: --with-64-bit-indices=1
+
+RESPONSIBLE AUTHORS
+==============
+Manuel Quezada
+Aron Ahmadia

source/figure_7.5/coarsen.py

+from pyclaw.solution import Solution
+from petclaw.io.petsc import read_petsc
+from petclaw.io.petsc import write_petsc
+from petsc4py import PETSc
+import petclaw
+import numpy as np
+import os
+
+def read_original_solution(frame,file_prefix='claw_p',path='./experiment/_output/_p'):
+    sol = Solution()
+    read_petsc(sol,frame=frame,path=path,file_prefix=file_prefix)
+    return sol
+
+def coarsen_xhalf(sol):
+    # Interpolation
+    sol.state.q_da.setBoundaryType(PETSc.DA.BoundaryType.NONE)
+    sol.state.q_da.setRefinementFactor(refine_x=2,refine_y=2)
+    qda_coarsen = sol.state.q_da.coarsen()
+    gqVec_coarsen = qda_coarsen.createGlobalVec()
+    qda_coarsen.setInterpolationType(interp_type=0)
+    mat,scale = qda_coarsen.getInterpolation(sol.state.q_da)
+    mat.multTranspose(sol.state.gqVec,gqVec_coarsen)
+    gqVec_coarsen = scale*gqVec_coarsen
+    qda_coarsen.setBoundaryType(PETSc.DA.BoundaryType.PERIODIC)
+
+    # get coarsened ng
+    ranges = qda_coarsen.getRanges()
+    ng = np.zeros(2)
+    for i,nrange in enumerate(ranges):
+        ng[i]=nrange[1]-nrange[0]
+
+    # reshape the coarsened vector
+    q_dim = [sol.state.meqn]
+    q_dim.extend(ng)
+    q_coarsen = gqVec_coarsen.getArray().copy().reshape(q_dim,order = 'F')
+    
+    # Create coarsen grid        
+    lower = sol.grid.lower; upper = sol.grid.upper
+    mx = qda_coarsen.getSizes()[0]; my = qda_coarsen.getSizes()[1]
+    x = petclaw.Dimension('x',lower[0],upper[0],mx)
+    y = petclaw.Dimension('y',lower[1],upper[1],my)
+    grid = petclaw.Grid([x,y])
+
+    # Creat coarsened state
+    state = petclaw.State(grid)
+    state.meqn = sol.state.meqn
+    state.t = sol.state.t
+    state.q = q_coarsen
+
+    return petclaw.Solution(state)
+
+if __name__== "__main__":
+    path='./experiment/_output/'
+    first_frame = 0
+    last_frame = 200
+    
+    coarsening_level = 3
+    
+    for frame in range(first_frame,last_frame+1):
+        print 'coarsening frame '+str(frame)
+        if frame<10:
+            filenumber=str(frame)
+            filenumber='000'+filenumber
+        elif frame<100:
+            filenumber=str(frame)
+            filenumber='00'+filenumber
+        elif frame<1000:
+            filenumber=str(frame)
+            filenumber='0'+filenumber
+        elif frame <10000:
+            filenumber=str(frame)
+
+        filename = ["claw_p.pkl%s" %filenumber,"claw_p.ptc%s" %filenumber,"claw_p.ptc%s.info" %filenumber]
+
+        original_sol = Solution()
+        read_petsc(original_sol,frame=frame,path=path,file_prefix='claw_p')
+        coarsened_sol = original_sol
+
+        for i in xrange(coarsening_level):
+            coarsened_sol = coarsen_xhalf(coarsened_sol)
+
+        # Write coarsened solution
+        #if not os.path.exists('./_output/_p/coarsened/'): 
+        #    os.mkdir('./_output/_p/coarsened/')
+        write_petsc(coarsened_sol,frame,'./',file_prefix='claw_p_coarsened')

source/figure_7.5/experiment/README.txt

+ABOUT THIS EXPERIMENT
+================
+Figure 7.5 shows the solution of the p-system in a checkerboard medium, demonstrating 2D cylindrical solitary wave formation
+
+The experimental output data is archived in compressed form (bzip2).
+
+Environment Software (IBM BlueGene/P) 
+==========================================
+- IBM XL Fortran 11.1 
+- IBM BlueGene/P System Driver V1R4M2_200_2010-100508P
+
+PyClaw Software 
+==========================================
+- Python 2.6
+- PETSc 3.1.8 (compiled with 64-bit integer support)
+- petsc4py 1.2
+- numpy 1.6.1 (patched on Shaheen to use IBM XL compilers)
+- http://github.com/clawpack/pyclaw/commit/XXX 
+   (use "git reset --hard XXX" in a PyClaw repository to switch to this version of
+   the code)
+- Riemann repository 
+- prod.sh script (included in this repository)
+
+RESPONSIBLE AUTHORS
+==============
+Manuel Quezada
+Aron Ahmadia