Changeset 8469

Timestamp:

Jul 19, 2012, 4:17:35 PM (13 years ago)

Author:

steve

Message:

Playing with viscosity operator

Location:

trunk/anuga_core/source/anuga/operators

Files:

• kinematic_viscosity_operator.py (modified) (9 diffs)
• test_kinematic_viscosity_operator.py (modified) (4 diffs)

trunk/anuga_core/source/anuga/operators/kinematic_viscosity_operator.py

@@ -13 +13 @@
 
 
-class Kinematic_Viscosity_Operator(Operator):
+class Kinematic_viscosity_operator(Operator):
     """
     Class for setting up structures and matrices for kinematic viscosity differential
@@ -27 +27 @@
     """
 
-    def __init__(self, domain, diffusivity=None, use_triangle_areas=True, verbose=False):
+    def __init__(self, domain, diffusivity='height', use_triangle_areas=True, verbose=False):
 
         if verbose: log.critical('Kinematic Viscosity: Beginning Initialisation')
@@ -51 +51 @@
             self.diffusivity.set_boundary_values(diffusivity)
 
+        if isinstance(diffusivity, str):
+            self.diffusivity = self.domain.get_quantity(diffusivity)
+            
     
         assert isinstance(self.diffusivity, Quantity)
@@ -128 +131 @@
 
         # diffusivity
-        h = domain.quantities['height']
-
-        #d = self.diffusivity
+        d = self.diffusivity
+
+
+        assert num.all(d.centroid_values >= 0.0)
 
         #d.set_values(num.where(h.centroid_values<1.0, 0.0, 1.0), location= 'centroids')
@@ -144 +148 @@
 
         #Update operator using current height
-        self.update_elliptic_matrix(h)
-
-        (u, self.u_stats) = self.parabolic_solve(u, u, h, u_out=u, update_matrix=False, output_stats=True)
-
-        (v, self.v_stats) = self.parabolic_solve(v, v, h, u_out=v, update_matrix=False, output_stats=True)
+        self.update_elliptic_matrix(d)
+
+        (u, self.u_stats) = self.parabolic_solve(u, u, d, u_out=u, update_matrix=False, output_stats=True)
+
+        (v, self.v_stats) = self.parabolic_solve(v, v, d, u_out=v, update_matrix=False, output_stats=True)
 
         # Update the conserved quantities
@@ -155 +159 @@
 
         self.dt = 0.0
-        #print self.timestepping_statistics()
+
+
 
     def statistics(self):
@@ -165 +170 @@
     def timestepping_statistics(self):
 
-        message = '    Kinematic Viscosity Operator: '
+        from anuga import indent
+
+        message = indent+'Kinematic Viscosity Operator: \n'
         if self.u_stats is not None:
-            message  += ' u iterations %.5g, ' % self.u_stats.iter
+            message  += indent + indent + 'u: ' + self.u_stats.__str__() +'\n'
 
         if self.v_stats is not None:
-            message += ' v iterations %.5g, ' % self.v_stats.iter
+            message  += indent + indent + 'v: ' + self.v_stats.__str__()
 
         return message
@@ -482 +489 @@
 
     def parabolic_solve(self, u_in, b, a = None, u_out = None, update_matrix=True, \
-                       imax=10000, tol=1.0e-8, atol=1.0e-8,
-                       iprint=None, output_stats=False):
+                       output_stats=False, use_dt_tol=True, iprint=None, imax=10000):
         """
         Solve for u in the equation
@@ -505 +511 @@
         """
 
+
+        if use_dt_tol:
+            tol  = min(self.dt,0.01)
+            atol = min(self.dt,0.01)
+        else:
+            tol  =  1.0e-5
+            atol = 1.0e-5
+
+
+
+        
         if u_out == None:
             u_out = Quantity(self.domain)

trunk/anuga_core/source/anuga/operators/test_kinematic_viscosity_operator.py

@@ -564 +564 @@
 
 
-        u_out = operator.parabolic_solve(u_in, b, a, iprint=1)
+        u_out = operator.parabolic_solve(u_in, b, a, iprint=1, use_dt_tol=False)
 
         assert num.allclose(u_out.centroid_values, U_mod[:n])
@@ -671 +671 @@
         kv.update_elliptic_matrix(h)
 
