Changeset 9364


Ignore:
Timestamp:
Nov 7, 2014, 2:58:40 PM (10 years ago)
Author:
davies
Message:

Refactoring plot_utils sww reading to use less memory for small sets of timeSlices

Location:
trunk/anuga_core
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/utilities/plot_utils.py

    r9360 r9364  
    172172           
    173173    """
    174     var=copy.copy(varIn) # avoid python pointer issues
     174
    175175    if (len(varIn.shape)==2):
    176         # Get particular time-slices, unless the variable is constant
    177         # (e.g. elevation is often constant)
     176        # There are multiple time-slices
    178177        if timeSlices is 'max':
    179178            # Extract the maxima over time, assuming there are multiple
    180179            # time-slices, and ensure the var is still a 2D array
    181180            if( not absMax):
    182                 var=var.max(axis=0)
     181                var = (varIn[:]).max(axis=0, keepdims=True)
    183182            else:
    184183                # For variables xmom,ymom,xvel,yvel we want the 'maximum-absolute-value'
    185                 # We could do this everywhere, but I assume the loop is a bit slower
    186                 varInds=abs(var).argmax(axis=0)
    187                 varNew=varInds*0.
     184                varInds = abs(varIn[:]).argmax(axis=0)
     185                varNew = varInds*0.
    188186                for i in range(len(varInds)):
    189                     varNew[i] = var[varInds[i],i]
    190                 #var=[var[varInds[i],i] for i in varInds]
    191                 var=varNew
    192             var=var.reshape((1,len(var)))
     187                    varNew[i] = varIn[varInds[i],i]
     188                var = varNew
     189                var=var.reshape((1,len(var)))
    193190        else:
    194             var=var[timeSlices,:]
     191            var=varIn[timeSlices,:]
     192            var.reshape((len(timeSlices), varIn.shape[1]))
     193    else:
     194        # There is 1 time slice only
     195        var = varIn[:]
    195196   
    196197    return var
    … …  
    219220             everything else is stored every time step for vertices
    220221    """
    221 
    222     # Import modules
    223 
    224 
    225222
    226223    # Open ncdf connection
    … …  
    258255    y=fid.variables['y'][:]
    259256
    260     stage=getInds(fid.variables['stage'][:], timeSlices=inds)
    261     elev=getInds(fid.variables['elevation'][:], timeSlices=inds)
     257    stage=getInds(fid.variables['stage'], timeSlices=inds)
     258    elev=getInds(fid.variables['elevation'], timeSlices=inds)
    262259
    263260    # Simple approach for volumes
    … …  
    266263    # Friction if it exists
    267264    if(fid.variables.has_key('friction')):
    268         friction=getInds(fid.variables['friction'][:],timeSlices=inds)
     265        friction=getInds(fid.variables['friction'],timeSlices=inds)
    269266    else:
    270267        # Set friction to nan if it is not stored
    271         friction=elev*0.+numpy.nan
    272    
    273     #@ Here we get 'all' of height / xmom /ymom
    274     #@ This could be done using less memory/computation in
    275     #@  the case of multiple time-slices
    276 
     268        friction = elev*0.+numpy.nan
     269
     270    # Trick to treat the case where inds == 'max'
     271    inds2 = copy.copy(inds)
     272    if inds == 'max':
     273        inds2 = range(len(fid.variables['time']))
     274   
     275    # Get height
    277276    if(fid.variables.has_key('height')):
    278         heightAll=fid.variables['height'][:]
     277        height = fid.variables['height'][inds2]
    279278    else:
    280279        # Back calculate height if it is not stored
    281         heightAll=fid.variables['stage'][:]
    282         if(len(heightAll.shape)==len(elev.shape)):
    283             heightAll=heightAll-elev
     280        height = fid.variables['stage'][inds2]+0.
     281        if(len(elev.shape)==2):
     282            height = height-elev
    284283        else:
    285             for i in range(heightAll.shape[0]):
    286                 heightAll[i,:]=heightAll[i,:]-elev
    287     heightAll=heightAll*(heightAll>0.) # Height could be negative for tsunami algorithm
    288     # Need xmom,ymom for all timesteps
    289     xmomAll=fid.variables['xmomentum'][:]
    290     ymomAll=fid.variables['ymomentum'][:]
    291 
    292     height=getInds(heightAll, timeSlices=inds)
    293     # For momenta, we want maximum-absolute-value events
    294     xmom=getInds(xmomAll, timeSlices=inds, absMax=True)
    295     ymom=getInds(ymomAll, timeSlices=inds, absMax=True)
    296 
    297     # velocity requires some intermediate calculation in general
    298     hInv=1./(heightAll+1.0e-12)
    299     tmp = xmomAll*hInv*(heightAll>minimum_allowed_height)
    300     xvel=getInds(tmp,timeSlices=inds, absMax=True)
    301     tmp = ymomAll*hInv*(heightAll>minimum_allowed_height)
    302     yvel=getInds(tmp,timeSlices=inds, absMax=True)
    303     tmp = (xmomAll**2+ymomAll**2)**0.5*hInv*(heightAll>minimum_allowed_height)
    304     vel=getInds(tmp, timeSlices=inds) # Vel is always >= 0.
     284            for i in range(height.shape[0]):
     285                height[i,:] = height[i,:]-elev
     286    height = height*(height>0.)
     287
     288    # Get xmom
     289    xmom = fid.variables['xmomentum'][inds2]
     290    ymom = fid.variables['ymomentum'][inds2]
     291
     292    # Get vel
     293    h_inv = 1.0/(height+1.0e-12)
     294    hWet = (height > minimum_allowed_height)
     295    xvel = xmom*h_inv*hWet
     296    yvel = ymom*h_inv*hWet
     297    vel = (xmom**2 + ymom**2)**0.5*h_inv*hWet
     298
     299    if inds == 'max':
     300        height = height.max(axis=0, keepdims=True)
     301        vel = vel.max(axis=0, keepdims=True)
     302        xvel = getInds(xvel, timeSlices=inds,absMax=True)
     303        yvel = getInds(yvel, timeSlices=inds,absMax=True)
     304        xmom = getInds(xmom, timeSlices=inds,absMax=True)
     305        ymom = getInds(ymom, timeSlices=inds,absMax=True)
    305306
    306307    fid.close()
    … …  
    342343                         minimum_allowed_height=minimum_allowed_height,\
    343344                         verbose=verbose)
     345
     346def _getCentVar(fid, varkey_c, inds, absMax=False,  vols = None):
     347    """
     348        Convenience function used to get centroid variables from netCDF
     349        file connection fid
     350
     351        The default arguments fid, vols0, vols1, vols2 exist in the
     352        _get_centroid_values function where this is used
     353
     354    """
     355    vols0 = vols[:,0]
     356    vols1 = vols[:,1]
     357    vols2 = vols[:,2]
     358
     359    if(fid.variables.has_key(varkey_c)==False):
     360        # It looks like centroid values are not stored
     361        # In this case, compute centroid values from vertex values
     362
     363        newkey=varkey_c.replace('_c','')
     364        if inds is not 'max':
     365            # Relatively efficient treatment is possible
     366            var_cent = fid.variables[newkey]
     367            if (len(var_cent.shape)>1):
     368                # array contain time slices
     369                var_cent = fid.variables[newkey][inds]
     370                var_cent = (var_cent[:,vols0]+var_cent[:,vols1]+var_cent[:,vols2])/3.0
     371            else:
     372                var_cent = fid.variables[newkey][:]
     373                var_cent = (var_cent[vols0]+var_cent[vols1]+var_cent[vols2])/3.0
     374        else:
     375            # Requires reading all the data
     376            tmp = fid.variables[newkey][:]
     377            try: # array contain time slices
     378                tmp=(tmp[:,vols0]+tmp[:,vols1]+tmp[:,vols2])/3.0
     379            except:
     380                tmp=(tmp[vols0]+tmp[vols1]+tmp[vols2])/3.0
     381            var_cent=getInds(tmp, timeSlices=inds, absMax=absMax)
     382    else:
     383        if inds is not 'max':
     384            if(len(fid.variables[varkey_c].shape)>1):
     385                var_cent = fid.variables[varkey_c][inds]
     386            else:
     387                var_cent = fid.variables[varkey_c][:]
     388        else:
     389            var_cent=getInds(fid.variables[varkey_c][:], timeSlices=inds, absMax=absMax)
     390    return var_cent
     391
    344392                                 
    345393
    … …  
    444492    y_cent=(y[vols0]+y[vols1]+y[vols2])/3.0
    445493
    446     def getCentVar(varkey_c, timeSlices=inds, absMax=False):
    447         """
    448             Convenience function, assumes knowedge of 'timeSlices' and vols0,1,2
    449         """
    450         if(fid.variables.has_key(varkey_c)==False):
    451             # It looks like centroid values are not stored
    452             # In this case, compute centroid values from vertex values
    453            
    454             newkey=varkey_c.replace('_c','')
    455             tmp = fid.variables[newkey][:]
    456             try: # array contain time slides
    457                 tmp=(tmp[:,vols0]+tmp[:,vols1]+tmp[:,vols2])/3.0
    458             except:
    459                 tmp=(tmp[vols0]+tmp[vols1]+tmp[vols2])/3.0
    460             var_cent=getInds(tmp, timeSlices=timeSlices, absMax=absMax)
    461         else:
    462             var_cent=getInds(fid.variables[varkey_c][:], timeSlices=timeSlices, absMax=absMax)
    463         return var_cent
    464 
    465494    # Stage and height and elevation
    466     stage_cent=getCentVar('stage_c', timeSlices=inds)
    467     elev_cent=getCentVar('elevation_c', timeSlices=inds)
     495    stage_cent = _getCentVar(fid, 'stage_c', inds=inds, vols=vols)
     496    elev_cent = _getCentVar(fid, 'elevation_c', inds=inds, vols=vols)
    468497
    469498    if(len(elev_cent.shape)==2):
    … …  
    471500        elev_cent=elev_cent[0,:]
    472501
    473     height_cent=stage_cent*0.
    474     for i in range(stage_cent.shape[0]):
    475         height_cent[i,:]=stage_cent[i,:]-elev_cent
    476 
    477502    # Friction might not be stored at all
    478503    try:
    479         friction_cent=getCentVar('friction_c')
     504        friction_cent = _getCentVar(fid, 'friction_c', inds=inds, vols=vols)
    480505    except:
    481506        friction_cent=elev_cent*0.+numpy.nan
    482 
    483     if(fid.variables.has_key('xmomentum_c')):
    484         # Assume that both xmomentum,ymomentum are stored at centroids
    485         # Because velocity is back computed, and we might want maxima,
    486         # we get all data for convenience
    487         xmomC=getCentVar('xmomentum_c', timeSlices=range(nts))
    488         ymomC=getCentVar('ymomentum_c', timeSlices=range(nts))
    489 
    490         # height might not be stored
    491         try:
    492             hC = getCentVar('height_c', timeSlices=range(nts))
    493         except:
    494             # Compute from stage
    495             hC = getCentVar('stage_c', timeSlices=range(nts))
    496             for i in range(hC.shape[0]):
    497                 hC[i,:]=hC[i,:]-elev_cent
    498        
    499         hInv=1.0/(hC+1.0e-06)
    500         xmom_cent = getInds(xmomC*(hC>minimum_allowed_height), timeSlices=inds,absMax=True)
    501         xvel_cent = getInds(xmomC*hInv*(hC>minimum_allowed_height), timeSlices=inds, absMax=True)
    502 
    503         ymom_cent = getInds(ymomC*(hC>minimum_allowed_height), timeSlices=inds,absMax=True)
    504         yvel_cent = getInds(ymomC*hInv*(hC>minimum_allowed_height), timeSlices=inds, absMax=True)
    505 
    506         tmp = (xmomC**2 + ymomC**2)**0.5*hInv*(hC>minimum_allowed_height)
    507         vel_cent=getInds(tmp, timeSlices=inds)
    508        
    509     else:
    510         #@ COMPUTE CENTROIDS FROM VERTEX VALUES
    511         #@
    512         #@ Here we get 'all' of height / xmom /ymom
    513         #@ This could be done using less memory/computation in
    514         #@  the case of multiple time-slices
    515         if(fid.variables.has_key('height')):
    516             heightAll=fid.variables['height'][:]
     507   
     508    # Trick to treat the case where inds == 'max'
     509    inds2 = copy.copy(inds)
     510    if inds == 'max':
     511        inds2 = range(len(fid.variables['time']))
     512   
     513    # height
     514    height_cent= stage_cent + 0.
     515    for i in range(stage_cent.shape[0]):
     516        height_cent[i,:] = stage_cent[i,:] - elev_cent
     517
     518    if fid.variables.has_key('xmomentum_c'):
     519        # Momenta
     520        xmom_cent = fid.variables['xmomentum_c'][inds2]
     521        ymom_cent = fid.variables['ymomentum_c'][inds2]
     522
     523        # Height -- need to do this again incase inds == 'max'
     524        if fid.variables.has_key('height_c'):
     525            height_c_tmp = fid.variables['height_c'][inds2]
    517526        else:
    518             # Back calculate height if it is not stored
    519             heightAll=fid.variables['stage'][:]
    520             elev = fid.variables['elevation'][:]
    521             if(len(heightAll.shape)==len(elev.shape)):
    522                 heightAll=heightAll-elev
    523             else:
    524                 for i in range(heightAll.shape[0]):
    525                     heightAll[i,:]=heightAll[i,:]-elev
    526         heightAll=heightAll*(heightAll>0.) # Height could be negative for tsunami algorithm
    527         # Need xmom,ymom for all timesteps
    528         xmomAll=fid.variables['xmomentum'][:]
    529         ymomAll=fid.variables['ymomentum'][:]
    530        
     527            height_c_tmp = fid.variables['stage_c'][inds2]+0.
     528            for i in range(height_c_tmp.shape[0]):
     529                height_c_tmp[i,:] = height_c_tmp[i,:] - elev_cent
     530        # Vel
     531        hInv = 1.0/(height_c_tmp + 1.0e-12)
     532        hWet = (height_c_tmp > minimum_allowed_height)
     533        xvel_cent = xmom_cent*hInv*hWet
     534        yvel_cent = ymom_cent*hInv*hWet
     535
     536    else:
     537        # Get important vertex variables
     538        xmom_v = fid.variables['xmomentum'][inds2]
     539        ymom_v = fid.variables['ymomentum'][inds2]
     540        stage_v = fid.variables['stage'][inds2]
     541        elev_v = fid.variables['elevation']
     542        # Fix elevation + get height at vertices
     543        if (len(elev_v.shape)>1):
     544            elev_v = elev_v[inds2]
     545            height_v = stage_v - elev_v
     546        else:
     547            elev_v = elev_v[:]
     548            height_v = stage_v + 0.
     549            for i in range(stage_v.shape[0]):
     550                height_v[i,:] = stage_v[i,:] - elev_v
     551
     552        # Height at centroids       
     553        height_c_tmp = (height_v[:, vols0] + height_v[:,vols1] + height_v[:,vols2])/3.0
     554       
     555        # Compute xmom/xvel/ymom/yvel
    531556        if velocity_extrapolation:
    532             # Compute velocity from vertex velocities, then back-compute
    533             # momentum from that
    534             hInv=1.0/(heightAll+1.0-06)
    535             tmp = xmomAll*hInv*(heightAll>minimum_allowed_height)
    536             xvel=(tmp[:,vols0]+tmp[:,vols1]+tmp[:,vols2])/3.0
    537             htc = (heightAll[:,vols0] + heightAll[:,vols1] + heightAll[:,vols2])/3.0
    538             xvel_cent=getInds(xvel, timeSlices=inds, absMax=True)
    539             xmom_cent=getInds(xvel*htc, timeSlices=inds, absMax=True)
    540 
    541             tmp = ymomAll*hInv*(heightAll>minimum_allowed_height)
    542             yvel=(tmp[:,vols0]+tmp[:,vols1]+tmp[:,vols2])/3.0
    543             yvel_cent=getInds(yvel, timeSlices=inds, absMax=True)
    544             ymom_cent=getInds(yvel*htc, timeSlices=inds, absMax=True)
    545            
    546             vel_cent=getInds( (xvel**2+yvel**2)**0.5, timeSlices=inds)
    547            
     557
     558            xvel_v = xmom_v*0.
     559            yvel_v = ymom_v*0.
     560
     561            hInv = 1.0/(height_v+1.0e-12)
     562            hWet = (height_v > minimum_allowed_height)
     563
     564            xvel_v = xmom_v*hInv*hWet
     565            yvel_v = ymom_v*hInv*hWet
     566
     567            # Final xmom/ymom centroid values
     568            xvel_cent = (xvel_v[:, vols0] + xvel_v[:,vols1] + xvel_v[:,vols2])/3.0
     569            xmom_cent = xvel_cent*height_c_tmp
     570            yvel_cent = (yvel_v[:, vols0] + yvel_v[:,vols1] + yvel_v[:,vols2])/3.0
     571            ymom_cent = yvel_cent*height_c_tmp
     572
    548573        else:
    549             # Compute momenta from vertex momenta, then back compute velocity from that
    550             tmp=xmomAll*(heightAll>minimum_allowed_height)
    551             htc = (heightAll[:,vols0] + heightAll[:,vols1] + heightAll[:,vols2])/3.0
    552             hInv=1./(htc+1.0e-06)
    553             xmom=(tmp[:,vols0]+tmp[:,vols1]+tmp[:,vols2])/3.0
    554             xmom_cent=getInds(xmom,  timeSlices=inds, absMax=True)
    555             xvel_cent=getInds(xmom*hInv, timeSlices=inds, absMax=True)
    556 
    557             tmp=ymomAll*(heightAll>minimum_allowed_height)
    558             ymom=(tmp[:,vols0]+tmp[:,vols1]+tmp[:,vols2])/3.0
    559             ymom_cent=getInds(ymom,  timeSlices=inds, absMax=True)
    560             yvel_cent=getInds(ymom*hInv, timeSlices=inds, absMax=True)
    561             vel_cent = getInds( (xmom**2+ymom**2)**0.5*hInv, timeSlices=inds)
     574            hInv = 1.0/(height_c_tmp + 1.0e-12)
     575            hWet = (height_c_tmp > minimum_allowed_height)
     576            xmom_v =  fid.variables['xmomentum'][inds2]
     577            xmom_cent = (xmom_v[:,vols0] + xmom_v[:,vols1] + xmom_v[:,vols2])/3.0
     578            xvel_cent = xmom_cent*hInv*hWet
     579            ymom_v =  fid.variables['ymomentum'][inds2]
     580            ymom_cent = (ymom_v[:,vols0] + ymom_v[:,vols1] + ymom_v[:,vols2])/3.0
     581            yvel_cent = ymom_cent*hInv*hWet
     582
     583    # Velocity
     584    vel_cent = (xvel_cent**2 + yvel_cent**2)**0.5
     585
     586    if inds == 'max':
     587        vel_cent = vel_cent.max(axis=0,keepdims=True)
     588        #vel_cent = getInds(vel_cent, timeSlices=inds)
     589        xmom_cent = getInds(xmom_cent, timeSlices=inds,absMax=True)
     590        ymom_cent = getInds(ymom_cent, timeSlices=inds,absMax=True)
     591        xvel_cent = getInds(xvel_cent, timeSlices=inds,absMax=True)
     592        yvel_cent = getInds(yvel_cent, timeSlices=inds,absMax=True)
    562593
    563594    fid.close()
    … …  
    568599
    569600
    570 def animate_1D(time, var, x, ylab=' '): #, x=range(var.shape[1]), vmin=var.min(), vmax=var.max()):
    571     # Input: time = one-dimensional time vector;
    572     #        var =  array with first dimension = len(time) ;
    573     #        x = (optional) vector width dimension equal to var.shape[1];
     601def animate_1D(time, var, x, ylab=' '):
     602    """Animate a 2d array with a sequence of 1d plots
     603
     604     Input: time = one-dimensional time vector;
     605            var =  array with first dimension = len(time) ;
     606            x = (optional) vector width dimension equal to var.shape[1];
     607            ylab = ylabel for plot
     608    """
    574609   
    575610    import pylab
    … …  
    921956        #    myTimeStep=[len(p2.time)-1]
    922957
     958        myTimeStep_Orig = myTimeStep
    923959        # Now, myTimeStep just holds indices we want to plot in p2
    924960        if(myTimeStep!='max'):
    … …  
    10031039       
    10041040    # Loop over all output quantities and produce the output
    1005     for myTSi in myTimeStep:
     1041    for myTSindex, myTSi in enumerate(myTimeStep):
    10061042        if(verbose):
    10071043            print 'Reduction = ', myTSi
    … …  
    10351071                    timestepString='max'
    10361072                else:
    1037                     timestepString=str(round(p2.time[myTS]))
     1073                    timestepString=str(myTimeStep_Orig[myTSindex])+'_Time_'+str(round(p2.time[myTS]))
    10381074            elif(myTS=='pointData'):
    10391075                gridq=myInterpFun(xyzPoints[:,2])
    … …  
    10751111
    10761112    x0=p.xllcorner
    1077     x1=p.yllcorner
     1113    y0=p.yllcorner
    10781114
    10791115    # Make vertices for PolyCollection Object

  • trunk/anuga_core/source/anuga/utilities/test_plot_utils.py

    r9303 r9364  
    88
    99flow_algorithms = ['DE0', 'DE1', '1_5', '2_0', 'tsunami']
    10 verbose=False
     10verbose = False
    1111
    1212class Test_plot_utils(unittest.TestCase):
    … …  
    142142            Is called by a test function for various flow algorithms
    143143        """
    144 
    145144        p=util.get_output('test_plot_utils.sww')
    146145        p2=util.get_centroids(p,velocity_extrapolation=ve)
    … …  
    206205        p_m=util.get_output('test_plot_utils.sww', timeSlices='max')
    207206        pc_m=util.get_centroids(p_m,velocity_extrapolation=ve)
    208 
     207       
    209208        assert(p_m.time==p.time.max())
    210209        assert(pc_m.time==p2.time.max())

  • trunk/anuga_core/validation_tests/analytical_exact/parabolic_basin/plot_results.py

    r9345 r9364  
    1414# Get the sww file
    1515p_st = util.get_output('parabola.sww')
    16 p2_st=util.get_centroids(p_st) #(p_st, velocity_extrapolation=True)
     16p2_st=util.get_centroids(p_st, velocity_extrapolation=True) #(p_st, velocity_extrapolation=True)
    1717
    1818index_origin = 3978
Note: See TracChangeset for help on using the changeset viewer.