""" 
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='amrclaw'):
#------------------------------
    
    """ 
    Define the parameters used for running Clawpack.

    INPUT:
        claw_pkg expected to be "amrclaw" for this setrun.

    OUTPUT:
        rundata - object of class ClawRunData 
    
    """ 
    
    from clawpack.clawutil import data 
    
    
    assert claw_pkg.lower() == 'amrclaw',  "Expected claw_pkg = 'amrclaw'"

    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')

    # Gamma EOS parameter for three materials: air, plastic and water
    probdata.add_param('gammagas',     1.4,  'gamma for ideal gas')
    probdata.add_param('gammaplas',    1.1, 'gamma est. for polystirene')
    probdata.add_param('gammawat',     7.15,  'gamma for water')

    # Pinf parameter for Tammann EOS for three materials: air, plastic and water
    # pinfplas Calculated with c^2=gamma*(p+pinf)/rho to make c =2240m/s (polystyrene speed of sound), 
    # (c= 1484 in water). Values from water obtained from kirsten's paper in the references
    probdata.add_param('pinfgas',     0.0,  'pinf for stiffend gas/plastic')
    probdata.add_param('pinfplas',    4789425947.72,  'pinf for stiffend gas/plastic') 
    probdata.add_param('pinfwat',      300000000.0,  'pinf for stiffend water')
    
    # Density at rest and at room temperature of: air, plastic and water
    probdata.add_param('rhog',     1.0,  'air density in kg/m^3')
    probdata.add_param('rhop',     1050.0,  'polystirene density in kg/m^3')
    probdata.add_param('rhow',     1000.0,  'water density in kg/m^3')

	# omegas, omeplas and omewat are not used in this code, but they can be used to extend the code 
    # to a more general EOS with an extra parameter, like de Van der Waals EOS.
    probdata.add_param('omegas',     0.0,  'omega (specific excluded volume) for stiffend gas/plastic')
    probdata.add_param('omeplas',     0.0,  'omega (specific excluded volume) for stiffend gas/plastic')
    probdata.add_param('omewat',     0.0,  'omega (specific excluded volume) for stiffend water')
    
    # Parameters for mapped grid (if changed, they also need to be changed in mapc2p.py, also
    # double check the mapc2p.py test works correctly with the new parameters) 
    probdata.add_param('rsqr',     0.039999,  'radius of outer square for mapped grid')
    probdata.add_param('rout',     0.015,  'radius of outer circle of ring inclusion for mapped grid')
    probdata.add_param('rinn',     0.010,  'radius of inner circle of ring inclusion for mapped grid')
    
    #------------------------------------------------------------------
    # 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] = -0.05 #-0.055          # xlower
    clawdata.upper[0] = 0.05 #0.055          # xupper
    clawdata.lower[1] = 0.000000e+00          # ylower
    clawdata.upper[1] = 0.05000e+00          # yupper
    
    # Number of grid cells:
    # For original mapped grid, used multiples of 20x10, so the interface is aligned
    clawdata.num_cells[0] = 40 #40 #40 #56-mymapping #40-randys     # mx
    clawdata.num_cells[1] = 20 #20 #14 #17-mymapping #14-randys      # my
    

    # ---------------
    # Size of system:
    # ---------------

    # Number of equations in the system:
    clawdata.num_eqn = 4

    # Number of auxiliary variables in the aux array (initialized in setaux)
    clawdata.num_aux = 7
    
    # 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.chkNNNNN' specified below should be in 
    # the OUTDIR indicated in Makefile.

    clawdata.restart = False               # True to restart from prior results
    clawdata.restart_file = 'fort.chk00006'  # 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 = 150
        clawdata.tfinal = 0.0002 #0.00025 #0.00015 #0.00015
        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 = 1
        clawdata.total_steps = 4
        clawdata.output_t0 = True  # output at initial (or restart) time?
        

    clawdata.output_format == 'ascii'      # 'ascii' or 'netcdf' 

    clawdata.output_q_components = 'all'   # could be list such as [True,True]
    clawdata.output_aux_components = 'all'  # could be list
    clawdata.output_aux_onlyonce = False    # 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.
    clawdata.verbosity = 0
    
    

    # --------------
    # 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-07
    
    # 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.5 #0.600000
    # max Courant number to allow without retaking step with a smaller dt:
    clawdata.cfl_max = 0.6 #0.700000
    
    # Maximum number of time steps to allow between output times:
    clawdata.steps_max = 500000


    # ------------------
    # Method to be used:
    # ------------------

    # Order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
    clawdata.order = 2
    
    # Use dimensional splitting? 
    clawdata.dimensional_split = 'unsplit' # 'godunov' #'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 = 1 #2
    
    
    # Number of waves in the Riemann solution:
    clawdata.num_waves = 3
    
    # 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 (see inline limiter.f to
    # enable or disable modified minmod (TVD, only first order, "add viscosity")
    clawdata.limiter = [1, 1, 1] #[1, 1, 1]
    
    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 = 1
    
    
    # --------------------
    # 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] = 'user'   # at xlower
    clawdata.bc_upper[0] = 'extrap'   # at xupper

    clawdata.bc_lower[1] = 'extrap'   # at ylower
    clawdata.bc_upper[1] = 'extrap'   # at yupper
                  
                  
    # --------------
    # Checkpointing:
    # --------------

    # Specify when checkpoint files should be created that can be
    # used to restart a computation.

    clawdata.checkpt_style = 0

    if clawdata.checkpt_style == 0:
      # Do not checkpoint at all
      pass

    elif clawdata.checkpt_style == 1:
      # Checkpoint only at tfinal.
      pass

    elif clawdata.checkpt_style == 2:
      # Specify a list of checkpoint times.
      ckouttime0 = 50.00/1000000.0
      ckouttime1 = 63.33/1000000.0
      clawdata.checkpt_times = [ckouttime0, ckouttime1]

    elif clawdata.checkpt_style == 3:
      # Checkpoint every checkpt_interval timesteps (on Level 1)
      # and at the final time.
      clawdata.checkpt_interval = 5

       
    # ---------------
    # Gauges:
    # ---------------
    gauges = rundata.gaugedata.gauges
    # for gauges append lines of the form  [gaugeno, x, y, t1, t2]
    gauges.append([0, -0.01, 0, 0., 1e9])
    gauges.append([1, -0.01, 0.005, 0., 1e9])
    gauges.append([2, -0.01, 0.01, 0., 1e9])
    gauges.append([3, 0.0, 0, 0., 1e9])
    gauges.append([4, 0.0, 0.005, 0., 1e9])
    gauges.append([5, 0.0, 0.01, 0., 1e9])
    gauges.append([6, 0.01, 0, 0., 1e9])
    gauges.append([7, 0.01, 0.005, 0., 1e9])
    gauges.append([8, 0.01, 0.01, 0., 1e9])

    # ---------------
    # AMR parameters:
    # ---------------
    amrdata = rundata.amrdata

    # max number of refinement levels:
    amrdata.amr_levels_max = 4

    # List of refinement ratios at each level (length at least amr_level_max-1)
    amrdata.refinement_ratios_x = [2, 2, 2, 2, 2]
    amrdata.refinement_ratios_y = [2, 2, 2, 2, 2]
    amrdata.refinement_ratios_t = [2, 2, 2, 2, 2]


    # Specify type of each aux variable in amrdata.auxtype.
    # This must be a list of length num_aux, each element of which is one
    # of:
    #   'center',  'capacity', 'xleft', or 'yleft'  (see documentation).
    amrdata.aux_type = ['center','center','center','center','center','center','center']


    # Flag for refinement based on Richardson error estimater:
    amrdata.flag_richardson = False    # use Richardson?
    amrdata.flag_richardson_tol = 1.000000e-00  # Richardson tolerance
    
    # Flag for refinement using routine flag2refine:
    amrdata.flag2refine = True     # use this?
    amrdata.flag2refine_tol = 5000.0 #10000.0 #10000.0 #10.000000 #100000.0 #5.000000e-02  # tolerance used in this routine (100000000.000000)
    # User can modify flag2refine to change the criterion for flagging.
    # Default: check maximum absolute difference of first component of q
    # between a cell and each of its neighbors.

    # steps to take on each level L between regriddings of level L+1:
    amrdata.regrid_interval = 2       

    # width of buffer zone around flagged points:
    # (typically the same as regrid_interval so waves don't escape):
    amrdata.regrid_buffer_width  = 2 #3

    # clustering alg. cutoff for (# flagged pts) / (total # of cells
    # refined)
    # (closer to 1.0 => more small grids may be needed to cover flagged
    # cells)
    amrdata.clustering_cutoff = 0.700000

    # print info about each regridding up to this level:
    amrdata.verbosity_regrid = 0      


    # -------------------
    # Refinement Regions:
    # -------------------
    regions = rundata.regiondata.regions
    
    # Remove initial spurious wave from interface coupling by not refining until t_0
    t_0 = 1.67/1000000.0
    
    # NOT REQUIRED IF MAPPING NOT PRESENT
    
    # All of the water region (for level 2 with mapping)
    #regions.append([2,2,0,1e9,-0.0155,0.0155, 0.0, 0.0155])
    
    # All of the water region (for level 3 with mapping)
    #regions.append([3,3,0,1e9,-0.0155,0.0155, 0.0, 0.0155])
    
    # Regions along interface (for level 3 with mapping)
    #regions.append([3,3,0,1e9,-0.0155,-0.0145, 0.0, 0.0155])
    #regions.append([3,3,0,1e9,0.0145,0.0155, 0.0, 0.0155])
    #regions.append([3,3,0,1e9,-0.0155, 0.0155, 0.0145, 0.0155])
    
    # Regions along interface (for level 4 with mapping)
    #regions.append([4,4,0,1e9,-0.0155,0.0155, 0.0, 0.0155])
    
    # Regions along interface (for level 5 with mapping)
    #regions.append([5,5,0,1e9,-0.0155,0.0155, 0.0, 0.0155])
    
    # Regions along interface (for level 5 with mapping)
    #regions.append([5,5,0,1e9,-0.02,0.02, 0.0, 0.02])
    
    # Regions along interface (for level 6 with mapping)
    #regions.append([6,6,0,1e9,-0.02,0.02, 0.0, 0.0155])
    
    # Regions along interface (for level 6 without mapping)
    #regions.append([5,6,t_0,1e9,-0.025,0.025, 0.0, 0.025])
    
    # Along one corner for Schlieren
    ##regions.append([6,6,t_0,1e9,-0.02,0.005, 0.005, 0.02])
    ##regions.append([6,6,t_0,1e9,-0.02,0.02, -0.02, 0.02])
    ##regions.append([6,6,t_0,1e9,-0.02,0.02, -0.02, 0.02])
    ##regions.append([5,5,t_0,1e9,-0.01,0.02, 0.01, 0.02])
    #regions.append([5,5,t_0,1e9,-0.02,-0.01, 0.01, 0.02])
    
    #OTHER COMBINATIONS
    # Interface corners
    #regions.append([3,3,0,1e9,-0.0155,-0.0145, 0.0145, 0.0155])
    #regions.append([3,3,0,1e9,0.0145,0.0155, 0.0145, 0.0155])

    #regions.append([2,3,4e-6,1e9,-0.03,-0.029, 0.0, 0.030])
    #regions.append([3,3,0,1e9,-0.02,0.02, 0.0, 0.02])
    # to specify regions of refinement append lines of the form
    #  [minlevel,maxlevel,t1,t2,x1,x2,y1,y2]


    #  ----- For developers ----- 
    # Toggle debugging print statements:
    amrdata.dprint = False      # print domain flags
    amrdata.eprint = False      # print err est flags