Changeset 3850


Ignore:
Timestamp:
Oct 24, 2006, 4:56:49 PM (19 years ago)
Author:
ole
Message:

Introduced a covariance measure to verify Okushiri timeseries plus
some incidental work

Files:
11 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r3846 r3850  
    641641        """
    642642        if name.endswith('.sww'):
    643             name = name[-4:]
     643            name = name[:-4]
    644644           
    645645        self.simulation_name = name

  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r3700 r3850  
    235235
    236236    #Get variables
    237     if verbose: print 'Get variables'   
     237    #if verbose: print 'Get variables'   
    238238    time = fid.variables['time'][:]   
    239239
    … …  
    318318                                  triangles,
    319319                                  interpolation_points,
    320                                   verbose = verbose)
     320                                  verbose=verbose)
    321321
    322322

  • anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

    r3689 r3850  
    7373              a mesh origin, since geospatial has its own mesh origin.
    7474        """
    75        
     75
    7676        #Convert input to Numeric arrays
    7777        triangles = ensure_numeric(triangles, Int)
    … …  
    7979                                             geo_reference = mesh_origin)
    8080
    81         #Don't pass geo_reference to mesh.  It doesn't work.
     81        #Don't pass geo_reference to mesh.  It doesn't work. (Still??)
     82       
     83        if verbose: print 'FitInterpolate: Building mesh'       
    8284        self.mesh = Mesh(vertex_coordinates, triangles)
    8385        self.mesh.check_integrity()
     86       
     87        if verbose: print 'FitInterpolate: Building quad tree'
    8488        self.root = build_quadtree(self.mesh,
    8589                                   max_points_per_cell = max_vertices_per_cell)

  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r3768 r3850  
    405405                 time,
    406406                 quantities,
    407                  quantity_names = None, 
    408                  vertex_coordinates = None,
    409                  triangles = None,
    410                  interpolation_points = None,
    411                  verbose = False):
     407                 quantity_names=None, 
     408                 vertex_coordinates=None,
     409                 triangles=None,
     410                 interpolation_points=None,
     411                 verbose=False):
    412412        """Initialise object and build spatial interpolation if required
    413413        """
    … …  
    417417
    418418
    419         #from anuga.abstract_2d_finite_volumes.util import mean, ensure_numeric
    420419        from anuga.config import time_format
    421420        import types
    422421
    423422
    424         #Check temporal info
     423        # Check temporal info
    425424        time = ensure_numeric(time)       
    426425        msg = 'Time must be a monotonuosly '
    … …  
    429428
    430429
    431         #Check if quantities is a single array only
     430        # Check if quantities is a single array only
    432431        if type(quantities) != types.DictType:
    433432            quantities = ensure_numeric(quantities)
    … …  
    438437
    439438
    440         #Use keys if no names are specified
     439        # Use keys if no names are specified
    441440        if quantity_names is None:
    442441            quantity_names = quantities.keys()
    443442
    444443
    445         #Check spatial info
     444        # Check spatial info
    446445        if vertex_coordinates is None:
    447446            self.spatial = False
    … …  
    454453
    455454             
    456         #Save for use with statistics
    457 
     455        # Save for use with statistics
    458456        self.quantities_range = {}
    459457        for name in quantity_names:
    … …  
    462460       
    463461        self.quantity_names = quantity_names       
    464         #self.quantities = quantities  #Takes too much memory     
    465462        self.vertex_coordinates = vertex_coordinates
    466463        self.interpolation_points = interpolation_points
    … …  
    470467       
    471468           
    472         #Precomputed spatial interpolation if requested
     469        # Precomputed spatial interpolation if requested
    473470        if interpolation_points is not None:
    474471            if self.spatial is False:
    … …  
    491488                self.precomputed_values[name] = zeros((p, m), Float)
    492489
    493             #Build interpolator
     490            # Build interpolator
     491            if verbose: print 'Building interpolation matrix'           
    494492            interpol = Interpolate(vertex_coordinates,
    495493                                   triangles,
    … …  
    497495                                   #self.interpolation_points,
    498496                                   #alpha = 0,
    499                                    verbose = verbose)
     497                                   verbose=verbose)
    500498
    501499            if verbose: print 'Interpolate'

  • anuga_core/source/anuga/fit_interpolate/search_functions.py

    r2201 r3850  
    88from Numeric import dot
    99
     10from numerical_tools import get_machine_precision
    1011
    1112def search_tree_of_vertices(root, mesh, x):
    … …  
    8384            sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
    8485
    85             #FIXME: Maybe move out to test or something
    86             epsilon = 1.0e-6
    87             assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
     86            # Integrity check - machine precision is too hard, but at
     87            # least check that we are within a couple of orders of magnitude
     88            epsilon = get_machine_precision()
     89            msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, 100*eps = %.15e'\
     90                  %(abs(sigma0+sigma1+sigma2-1), 100*epsilon)
     91            assert abs(sigma0 + sigma1 + sigma2 - 1.0) < 100*epsilon, msg
    8892
    8993            #Check that this triangle contains the data point
    9094           
    91             #Sigmas can get negative within
    92             #machine precision on some machines (e.g nautilus)
    93             #Hence the small eps
    94             eps = 1.0e-15
    95             if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
     95            # Sigmas are allowed to get negative within
     96            # machine precision on some machines (e.g nautilus)
     97            if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
    9698                element_found = True
    9799                break

  • anuga_core/source/anuga/load_mesh/loadASCII.py

    r3603 r3850  
    7676#  Needs to be defined
    7777##FIXME (DSG-DSG) if the ascii file isn't .xya give a better error message.
     78
     79# FIXME (Ole): Has this stuff been superseded by geospatial data?
     80# Much of this code is also there
    7881
    7982

  • anuga_core/source/anuga/utilities/numerical_tools.py

    r3849 r3850  
    1010#(this should move to somewhere central)
    1111#try:
    12 #    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float, arange   
     12#    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt,
     13# searchsorted, sort, concatenate, Float, arange   
    1314#except:
    1415#    #print 'Could not find scipy - using Numeric'
    15 #    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float, arange
    16 from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float, arange   
     16#    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt,
     17#searchsorted, sort, concatenate, Float, arange
     18
     19from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt,\
     20     searchsorted, sort, concatenate, Float, arange   
    1721
    1822# Getting an infinite number to use when using Numeric
    … …  
    3135    """
    3236
    33     error_msg = 'Input to acos is outside allowed domain [-1.0, 1.0]. I got %.12f' %x
     37    error_msg = 'Input to acos is outside allowed domain [-1.0, 1.0].'+\
     38                'I got %.12f' %x
    3439    warning_msg = 'Changing argument to acos from %.18f to %.1f' %(x, sign(x))
    3540
    36     # FIXME(Ole): Need function to compute machine precision as this is also
    37     # used in fit_interpolate/search_functions.py
    38    
    39    
    40     eps = 1.0e-15  # Machine precision
     41    eps = get_machine_precision() # Machine precision
    4142    if x < -1.0:
    4243        if x < -1.0 - eps:
    … …  
    5051            raise ValueError, errmsg
    5152        else:
    52             print 'NOTE: changing argument to acos from %.18f to 1.0' %x           
     53            print 'NOTE: changing argument to acos from %.18f to 1.0' %x
    5354            x = 1.0
    5455
    … …  
    148149        y = x
    149150
     151    x = ensure_numeric(x)
     152    y = ensure_numeric(y)       
    150153    assert(len(x)==len(y))
    151     N = len(x)
    152  
     154
     155    N = len(x)
    153156    cx = x - mean(x) 
    154157    cy = y - mean(y) 

  • anuga_validation/okushiri_2005/create_okushiri.py

    r3845 r3850  
    1616    """
    1717
    18     print 'Preparing time boundary from %s' %filename
     18    textversion = filename[:-4] + '.txt'
     19
     20    print 'Preparing time boundary from %s' %textversion
    1921    from Numeric import array
    2022
    21     fid = open(filename)
     23    fid = open(textversion)
    2224
    2325    #Skip first line
    … …  
    4345    from Scientific.IO.NetCDF import NetCDFFile
    4446
    45     outfile = filename[:-4] + '.tms'
    46     print 'Writing to', outfile
    47     fid = NetCDFFile(outfile, 'w')
     47    print 'Writing to', filename
     48    fid = NetCDFFile(filename, 'w')
    4849
    4950    fid.institution = 'Geoscience Australia'
    … …  
    6970if __name__ == "__main__":
    7071
     72
    7173    #Preparing points
    7274    #(Only when using original Benchmark_2_Bathymetry.txt file)
     75    # FIXME: Replace by using geospatial data
     76    #
    7377    #from anuga.abstract_2d_finite_volumes.data_manager import xya2pts
    7478    #xya2pts(project.bathymetry_filename, verbose = True,

  • anuga_validation/okushiri_2005/extract_timeseries.py

    r3845 r3850  
    55"""
    66
     7from Numeric import allclose
     8from Scientific.IO.NetCDF import NetCDFFile
     9
     10from anuga.abstract_2d_finite_volumes.util import file_function
     11from anuga.utilities.numerical_tools import\
     12     ensure_numeric, cov, get_machine_precision
     13
     14import project
     15
     16try:
     17    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
     18except:
     19    plotting = False
     20else:
     21    plotting = True
    722
    823
    9 from anuga.caching import cache
    10 
    11 
    12 
    13 
    14 gauges = [[0.000, 1.696]]  #Boundary
    15 gauges += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
    16 
    17 gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
    18 
     24#-------------------------
     25# Basic data
     26#-------------------------
    1927
    2028finaltime = 22.5
    2129timestep = 0.05
    2230
    23 #Read reference data
     31gauge_locations = [[0.000, 1.696]] # Boundary gauge
     32gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
     33gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
    2434
    25 #Input wave
    26 filename = 'Benchmark_2_input.tms'
    27 print 'Reading', filename
    28 from Scientific.IO.NetCDF import NetCDFFile
    29 from Numeric import allclose
    30 fid = NetCDFFile(filename, 'r')#
     35validation_data = {}
     36for key in gauge_names:
     37    validation_data[key] = []
    3138
    32 input_time = fid.variables['time'][:]
    33 input_stage = fid.variables['stage'][:]
     39   
     40expected_covariances = {'Boundary': 5.288392008865989e-05,
     41                        'ch5': 1.166748190444681e-04,
     42                        'ch7': 1.121816242516758e-04,
     43                        'ch9': 1.249543278366778e-04}
    3444
    3545
    36 #gauges
     46#-------------------------
     47# Read validation dataa
     48#-------------------------
     49
     50print 'Reading', project.boundary_filename
     51fid = NetCDFFile(project.boundary_filename, 'r')
     52input_time = fid.variables['time'][:]
     53validation_data['Boundary'] = fid.variables['stage'][:]
     54
    3755reference_time = []
    38 ch5 = []
    39 ch7 = []
    40 ch9 = []
    41 filename = 'output_ch5-7-9.txt'
    42 fid = open(filename)
     56fid = open(project.validation_filename)
    4357lines = fid.readlines()
    4458fid.close()
     59
    4560for i, line in enumerate(lines[1:]):
    4661    if i == len(input_time): break
    47 
     62   
    4863    fields = line.split()
    4964
    50     reference_time.append(float(fields[0]))
    51     ch5.append(float(fields[1])/100)   #cm2m
    52     ch7.append(float(fields[2])/100)   #cm2m
    53     ch9.append(float(fields[3])/100)   #cm2m
     65    reference_time.append(float(fields[0]))    # Record reference time
     66    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
     67        value = float(fields[1:][j])           # Omit time
     68        validation_data[key].append(value/100) # Convert cm2m
    5469
    5570
    56 from anuga.abstract_2d_finite_volumes.util import file_function
    57 from anuga.utilities.numerical_tools import ensure_numeric
    58 gauge_values = [ensure_numeric(input_stage),
    59                 ensure_numeric(ch5),
    60                 ensure_numeric(ch7),
    61                 ensure_numeric(ch9)] #Reference values
    62 
    63 
    64 
    65 #Read model output
    66 #filename = 'output.sww'
    67 filename = 'lwru2.sww'
    68 #filename = 'lwru2_.05s.sww'
    69 
    70 #f = file_function(filename,
    71 #                  quantities = 'stage',
    72 #                  interpolation_points = gauges,
    73 #                  verbose = True)
    74 
    75 f = cache(file_function,filename,
    76           {'quantities': 'stage',
    77            'interpolation_points': gauges,
    78            'verbose': True},
    79           #evaluate = True,
    80           dependencies = [filename],
    81           verbose = True)
    82 
    83 
    84 
    85 
    86 #Checks
    87 #print reference_time
    88 #print input_time
     71# Checks
    8972assert reference_time[0] == 0.0
    9073assert reference_time[-1] == finaltime
    9174assert allclose(reference_time, input_time)
    9275
     76for key in gauge_names:
     77    validation_data[key] = ensure_numeric(validation_data[key])
    9378
    9479
    95 #Validation
     80#--------------------------------------------------
     81# Read and interpolate model output
     82#--------------------------------------------------
     83#f = file_function('okushiri_new_limiters.sww',
     84#f = file_function('okushiri_as2005_with_mxspd=0.1.sww',
     85f = file_function(project.output_filename,
     86                  quantities='stage',
     87                  interpolation_points=gauge_locations,
     88                  use_cache=True,
     89                  verbose=True)
    9690
    9791
     92#--------------------------------------------------
     93# Compare model output to validation data
     94#--------------------------------------------------
     95
     96eps = get_machine_precision()
    9897for k, name in enumerate(gauge_names):
    9998    sqsum = 0
    10099    denom = 0
    101100    model = []
     101    print
    102102    print 'Validating ' + name
     103    observed_timeseries = validation_data[name]
    103104    for i, t in enumerate(reference_time):
    104         ref = gauge_values[k][i]
    105         val = f(t, point_id = k)[0]
    106         model.append(val)
     105        model.append(f(t, point_id=k)[0])
    107106
    108         sqsum += (ref - val)**2
    109         denom += ref**2
     107    res = cov(observed_timeseries, model)     
     108    print 'Covariance = %.18e' %res
    110109
    111     print sqsum
    112     print sqsum/denom
     110    if res < expected_covariances[name]-eps:
     111        print 'Result is better than expected by: %.18e'\
     112              %(res-expected_covariances[name])
     113        print 'Result = %.18e' %res
     114        print 'Expect = %.18e' %expected_covariances[name]   
     115    elif res > expected_covariances[name]+eps:
     116        print 'FAIL: %.18e' %res
     117    else:
     118        pass
    113119
    114     from pylab import *
    115     ion()
    116     hold(False)
    117     plot(reference_time, gauge_values[k], 'r.-',
    118          reference_time, model, 'k-')
    119     title('Gauge %s' %name)
    120     xlabel('time(s)')
    121     ylabel('stage (m)')   
    122     legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
    123     savefig(name, dpi = 300)
    124120
    125     raw_input('Next')
    126 show()
     121    if plotting is True:
     122        ion()
     123        hold(False)
     124   
     125        plot(reference_time, validation_data[name], 'r.-',
     126             reference_time, model, 'k-')
     127        title('Gauge %s' %name)
     128        xlabel('time(s)')
     129        ylabel('stage (m)')   
     130        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
     131        savefig(name, dpi = 300)       
     132
     133        raw_input('Next')
    127134
    128135
    129136
    130 #from pylab import *
    131 #plot(time, stage)
    132 #show()

  • anuga_validation/okushiri_2005/project.py

    r3845 r3850  
    22"""
    33
    4 boundary_filename = 'Benchmark_2_input.txt'
    5 bathymetry_filename = 'Benchmark_2_Bathymetry.txt'
     4# Given boundary wave
     5boundary_filename = 'Benchmark_2_input.tms'
     6
     7# Observed timeseries
     8validation_filename = 'output_ch5-7-9.txt'
     9
     10# A simple conversion from txt file to the more compact NetCDF format
     11bathymetry_filename = 'Benchmark_2_Bathymetry.pts'
     12
     13# Triangular mesh
    614mesh_filename = 'Benchmark_2.msh'
    715
    8 output_filename = 'okushiri.sww'
     16# Model output
     17output_filename = 'okushiri_old_limiters.sww'
    918
    1019

  • anuga_validation/okushiri_2005/run_okushiri.py

    r3845 r3850  
    4242domain.set_quantity('stage', 0.0)
    4343domain.set_quantity('elevation',
    44                     filename = project.bathymetry_filename[:-4] + '.pts',
    45                     alpha = 0.02,                   
    46                     verbose = True,
    47                     use_cache = True)
     44                    filename=project.bathymetry_filename,
     45                    alpha=0.02,                   
     46                    verbose=True,
     47                    use_cache=True)
    4848
    4949
    … …  
    5151# Domain parameters
    5252#-------------------------
    53 domain.set_name(project.output_filename)
    54 domain.set_default_order(2)
    55 domain.set_minimum_storable_height(0.001)
    56 domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
    57 
    58 # Set old (pre Sep 2006) defaults for limiters
    59 #domain.beta_w      = 0.9
    60 #domain.beta_w_dry  = 0.9
    61 #domain.beta_uh     = 0.9
    62 #domain.beta_uh_dry = 0.9
    63 #domain.beta_vh     = 0.9
    64 #domain.beta_vh_dry = 0.9
    65 
    66 #domain.check_integrity()
    67 #print domain.statistics()
    68 
     53domain.set_name(project.output_filename)  # Name of output sww file
     54domain.set_default_order(2)               # Apply second order scheme
     55domain.set_all_limiters(0.9)              # Maximal second order scheme (old lim)
     56domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
     57domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
    6958
    7059
    … …  
    7463
    7564# Create boundary function from timeseries provided in file
    76 function = file_function(project.boundary_filename[:-4] + '.tms',
     65function = file_function(project.boundary_filename,
    7766                         domain, verbose=True)
    7867
Note: See TracChangeset for help on using the changeset viewer.