def test_tr(grids=[1800,3600,7200],frame=2):
    reload(setrun)
    errs = []

    for i,mx in enumerate(grids):
        print i,mx
        rundata=setrun.setrun()

        rundata.probdata.trtime=600.
        rundata.probdata.t1= 10.
        rundata.probdata.a1= 0.1
        rundata.probdata.tw1=10.

        rundata.clawdata.mx=mx
        rundata.clawdata.tfinal=1200.
        rundata.clawdata.nout=60
        rundata.write()
        runclaw(xclawcmd='xsclaw',outdir='./_output')

        exact,approx,xc = timereverse(frame=frame)

        Linf_err = max(abs(approx-exact))
        print mx, Linf_err

        errs.append(Linf_err)

    return grids,errs

    def __init__(self, fault, name=""): 
        r"""
        Initialize a FaultJob object.
        
        See :class:`FaultJob` for full documentation
        
        """ 

        super(FaultJob, self).__init__()

        self.type = "tsunami"
        self.name = "final-tohoku-inversions"
        self.prefix = "fault_%s" % name
        self.executable = 'xgeoclaw'

        # Data objects
        import setrun
        self.rundata = setrun.setrun()

        # No variable friction for the time being
        self.rundata.friction_data.variable_friction = False

        # Replace dtopo file with our own
        self.dtopo_path = 'fault_%s.tt3' % name
        self.rundata.dtopo_data.dtopofiles = [[3, 4, 4, self.dtopo_path]]

        self.fault = fault

    def __init__(self, index, values, 
                       source_path='clawpack/geoclaw/examples/tsunami/chile2010'):
        r""""""

        super(LevelsJob, self).__init__()

        self.source_path = os.path.expandvars(source_path)

        self.type = "Run_Data"
        self.name = "AMR-Levels-%s" % \
                        os.path.split(os.path.splitext(self.source_path)[0])[-1]
        self.prefix = "sweep_%s" % (str(index))
        self.executable = 'xgeoclaw'

        # Data objects
        import setrun
        self.rundata = setrun.setrun()
        
        values = [int(i) for i in values]
        # Geo Claw Number of Cells Data
        self.rundata.clawdata.num_cells = [values[0],values[0]]

        # AMR Data max levels and refinement ratios
        self.rundata.amrdata.amr_levels_max = values[1]
        self.rundata.amrdata.refinement_ratios_x = values[2:]
        self.rundata.amrdata.refinement_ratios_y = values[2:]
        self.rundata.amrdata.refinement_ratios_t = values[2:]

    def __init__(self, index, values, 
                       source_path='$SRC/tohoku2011-paper1/sources/Ammon.txydz',
                       contours=(numpy.infty,0.0,-200.0,-numpy.infty)):
        r""""""

        super(FrictionJob, self).__init__()

        self.source_path = os.path.expandvars(source_path)

        self.type = "tohoku"
        self.name = "friction-sensitivity-%s" % \
                        os.path.split(os.path.splitext(self.source_path)[0])[-1]
        self.prefix = "fric_%s" % (str(index))
        self.executable = 'xgeoclaw'

        # Data objects
        import setrun
        self.rundata = setrun.setrun()

        # Set variable friction
        self.rundata.friction_data.variable_friction = True

        # Region based friction
        self.rundata.friction_data.friction_regions = [[
                       self.rundata.clawdata.lower, self.rundata.clawdata.upper,
                       contours, values]]

        # Set earthquake source model
        dtopo_data = self.rundata.dtopo_data
        dtopo_data.dtopofiles.append([1, 4, 4, self.source_path])

    def runTest(self):

        # Write out data files
        orig_path = os.getcwd()
        if sys.modules.has_key('setrun'):
            del(sys.modules['setrun'])
        sys.path.insert(0, self.test_path)
        import setrun
        rundata = setrun.setrun()
        os.chdir(self.temp_path)
        rundata.write()
        os.chdir(orig_path)
        sys.path.pop(0)

        # Run code
        runclaw_cmd = " ".join((
                            "cd %s ;" % self.temp_path,
                            "python",
                            "$CLAW/clawutil/src/python/clawutil/runclaw.py",
                            "xgeoclaw",
                            self.temp_path,
                            "True",
                            "False",
                            self.temp_path))
        subprocess.check_call(runclaw_cmd, stdout=self.stdout, 
                                           stderr=self.stderr,
                                           shell=True)
        self.stdout.flush()
        self.stderr.flush()

 def __init__(self):
     super(TestML1D,self).__init__()
     
     import setrun
     
     self.run_data = setrun.setrun()
     self.ml_data = setrun.set_multilayer_data()
     
     self.type = "ml_1d"

