source: trunk/anuga_profile/profile_structure.py @ 8453

""" This is a script for profiling the code. It uses the structures code.

This example includes a Model Topography that shows a TYPICAL Headwall Configuration

The aim is to change the Culvert Routine to Model more precisely the abstraction
from a vertical face.

The inflow must include the impact of Approach velocity.
Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style
abstraction

"""
print 'Starting.... Importing Modules...'
#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
import anuga

from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
from anuga.shallow_water.shallow_water_domain import Domain
from anuga.shallow_water.forcing import Rainfall, Inflow

from anuga.culvert_flows.culvert_class import Culvert_flow
from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
from math import pi, pow, sqrt

import numpy as num


#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
print 'Setting up domain'

length = 200. #x-Dir
width = 200.  #y-dir

dx = dy = 2.5          # Resolution: Length of subdivisions on both axes

points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
                                               len1=length, len2=width)
domain = Domain(points, vertices, boundary)   
domain.set_name('profile_structure') # Output name
domain.set_default_order(2)
domain.H0 = 0.01
domain.tight_slope_limiters = 1

print 'Size', len(domain)

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------

def topography(x, y):
    """Set up a weir
    
    A culvert will connect either side
    """
    # General Slope of Topography
    z=10.0-x/100.0  # % Longitudinal Slope
    
    #       NOW Add bits and Pieces to topography
    bank_hgt=10.0
    bridge_width = 50.0
    bank_width = 10.0
    
    us_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
    us_start_x = 10.0
    top_start_y = 50.0
    us_slope = 3.0  #Horiz : Vertic
    ds_slope = 3.0
    ds_apron_skew = 1.0 # 1.0 = 1 Length: 1 Width, 2.0 = 2 Length : 1 Width
    centre_line_y= top_start_y+bridge_width/2.0

    # CALCULATE PARAMETERS TO FORM THE EMBANKMENT
    us_slope_length = bank_hgt*us_slope
    us_end_x =us_start_x + us_slope_length
    us_toe_start_y =top_start_y - us_slope_length / us_apron_skew
    us_toe_end_y = top_start_y + bridge_width + us_slope_length / us_apron_skew

    top_end_y = top_start_y + bridge_width
    ds_slope_length = bank_hgt*ds_slope
    ds_start_x = us_end_x + bank_width
    ds_end_x = ds_start_x + ds_slope_length

    ds_toe_start_y =top_start_y - ds_slope_length / ds_apron_skew
    ds_toe_end_y = top_start_y + bridge_width + ds_slope_length / ds_apron_skew


    N = len(x)
    for i in range(N):

       # Sloping Embankment Across Channel
        if us_start_x < x[i] < us_end_x +0.1:   # For UPSLOPE on the Upstream FACE
            if us_toe_start_y +(x[i] - us_start_x)/us_apron_skew < y[i] < us_toe_end_y - (x[i] - us_start_x)/us_apron_skew:
               z[i]=z[i] # Flat Apron
            else:
               z[i] += z[i] + (x[i] - us_start_x)/us_slope    # Set elevation for Sloping Segment  U/S Face
        if us_end_x < x[i] < ds_start_x + 0.1: # FOR The TOP of BANK Segment
           if top_start_y < y[i] < top_end_y:
               z[i]=z[i] # Flat Apron
           else:
               z[i] +=  z[i]+bank_hgt        # Flat Crest of Embankment
        if ds_start_x < x[i] < ds_end_x: # DOWN SDLOPE Segment on Downstream face
            if  top_start_y-(x[i]-ds_start_x)/ds_apron_skew <  y[i]  < top_end_y + (x[i]-ds_start_x)/ds_apron_skew: # Cut Out Segment for Culvert FACE
               z[i]=z[i] # Flat Apron
            else:
               z[i] += z[i]+bank_hgt-(x[i] -ds_start_x)/ds_slope       # Sloping D/S Face
           
        

    return z

print 'Setting Quantities....'
domain.set_quantity('elevation', topography)  # Use function for elevation
domain.set_quantity('friction', 0.01)         # Constant friction 
domain.set_quantity('stage',
                    expression='elevation')   # Dry initial condition

#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
print 'Setting Boundary Conditions'
Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!

Bo = anuga.Dirichlet_boundary([-5.0, 0, 0])           # Outflow water at -5.0
Bd = anuga.Dirichlet_boundary([0,0,0])                # Outflow water at 0.0

Btus = anuga.Dirichlet_boundary([18.0, 0, 0])           # Outflow water at 5.0
Btds = anuga.Dirichlet_boundary([0.0, 0, 0])           # Outflow water at 1.0
domain.set_boundary({'left': Btus, 'right': Btds, 'top': Br, 'bottom': Br})

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------

def ev():

    for t in domain.evolve(yieldstep = 1, finaltime = 20):
        print domain.timestepping_statistics()

import cProfile
import time
import os
import platform

if not os.path.isdir('profile_structure'):
    os.mkdir('profile_structure')

t = time.ctime().replace(' ', '_')
t = t.replace(':', '-')
uname_ = platform.uname()

os.mkdir(os.path.join('profile_structure', t))

f = open(os.path.join('profile_structure', t, 'uname'), 'w')
f.write(str(uname_))
f.close()

cProfile.run('ev()', os.path.join('profile_structure', t, 'profile_structure.prf'))

#import pstats
#p = pstats.Stats('profile_structure.prf')
#p.strip_dirs().sort_stats(-1).print_stats()
#p.sort_stats('cumulative').print_stats(10)