source: trunk/anuga_core/source/anuga_parallel/print_stats.py @ 7841

Last change on this file since 7841 was 5242, checked in by steve, 16 years ago
File size: 5.7 KB
Line 
1#!/usr/bin/env python
2#########################################################
3#
4#
5# Calculate and print the norms of the domain
6#
7#
8#  The routines defined here are intended for debugging
9# use. They print the norms of the quantities in the
10# domain. As opposed to the definitions given
11# in utiltites.norms these calculations take a
12# parallel domain into account.
13#
14#  Authors: Linda Stals and April 2006
15#
16#
17#########################################################
18
19import sys
20
21import pypar
22
23from Numeric import array, Int8, zeros, ones, take, nonzero, Float
24from anuga.utilities.norms import l1_norm, l2_norm, linf_norm
25
26#########################################################
27#
28# Find out which triangles are full triangles (only these
29# triangles should be included in the norm calculations)
30#
31# *)
32#
33# -------------------------------------------------------
34#
35# *) A 1-D array, tri_full_flag, is returned.
36#
37# *) The size of tri_full_flag is the same as the number
38# of vertices in the domain
39#
40# *) If tri_full_flag[i] = 1, then triangle number i is
41# a full triangle, if tri_full_flag[i] = 0  the triangle
42# is a ghost triangle
43#
44#########################################################
45
46def build_full_flag(domain, ghost_recv_dict):
47
48    tri_full_flag = ones(len(domain.get_triangles()), Int8)
49    for i in ghost_recv_dict.keys():
50        for id in ghost_recv_dict[i][0]:
51            tri_full_flag[id] = 0
52       
53
54    return tri_full_flag
55
56#########################################################
57#
58# Print the l1 norm of the given quantity
59#
60# *) The quantity is an array containing three columns
61#
62# -------------------------------------------------------
63#
64# *) The l1 norm is calculated along each axis
65#
66# *) The l1 norm is printed
67#
68# *) Processor 0 prints the results
69#
70#
71#########################################################
72def print_l1_stats(full_edge):
73
74    numprocs = pypar.size()
75    myid = pypar.rank()
76   
77    tri_norm = zeros(3, Float)
78    recv_norm = zeros(3, Float)
79    tri_norm[0] = l1_norm(full_edge[:, 0])
80    tri_norm[1] = l1_norm(full_edge[:, 1])
81    tri_norm[2] = l1_norm(full_edge[:, 2])
82    if myid == 0:
83        for p in range(numprocs-1):
84            pypar.receive(p+1, recv_norm)
85            tri_norm[0] = tri_norm[0]+recv_norm[0]
86            tri_norm[1] = tri_norm[1]+recv_norm[1]
87            tri_norm[2] = tri_norm[2]+recv_norm[2]
88        print 'l1_norm along each axis : [', tri_norm[0],', ', tri_norm[1], ', ', tri_norm[2], ']'
89
90    else:
91        pypar.send(tri_norm, 0)
92
93#########################################################
94#
95# Print the l2 norm of the given quantity
96#
97# *) The quantity is an array containing three columns
98#
99# -------------------------------------------------------
100#
101# *) The l2 norm is calculated along each axis
102#
103# *) The l2 norm is printed
104#
105# *) Processor 0 prints the results
106#
107#
108#########################################################
109def print_l2_stats(full_edge):
110
111    numprocs = pypar.size()
112    myid = pypar.rank()
113   
114    tri_norm = zeros(3, Float)
115    recv_norm = zeros(3, Float)
116    tri_norm[0] = pow(l2_norm(full_edge[:, 0]), 2)
117    tri_norm[1] = pow(l2_norm(full_edge[:, 1]), 2)
118    tri_norm[2] = pow(l2_norm(full_edge[:, 2]), 2)
119    if myid == 0:
120        for p in range(numprocs-1):
121            pypar.receive(p+1, recv_norm)
122            tri_norm[0] = tri_norm[0]+recv_norm[0]
123            tri_norm[1] = tri_norm[1]+recv_norm[1]
124            tri_norm[2] = tri_norm[2]+recv_norm[2]
125        print 'l2_norm along each axis : [', pow(tri_norm[0], 0.5),', ', pow(tri_norm[1], 0.5), \
126              ', ', pow(tri_norm[2], 0.5), ']'
127    else:
128        pypar.send(tri_norm, 0)
129
130
131#########################################################
132#
133# Print the linf norm of the given quantity
134#
135# *) The quantity is an array containing three columns
136#
137# -------------------------------------------------------
138#
139# *) The linf norm is calculated along each axis
140#
141# *) The linf norm is printed
142#
143# *) Processor 0 prints the results
144#
145#
146#########################################################
147def print_linf_stats(full_edge):
148
149    numprocs = pypar.size()
150    myid = pypar.rank()
151   
152    tri_norm = zeros(3, Float)
153    recv_norm = zeros(3, Float)
154       
155    tri_norm[0] = linf_norm(full_edge[:, 0])
156    tri_norm[1] = linf_norm(full_edge[:, 1])
157    tri_norm[2] = linf_norm(full_edge[:, 2])
158    if myid == 0:
159        for p in range(numprocs-1):
160            pypar.receive(p+1, recv_norm)
161            tri_norm[0] = max(tri_norm[0], recv_norm[0])
162            tri_norm[1] = max(tri_norm[1], recv_norm[1])
163            tri_norm[2] = max(tri_norm[2], recv_norm[2])
164        print 'linf_norm along each axis : [', tri_norm[0],', ', tri_norm[1], ', ', tri_norm[2], ']'
165    else:
166        pypar.send(tri_norm, 0)
167
168       
169#########################################################
170#
171# Print the norms of the quantites assigned to the domain
172# (this is useful for checking the numerical results
173# in the parallel computation)
174#
175# *) tri_full_flag states which of the triangles are
176# full triangles
177#
178#
179# -------------------------------------------------------
180#
181# *) For each quantity, the l1, l2 and linf norms are
182# printed
183#
184# *) The size of tri_full_flag is the same as the number
185# of vertices in the domain
186#
187# *) Only the full triangles are used in the norm
188# calculations
189#
190# *) Processor 0 prints the results
191#
192#########################################################
193def print_test_stats(domain, tri_full_flag):
194
195    myid = pypar.rank()
196
197    for k in domain.quantities.keys():
198        TestStage = domain.quantities[k]
199        if myid == 0:
200            print " ===== ", k, " ===== "
201        full_edge = take(TestStage.edge_values, nonzero(tri_full_flag))
202        print_l1_stats(full_edge)
203        print_l2_stats(full_edge)
204        print_linf_stats(full_edge)
205
Note: See TracBrowser for help on using the repository browser.