def getent(K,rho):
  rundata=setrun.setrun()
  rundata.probdata.K_B=K
  rundata.probdata.rho_B=rho
  print K,rho
  rundata.write()
  runclaw(outdir='./_output')

  ent=entropy(rundata.clawdata.nout)
  return ent[-1]

    def __init__(self, levels=3):
        super(SwirlTest, self).__init__()

        self.type = "AMRClaw"
        self.name = "swirl"
        self.setplot = "setplot"

        # Data objects
        import setrun
        self.rundata = setrun.setrun()

        self.rundata.amrdata.amr_levels_max = levels

    def load_rundata(self):
        r"""(Re)load setrun module and create *rundata* object


        """

        if sys.modules.has_key("setrun"):
            del (sys.modules["setrun"])
        sys.path.insert(0, self.test_path)
        import setrun

        self.rundata = setrun.setrun()
        sys.path.pop(0)

def gauges2kml(rundata=None, fname='gauges.kml'):

    """
    Read in the gauge locations from setrun.py and create a kml point for each.
    """


    if rundata is None:
        try:
            import setrun
            reload(setrun)
            rundata = setrun.setrun()
        except:
            raise IOError("*** cannot execute setrun file")

    elev = 0.
    kml_text = kml_header()
    

    gauges = rundata.gaugedata.gauges
    if len(gauges)==0:
        print "No gauges found in setrun.py"


    for rnum,gauge in enumerate(gauges):
        t1,t2 = gauge[3:5]
        x1,y1 = gauge[1:3]
        gaugeno = gauge[0]
        print "Gauge %i: %10.6f  %10.6f  \n" % (gaugeno,x1,y1) \
                + "  t1 = %10.1f,  t2 = %10.1f" % (t1,t2)
        mapping = {}
        mapping['gaugeno'] = gaugeno
        mapping['t1'] = t1
        mapping['t2'] = t2
        mapping['x1'] = x1
        mapping['y1'] = y1
        mapping['elev'] = elev
        mapping['name'] = 'Gauge %i' % rnum
        description = "  t1 = %g, t2 = %g\n" % (t1,t2) \
            + "  x1 = %g, y1 = %g\n" % (x1,y1)
        mapping['desc'] = description

        gauge_text = kml_gauge(mapping)
        kml_text = kml_text + gauge_text
    kml_text = kml_text + kml_footer()
    kml_file = open(fname,'w')
    kml_file.write(kml_text)
    kml_file.close()
    print "Created ",fname

