Commit 76d3b52d authored by Dion Haefner's avatar Dion Haefner

updated references - all tests should pass now

parent eb60c8c1
......@@ -10,28 +10,6 @@ from dorie.testtools.utilities.decorators import evaluation_script
from dorie.testtools.utilities import richards_equation
from dorie.parfield.parameter_file import read_parameter_file
def average_2d_mesh(xy, z, connectivity):
"""
Averages an array z defined on a triangular or rectangular mesh with vertices
xy and a certain connectivity.
"""
npt, dim = xy.shape
ncells, dim1 = connectivity.shape
assert ncells == len(z), "Shape mismatch (got: {}, {})".format(connectivity.shape, z.shape)
assert dim == 2, "Only two-dimensional geometry supported"
if dim1 == dim+1:
areafun = lambda t: .5 * abs(t[:,0,0] * t[:,1,1] - t[:,0,1] * t[:,1,0])
elif dim1 == dim+2:
areafun = lambda t: abs(t[:,0,0] * t[:,1,1])
else:
raise ValueError("Only triangular and rectangular geometry supported")
corners = xy[connectivity]
t = corners[:,1:,:] - corners[:,0,:][:,np.newaxis,:]
area = areafun(t)
zsum = np.sum(area * z)
areasum = np.sum(area)
return zsum / areasum
@evaluation_script
def evaluate(iniinfo,runtime):
pvd = iniinfo["output.fileName"] + ".pvd"
......@@ -79,6 +57,7 @@ def evaluate(iniinfo,runtime):
tri, connectivity = grid.triangulation()
data = grid.data()
iy = np.argsort(y)
iy_inv = np.arange(len(y))[iy]
y = y[iy]
k0_dorie = data["K_0"].flatten()[iy]
......@@ -103,7 +82,8 @@ def evaluate(iniinfo,runtime):
r_y = np.array(r_y)
head = data["head"].flatten()[iy]
res_head = r_y - head
l2_head = np.sqrt(average_2d_mesh(tri.T,np.square(res_head),connectivity))
average = lambda v: grid.integrate(v) / grid.integrate(1)
l2_head = np.sqrt(average(np.square(res_head[iy_inv])))
print("head L2 error: {:.2e}".format(float(l2_head)))
plt.figure()
......@@ -120,7 +100,7 @@ def evaluate(iniinfo,runtime):
flux = data["flux"][iy,:]
res_flux = flux - np.array([0.,influx,0.])[np.newaxis,:]
l2_flux = [np.sqrt(average_2d_mesh(tri.T,np.square(r),connectivity)) for r in res_flux.T]
l2_flux = [np.sqrt(average(np.square(r[iy_inv]))) for r in res_flux.T]
print("flux L2 errors: {:.2e}, {:.2e}, {:.2e}".format(*l2_flux))
plt.figure()
......
......@@ -18,15 +18,18 @@ def evaluate(iniinfo,runtime):
except KeyError:
raise ValueError("The reference evaluator assumes the key _reference.path to be existent in the inifile")
if os.path.isfile(reference):
vtkfile_run = vtkfile.VTKFile(vtks[-1].path)
vtkfile_run.read()
vtkfile_ref = vtkfile.VTKFile(reference)
vtkfile_ref.read()
# demand equal grids when not using adaptivity
strict_compare = text_to_bool(iniinfo["adaptivity.useAdaptivity"])
ret = fuzzy_compare_grid.compare(vtkfile_run.grid, vtkfile_ref.grid, strict=strict_compare)
else:
if not os.path.isfile(reference):
raise RuntimeError("Reference file {} does not exist".format(reference))
vtkfile_run = vtkfile.VTKFile(vtks[-1].path)
vtkfile_run.read()
vtkfile_ref = vtkfile.VTKFile(reference)
vtkfile_ref.read()
print("Operating on files: {}, {}".format(vtks[-1].path,reference))
# demand equal grids when not using adaptivity
strict_compare = not text_to_bool(iniinfo["adaptivity.useAdaptivity"])
ret = fuzzy_compare_grid.compare(vtkfile_run.grid, vtkfile_ref.grid, strict=strict_compare)
return ret
......@@ -3,66 +3,104 @@ from scipy import interpolate
from dorie.utilities.grid import BaseGrid
def compare(grid1,grid2,abstol=1E-6,reltol=1E-3,strict=False):
if not isinstance(grid1,BaseGrid) or not isinstance(grid2,BaseGrid):
raise TypeError("Inputs must be of type BaseGrid (got: {} and {})"\
.format(type(grid1),type(grid2)))
def sort_grid(grid):
new_shape = (grid.cellCenters().shape[0],-1)
cells = grid.cellCenters().reshape(*new_shape)
ix = np.lexsort(cells)
sorted_cells = np.array([c[ix] for c in cells])
data = grid.data()
for key in data:
data[key] = data[key].flatten()
return sorted_cells, data
# sort grid to correct for possibly different order of grid elements
cells1, data1 = sort_grid(grid1)
cells2, data2 = sort_grid(grid2)
interpolate_grid = False
# check if number of grid elements is equal
if not cells1.shape == cells2.shape:
def _equal_grid_shapes(cells1,cells2):
shapes_equal = cells1.shape == cells2.shape
if not shapes_equal:
print("Grid shapes do not match: {}, {}".format(cells1.shape,cells2.shape))
if strict:
return False
else:
interpolate_grid = True
return shapes_equal
# check if coordinates are equal
if not np.allclose(cells1,cells2,atol=abstol,rtol=reltol):
def _equal_coordinates(cells1,cells2,abstol,reltol):
if cells1.shape != cells2.shape:
return False
coords_equal = np.allclose(cells1,cells2,atol=abstol,rtol=reltol)
if not coords_equal:
print("Grids do not match")
print(cells1)
print(cells2)
idiff = np.argmax(np.sqrt(np.sum((cells1-cells2)**2, axis=0)))
print("\tLargest difference: got {} and {}".format(cells1[:,idiff],cells2[:,idiff]))
if strict:
return False
else:
interpolate_grid = True
return coords_equal
# check if the same variables are contained
if set(data1.keys()) != set(data2.keys()):
def _equal_variable_keys(grid1,grid2):
data1, data2 = [g.data() for g in (grid1,grid2)]
keys_equal = set(data1.keys()) == set(data2.keys())
if not keys_equal:
print("Variable keys do not match")
return False
return keys_equal
def _sort_cells(grid):
cells = grid.cellCenters()
ix = np.lexsort(cells.T)
sorted_cells = cells[ix,...]
return sorted_cells
def _sort_grid(grid):
cells = grid.cellCenters()
ix = np.lexsort(cells.T)
sorted_cells = cells[ix,...]
sorted_data = {key: val[ix] for key,val in grid.data().items()}
return sorted_cells, sorted_data
def _equal_grids(grid1,grid2,abstol,reltol):
cells1 = _sort_cells(grid1)
cells2 = _sort_cells(grid2)
return _equal_grid_shapes(cells1,cells2) and _equal_coordinates(cells1,cells2,abstol,reltol)
def _calc_residuals(grid1,grid2,interpolate_grid):
cells1 = grid1.cellCenters()
cells2 = grid2.cellCenters()
data1 = grid1.data()
data2 = grid2.data()
# interpolate (nearest neighbor) if desired
if interpolate_grid:
for key in data1:
residuals = {}
for key in data1.keys():
print("Interpolating variable {} ...".format(key))
interpolator = interpolate.NearestNDInterpolator(cells1.T,data1[key])
data1[key] = interpolator(cells2.T)
interpolator = interpolate.NearestNDInterpolator(cells1,data1[key])
residuals[key] = interpolator(cells2) - data2[key]
else:
_, sorted_data1 = _sort_grid(grid1)
_, sorted_data2 = _sort_grid(grid2)
residuals = {var: v1 - v2 for var, v1, v2 in zip(data1.keys(),sorted_data1.values(),sorted_data2.values())}
return residuals, grid2
# compare variable data for each variable
def _compare_max_norm(residuals, atol, rtol):
passed = True
for key in data1:
if not np.allclose(data1[key], data2[key], atol=abstol, rtol=reltol):
print("Data in variable {} does not match".format(key))
idiff = np.argmax(np.abs(data1[key] - data2[key]))
print("\tLargest difference: got {} and {}".format(data1[key][idiff],data2[key][idiff]))
for var, res in residuals.items():
if not np.allclose(res, 0., atol=atol, rtol=rtol):
print(var, res)
print("Data in variable {} does not match".format(var))
print("\tLargest difference: {}".format(np.max(np.abs(res))))
passed = False
return passed
def _compare_l2_norm(residuals, grid, abstol):
passed = True
for var, res in residuals.items():
l2 = np.sqrt(grid.integrate(res**2) / grid.integrate(1))
if np.all(l2 > abstol):
print("Data in variable {} does not match".format(var))
print("\tL2 error: {}".format(l2))
passed = False
return passed
def compare(grid1,grid2,abstol=1E-6,reltol=1E-3,strict=False):
if not isinstance(grid1,BaseGrid) or not isinstance(grid2,BaseGrid):
raise TypeError("Inputs must be of type BaseGrid (got: {} and {})"\
.format(type(grid1),type(grid2)))
interpolate_grid = False
if not _equal_grids(grid1,grid2,abstol,reltol):
if strict:
return False
else:
interpolate_grid = True
if not _equal_variable_keys(grid1,grid2):
return False
residuals, res_grid = _calc_residuals(grid1,grid2,interpolate_grid)
if strict:
return _compare_max_norm(residuals, abstol, reltol)
else:
return _compare_l2_norm(residuals, res_grid, abstol)
......@@ -132,6 +132,39 @@ class UnstructuredVTKGrid(BaseGrid):
"""
return self.points.T, self.connectivity.reshape(self._numCells,self.ppc)
def integrate(self,v):
"""
Integrates an array v defined on each cell center across the grid.
:param v: Input array of shape (ncells, dim) or scalar
"""
if np.asarray(v).size > 1 and self._numCells != np.asarray(v).shape[0]:
raise ValueError("Shape mismatch (got: {}, {})".format(self._numCells, np.asarray(v).shape[0]))
if self.dim == 2:
if self.ppc == 3:
areafun = lambda t: .5 * abs(t[:,0,0] * t[:,1,1] - t[:,0,1] * t[:,1,0])
elif self.ppc == 4:
areafun = lambda t: abs(t[:,0,0] * t[:,1,1])
else:
raise ValueError("Only triangular and rectangular geometry supported")
elif self.dim == 3:
if self.ppc == 4:
areafun = lambda t: 1./6. * abs(np.sum(t[:,0,:] * np.cross(t[:,1,:],t[:,2,:]),axis=-1))
elif self.ppc == 6:
areafun = lambda t: abs(t[:,0,0] * t[:,1,1] * t[:,2,2])
else:
raise ValueError("Only tetrahedral and cubic geometry supported")
else:
raise ValueError("Only two and three-dimensional grids supported")
corners = self.points[self.connectivity.reshape(self._numCells,self.ppc)]
t = corners[:,1:,:] - corners[:,0,:][:,np.newaxis,:]
cell_area = areafun(t)
broadcast_shape = (self._numCells,) + (1,) * (np.asarray(v).ndim-1)
integral = np.sum(cell_area.reshape(broadcast_shape) * v, axis=0)
return np.squeeze(integral)
def _calcCellCenters(self):
"""
:returns: Coordinates of cell centers
......
......@@ -27,7 +27,7 @@ parameters.arrayFile = "{_asset_path}/parfields/sand.h5", "{_asset_path}/parfiel
[_ode]
flux = -3e-6
head_abstol = 1E-5, 5E-4 | expand gridtype
flux_abstol = 1E-8, 2E-7 | expand gridtype
flux_abstol = 1E-8, 5E-7 | expand gridtype
[NewtonParameters]
AbsoluteLimit = 1E-8
......
......@@ -27,8 +27,8 @@ parameters.interpolation = linear
[_ode]
flux = -3e-6
head_abstol = 1E-5, 1E-4 | expand gridtype
flux_abstol = 1E-8, 5E-8 | expand gridtype
head_abstol = 1E-5, 5E-4 | expand gridtype
flux_abstol = 1E-8, 2E-7 | expand gridtype
[NewtonParameters]
AbsoluteLimit = 1E-8
......
......@@ -5,6 +5,7 @@ _test_command = run
_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_evaluation = reference
_reference.path = {_asset_path}/references/2d/
_reference.abstol = 1e-8, 0.01 | expand adaptivity
output.fileName = reference_2d | unique name
output.outputPath = reference_2d | unique name
......
......@@ -6,6 +6,7 @@ _asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_evaluation = reference
_reference.path = {_asset_path}/references/evaporation/
_reference.abstol = 1e-8, 0.05 | expand adaptivity
output.fileName = reference_evaporation | unique name
output.outputPath = reference_evaporation | unique name
......
......@@ -4,7 +4,9 @@ __name = reference_interpolators
_test_command = run
_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_evaluation = reference
_reference.path = {_asset_path}/references/interpolators/
_reference.abstol = 0.05
output.fileName = reference_interpolators | unique name
output.outputPath = reference_interpolators | unique name
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment