source: trunk/anuga_core/source/anuga/operators/test_kinematic_viscosity.py @ 8152

from anuga import Domain
from anuga import Quantity
from anuga import Dirichlet_boundary
from kinematic_viscosity import Kinematic_Viscosity

import numpy as num
from math import sqrt
import unittest

class Test_Kinematic_Viscosity(unittest.TestCase):
    def setUp(self):
        pass

    def tearDown(self):
        pass
    
    #First test operator class (1 triangle)
    def operator1(self):
        points = num.array([[0.0,0.0],[1.0,0.0],[0.0,1.0]])
        elements = num.array([[0,1,2]])
        boundary_map = {}
        boundary_map[(0,0)] = 'edge0'
        boundary_map[(0,1)] = 'edge1'
        boundary_map[(0,2)] = 'edge2'
        domain = Domain(coordinates=points,vertices=elements,boundary=boundary_map)

        D0 = Dirichlet_boundary([1,0,3])
        D1 = Dirichlet_boundary([2,1,0])
        D2 = Dirichlet_boundary([3,1,2])
        domain.set_boundary({'edge0': D0, 'edge1': D1, 'edge2': D2})

        domain.set_quantity('stage', lambda x,y : x+2*y )
        domain.set_quantity('elevation', lambda x,y : 3*x+5*y )



        #print domain.quantities['stage'].vertex_values

        #print domain.quantities['stage'].edge_values

        domain.update_boundary()


        #print domain.quantities['stage'].boundary_values
        
        return Kinematic_Viscosity(domain)

    #Second test operator class (2 triangles)
    def operator2(self):
        points = num.array([[0.0,0.0],[1.0,0.0],[1.0,1.0],[0.0,1.0]])
        elements = num.array([[0,1,3],[1,2,3]])
        boundary_map = {}
        boundary_map[(0,1)] = 'edge0'
        boundary_map[(0,2)] = 'edge1'
        boundary_map[(1,0)] = 'edge2'
        boundary_map[(1,2)] = 'edge3'
        domain = Domain(coordinates=points,vertices=elements,boundary=boundary_map)

        D0 = Dirichlet_boundary([1,1,2])
        D1 = Dirichlet_boundary([1,2,2])
        D2 = Dirichlet_boundary([1,1,0])
        D3 = Dirichlet_boundary([1,2,1])

        domain.set_boundary({'edge0': D0, 'edge1': D1, 'edge2': D2, 'edge3': D3})
        domain.update_boundary()



        return Kinematic_Viscosity(domain)

    def test_enumerate_boundary(self):
        operator1 = self.operator1()
        boundary_enumeration = operator1.domain.boundary_enumeration

        assert boundary_enumeration[(0,0)] == 0
        assert boundary_enumeration[(0,1)] == 1
        assert boundary_enumeration[(0,2)] == 2

        operator2 = self.operator2()
        boundary_enumeration = operator2.domain.boundary_enumeration


        assert boundary_enumeration[(0,1)] == 0
        assert boundary_enumeration[(0,2)] == 1
        assert boundary_enumeration[(1,0)] == 2
        assert boundary_enumeration[(1,2)] == 3

    def test_geo_structure(self):
        operator1 = self.operator1()
        indices = operator1.geo_structure_indices
        values = operator1.geo_structure_values

        assert num.allclose(indices, num.array([[1, 2, 3]]))
        assert num.allclose(values, num.array([[-6.0, -6.0/sqrt(5), -6.0/sqrt(5)]]))

        operator2 = self.operator2()
        indices = operator2.geo_structure_indices
        values = operator2.geo_structure_values
        assert num.allclose(indices, num.array([[1,2,3],[4,0,5]]))
        assert num.allclose(values, num.array([[-3.0,-6.0/sqrt(5),-6.0/sqrt(5)],[-6.0/sqrt(5),-3.0,-6.0/sqrt(5)]]))

    def test_elliptic_matrix_one_triangle(self):

        operator = self.operator1()
        domain = operator.domain
        
        operator.update_elliptic_matrix()

        A = operator.elliptic_matrix

        assert num.allclose(A.todense(), num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)]))

        diffusivity = operator.diffusivity
        diffusivity.set_values(10.0)
        diffusivity.set_boundary_values(10.0)
        
        operator.update_elliptic_matrix()

        assert num.allclose(A.todense(), 10*num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)]))
    

    def test_elliptic_matrix_two_triangles(self):


        operator = self.operator2()

        domain = operator.domain
        diffusivity = operator.diffusivity

        A = operator.elliptic_matrix

        diffusivity.set_values(1.0)
        diffusivity.set_boundary_values(1.0)
        operator.update_elliptic_matrix()

        A0 = num.array([[-3.0,3.0,0.0,0.0,0.0,0.0],
                        [0.0,-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0]])
        A1 = num.array([[-6.0/sqrt(5.0),0.0,6.0/sqrt(5.0),0.0,0.0,0.0],\
                        [3.0,-3.0,0.0,0.0,0.0,0.0]])
        A2 = num.array([[-6.0/sqrt(5.0),0.0,0.0,6.0/sqrt(5.0),0.0,0.0],\
                        [0.0, -6.0/sqrt(5.0), 0.0, 0.0, 0.0, 6.0/sqrt(5.0)]])


        assert num.allclose(A.todense(), A0+A1+A2)

        diffusivity.set_values([2.0, 1.0], location = 'centroids')
        diffusivity.set_boundary_values(1.0)
        operator.update_elliptic_matrix()

        A = operator.elliptic_matrix
        

        assert num.allclose(A.todense()[0,:], 1.5*A0[0,:]+1.5*A1[0,:]+1.5*A2[0,:])
        assert num.allclose(A.todense()[1,:], A0[1,:]+1.5*A1[1,:]+A2[1,:])

        diffusivity.set_values([-2.0, -2.0], location = 'centroids')
        diffusivity.set_boundary_values(1.0)
        operator.update_elliptic_matrix()

        assert num.allclose(A.todense()[0,:], -2*A0[0,:]-0.5*A1[0,:]-0.5*A2[0,:])
        assert num.allclose(A.todense()[1,:], -0.5*A0[1,:]-2*A1[1,:]-0.5*A2[1,:])

    def test_elliptic_multiply_include_boundary_one_triangle(self):
        operator = self.operator1()
        operator.set_triangle_areas(False)

        print operator.apply_triangle_areas
        
        n = operator.n

        q_in = Quantity(operator.domain)
        q_in.set_values(1.0)
        q_in.set_boundary_values(1.0)
        operator.update_elliptic_matrix()

        
        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])


        q_1 = operator.elliptic_multiply(q_in)

        q_2 = operator.elliptic_multiply(q_in, quantity_out = q_in)

        assert id(q_in) == id(q_2)

        assert num.allclose(q_1.centroid_values,q_2.centroid_values)

        assert num.allclose( num.zeros((n,), num.float), q_1.centroid_values )

        #Now have different boundary values

        q_in.set_values(1.0)
        q_in.set_boundary_values(0.0)
        operator.update_elliptic_matrix()


        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])


        q_1 = operator.elliptic_multiply(q_in)

        assert num.allclose( [-6.0-12.0/sqrt(5)], q_1.centroid_values )
        
    def test_elliptic_multiply_exclude_boundary_one_triangle(self):
        operator = self.operator1()
        operator.set_triangle_areas(False)

        print operator.apply_triangle_areas
        #n = operator.n

        q_in = Quantity(operator.domain)
        q_in.set_values(1.0)
        q_in.set_boundary_values(1.0)
        operator.update_elliptic_matrix()


        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])


        q_1 = operator.elliptic_multiply(q_in, include_boundary=False)


        assert num.allclose( [-6.0-12.0/sqrt(5)], q_1.centroid_values )

    def test_elliptic_multiply_include_boundary_one_triangle(self):
        operator = self.operator1()
        operator.set_triangle_areas(True)

        n = operator.n

        q_in = Quantity(operator.domain)
        q_in.set_values(1.0)
        q_in.set_boundary_values(1.0)
        operator.update_elliptic_matrix()


        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])


        q_1 = operator.elliptic_multiply(q_in)

        q_2 = operator.elliptic_multiply(q_in, quantity_out = q_in)

        assert id(q_in) == id(q_2)

        assert num.allclose(q_1.centroid_values,q_2.centroid_values)

        assert num.allclose( num.zeros((n,), num.float), q_1.centroid_values )

        #Now have different boundary values

        q_in.set_values(1.0)
        q_in.set_boundary_values(0.0)
        operator.update_elliptic_matrix()


        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])


        q_1 = operator.elliptic_multiply(q_in)

        assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values )

    def test_elliptic_multiply_exclude_boundary_one_triangle(self):
        operator = self.operator1()
        operator.set_triangle_areas(True)

        q_in = Quantity(operator.domain)
        q_in.set_values(1.0)
        q_in.set_boundary_values(1.0)
        operator.update_elliptic_matrix()


        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])


        q_1 = operator.elliptic_multiply(q_in, include_boundary=False)


        assert num.allclose( [-12.0-24.0/sqrt(5)], q_1.centroid_values )



    def test_mul(self):
        operator = self.operator1()

        q = Quantity(operator.domain)
        q.set_values(2.0)
        #q boundary_values should equal 0.0

        operator.build_elliptic_boundary_term()
        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
        V1 = num.array([2.0]) #(uh)=2
        U1 = num.array([[2.0],[0.0],[0.0],[0.0]])
        assert num.allclose(operator * q, 2*num.array(num.mat(A)*num.mat(U1)).reshape(1,))

    def test_cg_solve(self):
        #cf self.test_mul()
        operator1 = self.operator1()
        operator1.apply_stage_heights(num.array([[1.0]])) #h=1
        operator1.set_qty_considered('u')
        V = num.array([2.0]) #h=1, (uh)=2
        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
        U = num.array([[2.0,2.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
        test = 2*num.mat(A)*num.mat(U[:,0].reshape(4,1))
        X = operator1.cg_solve(num.array(test).reshape(1,))
        assert num.allclose(V, X)

    def test_parabolic_solve(self):
        operator1 = self.operator1()
        operator1.apply_stage_heights(num.array([[1.0]])) #h=1
        A = num.array([-6.0-12.0/sqrt(5), 6.0,  6.0/sqrt(5), 6.0/sqrt(5)])
        U = num.array([[2.0,1.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
        u = num.array([[2.0,1.0]])
        U_new = operator1.parabolic_solver(u)
        U_mod = num.array([[0.0,0.0],[2.0,1.0],[1.0,2.0],[1.0,0.0]])
        U_mod[0,:] = U_new
        assert num.allclose(U_new - operator1.dt * 2 * num.mat(A)*num.mat(U_mod), U[0,:])

################################################################################

if __name__ == "__main__":
    suite = unittest.makeSuite(Test_Kinematic_Viscosity, 'test')
    runner = unittest.TextTestRunner()
    runner.run(suite)