def run_sw(hl=3.,ul=0.,hr=1.,ur=0.):
	os.system("make .exe")
	import setrun
	rundata = setrun.setrun()  # initialize most run-time variables for clawpack
	outdir = '_output'

	# IC for shallow water code
	rundata.probdata.hl = hl
	rundata.probdata.ul = ul
	rundata.probdata.hr = hr
	rundata.probdata.ur = ur
	rundata.write()

	os.system("make output > output.txt")
	os.system("make plots >> output.txt")

    def __init__(self, test = "A"):
        super(BenchmarkBaseTest,self).__init__()
        
        import setrun
        self.run_data = setrun.setrun()

        self.name = "TestOutput"
        self.test = test
       
        if test == "A":
            self.run_data.clawdata.nout = 22
        elif test == "C":
            self.run_data.clawdata.nout = 28

    # Convert angle to degrees for the label
        self.prefix = "Case%s" % test

    def __init__(self, storm, storm_number):
        r"""
        Initialize Habanero storm surge job

        See :class:`StormJob` for full documentation
        """

        super(StormJob, self).__init__()

        self.storm_number = storm_number
        self.storm = storm

        # Habanero queue settings
        self.omp_num_threads = 24
        self.time = "1:00:00"
        self.queue = ""

        # General job info
        self.type = "surge"
        self.name = "global_1"
        self.prefix = "storm_%s" % self.storm_number
        self.executable = "xgeoclaw"

        # Modify run data
        import setrun
        self.rundata = setrun.setrun()

        # Modify output times
        self.rundata.clawdata.output_style = 2
        recurrence = 6.0
        tfinal = (storm.t[-1] - storm.t[0]).total_seconds()
        N = int(tfinal / (recurrence * 60**2))
        self.rundata.clawdata.output_times = [t for t in
                 numpy.arange(0.0, N * recurrence * 60**2, recurrence * 60**2)]
        self.rundata.clawdata.output_times.append(tfinal)

        # Modify storm data
        surge_data = self.rundata.surge_data
        base_path = os.path.expandvars(os.path.join("$DATA_PATH", "storms",
                                                    "global", "storms"))
        surge_data.storm_file = os.path.join(base_path,
                                             'storm_%s.storm'
                                             % (str(i).zfill(5)))
        self.storm.time_offset = storm.t[0]

        print("Writing out GeoClaw storms...")
        self.storm.write(self.rundata.surge_data.storm_file)

    def __init__(self, ensemble, storm_num, base_path='./', storms_path='./'):

        super(StormJob, self).__init__()

        self.type = ""
        self.name = ""
        self.prefix = str(storm_num).zfill(5)
        self.storm_num = storm_num
        self.executable = "xgeoclaw"

        # Create base data object
        import setrun
        self.rundata = setrun.setrun()

        # Storm specific data
        self.ensemble = ensemble
        self.storm_file_path = os.path.join(storms_path
                               os.path.abspath("./%s.storm" % storm_num))

        # Set storm file
        self.rundata.storm_data.storm_type = 1
        self.rundata.storm_data.storm_file = self.storm_file_path

    def __init__(self, angle=0.0, bathy_angle=0.0, location=[-0.1, 0.0], 
                       bathy_location=0.15):

        super(PlaneWaveTest, self).__init__()

        self.executable = 'xgeoclaw'
        self.type = "multilayer"
        self.name = "planewave"

        # Convert angle to degrees for the label
        self.prefix = "ml_2d_ia%s_ba%s" % (int(angle * 180.0 / numpy.pi),
            int(bathy_angle * 180.0 / numpy.pi))

        # Data objects
        self.rundata = setrun.setrun()

        self.rundata.qinit_data.angle = angle
        self.rundata.qinit_data.init_location = location

        self.bathy_location = bathy_location
        self.bathy_angle = bathy_angle

        # Add gauges down perpendicular
        self.rundata.gaugedata.gauges = []
        gauge_locations = [-0.1,0.0,0.1,0.2,0.3]
        for (i,x_c) in enumerate(gauge_locations):
            # y0 = (self.run_data.clawdata.yupper - self.run_data.clawdata.ylower) / 2.0
            # x_p,y_p = transform_c2p(x_c,0.0,location[0],location[1],angle)
            x_p = x_c * numpy.cos(angle)
            y_p = x_c * numpy.sin(angle)
            # print "+=====+"
            # print x_c,0.0
            # print x_p,y_p
            if (self.rundata.clawdata.lower[0] < x_p < self.rundata.clawdata.upper[0] and
                    self.rundata.clawdata.lower[1] < y_p < self.rundata.clawdata.upper[1]):
                self.rundata.gaugedata.gauges.append([i, x_p, y_p, 0.0, 1e10])

    def __init__(self, deformation_file, topo_type=None):

        super(FaultJob, self).__init__()

        self.type = "compsyn"
        self.name = "guerrero_gap"
        self.executable = 'xgeoclaw'

        # Given an xyzt file, use this to run simulation
        self.deformation_file = os.path.abspath(deformation_file)
        self.prefix = os.path.basename(deformation_file).split('.')[0]

        # Data objects
        import setrun
        self.rundata = setrun.setrun()

        # Add deformation file
        self.rundata.dtopo_data.dtopofiles = []
        if topo_type is None:
            # Try to look at suffix for type
            extension = os.path.splitext(deformation_file)[1][1:]
            if extension[:2] == "tt":
                topo_type = int(extension[2])
            elif extension == 'xyz':
                topo_type = 1
            else:
                # Default to 3
                topo_type = 3
        self.rundata.dtopo_data.dtopofiles.append(
                                          [topo_type,3,5,self.deformation_file])
        # Add earthquake to regions list, force refinement for 1 minute
        regions = self.rundata.regiondata.regions
        regions.append([7, 7, 0.0, 60e0, -99.463937141374046, 
                                         -98.485815563597853, 
                                          16.622495995962375, 
                                          17.490586994378546 ])

