From 55adef5b1677979da377bd6d67794fa536e5f2b8 Mon Sep 17 00:00:00 2001 From: Connor Jordan Date: Tue, 30 Jul 2024 15:53:35 +0100 Subject: [PATCH 1/2] Replace .at() interpolation in 2D callbacks - .at() raises an error when the number of processors used is large, use a VertexOnlyMesh instead for interpolation --- thetis/callback.py | 56 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 52 insertions(+), 4 deletions(-) diff --git a/thetis/callback.py b/thetis/callback.py index 52195e926..d92845dda 100644 --- a/thetis/callback.py +++ b/thetis/callback.py @@ -529,6 +529,20 @@ def __init__(self, solver_obj, self.field_names = field_names self._name = name + # create VertexOnlyMesh and function spaces for interpolation + self.mesh = solver_obj.fields[field_names[0]].function_space().mesh() + self.vom = VertexOnlyMesh(self.mesh, detector_locations, redundant=True) + self.function_spaces = {} + for field_name in field_names: + field = solver_obj.fields[field_name] + if isinstance(field.function_space().ufl_element(), VectorElement): + P0DG = VectorFunctionSpace(self.vom, "DG", 0) + P0DG_input_ordering = VectorFunctionSpace(self.vom.input_ordering, "DG", 0) + else: + P0DG = FunctionSpace(self.vom, "DG", 0) + P0DG_input_ordering = FunctionSpace(self.vom.input_ordering, "DG", 0) + self.function_spaces[field_name] = (P0DG, P0DG_input_ordering) + @property def name(self): return self._name @@ -539,7 +553,8 @@ def variable_names(self): def _values_per_field(self, values): """ - Given all values evaulated in a detector location, return the values per field""" + Given all values evaulated in a detector location, return the values per field + """ i = 0 result = [] for dim in self.field_dims: @@ -554,7 +569,13 @@ def message_str(self, *args): for name, values in zip(self.detector_names, args)) def _evaluate_field(self, field_name): - return self.solver_obj.fields[field_name](self.detector_locations) + field = self.solver_obj.fields[field_name] + P0DG, P0DG_input_ordering = self.function_spaces[field_name] + f_at_points = Function(P0DG) + f_at_input_points = Function(P0DG_input_ordering) + f_at_points.interpolate(field) + f_at_input_points.interpolate(f_at_points) + return f_at_input_points.dat.data_ro def __call__(self): """ @@ -673,6 +694,7 @@ def __init__(self, solver_obj, fieldnames, x, y, self.tolerance = tolerance self.eval_func = eval_func self._initialized = False + self.bathymetry_val = None @PETSc.Log.EventDecorator("thetis.TimeSeriesCallback2D._initialize") def _initialize(self): @@ -684,10 +706,31 @@ def _initialize(self): xyz = (self.x, self.y, self.z) if self.on_sphere else (self.x, self.y) self.xyz = numpy.array([xyz]) + # create VertexOnlyMesh and function spaces + self.mesh = self.solver_obj.fields[self.fieldnames[0]].function_space().mesh() + self.vom = VertexOnlyMesh(self.mesh, self.xyz, redundant=True) + self.function_spaces = {} + for field_name in self.fieldnames: + field = self.solver_obj.fields[field_name] + if isinstance(field.function_space().ufl_element(), VectorElement): + P0DG = VectorFunctionSpace(self.vom, "DG", 0) + P0DG_input_ordering = VectorFunctionSpace(self.vom.input_ordering, "DG", 0) + else: + P0DG = FunctionSpace(self.vom, "DG", 0) + P0DG_input_ordering = FunctionSpace(self.vom.input_ordering, "DG", 0) + self.function_spaces[field_name] = (P0DG, P0DG_input_ordering) + # test evaluation try: if self.eval_func is None: - self.solver_obj.fields.bathymetry_2d.at(self.xyz, tolerance=self.tolerance) + bathymetry_field = self.solver_obj.fields.bathymetry_2d + P0DG = FunctionSpace(self.vom, "DG", 0) + P0DG_input_ordering = FunctionSpace(self.vom.input_ordering, "DG", 0) + f_at_points = Function(P0DG) + f_at_input_points = Function(P0DG_input_ordering) + f_at_points.interpolate(bathymetry_field) + f_at_input_points.interpolate(f_at_points) + self.bathymetry_val = f_at_input_points.dat.data_ro[:] else: self.eval_func(self.solver_obj.fields.bathymetry_2d, self.xyz, tolerance=self.tolerance) except PointNotInDomainError as e: @@ -707,7 +750,12 @@ def __call__(self): try: field = self.solver_obj.fields[fieldname] if self.eval_func is None: - val = field.at(self.xyz, tolerance=self.tolerance) + P0DG, P0DG_input_ordering = self.function_spaces[fieldname] + f_at_points = Function(P0DG) + f_at_input_points = Function(P0DG_input_ordering) + f_at_points.interpolate(field) + f_at_input_points.interpolate(f_at_points) + val = f_at_input_points.dat.data_ro[:] else: val = self.eval_func(field, self.xyz, tolerance=self.tolerance) arr = numpy.array(val) From 32dd936535c52a5bf93f1fa3400dc1007c911619 Mon Sep 17 00:00:00 2001 From: Connor Jordan Date: Thu, 10 Oct 2024 09:49:43 +0100 Subject: [PATCH 2/2] Add helper function for VOM interpolation Functions --- thetis/callback.py | 51 +++++++++------------------------------------- thetis/utility.py | 35 +++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 41 deletions(-) diff --git a/thetis/callback.py b/thetis/callback.py index d92845dda..878dc65a7 100644 --- a/thetis/callback.py +++ b/thetis/callback.py @@ -529,19 +529,8 @@ def __init__(self, solver_obj, self.field_names = field_names self._name = name - # create VertexOnlyMesh and function spaces for interpolation - self.mesh = solver_obj.fields[field_names[0]].function_space().mesh() - self.vom = VertexOnlyMesh(self.mesh, detector_locations, redundant=True) - self.function_spaces = {} - for field_name in field_names: - field = solver_obj.fields[field_name] - if isinstance(field.function_space().ufl_element(), VectorElement): - P0DG = VectorFunctionSpace(self.vom, "DG", 0) - P0DG_input_ordering = VectorFunctionSpace(self.vom.input_ordering, "DG", 0) - else: - P0DG = FunctionSpace(self.vom, "DG", 0) - P0DG_input_ordering = FunctionSpace(self.vom.input_ordering, "DG", 0) - self.function_spaces[field_name] = (P0DG, P0DG_input_ordering) + # initialise interpolation functions using vom_interpolator_functions + self.interp_functions = vom_interpolator_functions(solver_obj, field_names, detector_locations) @property def name(self): @@ -570,12 +559,10 @@ def message_str(self, *args): def _evaluate_field(self, field_name): field = self.solver_obj.fields[field_name] - P0DG, P0DG_input_ordering = self.function_spaces[field_name] - f_at_points = Function(P0DG) - f_at_input_points = Function(P0DG_input_ordering) + f_at_points, f_at_input_points = self.interp_functions[field_name] f_at_points.interpolate(field) f_at_input_points.interpolate(f_at_points) - return f_at_input_points.dat.data_ro + return f_at_input_points.dat.data_ro[:] def __call__(self): """ @@ -694,7 +681,6 @@ def __init__(self, solver_obj, fieldnames, x, y, self.tolerance = tolerance self.eval_func = eval_func self._initialized = False - self.bathymetry_val = None @PETSc.Log.EventDecorator("thetis.TimeSeriesCallback2D._initialize") def _initialize(self): @@ -706,31 +692,16 @@ def _initialize(self): xyz = (self.x, self.y, self.z) if self.on_sphere else (self.x, self.y) self.xyz = numpy.array([xyz]) - # create VertexOnlyMesh and function spaces - self.mesh = self.solver_obj.fields[self.fieldnames[0]].function_space().mesh() - self.vom = VertexOnlyMesh(self.mesh, self.xyz, redundant=True) - self.function_spaces = {} - for field_name in self.fieldnames: - field = self.solver_obj.fields[field_name] - if isinstance(field.function_space().ufl_element(), VectorElement): - P0DG = VectorFunctionSpace(self.vom, "DG", 0) - P0DG_input_ordering = VectorFunctionSpace(self.vom.input_ordering, "DG", 0) - else: - P0DG = FunctionSpace(self.vom, "DG", 0) - P0DG_input_ordering = FunctionSpace(self.vom.input_ordering, "DG", 0) - self.function_spaces[field_name] = (P0DG, P0DG_input_ordering) + # initialise interpolation functions using vom_interpolator_functions + self.interp_functions = vom_interpolator_functions(self.solver_obj, self.fieldnames, self.xyz) # test evaluation try: if self.eval_func is None: - bathymetry_field = self.solver_obj.fields.bathymetry_2d - P0DG = FunctionSpace(self.vom, "DG", 0) - P0DG_input_ordering = FunctionSpace(self.vom.input_ordering, "DG", 0) - f_at_points = Function(P0DG) - f_at_input_points = Function(P0DG_input_ordering) - f_at_points.interpolate(bathymetry_field) + field = self.solver_obj.fields[self.fieldnames[0]] + f_at_points, f_at_input_points = self.interp_functions[self.fieldnames[0]] + f_at_points.interpolate(field) f_at_input_points.interpolate(f_at_points) - self.bathymetry_val = f_at_input_points.dat.data_ro[:] else: self.eval_func(self.solver_obj.fields.bathymetry_2d, self.xyz, tolerance=self.tolerance) except PointNotInDomainError as e: @@ -750,9 +721,7 @@ def __call__(self): try: field = self.solver_obj.fields[fieldname] if self.eval_func is None: - P0DG, P0DG_input_ordering = self.function_spaces[fieldname] - f_at_points = Function(P0DG) - f_at_input_points = Function(P0DG_input_ordering) + f_at_points, f_at_input_points = self.interp_functions[fieldname] f_at_points.interpolate(field) f_at_input_points.interpolate(f_at_points) val = f_at_input_points.dat.data_ro[:] diff --git a/thetis/utility.py b/thetis/utility.py index d822f432e..f5ea1898f 100755 --- a/thetis/utility.py +++ b/thetis/utility.py @@ -1154,3 +1154,38 @@ def form2indicator(F): }, ) return indicator + + +@PETSc.Log.EventDecorator("thetis.vom_interpolator_functions") +def vom_interpolator_functions(solver_obj, field_names, locations): + r""" + Creates function spaces and associated Functions for interpolation + on a VertexOnlyMesh (VOM) and returns them for reuse. + + :arg solver_obj: Thetis solver object + :arg field_names: List of field names to create functions for. + :arg locations: List of locations for interpolation. + :return: A dictionary mapping field names to a tuple of (f_at_points, f_at_input_points) + which are Functions for interpolation. + """ + vom = VertexOnlyMesh(solver_obj.mesh2d, locations, redundant=True) + + functions_dict = {} + + for field_name in field_names: + field = solver_obj.fields[field_name] + + if isinstance(field.function_space().ufl_element(), VectorElement): + P0DG = VectorFunctionSpace(vom, "DG", 0) + P0DG_input_ordering = VectorFunctionSpace(vom.input_ordering, "DG", 0) + else: + P0DG = FunctionSpace(vom, "DG", 0) + P0DG_input_ordering = FunctionSpace(vom.input_ordering, "DG", 0) + + f_at_points = Function(P0DG) + f_at_input_points = Function(P0DG_input_ordering) + + # Store the Functions in the dictionary keyed by field name + functions_dict[field_name] = (f_at_points, f_at_input_points) + + return functions_dict