source: trunk/anuga_validation/validation_tests/validation.py @ 8538

import sys
import anuga
from math import cos
from numpy import zeros, float
from time import localtime, strftime, gmtime


class Validation:
    """
    """

    def __init__(self, output_file, output_dir):
        """
        """
        self.output_file = output_file
        self.output_dir = output_dir
        
        
    def set_up_domain(self, dx, dy, length, width):
        """
        """
        points, vertices, boundary = anuga.rectangular_cross(int(length/dx), int(width/dy), length, width, (0.0, -width/2))

        self.domain = anuga.Domain(points, vertices, boundary) 

        self.domain.set_name(self.output_file)                
        self.domain.set_datadir(self.output_dir)  
        
    def set_up_algorithm(self, 
                         timestepping_method, 
                         default_order,
                         use_edge_limiter,
                         tight_slope_limiters,
                         use_centroid_velocities,
                         CFL,
                         beta_w,
                         beta_uh,
                         beta_vh):
        """
        """
        self.domain.set_timestepping_method(timestepping_method)
        self.domain.set_default_order(default_order)

        self.domain.use_edge_limiter = use_edge_limiter
        self.domain.tight_slope_limiters = tight_slope_limiters
        self.domain.use_centroid_velocities = use_centroid_velocities

        self.domain.CFL = CFL

        self.domain.beta_w = beta_w
        self.domain.beta_uh = beta_uh
        self.domain.beta_vh = beta_vh
        
    def set_initial_conditions(self,
                               elevation,
                               friction,
                               h0,
                               h1):
        """
        """
        self.domain.set_quantity('elevation', elevation)
        self.domain.set_quantity('friction', friction)

        def height(x,y):
            z = zeros(len(x), float)
            for i in range(len(x)):
                if x[i]<=50000.0:
                    z[i] = h0
                else:
                    z[i] = h1
            return z

        self.domain.set_quantity('stage', height)
        
    def set_boundary_conditions(self):
        """
        """
        from math import sin, pi, exp
        Br = anuga.Reflective_boundary(self.domain)      # Solid reflective wall
        Bt = anuga.Transmissive_boundary(self.domain)    # Continue all values on boundary 
        Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values

        # Associate boundary tags with boundary objects
        self.domain.set_boundary({'left': Bt, 'right': Bt, 'top': Br, 'bottom': Br})

    def evolve(self, yieldstep=None, finaltime=None, duration=None):
        """
        """
        for t in self.domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
            print t
        
        
if __name__ == "__main__":
    
    time = strftime('%Y%m%d_%H%M%S', localtime())

    output_dir = 'dam_break_' + time
    output_file = 'dam_break'

    #anuga.copy_code_files(output_dir,__file__)
    
    valid = Validation(output_file, output_dir)
    
    dx = 1000.
    dy = dx
    length = 100000.
    width = 10*dx
    
    valid.set_up_domain(dx, dy, length, width)

    CFL = 1.0
    beta_w = 0.6
    beta_uh = 0.6
    beta_vh = 0.6

    valid.set_up_algorithm(timestepping_method = 'rk2',
                           default_order =2,
                           use_edge_limiter = True,
                           tight_slope_limiters = True,
                           use_centroid_velocities = False,
                           CFL = CFL,
                           beta_w = beta_w,
                           beta_uh = beta_uh,
                           beta_vh = beta_vh)

    h0 = 10.0
    h1 = 1.0
    
    valid.set_initial_conditions(0, 0, h0, h1)
    valid.set_boundary_conditions()
    valid.evolve(yieldstep = 100.0, finaltime = 60*60.)