Changeset 9221


Ignore:
Timestamp:
Jun 26, 2014, 10:58:41 PM (10 years ago)
Author:
steve
Message:

Updating some validation test and added get_total_volume to domain class (parallel and sequential)

Location:
trunk/anuga_core
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r9218 r9221  
    558558
    559559    def set_DE2_defaults(self):
    560         """Set up the defaults for running the flow_algorithm "DE1"
     560        """Set up the defaults for running the flow_algorithm "DE2"
    561561           A 'discontinuous elevation' method
    562562        """
    563         self.set_CFL(0.95)
     563        self.set_CFL(1.0)
    564564        self.set_use_kinematic_viscosity(False)
    565565        #self.timestepping_method='rk2'#'rk3'#'euler'#'rk2'
    … …  
    625625
    626626
     627    def set_DE3_defaults(self):
     628        """Set up the defaults for running the flow_algorithm "DE3"
     629           A 'discontinuous elevation' method
     630        """
     631        self.set_CFL(0.9)
     632        self.set_use_kinematic_viscosity(False)
     633        #self.timestepping_method='rk2'#'rk3'#'euler'#'rk2'
     634        self.set_timestepping_method(1)
     635       
     636        self.set_using_discontinuous_elevation(True)
     637        self.set_compute_fluxes_method('DE')
     638        self.set_distribute_to_vertices_and_edges_method('DE')
     639       
     640        # Don't place any restriction on the minimum storable height
     641        self.minimum_storable_height=-99999999999.0
     642        self.minimum_allowed_height=1.0e-12
     643
     644        self.use_edge_limiter=True
     645        self.set_default_order(2)
     646        self.set_extrapolate_velocity()
     647
     648        self.beta_w=0.7
     649        self.beta_w_dry=0.1
     650        self.beta_uh=0.7
     651        self.beta_uh_dry=0.1
     652        self.beta_vh=0.7
     653        self.beta_vh_dry=0.1
     654       
     655
     656        #self.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2,
     657        #         'ymomentum': 2, 'elevation': 2, 'height':2})
     658        #self.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2,
     659        #         'ymomentum': 2, 'elevation': 1})
     660        self.set_store_centroids(True)
     661
     662        self.optimise_dry_cells=False
     663
     664        # We need the edge_coordinates for the extrapolation
     665        self.edge_coordinates=self.get_edge_midpoint_coordinates()
     666
     667        # By default vertex values are NOT stored uniquely
     668        # for storage efficiency. We may override this (but not so important since
     669        # centroids are stored anyway
     670        # self.set_store_vertices_smoothly(False)
     671
     672        self.maximum_allowed_speed=0.0
     673
     674        ## FIXME: Should implement tracking of boundary fluxes
     675        ## Keep track of the fluxes through the boundaries
     676        self.boundary_flux_integral=num.ndarray(1)
     677        self.boundary_flux_integral[0]=0.
     678        self.boundary_flux_sum=num.ndarray(1)
     679        self.boundary_flux_sum[0]=0.
     680
     681        self.call=1 # Integer counting how many times we call compute_fluxes_central
     682
     683        if self.processor == 0 and self.verbose:
     684            print '##########################################################################'
     685            print '#'
     686            print '# Using discontinuous elevation solver DE3'
     687            print '#'
     688            print '# A slightly less diffusive version than DE0, uses euler timestepping'
     689            print '#'
     690            print '# Make sure you use centroid values when reporting on important output quantities'
     691            print '#'
     692            print '##########################################################################'
     693
     694
    627695
    628696    def update_special_conditions(self):
    … …  
    780848           DE1
    781849           DE2
     850           DE3
    782851        """
    783852
    … …  
    787856            flag = str(float(str(flag))).replace(".","_")
    788857
    789         flow_algorithms = ['1_0', '1_5', '1_75', '2_0', '2_0_limited', '2_5', 'tsunami', 'yusuke', 'DE0', 'DE1', 'DE2']
     858        flow_algorithms = ['1_0', '1_5', '1_75', '2_0', '2_0_limited', '2_5', \
     859                           'tsunami', 'yusuke', 'DE0', 'DE1', 'DE2', 'DE3']
    790860
    791861        if flag in flow_algorithms:
    … …  
    903973        if self.flow_algorithm == 'DE2':
    904974            self.set_DE2_defaults()
     975           
     976        if self.flow_algorithm == 'DE3':
     977            self.set_DE3_defaults()
    905978
    906979    def get_flow_algorithm(self):
    … …  
    908981
    909982        Currently
    910            1_0, 1_5, 1_75 2_0, 2_5, tsunami, DE0, DE1, DE2
     983           1_0, 1_5, 1_75 2_0, 2_5, tsunami, DE0, DE1, DE2, DE3
    911984        """
    912985
    … …  
    12201293        return self.get_quantity('elevation').\
    12211294                   get_maximum_location(indices=wet_elements)
     1295
     1296
     1297    def get_total_volume(self):
     1298   
     1299        from anuga import myid, numprocs, send, receive, barrier
     1300
     1301        if self.get_using_discontinuous_elevation():   
     1302            Height = self.quantities['height']
     1303            volume = Height.get_integral()
     1304        else:
     1305            Stage = self.quantities['stage']
     1306            Elev =  self.quantities['elevation']
     1307            Height = Stage-Elev
     1308            volume = Height.get_integral()
     1309           
     1310        if myid == 0:
     1311            total_volume = volume
     1312            for i in range(1,numprocs):
     1313                remote_volume = receive(i)
     1314                total_volume = total_volume + remote_volume
     1315        else:
     1316            send(volume,0)
     1317       
     1318        barrier()
     1319   
     1320        if myid == 0:
     1321            for i in range(1,numprocs):
     1322                send(total_volume,i)
     1323        else:
     1324            total_volume = receive(0)
     1325       
     1326        return total_volume
    12221327
    12231328
    … …  
    17291834                print 'Cumulative mass protection: '+str(mass_error)+' m^3 '
    17301835
    1731         elif self.flow_algorithm == 'DE1'  or self.flow_algorithm == 'DE0':
     1836        elif self.flow_algorithm.startswith('DE'):
    17321837
    17331838            from swDE1_domain_ext import protect
    … …  
    18311936           
    18321937            if len(negative_ids)>0:
    1833                 #print 'NEGATIVE INDICES'
     1938                import warnings
     1939                warnings.warn('Negative cells being set to zero depth, possible loss of conservation')
    18341940                Stage.centroid_values[negative_ids] = Elev.centroid_values[negative_ids]
    18351941                Xmom.centroid_values[negative_ids] = 0.0

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

    r9174 r9221  
    4343
    4444#--------------------
    45 # Make plot animation
     45# Make plots
    4646#--------------------
    4747# Find an y value close to y==50
    … …  
    6767#pyplot.ylim([0.0,0.05])
    6868pyplot.xlabel('Xposition m')
    69 pyplot.ylabel('Depth m')
    70 pyplot.title('Velocity down the slope (along y=50.)')
     69pyplot.ylabel('Velocity m/s')
     70pyplot.title('X Velocity down the slope (along y=50.)')
    7171pyplot.legend(loc='best')
    7272pyplot.savefig('xvelocity_x.png')
    … …  
    9898pyplot.ylim([0.0,0.35])
    9999pyplot.xlabel('Yposition along the line x=50')
    100 pyplot.ylabel('Xvelocity m/s')
    101 pyplot.title('Final Xvelocity around the line x=50.')
     100pyplot.ylabel('Velocity m/s')
     101pyplot.title('X Velocity around the line x=50.')
    102102pyplot.legend(loc='best')
    103103pyplot.savefig('x_velocity.png')
    … …  
    108108#pyplot.ylim([0.0,0.3])
    109109pyplot.xlabel('Yposition along the line x=50')
    110 pyplot.ylabel('Yvelocity')
    111 pyplot.title('Final Yvelocity around the line x=50.')
     110pyplot.ylabel('Velocity m/s')
     111pyplot.title('Y Velocity around the line x=50.')
    112112pyplot.legend(loc='best')
    113113pyplot.savefig('y_velocity.png')

  • trunk/anuga_core/validation_tests/analytical_exact/rundown_mild_slope_coarse/results.tex

    r9038 r9221  
    4141\begin{figure}
    4242\begin{center}
    43 \includegraphics[width=0.8\textwidth]{final_depth.png}
     43\includegraphics[width=0.8\textwidth]{depth_x.png}
    4444\caption{Depth in the downstream direction}
    4545\label{fig:depthdownchan}
     46\end{center}
     47\end{figure}
     48
     49\begin{figure}
     50\begin{center}
     51\includegraphics[width=0.8\textwidth]{xvelocity_x.png}
     52\caption{X velocity in the downstream direction}
     53\label{fig:xveldownchan}
     54\end{center}
     55\end{figure}
     56
     57\begin{figure}
     58\begin{center}
     59\includegraphics[width=0.8\textwidth]{depth_y.png}
     60\caption{Depth in the downstream direction}
     61\label{fig:depthacrosschan}
    4662\end{center}
    4763\end{figure}

  • trunk/anuga_core/validation_tests/analytical_exact/trapezoidal_channel/numerical_channel_floodplain.py

    r9217 r9221  
    1010import anuga
    1111import numpy
    12 from anuga import myid, finalize, distribute
     12from anuga import myid, finalize, distribute, numprocs, receive, send
    1313
    1414#------------------------------------------------------------------------------
    … …  
    2828          man_n, l0
    2929
    30 #l0 = 2.0
     30l0 = 2.0
    3131
    3232assert chan_width < floodplain_width, \
    … …  
    5454# do not overlap.
    5555channel_polygon = [ [floodplain_width/2. - chan_width/2., +l0],
    56                     [floodplain_width/2. - chan_width/2., floodplain_length-l0],
    57                     [floodplain_width/2. + chan_width/2., floodplain_length-l0],
     56                    [floodplain_width/2. - chan_width/2., floodplain_length - l0],
     57                    [floodplain_width/2. + chan_width/2., floodplain_length - l0],
    5858                    [floodplain_width/2. + chan_width/2., +l0]
    5959                    ]
    … …  
    127127# Define inlet operator
    128128#---------------------------------------------------------------------
    129 flow_in_yval=0.0
     129flow_in_yval=5.0
    130130line1 = [ [floodplain_width/2. - chan_width/2., flow_in_yval],\
    131131          [floodplain_width/2. + chan_width/2., flow_in_yval] ]
    … …  
    173173    parameter_file.close()
    174174
     175
     176
     177
    175178#------------------------------------------------------------------------------
    176179# Evolve system through time
    … …  
    178181for t in domain.evolve(yieldstep=10.0, finaltime=1000.0):
    179182    if myid == 0 and verbose: print domain.timestepping_statistics()
     183   
     184#     if numpy.allclose(t,0.0):
     185#         exact_volume = domain.get_total_volume()
     186#     else:
     187#         exact_volume = exact_volume + Qin*10.0
     188#         
     189#     total_volume= domain.get_total_volume()
     190   
     191   
     192    #if myid == 0 and verbose: print anuga.indent,'total volume - exact_volume ', total_volume - exact_volume   
    180193
    181194domain.sww_merge(delete_old=True)

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

    r9216 r9221  
    1414print numpy.any(v)
    1515# Numerical results along a central channel 'slice'
    16 index=p.stage.shape[0]-1
     16index= -1
    1717V1 = p.stage[index,v] - p.elev[v]
    1818V2 = p.yvel[index,v]
    … …  
    4949dc_analytical = scipy.optimize.fmin(minme, x0=1.0)[0]
    5050
    51 
     51print 'dc_analytic ',dc_analytical
    5252
    5353##################################

  • trunk/anuga_core/validation_tests/analytical_exact/trapezoidal_channel/results.tex

    r9216 r9221  
    2828
    2929Figure~\ref{fig:hydrographs} show the hydrographs through various cross sections showing the flows limiting to the
    30 expected inflow $Q$.
     30expected inflow $Q$ (For coarser grids there is a discrepency between the expected and calculated limiting
     31hydrograph due to the error in calculating the hydrograph). It is also noted that the transient flow is quite different
     32for different grid sizes. We theorize that the coarser grids produce a rougher bed which slows down the flow.
    3133
    3234\begin{figure}
Note: See TracChangeset for help on using the changeset viewer.