def run_code_or_restart():

    import time
    tm = time.localtime()
    year = str(tm[0]).zfill(4)
    month = str(tm[1]).zfill(2)
    day = str(tm[2]).zfill(2)
    hour = str(tm[3]).zfill(2)
    minute = str(tm[4]).zfill(2)
    second = str(tm[5]).zfill(2)
    timestamp = '%s-%s-%s-%s%s%s'  % (year,month,day,hour,minute,second)

    finished, latest, t_latest = examine_outdir(outdir)

    if finished:
        print("Code has finished running, remove %s to run again" % outdir)
        return

    restart = (latest is not None)

    fname_output = 'run_output.txt'
    fname_errors = 'run_errors.txt'

    if restart:
        print("Will attempt to restart using checkpoint file %s at t = %s" \
                % (latest, t_latest))
        print("Appending output stream to %s" % fname_output)
        access = 'a'
    else:
        print("Will run code -- no restart")
        print("Writing output stream to %s" % fname_output)
        access = 'w'

    fout = open(fname_output, access)
    ferr = open(fname_errors, access)

    if restart:
        fout.flush()
        fout.write("\n=========== RESTART =============\n" + \
                "Local time: %s\n" % timestamp + \
                "Will attempt to restart using checkpoint file %s at t = %s\n" \
                % (latest, t_latest))
        fout.flush()
        make_args = ['make','output','RESTART=True']
    else:
        make_args = ['make','output']

    #if restart:
    #    No longer need to do this since new restart now adds to gauge*.txt files
    #    fortgauge = os.path.join(outdir,'fort.gauge')
    #    fortgauge2 = os.path.join(outdir,'fort.gauge_%s' % timestamp)
    #    os.system("mv %s %s" % (fortgauge,fortgauge2))
    #    fout.write("Moving %s to %s \n" % (fortgauge,fortgauge2))
    #    fout.flush()

    rundata = setrun('amrclaw')
    rundata.clawdata.restart = restart
    rundata.clawdata.restart_file = 'fort.chk' + str(latest)
    if restart:
        rundata.clawdata.output_t0 = False  # to avoid plotting at restart times
    rundata.write()
    
    job = subprocess.Popen(make_args, stdout=fout,stderr=ferr,env=env)
    return_code = job.wait()
    
    if return_code == 0:
        print("Successful run\n")
    else:
        print("Problem running code\n")
    print("See %s and %s" % (fname_output,fname_errors))
    
    fout.close()
    ferr.close()

#! /usr/bin/python

"""
Run several tests with different parameters, in this case different
values of d.
"""

from setrun import setrun
from setplot import setplot
from pyclaw.runclaw import runclaw
from pyclaw.plotters.plotclaw import plotclaw

for n in [1,2,4]:
    for d in [0.061, 0.080, 0.100, 0.120, 0.140, 0.149, 0.189]:
        rundata = setrun(d_param=d)
        mx = 72*n
        my = 18*n
        rundata.clawdata.mx = mx
        rundata.clawdata.my = my
        rundata.write()

        runclaw(xclawcmd = "xgeoclaw", outdir="_output_%s_my%s" % (d,my))

        setplot.d = d
        setplot.my = my
        plotclaw(outdir="_output_%s_my%s" % (d,my), \
                 plotdir="_plots_%s_my%s" % (d,my),\
                 setplot=setplot)

def gauges2kml(rundata=None, fname='gauges.kml', verbose=True):

    """

    Create a KML marker for each gauge specified for a GeoClaw run.

    :Inputs:

      - *rundata* - an object of class *ClawRunData* or None

        If *rundata==None*, try to create based on executing function *setrun*
        from the `setrun.py` file in the current directory.

      - *fname* (str) - resulting kml file.

      - *verbose* (bool) - If *True*, print out info about each region found


    :Example:

        >>> from clawpack.geoclaw import kmltools
        >>> kmltools.gauges2kml()

    is equivalent to:

        >>> from clawpack.geoclaw import kmltools
        >>> from setrun import setrun
        >>> rundata = setrun()
        >>> kmltools.gauges2kml(rundata)

    By default this creates a file named *gauges.kml* that can be opened in
    Google Earth.

    """


    if rundata is None:
        try:
            import setrun
            reload(setrun)
            rundata = setrun.setrun()
        except:
            raise IOError("*** cannot execute setrun file")

    elev = 0.
    kml_text = kml_header(fname)


    gauges = rundata.gaugedata.gauges
    if len(gauges)==0 and verbose:
        print("No gauges found in setrun.py")


    for rnum,gauge in enumerate(gauges):
        t1,t2 = gauge[3:5]
        x1,y1 = gauge[1:3]
        gaugeno = gauge[0]
        if verbose:
            print("Gauge %i: %10.6f  %10.6f  \n" % (gaugeno,x1,y1) \
                    + "  t1 = %10.1f,  t2 = %10.1f" % (t1,t2))
        mapping = {}
        mapping['gaugeno'] = gaugeno
        mapping['t1'] = t1
        mapping['t2'] = t2
        mapping['x1'] = x1
        mapping['y1'] = y1
        mapping['elev'] = elev
        mapping['name'] = 'Gauge %i' % rnum
        description = "  t1 = %g, t2 = %g\n" % (t1,t2) \
            + "  x1 = %g, y1 = %g\n" % (x1,y1)
        mapping['desc'] = description

        gauge_text = kml_gauge(mapping)
        kml_text = kml_text + gauge_text
    kml_text = kml_text + kml_footer()
    kml_file = open(fname,'w')
    kml_file.write(kml_text)
    kml_file.close()
    if verbose:
        print("Created ",fname)

