ANUGA Kinematic Viscosity Module
Lindon Roberts, 2010

Summary
============================================================================
This module applies kinematic viscosity smoothing to the velocities in a 
domain, given by the parabolic equations
  du/dt = (h*u_x)_x + (h*u_y)_y
  dv/dt = (h*v_x)_x + (h*v_y)_y
for u(x,y) and v(x,y) being water velocities, and h(x,y) is the stage height.
For more details on the implementation, see the report. In particular, the
code implements a pure elliptic solver. Both solvers are implicit.

Files Included
============================================================================
- kinematic_viscosity.py: The main python module.
- kinematic_viscosity_ext.c: The C extension to the python module.
- test_kinematic_viscosity.py: A unit test for the python module.
- util_ext.*: Copies of the files in anuga.utilities (for kinematic_viscosity_ext.c).
- sparse.py: A modified version of anuga.utilities.sparse (which will supersede it).
- test_sparse.py: A unit test for the sparse module.
- make.bat: A makefile script for MinGW on Windows (this may have to be modified).

Usage
============================================================================
Some sample code can be found in the unit tests. In particular, one 
requires an anuga.abstract_2d_finite_volumes.domain domain, a numpy 
vector h with h.shape == (n,) of stage heights and momentum numpy matrix 
p with p == [(hu) | (hv)] where u and v are x- and y-velocity and
p.shape == (n,2).

The stage height vector is the stage heights for the current time step.

To initialise the kinematic viscosity operator, one uses
kv_operator = kinematic_viscosity.Kinematic_Viscosity_Operator(domain)
with optional parameters

1. apply_triangle_areas (default True)
Whether or not to include the factor of 1/|Triangle area| in the numerical 
approximation (see the report for more details).

2. verbose (default False)
Whether to include detailed logs for anuga.log.

Parabolic Solver:
We use the commands
> kv_operator.apply_stage_heights(h)
> p_new = kv_operator.parabolic_solver(p)

Elliptic Solver:
This solves
  (h*u_x)_x + (h*u_y)_y = f_1
  (h*v_y)_y + (h*v_y)_y = f_2
We also need a right-hand side numpy matrix f with f.shape == (n,2)
and f == [f_1 | f_2]. Then
> kv_operator.apply_stage_heights(h)
> p = kv_operator.cg_solve(f)
You can also solve for in only one dimension with
> kv_operator.apply_stage_heights(h)
> kv_operator.set_qty_considered('u')
> p_x = kv_operator.cg_solve(f1)
(and similarly for v)

============================================================================

Last Modified: 25 October 2010