| 1 | import matplotlib |
|---|
| 2 | import matplotlib.pyplot as plt |
|---|
| 3 | import matplotlib.tri as tri |
|---|
| 4 | from pylab import * |
|---|
| 5 | |
|---|
| 6 | # The abstract Python-MPI interface |
|---|
| 7 | from anuga_parallel.parallel_abstraction import size, rank, get_processor_name |
|---|
| 8 | from anuga_parallel.parallel_abstraction import finalize, send, receive |
|---|
| 9 | from anuga_parallel.parallel_abstraction import pypar_available, barrier |
|---|
| 10 | |
|---|
| 11 | import random |
|---|
| 12 | |
|---|
| 13 | import numpy as num |
|---|
| 14 | |
|---|
| 15 | ''' |
|---|
| 16 | Adds random noise to stage quantities (centroid) in the domain. |
|---|
| 17 | domain - domain where noise is to be added |
|---|
| 18 | samples - number of centroids to add noise to |
|---|
| 19 | dthr - minimum depth threshold, stages corresponding to depths less than this value are not modified |
|---|
| 20 | mag - (relative) magnitude of the noise, uniformly distributed in [-mag, mag] |
|---|
| 21 | ''' |
|---|
| 22 | def addRandomNoiseToStage(domain, samples = 50, dthr = 0.1, mag = 1E-15): |
|---|
| 23 | |
|---|
| 24 | # Calculate depth |
|---|
| 25 | depth = domain.get_quantity('stage').centroid_values[domain.tri_full_flag == 1] \ |
|---|
| 26 | - domain.get_quantity('elevation').centroid_values[domain.tri_full_flag == 1] |
|---|
| 27 | |
|---|
| 28 | stage = domain.get_quantity('stage') |
|---|
| 29 | |
|---|
| 30 | # Sample centroid values with depth greater than dthr |
|---|
| 31 | n_tri = int(num.sum(domain.tri_full_flag)) |
|---|
| 32 | tri_index = num.array(range(0,n_tri)) |
|---|
| 33 | submerged = tri_index[depth > dthr] |
|---|
| 34 | submerged = random.sample(submerged, min(samples, len(submerged))) |
|---|
| 35 | |
|---|
| 36 | # Add random noise to sampled centroid values |
|---|
| 37 | if len(submerged) > 0: |
|---|
| 38 | for i in range(len(submerged)): |
|---|
| 39 | new_val = (1.0 + random.uniform(-1.0, 1.0)*mag)*stage.centroid_values[submerged[i]] |
|---|
| 40 | domain.get_quantity('stage').centroid_values.put(submerged[i], new_val) |
|---|
| 41 | |
|---|
| 42 | ''' |
|---|
| 43 | Plot errors in centroid quantities, triangles with errors are denoted in blue, while |
|---|
| 44 | triangles without errors are denoted in green. |
|---|
| 45 | |
|---|
| 46 | domain - domain for which the error plots will be displayed |
|---|
| 47 | control_data - control data to compare quantity values against (must have serial indexing) |
|---|
| 48 | rthr - relative error threshold |
|---|
| 49 | athr - absolute error threshold |
|---|
| 50 | quantity - quantity e.g stage |
|---|
| 51 | filename - output filename |
|---|
| 52 | ''' |
|---|
| 53 | |
|---|
| 54 | def plotCentroidError(domain, control_data, rthr = 1E-7, athr = 1E-12, |
|---|
| 55 | quantity = 'stage', filename = 'centroid_error.png'): |
|---|
| 56 | |
|---|
| 57 | n_triangles = num.sum(domain.tri_full_flag) |
|---|
| 58 | |
|---|
| 59 | if size() > 1: |
|---|
| 60 | # If parallel, translate control data to parallel indexing |
|---|
| 61 | local_control_data = num.zeros(n_triangles) |
|---|
| 62 | inv_tri_map = domain.get_inv_tri_map() |
|---|
| 63 | |
|---|
| 64 | for i in range(n_triangles): |
|---|
| 65 | local_control_data[i] = control_data[inv_tri_map[(rank(), i)]] |
|---|
| 66 | else: |
|---|
| 67 | local_control_data = control_data |
|---|
| 68 | |
|---|
| 69 | # Evaluate absolute and relative difference between control and actual values |
|---|
| 70 | stage = domain.get_quantity(quantity) |
|---|
| 71 | actual_data = stage.centroid_values[:n_triangles] |
|---|
| 72 | adiff = num.fabs((actual_data - local_control_data)) |
|---|
| 73 | rdiff = adiff/num.fabs(local_control_data) |
|---|
| 74 | |
|---|
| 75 | # Compute masks for error (err_mask) and non-error (acc_mask) vertex indices based on thresholds |
|---|
| 76 | vertices = domain.get_vertex_coordinates() |
|---|
| 77 | err_mask = rdiff > rthr |
|---|
| 78 | err_mask[adiff <= athr] = False |
|---|
| 79 | err_mask = num.repeat(err_mask, 3) |
|---|
| 80 | |
|---|
| 81 | acc_mask = ~err_mask |
|---|
| 82 | inv_tri_map = domain.get_inv_tri_map() |
|---|
| 83 | |
|---|
| 84 | # Plot error and non-error triangle |
|---|
| 85 | if rank() == 0: |
|---|
| 86 | fx = {} |
|---|
| 87 | fy = {} |
|---|
| 88 | gx = {} |
|---|
| 89 | gy = {} |
|---|
| 90 | |
|---|
| 91 | fx[0] = vertices[acc_mask,0] |
|---|
| 92 | fy[0] = vertices[acc_mask,1] |
|---|
| 93 | gx[0] = vertices[err_mask,0] |
|---|
| 94 | gy[0] = vertices[err_mask,1] |
|---|
| 95 | |
|---|
| 96 | # Receive vertex indices of non-error triangles (fx, fy) and error triangles (gx, gy) |
|---|
| 97 | for i in range(1,size()): |
|---|
| 98 | fx[i] = receive(i) |
|---|
| 99 | fy[i] = receive(i) |
|---|
| 100 | gx[i] = receive(i) |
|---|
| 101 | gy[i] = receive(i) |
|---|
| 102 | |
|---|
| 103 | # Plot non-error triangles in green |
|---|
| 104 | for i in range(0,size()): |
|---|
| 105 | n = int(len(fx[i])/3) |
|---|
| 106 | triang = num.array(range(0,3*n)) |
|---|
| 107 | triang.shape = (n, 3) |
|---|
| 108 | |
|---|
| 109 | if len(fx[i]) > 0: |
|---|
| 110 | plt.triplot(fx[i], fy[i], triang, 'g-') |
|---|
| 111 | |
|---|
| 112 | # Plot error triangles in blue |
|---|
| 113 | for i in range(0,size()): |
|---|
| 114 | n = int(len(gx[i])/3) |
|---|
| 115 | |
|---|
| 116 | triang = num.array(range(0,3*n)) |
|---|
| 117 | triang.shape = (n, 3) |
|---|
| 118 | |
|---|
| 119 | if len(gx[i]) > 0: |
|---|
| 120 | plt.triplot(gx[i], gy[i], triang, 'b--') |
|---|
| 121 | |
|---|
| 122 | # Save plot |
|---|
| 123 | plt.savefig(filename) |
|---|
| 124 | |
|---|
| 125 | else: |
|---|
| 126 | # Send error and non-error vertex indices to Proc 0 |
|---|
| 127 | send(vertices[acc_mask,0], 0) |
|---|
| 128 | send(vertices[acc_mask,1], 0) |
|---|
| 129 | send(vertices[err_mask,0], 0) |
|---|
| 130 | send(vertices[err_mask,1], 0) |
|---|