def regions2kml(rundata=None,fname='regions.kml',verbose=True,combined=True):

    """
    Create a KML box for each AMR region specified for a GeoClaw run.

    :Inputs:

      - *rundata* - an object of class *ClawRunData* or None

        If *rundata==None*, try to create based on executing function *setrun*
        from the `setrun.py` file in the current directory.

      - *fname* (str) - resulting kml file.

      - *verbose* (bool) - If *True*, print out info about each region found

      - *combined* (bool) - If *True*, combine into single kml file with
        name given by *fname*.  This is the default. 
        If False, *fname* is ignored and individual files are created for
        each region with names are Domain.kml, Region00.kml, etc.
        These will show up separately in GoogleEarth so they can be turned
        on or off individually.

    First create a box for the entire domain (in red) and then a box
    for each region (in white).

    :Example:

        >>> from clawpack.geoclaw import kmltools
        >>> kmltools.regions2kml()

    is equivalent to:

        >>> from clawpack.geoclaw import kmltools
        >>> from setrun import setrun
        >>> rundata = setrun()
        >>> kmltools.regions2kml(rundata)

    By default this creates a file named *regions.kml* that can be opened in
    Google Earth.

    """

    from numpy import cos,pi,floor

    if rundata is None:
        try:
            import setrun
            reload(setrun)
            rundata = setrun.setrun()
        except:
            raise IOError("*** cannot execute setrun file")

    clawdata = rundata.clawdata
    x1,y1 = clawdata.lower[0:]
    x2,y2 = clawdata.upper[0:]
    description = "  x1 = %g, x2 = %g\n" % (x1,x2) \
            + "  y1 = %g, y2 = %g\n" % (y1,y2)

    mx,my = clawdata.num_cells[0:]
    dx = (x2-x1)/float(mx)
    dx_meters = dx*111e3*cos(pi*0.5*(y1+y2)/180.)
    dy = (y2-y1)/float(my)
    dy_meters = dy*111e3
    if verbose:
        print("Domain:   %10.6f  %10.6f  %10.6f  %10.6f" % (x1,x2,y1,y2))
    dx_deg,dx_min,dx_sec = deg2dms(dx)
    dy_deg,dy_min,dy_sec = deg2dms(dy)
    #print "Level 1 resolution:  dx = %g deg, %g min, %g sec = %g meters" \
    #   % (dx_deg,dx_min,dx_sec,dx_meters)
    levtext = "Level 1 resolution:  dy = %g deg, %g min, %g sec = %g meters\n" \
        % (dy_deg,dy_min,dy_sec,dy_meters)
    if verbose:
        print(levtext)
    description = description + levtext

    amr_levels_max = rundata.amrdata.amr_levels_max
    refinement_ratios_y = rundata.amrdata.refinement_ratios_y
    num_ref_ratios = len(refinement_ratios_y)
    if amr_levels_max > num_ref_ratios+1:
        raise IOError("*** Too few refinement ratios specified for " \
            + "amr_levels_max = %i" % amr_levels_max)
    dy_levels = (num_ref_ratios+1) * [dy]
    for k,r in enumerate(refinement_ratios_y):
        level = k+2
        dy = dy_levels[k] / r
        dy_levels[k+1] = dy
        dy_meters = dy*111e3
        dy_deg,dy_min,dy_sec = deg2dms(dy)
        levtext = "Level %s resolution:  dy = %g deg, %g min, %g sec = %g meters  (refined by %i)\n" \
                % (level,dy_deg,dy_min,dy_sec,dy_meters,r)
        if verbose:
            print(levtext)
        description = description + levtext

    if verbose:
        print("Allowing maximum of %i levels" % amr_levels_max)

    elev = 0.
    if not combined:
#.........这里部分代码省略.........