-        kv.parabolic_solve(u, u, h, u_out=u, update_matrix=False, iprint=1)
-
-        kv.parabolic_solve(v, v, h, u_out=v, update_matrix=False, iprint=1)
-
+        kv.parabolic_solve(u, u, h, u_out=u, update_matrix=False, iprint=1, use_dt_tol=False)
+
+        kv.parabolic_solve(v, v, h, u_out=v, update_matrix=False, iprint=1, use_dt_tol=False)
+
+
+        #print u.centroid_values
+        #print u.boundary_values
         assert num.allclose(u.centroid_values, num.ones_like(u.centroid_values), rtol=1.0e-1)
         assert num.allclose(u.boundary_values, num.ones_like(u.boundary_values))
@@ -762 +765 @@
         kv.update_elliptic_matrix(h)
 
-        kv.parabolic_solve(u, u, h, u_out=u, update_matrix=False, iprint=1)
-
-        kv.parabolic_solve(v, v, h, u_out=v, update_matrix=False, iprint=1)
+        kv.parabolic_solve(u, u, h, u_out=u, update_matrix=False, iprint=1, use_dt_tol=False)
+
+        kv.parabolic_solve(v, v, h, u_out=v, update_matrix=False, iprint=1, use_dt_tol=False)
 
 
@@ -793 +796 @@
         assert num.allclose(uh.centroid_values, u.centroid_values*h.centroid_values )
         assert num.allclose(vh.centroid_values, v.centroid_values*h.centroid_values )
+
+    def test_kinematic_operator_default(self):
+
+        from anuga import rectangular_cross_domain
+        from anuga import Reflective_boundary
+
+        m1 = 10
+        n1 = 10
+        domain = rectangular_cross_domain(m1,n1)
+
+        domain.set_flow_algorithm('2_0')
+
+        #
+        domain.set_quantity('elevation', expression='x')
+        domain.set_quantity('friction', 0.03)
+        domain.set_quantity('stage',expression='elevation + 2*(x-0.5)')
+        domain.set_quantity('xmomentum', expression='2*x+3*y')
+        domain.set_quantity('ymomentum', expression='5*x+7*y')
+
+        B = Reflective_boundary(domain)
+        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
+
+        kv = Kinematic_Viscosity_Operator(domain)
+
+
+        # let's make timestep large so that the final solution will look like
+        #the solution of hte elliptic problem. In this case u -> 1, v -> 2.
+
+
+        for t in domain.evolve(yieldstep = 1.0, finaltime = 10.0):
+            domain.write_time()
+            domain.print_operator_timestepping_statistics()
+
+#        dt = 1000.0
+#        kv.dt = dt
+#        n = kv.n
+#        nt = kv.tot_len
+#
+#        kv.update_elliptic_matrix(h)
+#
+#        kv.parabolic_solve(u, u, h, u_out=u, update_matrix=False, iprint=1)
+#
+#        kv.parabolic_solve(v, v, h, u_out=v, update_matrix=False, iprint=1)
+#
+#
+#        #print 'u'
+#        #print u.centroid_values
+#        #print u.boundary_values
+#
+#        #print num.where(h.centroid_values > 0.0, 1.0, 0.0)
+#
+#        assert num.allclose(u.centroid_values, num.where(h.centroid_values > 0.0, 1.0, 0.0), rtol=1.0e-1)
+#        assert num.allclose(u.boundary_values, num.ones_like(u.boundary_values))
+#
+#        assert num.allclose(v.centroid_values, num.where(h.centroid_values > 0.0, 2.0, 0.0), rtol=1.0e-1)
+#        assert num.allclose(v.boundary_values, 2.0*num.ones_like(v.boundary_values))
+#
+#
+#        domain.update_centroids_of_momentum_from_velocity()
+#
+#        domain.distribute_to_vertices_and_edges()
+#
+#        uh = domain.quantities['xmomentum']
+#        vh = domain.quantities['ymomentum']
+#
+#        #print 'uh'
+#        #print uh.centroid_values
+#        #print uh.boundary_values
+#
+#        assert num.allclose(uh.centroid_values, u.centroid_values*h.centroid_values )
+#        assert num.allclose(vh.centroid_values, v.centroid_values*h.centroid_values )