From 3e2d75af62af68f2587d71c458ab6695ea52fa20 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Tue, 13 Aug 2024 10:23:07 +0800 Subject: [PATCH 1/3] Update _utils.py update gather data fn: Uses a mask value that will work for most dtypes (except bool?) Automatically gets dtype from array/value being passed in --- src/underworld3/utilities/_utils.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/underworld3/utilities/_utils.py b/src/underworld3/utilities/_utils.py index cc68c73e9..c1fffc4d3 100755 --- a/src/underworld3/utilities/_utils.py +++ b/src/underworld3/utilities/_utils.py @@ -90,12 +90,14 @@ def mem_footprint(): return python_process.memory_info().rss // 1000000 -def gather_data(val, bcast=False, dtype="float64"): +def gather_data(val, masked_value=-999999, bcast=False): """ gather values on root (bcast=False) or all (bcast = True) processors Parameters: vals : Values to combine into a single array on the root or all processors + masked_value : value to use as mask value + bcast : broadcast array/value to all processors returns: val_global : combination of values form all processors @@ -106,12 +108,15 @@ def gather_data(val, bcast=False, dtype="float64"): rank = uw.mpi.rank size = uw.mpi.size + dtype = val.dtype + + ### make sure all data comes in the same order - with uw.mpi.call_pattern(pattern="sequential"): - if len(val > 0): - val_local = np.ascontiguousarray(val.copy()) - else: - val_local = np.array([np.nan], dtype=dtype) + # with uw.mpi.call_pattern(pattern="sequential"): + if len(val > 0): + val_local = np.ascontiguousarray(val.copy(), dtype=dtype) + else: + val_local = np.ascontiguousarray([masked_value], dtype=dtype) comm.barrier() @@ -134,7 +139,7 @@ def gather_data(val, bcast=False, dtype="float64"): if uw.mpi.rank == 0: ### remove rows with NaN - val_global = val_global[~np.isnan(val_global)] + val_global = val_global[val_global != masked_value] comm.barrier() From b37bc205fccba717015971afe08905f6ce7a32bc Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Thu, 15 Aug 2024 18:05:39 +0800 Subject: [PATCH 2/3] Update discretisation.py update for changes in gather_data fn --- src/underworld3/discretisation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/underworld3/discretisation.py b/src/underworld3/discretisation.py index a8fdfbef3..5bec47817 100644 --- a/src/underworld3/discretisation.py +++ b/src/underworld3/discretisation.py @@ -440,7 +440,7 @@ def view(self): else: i = 0 - ii = uw.utilities.gather_data(np.array([i]), dtype="int") + ii = uw.utilities.gather_data(np.array([i])) if uw.mpi.rank == 0: print( From d651d49d9a6050946aae12ad7676caafd23cee28 Mon Sep 17 00:00:00 2001 From: Ben Knight <55677727+bknight1@users.noreply.github.com> Date: Thu, 26 Sep 2024 11:04:37 +0800 Subject: [PATCH 3/3] add L2 norm function For benchmarking purposes --- src/underworld3/maths/functions.py | 50 +++++++++++++++++++++++++++++ src/underworld3/utilities/_utils.py | 4 +-- 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/src/underworld3/maths/functions.py b/src/underworld3/maths/functions.py index 599230cdd..720a4eda3 100644 --- a/src/underworld3/maths/functions.py +++ b/src/underworld3/maths/functions.py @@ -4,6 +4,7 @@ from sympy import sympify from typing import Optional, Callable from underworld3 import function +from underworld3.cython.petsc_maths import Integral def delta( @@ -15,3 +16,52 @@ def delta( delta_fn = sympy.exp(-(x**2) / (2 * epsilon**2)) / (epsilon * sqrt_2pi) return delta_fn + + + +def L2_norm(U_exp, A_exp, mesh): + ''' + Compute the L2 norm between the computed and analytical fields, both of which are sympy.Matrix objects. + + The L2-norm is defined as: + + .. math:: + L_2\text{-norm} = \left( \int_{\Omega} |u_{\text{numeric}}(x) - u_{\text{analytic}}(x)|^2 \, dx \right)^{1/2} + + Where: + .. math:: + - \( u_{\text{numeric}}(x) \) is the numerical solution. + - \( u_{\text{analytic}}(x) \) is the analytical solution. + - \( \Omega \) is the domain. + - The integral computes the total squared error across the entire domain, and the square root is applied after the integration to give the L2-norm. + + Parameters: + U_exp : sympy.Matrix representing a scalar (1x1) or vector field (Nx1 or 1xN) + A_exp : sympy.Matrix representing a scalar (1x1) or vector field (Nx1 or 1xN) + mesh : the mesh over which to integrate + + Returns: + L2_norm : the computed L2 norm + ''' + # Check if the inputs are matrices + if isinstance(U_exp, sympy.Matrix) and isinstance(A_exp, sympy.Matrix): + # Ensure that the dimensions of U_exp and A_exp match + assert U_exp.shape == A_exp.shape, "U_exp and A_exp must have the same shape." + + # Initialize the squared difference + squared_diff = 0 + + # Loop over the components of the matrices (scalar or vector case) + for i in range(U_exp.shape[0]): + squared_diff += (U_exp[i] - A_exp[i])**2 + + else: + raise TypeError("U_exp and A_exp must be sympy.Matrix objects.") + + # Perform the integral over the mesh + I = Integral(mesh, squared_diff) + + # Compute the L2 norm (sqrt of the integral result) + L2_norm = sympy.sqrt( I.evaluate() ) + + return L2_norm diff --git a/src/underworld3/utilities/_utils.py b/src/underworld3/utilities/_utils.py index c1fffc4d3..4c0c5675e 100755 --- a/src/underworld3/utilities/_utils.py +++ b/src/underworld3/utilities/_utils.py @@ -132,13 +132,13 @@ def gather_data(val, masked_value=-999999, bcast=False): comm.barrier() - ## gather x values, can't do them together + ## gather the data on the root proc comm.Gatherv(sendbuf=val_local, recvbuf=(val_global, sendcounts), root=0) comm.barrier() if uw.mpi.rank == 0: - ### remove rows with NaN + ### remove rows with mask value val_global = val_global[val_global != masked_value] comm.barrier()