def compute_error(degree, simulation_index, number_of_points):
# Number of grid points. This is assumed to be the same in the x and y directions.
nx = number_of_points
ny = number_of_points
# Number of halo nodes at each end
halo = degree/2
# Read in the simulation output
path = "./mms_%d_%d/mms_%d_%d_opsc_code/" % (degree, simulation_index, degree, simulation_index)
dump = glob.glob(path + "/mms_*.h5")
if not dump or len(dump) > 1:
print "Error: No dump file found, or more than one dump file found."
sys.exit(1)
f = h5py.File(dump[-1], 'r')
group = f["mms_%d_%d_block" % (degree, simulation_index)]
# Get the numerical solution field
phi = group["phi"].value
# Ignore the 2 halo nodes at the left (and bottom) end of the domain. Include one strip of halo points at the right (and top) of the domain to enforce periodicity in the solution field plot.
phi = phi[halo:nx+halo+1, halo:ny+halo+1]
print phi.shape
# Grid spacing. Note: The length of the domain is divided by nx (or ny) and not nx-1 (or ny-1) because of the periodicity. In total we have nx+1 points, but we only solve nx points; the (nx+1)-th point is set to the same value as the 0-th point to give a full period, to save computational effort.
dx = (2.0*pi)/(nx)
dy = (2.0*pi)/(ny)
# Coordinate arrays.
x = numpy.zeros((nx+1)*(ny+1)).reshape((nx+1, ny+1))
y = numpy.zeros((nx+1)*(ny+1)).reshape((nx+1, ny+1))
phi_analytical = numpy.zeros((nx+1)*(ny+1)).reshape((nx+1, ny+1))
# Compute the error by first interpolating the numerical and analytical results onto a much finer grid of points and computing the L2 norm of the absolute difference.
grid_points = []
grid_numerical = []
grid_analytical = []
target_grid_x, target_grid_y = numpy.mgrid[0:2*pi:32j, 0:2*pi:32j]
for i in range(0, nx+1):
for j in range(0, ny+1):
# Work out the x and y coordinates. Note the swapping of the 'j' and 'i' here.
x[i,j] = j*dx
y[i,j] = i*dy
grid_points.append([x[i,j], y[i,j]])
grid_numerical.append(phi[i,j])
phi_analytical[i,j] = sin(x[i,j])*cos(y[i,j])
grid_analytical.append(phi_analytical[i,j])
grid_points = numpy.array(grid_points)
grid_numerical = numpy.array(grid_numerical)
grid_analytical = numpy.array(grid_analytical)
interpolated_numerical = griddata(grid_points, grid_numerical, (target_grid_x, target_grid_y), method='nearest')
interpolated_analytical = griddata(grid_points, grid_analytical, (target_grid_x, target_grid_y), method='nearest')
# Only plot phi for the 6th order simulations.
if degree == 12:
plot_phi(simulation_index, phi, phi_analytical)
return numpy.linalg.norm(abs(interpolated_numerical - interpolated_analytical), ord=2)
评论列表
文章目录
