1 | from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain |
---|
2 | import anuga.abstract_2d_finite_volumes.mesh_factory as mesh_factory |
---|
3 | import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh |
---|
4 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Dirichlet_boundary |
---|
5 | from kinematic_viscosity import Kinematic_Viscosity_Operator |
---|
6 | import numpy as num |
---|
7 | import anuga.utilities.log as log |
---|
8 | |
---|
9 | |
---|
10 | print "Setting up" |
---|
11 | grid_rows = 40 |
---|
12 | grid_cols = 40 |
---|
13 | #Set up a rectangular grid |
---|
14 | points,elements,boundary = mesh_factory.rectangular_cross(grid_rows,grid_cols) |
---|
15 | mesh = neighbour_mesh.Mesh(points,elements) |
---|
16 | n = len(mesh) |
---|
17 | print "n =",n |
---|
18 | boundary_map = {} |
---|
19 | for (vol_id,edge) in boundary.keys(): |
---|
20 | x = mesh.get_edge_midpoint_coordinate(vol_id,edge)[0] |
---|
21 | y = mesh.get_edge_midpoint_coordinate(vol_id,edge)[1] |
---|
22 | #Define boundary conditions: u=x^2-y^2, v=x-y+2 |
---|
23 | boundary_map[(vol_id,edge)] = Dirichlet_boundary([1,x*x-y*y,x-y+2]) |
---|
24 | #Now define an operator |
---|
25 | domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map) |
---|
26 | KV_Operator = Kinematic_Viscosity_Operator(domain) |
---|
27 | print "Applying stage heights" |
---|
28 | #Set up the operator and RHS |
---|
29 | h = num.ones((n,),num.float) |
---|
30 | |
---|
31 | def solve(): |
---|
32 | KV_Operator.apply_stage_heights(h) |
---|
33 | rhs = num.zeros((n,2),num.float) #we will solve for u and v |
---|
34 | print "Solving" |
---|
35 | #Solve |
---|
36 | velocities = KV_Operator.cg_solve(rhs) |
---|
37 | new_velocities = KV_Operator.parabolic_solver(velocities) |
---|
38 | |
---|
39 | #Do the profiling |
---|
40 | import profile, pstats |
---|
41 | FN = 'kv_profile.dat' |
---|
42 | print "Starting solve..." |
---|
43 | profile.run('solve()', FN) |
---|
44 | S = pstats.Stats(FN) |
---|
45 | s = S.sort_stats('cumulative').print_stats(30) |
---|