[Anuga-user] how to ensure that the depth value is greater than or equal to 0 in the complex terrain?

[刘军志 <li...@lr...>]

Dear ANUGA users,



    We are currently using Anuga to simulate water flow motion. A binary installation package, i.e., Anuga-2.0-cp27-none-win32.whl, was used. We found that, if the terrain is complex, Anuga simulation would generate some negative depth values (i.e., stage<elevation). Anuga simulation would not generate negative depth values if the terrain is simple. 

    

    Our question is how to ensure that the depth value is greater than or equal to 0 in the complex terrain.




    Thank you very much for your help!
    Regards,


    

Jingchao Jiang

Hangzhou Dianzi University, China


The test code we used is as follows:

import unittest
import os.path
import sys

import numpy
import anuga

from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
from anuga.shallow_water.shallow_water_domain import Domain
from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.utilities.system_tools import get_pathname_from_package

from anuga.structures.inlet_operator import Inlet_operator
from anuga.file.netcdf import NetCDFFile
import numpy as num

domain_length = 40.0
domain_width = 50.0
elevation_0=10.0
elevation_1=10.0  
class Test_inlet_operator():
   
           
    def _create_domain(self,swwfilename,elevation_fun,d_length,
                            d_width,
                            dx,
                            dy):
        
        points, vertices, boundary = rectangular_cross(int(d_length/dx), int(d_width/dy),
                                                        len1=d_length, len2=d_width)
        domain = Domain(points, vertices, boundary)   
        domain.set_name(swwfilename)                 # Output name
        domain.set_store()
        domain.set_default_order(2)
        domain.H0 = 0.01
        domain.tight_slope_limiters = 1

        #print 'Size', len(domain)

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

      
        #print 'Setting Quantities....'
        domain.set_quantity('elevation', elevation_fun)  # Use function for elevation
        domain.set_quantity('stage',  expression='elevation')   # Use function for elevation

        Br = anuga.Reflective_boundary(domain)
        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
        
        return domain

    def test_inlet_constant_Q(self,swwfilename,elevation_fun):
        """test_inlet_Q
        
        This tests that the inlet operator adds the correct amount of water
        """

      

        domain = self._create_domain(swwfilename=swwfilename,elevation_fun=elevation_fun,d_length=domain_length,
                                     d_width=domain_width,
                                     dx = 1,
                                     dy = 1)

        vol0 = domain.compute_total_volume()

        finaltime = 30.0
        Q1 = 5.00
        Q2 = 10.0
        
        center1 = (1,1)
        radius1 = 1
        region1 = anuga.Region(domain, center=center1, radius=radius1)
        Inlet_operator(domain, region1, Q1, logging=False)
       
        center2 = (11,30)
        radius2 = 1
        region2 = anuga.Region(domain, center=center2, radius=radius2)
        Inlet_operator(domain, region2, Q2, logging=False)
        for t in domain.evolve(yieldstep = 1.0, finaltime = finaltime):
            #domain.write_time()
            #print domain.volumetric_balance_statistics()
            pass
 

        vol1 = domain.compute_total_volume()

        assert numpy.allclose((Q1+Q2)*finaltime, vol1-vol0, rtol=1.0e-8) 
        assert numpy.allclose((Q1+Q2)*finaltime, domain.fractional_step_volume_integral, rtol=1.0e-8) 

def simple_topography(x, y):
    """Set up a elevation
    """
    z = numpy.zeros(x.shape,dtype='d')
    z[:] = elevation_0
    numpy.putmask(z, x > domain_length/2, elevation_1)
    return z
            
def complex_topography(x,y):
    """Complex topography defined by a function of vectors x and y."""
    z = -x/10
    N = len(x)
    for i in range(N):
        # Step
        if 10 < x[i] < 12:
            z[i] += 0.4 - 0.05*y[i]
        # Constriction
        if 27 < x[i] < 29 and y[i] > 3:
            z[i] += 2
        # Pole
        if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
            z[i] += 2
    return z
def check_depth(swwfilename):
    frame_number=3
    fin = NetCDFFile(swwfilename+".sww", 'r')
    ele=num.ravel(num.array(fin.variables["elevation"][:], num.float))
    stage=num.array(fin.variables["stage"][frame_number], num.float)
    print("number of negative depths:")
    print(num.sum(stage-ele<0))
    print("-------------------------")
# =========================================================================
if __name__ == "__main__":

    test=Test_inlet_operator()
    test.test_inlet_constant_Q('Test_Outlet_Inlet1',simple_topography)
    check_depth("Test_Outlet_Inlet1")
    test.test_inlet_constant_Q('Test_Outlet_Inlet2',complex_topography)
    check_depth("Test_Outlet_Inlet2")