Changeset 8466


Ignore:
Timestamp:
Jul 11, 2012, 9:42:10 PM (13 years ago)
Author:
davies
Message:

Adding new 'okada_tsunami' routines, which are well tested
and can treat both simple and complex faults

Location:
trunk
Files:
9 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/shallow_water/eqf.py

    r6157 r8466  
    22# earthquake_tsunami function
    33#
    4 
     4#import okada
    55"""This function returns a callable object representing an initial water
    66   displacement generated by a submarine earthqauke.

  • trunk/anuga_core/source/anuga/shallow_water/eqf_v2.py

    r7769 r8466  
    3737    from math import sin, radians
    3838
     39    raise Exception, 'WARNING: (GD, 11/07/2012) there seem to be bugs in eqf_v2.py\
     40    , I suggest you use okada_tsunami'
     41
    3942    if domain is not None:
    4043        xllcorner = domain.geo_reference.get_xllcorner()
    … …  
    156159
    157160        for i in range(N-1):
    158             self.SRECTF(alp, xr[i]*.001, yr[i]*.001, depth*.001, zero, length,\
     161            # GD CHANGE
     162            self.SRECTF(alp, xr[i]*.001, yr[i]*.001, depth, zero, length,\
    159163                        zero, width, sd, cd, disl1, disl2, disl3)
     164            #self.SRECTF(alp, xr[i]*.001, yr[i]*.001, depth*0.001, zero, length,\
     165            #            zero, width, sd, cd, disl1, disl2, disl3)
    160166           
    161167            z2[i] = self.U3
    … …  
    281287
    282288        if Q != F0:
    283             TT = atan( radians( XI*ET/(Q*R) ))
     289            #TT = atan( radians( XI*ET/(Q*R) ))
     290            # GD CHANGE
     291            TT = atan( XI*ET/(Q*R) )
    284292        else:
    285293            TT = F0
    … …  
    299307        if CD == 0:
    300308            #C==============================     
    301             #C=====   INCLINED FAULT   =====     
     309            #C=====   VERTICAL FAULT   =====     
    302310            #C==============================
    303311            RD2=RD*RD
    … …  
    312320        else:
    313321            #C==============================   
    314             #C=====   VERTICAL FAULT   =====     
     322            #C=====   INCLINED FAULT   =====     
    315323            #C==============================
    316324            TD=SD/CD                                                         
    … …  
    319327                A5=F0
    320328            else:                                                   
    321                 A5= ALP*F2/CD*atan( radians((ET*(X+Q*CD)+X*(R+X)*SD) / \
     329                #A5= ALP*F2/CD*atan( radians((ET*(X+Q*CD)+X*(R+X)*SD) / \
     330                #                    (XI*(R+X)*CD) ))   
     331                # GD CHANGE
     332                A5= ALP*F2/CD*atan( ((ET*(X+Q*CD)+X*(R+X)*SD) / \
    322333                                    (XI*(R+X)*CD) ))   
    323334            A4= ALP/CD*(  log(RD) - SD*DLE )                                 

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

    r8450 r8466  
    493493                //Bedslope approx 1:
    494494                //bedslope_work = g*hc*(ql[0])*length-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length;
     495                bedslope_work = g*length*( hc*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
    495496                //
    496497                // Bedslope approx 2
    … …  
    505506
    506507                // Bedslope approx 3
    507                 bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
     508                //bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
    508509                //
    509510                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
    … …  
    535536                    // Bedslope approx 1:
    536537                    //bedslope_work = g*hc_n*(qr[0])*length-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
     538                    bedslope_work = g*length*(hc_n*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
    537539                    //
    538540                    // Bedslope approx 2:
    … …  
    547549                    //
    548550                    // Bedslope approx 3
    549                     bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
     551                    //bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
    550552                    //
    551553                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
    … …  
    10791081      dqv[1] = a*dxv1 + b*dyv1;
    10801082      dqv[2] = a*dxv2 + b*dyv2;
    1081      
     1083   
     1084      // Idea: New limiter, computed using only neighbouring centroid values and local centroid value
     1085      // Step1: Compute the value of stage at the centroid of triangle k ,
     1086      //       using only the neighbouring centroid stage values
     1087      //tmp = a*(x-x0) + b*(y-y0) + stage_centroid_values[k0];
     1088      // Step2: Now, compute a limiting factor, so that if gradients are multiplied by this,
     1089      //        then for a triangle with these limited gradients, based at the local centroid value
     1090      //        ,the re-constructed neighbour centroid values will not over shoot (if larger) or undershoot (if smaller)
     1091      //limiting_factor=1 - max_all_k0((tmp - stage_centroid_values[k])/(tmp - stage_centroid_values[k0]))
     1092      //limiting_factor=min(max(limiting_factor,0),1)
     1093      // Advantages: No limiting for linear problem. Worth a look anyway
     1094 
    10821095      // Now we want to find min and max of the centroid and the
    10831096      // vertices of the auxiliary triangle and compute jumps
    … …  
    10881101   
    10891102      // Limit the gradient
     1103      //if(hmin/hc<0.5){
    10901104      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1105      //}
    10911106      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
    10921107     

  • trunk/anuga_work/development/gareth/tests/dam_break/plotme.py

    r8446 r8466  
    99p2_st=util.get_centroids(p_st)
    1010
    11 p_dev = util.get_output('dam_break_20120517_003225/dam_break.sww', 0.001)
     11p_dev = util.get_output('dam_break_20120627_231618/dam_break.sww', 0.001)
    1212p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
    1313

  • trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

    r8446 r8466  
    7474# Associate boundary tags with boundary objects
    7575#----------------------------------------------
    76 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
     76domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
    7777
    7878#------------------------------

  • trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_2plot.py

    r8354 r8466  
    88import numpy
    99import scipy
    10 import pylab
    11 
    12 import util
     10from matplotlib import pyplot as pyplot
     11from anuga.utilities import plot_utils as util
    1312#--------------
    1413# Get variables
    … …  
    2625# Define variables of case study
    2726#-------------------------------
    28 #mann=0.01 # Manning's coef
    29 #bedslope=-0.01
    30 #fluxin=0.1 #The momentum flux at the upstream boundary
     27mann=0.03 # Manning's coef
     28bedslope=-0.1
     29fluxin=20./100. #The momentum flux at the upstream boundary ( = discharge / width)
    3130
    3231#---------------------------------------------
    … …  
    3635# 1) Friction_slope=bedslope , and
    3736# 2) flux = depth*velocity = flux_in
    38 #uana= ( mann**(-2.)*abs(bedslope)*fluxin**(4./3.) )**(3./10.) # Velocity
    39 #dana= fluxin/uana # Depth
     37uana= ( mann**(-2.)*abs(bedslope)*fluxin**(4./3.) )**(3./10.) # Velocity
     38dana= fluxin/uana # Depth
    4039
    4140#--------------------
    4241# Make plot animation
    4342#--------------------
    44 pylab.close() #If the plot is open, there will be problems
    45 pylab.ion()
     43pyplot.close() #If the plot is open, there will be problems
     44pyplot.ion()
    4645
    47 line, = pylab.plot( (p2.x[v].min(),p2.x[v].max()) ,( (p2.stage[:,v]-p2.elev[:,v]).max(),(p2.stage[:,v]-p2.elev[v]).min() ) )
    48 #line.set_label('Numerical')
    49 #pylab.plot( (0,100),(dana,dana), 'r',label='Analytical' )
    50 #pylab.legend()
    51 #pylab.plot(x[v],elev[v],'green')
     46line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,( (p2.stage[:,v]-p2.elev[:,v]).max(),(p2.stage[:,v]-p2.elev[v]).min() ) )
     47line.set_label('Numerical')
     48pyplot.plot( (0,100),(dana,dana), 'r',label='Analytical' )
     49pyplot.legend()
     50#pyplot.plot(x[v],elev[v],'green')
    5251for i in range(p2.xmom.shape[0]):
    5352    line.set_xdata(p2.x[v])
    5453    line.set_ydata(p2.stage[i,v]-p2.elev[v])
    55     pylab.draw()
    56     pylab.title('Flow depth: should be converging to steady uniform flow ' )
    57     pylab.xlabel('x')
    58     pylab.ylabel('depth (m)')
     54    pyplot.draw()
     55    pyplot.title('Flow depth: should be converging to steady uniform flow ' )
     56    pyplot.xlabel('Distance down-slope (m)')
     57    pyplot.ylabel('depth (m)')
    5958
    60 pylab.title('Flow depth at 400s -- there should be a central region with steady uniform flow ' )
    61 pylab.savefig('final_depth_v2.png')
     59pyplot.title('Flow depth at 400s -- there should be a central region with steady uniform flow ' )
     60pyplot.savefig('final_depth_v2.png')
    6261
    6362#--------------------------------------------
    6463# Compare x-velocity with analytical solution
    6564#--------------------------------------------
    66 #pylab.clf()
    67 #pylab.plot(x[v],xvel[200,v], label='Numerical')
    68 #pylab.plot((0,100),(uana,uana),'r',label='Analytical')
    69 #pylab.legend()
    70 #pylab.xlabel('x')
    71 #pylab.ylabel('Velocity m/s')
     65#pyplot.clf()
     66#pyplot.plot(x[v],xvel[200,v], label='Numerical')
     67#pyplot.plot((0,100),(uana,uana),'r',label='Analytical')
     68#pyplot.legend()
     69#pyplot.xlabel('x')
     70#pyplot.ylabel('Velocity m/s')
    7271#plottitle = "Final velocity in the x-direction -- some wiggles. \n Can this be improved with\
    7372# a different downstream boundary condition?" #\n Setting depth only at the downstream boundary does not help "
    74 #pylab.title(plottitle)
    75 #pylab.savefig('SU_x_velocity.png')
     73#pyplot.title(plottitle)
     74#pyplot.savefig('SU_x_velocity.png')
     75
     76
     77# Find an x value close to x==50
     78tmp=(abs(p2.x-50.)).argmin()
     79v=(p2.x==p2.x[tmp])
    7680
    7781#--------------------------------------------------
    7882# Check if y velocity is zero everywhere (it should be)
    7983#--------------------------------------------------
    80 pylab.clf()
    81 #v=(x==50.)
    82 #pylab.clf()
    83 pylab.plot(p2.y[v],p2.yvel[100,v], 'rs-')
    84 ##pylab.plot(y[v],yvel[200,v], 'bp--')
    85 pylab.xlabel('y')
    86 pylab.ylabel('Velocity m/s')
    87 pylab.title('Final velocity in the y-direction (i.e. perpendicular to the channel cross-section) \n at x=50. (similar for other x) ')
    88 pylab.savefig('y_velocity_v2.png')
     84pyplot.clf()
     85pyplot.plot(p2.y[v],p2.yvel[100,v],'o', label='Numerical')
     86pyplot.plot((0,100),(0.0, 0.0),label='Analytical')
     87pyplot.xlabel('Distance along the line x==50 (across the slope) (m)')
     88pyplot.ylabel('Velocity m/s')
     89pyplot.title('Final y-velocity near x=50.')
     90pyplot.legend()
     91pyplot.savefig('y_velocity_v2.png')
    8992
    9093
    91 pylab.clf()
    92 v=(p2.x==50.)
    93 #pylab.clf()
    94 pylab.plot(p2.y[v],p2.xvel[100,v], 'rs-')
    95 ##pylab.plot(y[v],yvel[200,v], 'bp--')
    96 pylab.xlabel('y')
    97 pylab.ylabel('Velocity m/s')
    98 pylab.title('Final velocity in the x-direction (downstream) \n at x=50. (similar for other x) ')
    99 pylab.savefig('x_velocity_v2.png')
     94pyplot.clf()
     95pyplot.plot(p2.y[v],p2.xvel[100,v], 'o', label='Numerical')
     96pyplot.plot((0,100),(uana,uana),label='Analytical')
     97pyplot.ylim([1.0,3.0])
     98pyplot.xlabel('Distance along the line x==50 (across the slope) (m)')
     99pyplot.ylabel('Velocity m/s')
     100pyplot.title('Final velocity in the x-direction (downstream) \n at x=50. (similar for other x) ')
     101pyplot.legend()
     102pyplot.savefig('x_velocity_v2.png')
     103pyplot.clf()
    100104
    101 pylab.clf()
    102 
Note: See TracChangeset for help on using the changeset viewer.