diff --git a/docs/advanced/index.md b/docs/advanced/index.md index 41655844f..e83b680e2 100644 --- a/docs/advanced/index.md +++ b/docs/advanced/index.md @@ -74,4 +74,5 @@ curved-boundary-conditions mesh-adaptation troubleshooting api-patterns +SWARM-INTEGRATION-STATISTICS ``` \ No newline at end of file diff --git a/docs/api/constitutive_models.md b/docs/api/constitutive_models.md index dcc27bf84..4eb974f6b 100644 --- a/docs/api/constitutive_models.md +++ b/docs/api/constitutive_models.md @@ -1,10 +1,5 @@ # Constitutive Models -```{eval-rst} -.. automodule:: underworld3.constitutive_models - :no-members: -``` - ## Viscous Flow Models ### ViscousFlowModel diff --git a/docs/api/coordinates.md b/docs/api/coordinates.md index 39f26b920..522cda5b1 100644 --- a/docs/api/coordinates.md +++ b/docs/api/coordinates.md @@ -1,10 +1,5 @@ # Coordinate Systems -```{eval-rst} -.. automodule:: underworld3.coordinates - :no-members: -``` - ## CoordinateSystem The base class for coordinate system handling on meshes. diff --git a/docs/api/discretisation.md b/docs/api/discretisation.md index d47e1c0f7..d2adcacbb 100644 --- a/docs/api/discretisation.md +++ b/docs/api/discretisation.md @@ -1,10 +1,5 @@ # Discretisation -```{eval-rst} -.. automodule:: underworld3.discretisation - :no-members: -``` - ## Mesh The computational mesh class that provides the spatial discretisation for finite element calculations. diff --git a/docs/api/function.md b/docs/api/function.md index 4b948370e..b3a052b37 100644 --- a/docs/api/function.md +++ b/docs/api/function.md @@ -1,16 +1,11 @@ # Function and Expressions -```{eval-rst} -.. automodule:: underworld3.function - :no-members: -``` - ## Expressions ### UWexpression ```{eval-rst} -.. autoclass:: underworld3.function.UWexpression +.. autoclass:: underworld3.function.expressions.UWexpression :members: :show-inheritance: ``` @@ -75,12 +70,6 @@ Factory function for creating UWQuantity objects with units. .. autofunction:: underworld3.function.global_evaluate ``` -### evalf - -```{eval-rst} -.. autofunction:: underworld3.function.evalf -``` - ### evaluate_gradient ```{eval-rst} diff --git a/docs/api/maths.md b/docs/api/maths.md index 106e45ba9..8b1635d1f 100644 --- a/docs/api/maths.md +++ b/docs/api/maths.md @@ -1,11 +1,5 @@ # Mathematical Operations -```{eval-rst} -.. automodule:: underworld3.maths - :members: - :show-inheritance: -``` - ## Integral Operations ```{eval-rst} diff --git a/docs/api/solvers.md b/docs/api/solvers.md index ad6add24f..6b5c464f3 100644 --- a/docs/api/solvers.md +++ b/docs/api/solvers.md @@ -1,10 +1,5 @@ # Solvers -```{eval-rst} -.. automodule:: underworld3.systems.solvers - :no-members: -``` - ## Stokes Flow ### SNES_Stokes diff --git a/docs/api/swarm.md b/docs/api/swarm.md index 8486fe160..8b9fd1b1d 100644 --- a/docs/api/swarm.md +++ b/docs/api/swarm.md @@ -1,10 +1,5 @@ # Particle Swarms -```{eval-rst} -.. automodule:: underworld3.swarm - :no-members: -``` - ## Swarm ```{eval-rst} @@ -21,14 +16,6 @@ :show-inheritance: ``` -## Population Control - -```{eval-rst} -.. autoclass:: underworld3.swarm.PopulationControl - :members: - :show-inheritance: -``` - ## Index Swarm Variable ```{eval-rst} diff --git a/docs/api/systems_ddt.md b/docs/api/systems_ddt.md index b402bc9dc..1bc7e9428 100644 --- a/docs/api/systems_ddt.md +++ b/docs/api/systems_ddt.md @@ -1,10 +1,5 @@ # Time Derivatives -```{eval-rst} -.. automodule:: underworld3.systems.ddt - :no-members: -``` - Time derivative operators approximate $D\phi/Dt$ or $DF/Dt$ for transient solvers. All operators share a common interface: ``update_pre_solve(dt)`` before each timestep, ``bdf()`` for the BDF approximation in the weak form, diff --git a/docs/beginner/quickstart.md b/docs/beginner/quickstart.md index e5dccc8c0..1a179618f 100644 --- a/docs/beginner/quickstart.md +++ b/docs/beginner/quickstart.md @@ -67,7 +67,7 @@ $$ where $t(N_{ref})$ is the run-time for a reference number of processors, $N_{ref}$, and $t(N)$ is the run-time when $N$ processors are used. In the ideal case, $N$ additional processors should contribute all of its resources in solving the problem and reduce the compute time by a factor of $N$ relative to the reference run time. For example, using $2 N_{ref}$ processors will ideally halve the run-time resulting to a speed-up = 2. -```{figure} media/UW3-StrongScalingSolvers.png +```{figure} ../media/UW3-StrongScalingSolvers.png :name: fig-strong-scaling Strong parallel-scaling tests run on Australia's peak computing system, [GADI, at the National Computational Infrastructure](https://nci.org.au/our-systems/hpc-systems?ref=underworldcode.org). This is a typical High Performance Computing facility with large numbers of dedicated, identical CPUs and fast communication links. diff --git a/docs/developer/UW3_Developers_MathematicalObjects.md b/docs/developer/UW3_Developers_MathematicalObjects.md index 2178b5ae2..0f4a323ab 100644 --- a/docs/developer/UW3_Developers_MathematicalObjects.md +++ b/docs/developer/UW3_Developers_MathematicalObjects.md @@ -574,8 +574,7 @@ def test_mathematical_operations(): ## Planned Enhancements ### 1. Operator Overloading -```python - +```text # Future possibilities laplacian = velocity @ nabla**2 # Operator syntax curl = nabla × velocity # Cross product notation diff --git a/docs/developer/container/README.md b/docs/developer/container/README.md index 6d7a77885..df60be618 100644 --- a/docs/developer/container/README.md +++ b/docs/developer/container/README.md @@ -1,3 +1,7 @@ +--- +orphan: true +--- + The `Containerfile` is used to build the container and can be used with docker or podman. At present only amd64 architecture is built, because vtk-osmesa isn't available for arm by default. This may change in future, or we could build vtk-osmesa from source (see https://docs.pyvista.org/extras/building_vtk.html) diff --git a/docs/developer/design/COORDINATE_MIGRATION_GUIDE.md b/docs/developer/design/COORDINATE_MIGRATION_GUIDE.md index d38b4c319..63d964f2f 100644 --- a/docs/developer/design/COORDINATE_MIGRATION_GUIDE.md +++ b/docs/developer/design/COORDINATE_MIGRATION_GUIDE.md @@ -7,18 +7,18 @@ ### Pattern Migration Cheat Sheet -```python -# ❌ DEPRECATED → ✅ RECOMMENDED -mesh.data → mesh.X.coords or mesh.points -mesh.data.shape → mesh.X.coords.shape -mesh.data[:, 0] → mesh.X.coords[:, 0] -mesh.data.min() → mesh.X.coords.min() -mesh.data.copy() → mesh.X.coords.copy() - -# ✅ UNCHANGED (these work exactly as before) -mesh.X[0] → mesh.X[0] (symbolic x-coordinate) -x, y = mesh.X → x, y = mesh.X (symbolic unpacking) -mesh.units → mesh.units (coordinate units) +```text +DEPRECATED RECOMMENDED +mesh.data -> mesh.X.coords or mesh.points +mesh.data.shape -> mesh.X.coords.shape +mesh.data[:, 0] -> mesh.X.coords[:, 0] +mesh.data.min() -> mesh.X.coords.min() +mesh.data.copy() -> mesh.X.coords.copy() + +UNCHANGED (these work exactly as before) +mesh.X[0] -> mesh.X[0] (symbolic x-coordinate) +x, y = mesh.X -> x, y = mesh.X (symbolic unpacking) +mesh.units -> mesh.units (coordinate units) ``` ## Migration Examples by Use Case diff --git a/docs/developer/design/WHY_UNITS_NOT_DIMENSIONALITY.md b/docs/developer/design/WHY_UNITS_NOT_DIMENSIONALITY.md index 001438206..aaeb8f1ae 100644 --- a/docs/developer/design/WHY_UNITS_NOT_DIMENSIONALITY.md +++ b/docs/developer/design/WHY_UNITS_NOT_DIMENSIONALITY.md @@ -20,7 +20,7 @@ - **Dimensionality**: ⚠️ Abstract physics concept ([length], [time], [mass]) ### What the API Shows: -```python +```text # User-friendly - clear what's happening var.units = "m" var.to("km") @@ -345,7 +345,7 @@ from underworld3.utilities.units_mixin import UnitAwareMixin - What does that even mean? ### 3. API Clarity ✅ -```python +```text # Clear and intuitive var.units = "m" var.to("km") diff --git a/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md b/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md index 5499bed6f..4070330f4 100644 --- a/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md +++ b/docs/developer/guides/HOW-TO-WRITE-UW3-SCRIPTS.md @@ -612,7 +612,7 @@ def test_swarm_functionality(): If you see generated C code with symbolic expressions instead of numbers: -```c +```text // ERROR symptom in generated code: out[0] = 1.0/{ \eta \hspace{ 0.0006pt } }; // Should be numeric! ``` diff --git a/docs/developer/subsystems/constitutive-models-anisotropy.md b/docs/developer/subsystems/constitutive-models-anisotropy.md index f43659277..994c24576 100644 --- a/docs/developer/subsystems/constitutive-models-anisotropy.md +++ b/docs/developer/subsystems/constitutive-models-anisotropy.md @@ -45,7 +45,7 @@ In `underworld` tensor rotation is provided for rank 2 and rank 4 tensors by `uw If we apply constraints to the deformation, we expect to reduce the number of independent material constants. Incompressibility should reduce the number of *independent* material constants. In an isotropic medium (or a medium with cubic symmetry), incompressibility eliminates the one volumetric material modulus (e.g. the bulk modulus). In general anisotropic media, it is not the case that changes in pressure result in uniform expansion or contraction, and an incompressibility constraint reduces the number of *independent* material constants. In the transversely isotropic case, there are five independent materia constants in general, reducing to four when the material is incompressible. -It is not a given that the stiffness matrix is trivial to construct / meaningful for incompressible anisotropy and there is some discussion here: https://rastgaragah.wordpress.com/2013/03/12/incompressibility-of-linearly-elastic-material/ (identifies the issue) and this is explained in more detail by {cite}`destradeIncompressibleLimitLinear2002` +It is not a given that the stiffness matrix is trivial to construct / meaningful for incompressible anisotropy and there is some discussion here: https://rastgaragah.wordpress.com/2013/03/12/incompressibility-of-linearly-elastic-material/ (identifies the issue) and this is explained in more detail by Destrade et al. (2002) ### Example @@ -192,7 +192,7 @@ display(C_IJm_MM_iso) ## Han & Wahr, 1997 (full transverse isotropic tensor) -In the {cite}`hanAnalysisAnisotropicMantle1997` paper, the expression for incompressible transverse-isotropy is as follows +In the Han & Wahr (1997) paper, the expression for incompressible transverse-isotropy is as follows $$\left[\begin{matrix}2 \eta_{0} + \mu_{0} & \mu_{0} & 0 & 0 & 0 & 0\\\mu_{0} & 2 \eta_{0} + \mu_{0} & 0 & 0 & 0 & 0\\0 & 0 & - 2 \Delta\eta + 2 \eta_{0} + \mu_{1} & 0 & 0 & 0\\0 & 0 & 0 & - 2 \Delta\eta + 2 \eta_{0} & 0 & 0\\0 & 0 & 0 & 0 & - 2 \Delta\eta + 2 \eta_{0} & 0\\0 & 0 & 0 & 0 & 0 & 2 \eta_{0}\end{matrix}\right]$$ @@ -249,7 +249,7 @@ display(C_IJm_HW_R_s2) ## Orthotropic medium -**Note** all the caveats above regarding incompressibility. The {cite}`browaeysDecompositionElasticTensor2004` elastic tensors have a bulk modulus term, so it is not completely obvious how to square the assumptions in the first two implementations with this set. +**Note** all the caveats above regarding incompressibility. The Browaeys & Chevrot (2004) elastic tensors have a bulk modulus term, so it is not completely obvious how to square the assumptions in the first two implementations with this set. The full formulation should look like this: diff --git a/src/underworld3/discretisation/__init__.py b/src/underworld3/discretisation/__init__.py index 782fc7dc7..30ce1d1a9 100644 --- a/src/underworld3/discretisation/__init__.py +++ b/src/underworld3/discretisation/__init__.py @@ -4,20 +4,13 @@ This module provides the core discretisation infrastructure for finite element computations in Underworld3. -Classes -------- -Mesh : class - Unstructured mesh with PETSc DMPlex backend. Supports various cell - types and coordinate systems. -MeshVariable : class - Field variable defined on mesh (nodal or cell-based storage). +**Mesh** -- Unstructured mesh with PETSc DMPlex backend. -Functions ---------- -checkpoint_xdmf : function - Save mesh and variables to XDMF format for visualization. -meshVariable_lookup_by_symbol : function - Find mesh variable by its symbolic representation. +**MeshVariable** -- Field variable defined on mesh (nodal or cell-based storage). + +**checkpoint_xdmf** -- Save mesh and variables to XDMF format. + +**meshVariable_lookup_by_symbol** -- Find mesh variable by its symbol. See Also -------- diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index cd8a717eb..5fb4f5279 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -1869,11 +1869,11 @@ def write_checkpoint( index: Optional[int] = 0, unique_id: Optional[bool] = False, ): - """Write data in a format that can be restored for restarting the simulation + """Write data in a format that can be restored for restarting the simulation. + The difference between this and the visualisation is 1) the parallel section needs to be stored to reload the data correctly, and 2) the visualisation information (vertex form of fields) - is not stored. This routines uses dmplex *VectorView and *VectorLoad functionality - + is not stored. This routine uses dmplex VectorView and VectorLoad functionality. """ # The mesh checkpoint is the same as the one required for visualisation diff --git a/src/underworld3/function/__init__.py b/src/underworld3/function/__init__.py index f18c64849..eb09d37e1 100644 --- a/src/underworld3/function/__init__.py +++ b/src/underworld3/function/__init__.py @@ -5,21 +5,16 @@ on meshes and swarms. It bridges SymPy symbolic mathematics with PETSc numerical evaluation. -Key Components --------------- -expression : class - User-facing symbolic expression wrapper (UWexpression). -evaluate : function - Evaluate expressions at mesh/swarm points. -global_evaluate : function - Parallel-safe evaluation gathering results across MPI ranks. -UnderworldFunction : class - Core function evaluation machinery. - -Unit Conversion ---------------- -convert_quantity_units, make_dimensionless, add_units - Utilities for handling physical units in expressions. +**expression** -- User-facing symbolic expression wrapper (UWexpression). + +**evaluate** -- Evaluate expressions at mesh/swarm points. + +**global_evaluate** -- Parallel-safe evaluation gathering results across MPI ranks. + +**UnderworldFunction** -- Core function evaluation machinery. + +**convert_quantity_units**, **make_dimensionless**, **add_units** -- +Utilities for handling physical units in expressions. The function module integrates with pint for optional unit-aware computations when the mesh has an associated unit registry. diff --git a/src/underworld3/function/_function.pyx b/src/underworld3/function/_function.pyx index 7ad5820de..7913208d1 100644 --- a/src/underworld3/function/_function.pyx +++ b/src/underworld3/function/_function.pyx @@ -133,18 +133,16 @@ class UnderworldAppliedFunctionDeriv(UnderworldAppliedFunction): class UnderworldFunction(sympy.Function): """ - This is a metaclass, so it returns programmatic class objects rather - than instances. This basically follows the pattern of the `sympy.Function` - metaclass, with two key differences: - 1. We set `UnderworldAppliedFunction` as the base class. This is really just a - dummy class (see its definition) which allows us to do calls such - as `isinstance(someobj, UnderworldAppliedFunction)` to test if a `sympy` - object is one we've defined. - 2. We grab a weakref of the owning meshvariable onto the *class* itself. Note - that it's important that it's recorded onto the class (instead of the instance), - as Sympy internally sometimes uses calls such as `type(obj)(obj.args)` to - replace objects with cloned instances, and therefore 'extra' info must be - recorded onto the class so that the clones are _complete_. + Metaclass that returns programmatic class objects rather than instances. + + This basically follows the pattern of the ``sympy.Function`` + metaclass, with two key differences. First, we set + ``UnderworldAppliedFunction`` as the base class, which allows + ``isinstance(someobj, UnderworldAppliedFunction)`` checks. + Second, we grab a weakref of the owning meshvariable onto the + class itself (not the instance), because SymPy internally uses + ``type(obj)(obj.args)`` to clone instances and extra info must + be on the class so that clones are complete. Consider the calling pattern @@ -157,16 +155,16 @@ class UnderworldFunction(sympy.Function): Parameters ---------- - name: + name : str The name of the function. - meshvar: + meshvar : MeshVariable The mesh variable corresponding to this function. - vtype: - The variable type (scalar,vector,etc). - component: + vtype : VarType + The variable type (scalar, vector, etc). + component : int or tuple For vector functions, this is the component of the vector. - For example, component `1` might correspond to `v_y`. - For tensors, the component is a tuple + For example, component ``1`` might correspond to ``v_y``. + For tensors, the component is a tuple. For scalars, this value is ignored. """ def __new__(cls, diff --git a/src/underworld3/function/functions_unit_system.py b/src/underworld3/function/functions_unit_system.py index 2705d4eee..fddc8d873 100644 --- a/src/underworld3/function/functions_unit_system.py +++ b/src/underworld3/function/functions_unit_system.py @@ -72,14 +72,11 @@ def evaluate( evalf : bool, optional Force numerical evaluation via sympy evalf (default: False) mode : str, optional - Evaluation mode controlling accuracy vs speed tradeoff: - - "default": Accurate evaluation. Projection for derivatives (O(h²)), - direct calculation otherwise. DMInterpolation inside mesh, RBF outside. - - "fast": Quick visualization mode. Clement gradient recovery for - derivatives (O(h), no solve), RBF interpolation everywhere. - - "projection": Always use L2 projection (even without derivatives), - DMInterpolation inside mesh, RBF outside. - Default: "default" + Evaluation mode controlling accuracy vs speed tradeoff. + Options: ``"default"`` (accurate, projection for derivatives), + ``"fast"`` (Clement gradient, RBF everywhere), + ``"projection"`` (always L2 projection). + Default: ``"default"`` data_layout : callable, optional Data layout specification (default: None) check_extrapolated : bool, optional @@ -95,10 +92,9 @@ def evaluate( Returns ------- UWQuantity, UnitAwareArray, or ndarray - - If non-dimensional scaling is active: plain ndarray (non-dimensional) - - If expression has units and result is scalar: UWQuantity - - If expression has units and result is array: UnitAwareArray - - If expression has no units: plain ndarray (as before) + If non-dimensional scaling is active, returns plain ndarray. + If expression has units, returns UWQuantity (scalar) or UnitAwareArray. + Otherwise returns plain ndarray. Notes ----- @@ -402,14 +398,11 @@ def global_evaluate( evalf : bool, optional Force numerical evaluation via sympy evalf (default: False) mode : str, optional - Evaluation mode controlling accuracy vs speed tradeoff: - - "default": Accurate evaluation. Projection for derivatives (O(h²)), - direct calculation otherwise. DMInterpolation inside mesh, RBF outside. - - "fast": Quick visualization mode. Clement gradient recovery for - derivatives (O(h), no solve), RBF interpolation everywhere. - - "projection": Always use L2 projection (even without derivatives), - DMInterpolation inside mesh, RBF outside. - Default: "default" + Evaluation mode controlling accuracy vs speed tradeoff. + Options: ``"default"`` (accurate, projection for derivatives), + ``"fast"`` (Clement gradient, RBF everywhere), + ``"projection"`` (always L2 projection). + Default: ``"default"`` data_layout : callable, optional Data layout specification (default: None) check_extrapolated : bool, optional @@ -425,8 +418,8 @@ def global_evaluate( Returns ------- UWQuantity, UnitAwareArray, or ndarray - - If non-dimensional scaling is active: plain ndarray (non-dimensional) - - Otherwise: result with appropriate unit tracking + If non-dimensional scaling is active, returns plain ndarray. + Otherwise returns result with appropriate unit tracking. Notes ----- diff --git a/src/underworld3/materials.py b/src/underworld3/materials.py index 4ba1ddf0d..b9a9be69f 100644 --- a/src/underworld3/materials.py +++ b/src/underworld3/materials.py @@ -294,10 +294,10 @@ def add_callback(self, callback: Callable): """ Add a callback function for material changes. - Parameters: - ----------- + Parameters + ---------- callback : callable - Function called as callback(event_type, *args) + Function called as callback(event_type, \*args) """ self._callbacks.append(callback) diff --git a/src/underworld3/maths/__init__.py b/src/underworld3/maths/__init__.py index 5e99b6f3c..ffa774e3a 100644 --- a/src/underworld3/maths/__init__.py +++ b/src/underworld3/maths/__init__.py @@ -16,16 +16,13 @@ vector_calculus_spherical : module Vector calculus operators in spherical coordinates. -Functions ---------- -delta_function : function - Regularized delta function for localized source terms. -L2_norm : function - L2 norm computation for fields. -Integral : class - Domain integration over mesh. -CellWiseIntegral : class - Cell-by-cell integration. +**delta_function** -- Regularized delta function for localized source terms. + +**L2_norm** -- L2 norm computation for fields. + +**Integral** -- Domain integration over mesh. + +**CellWiseIntegral** -- Cell-by-cell integration. See Also -------- diff --git a/src/underworld3/meshing/annulus.py b/src/underworld3/meshing/annulus.py index 9a9a861b0..1a3631d67 100644 --- a/src/underworld3/meshing/annulus.py +++ b/src/underworld3/meshing/annulus.py @@ -559,60 +559,45 @@ def SegmentofAnnulus( Generates a segment of an annulus using Gmsh. This function creates a 2D mesh of an annular segment defined by outer and inner radii, and the extent of the angle. The mesh can be customized with various parameters like cell size, element degree, and verbosity. - Parameters: - ----------- + Parameters + ---------- radiusOuter : float, optional The outer radius of the annular segment. Default is 1.0. - radiusInner : float, optional The inner radius of the annular segment. Default is 0.547. - angleExtent : float, optional The angular extent of the segment in degrees. Default is 45. - cellSize : float, optional - The target size for the mesh elements. This controls the density of the mesh. Default is 0.1. - + The target size for the mesh elements. Default is 0.1. centre : bool, optional - If True, the segment will be centered at the origin. If False, the segment is positioned based on the radii. Default is False. - + If True, the segment will be centered at the origin. Default is False. degree : int, optional - The polynomial degree of the finite elements used in the mesh. Default is 1. - + The polynomial degree of the finite elements. Default is 1. qdegree : int, optional - The quadrature degree for integration. Higher values may improve accuracy but increase computation time. Default is 2. - + The quadrature degree for integration. Default is 2. filename : str, optional - The name of the file where the mesh will be saved. If None, a default name is generated based on the parameters. Default is None. - + The name of the file where the mesh will be saved. Default is None. refinement : optional - Refinement level or method for the mesh. Used to increase the resolution of the mesh in certain regions. Default is None. - + Refinement level for the mesh. Default is None. gmsh_verbosity : int, optional - Controls the verbosity of Gmsh output. Set to 0 for minimal output, higher numbers for more detailed logs. Default is 0. - + Gmsh output verbosity (0=quiet). Default is 0. verbose : bool, optional - If True, the function prints additional information during execution. Default is False. + If True, print additional information. Default is False. - Returns: - -------- - None - The function generates and saves a mesh file according to the specified parameters. + Returns + ------- + Mesh + The generated annular segment mesh. - Example: + Examples -------- - mesh = uw.meshing.SegmentofAnnulus( - radiusOuter=2.0, - radiusInner=1.0, - angleExtent=90.0, - cellSize=0.05, - centre=True, - degree=2, - qdegree=3, - filename="custom_annulus_segment.msh", - gmsh_verbosity=1, - verbose=True - ) + >>> mesh = uw.meshing.SegmentofAnnulus( + ... radiusOuter=2.0, + ... radiusInner=1.0, + ... angleExtent=90.0, + ... cellSize=0.05, + ... centre=True, + ... ) """ # Convert unit-aware quantities to non-dimensional units radiusOuter = uw.scaling.non_dimensionalise(radiusOuter) diff --git a/src/underworld3/meshing/faults.py b/src/underworld3/meshing/faults.py index d59ade94d..17f4f16f7 100644 --- a/src/underworld3/meshing/faults.py +++ b/src/underworld3/meshing/faults.py @@ -477,18 +477,27 @@ def compute_distance_field( computation. The returned field contains the absolute distance to the nearest fault surface at each mesh point. - Args: - mesh: The mesh to compute distances on - distance_var: Optional existing MeshVariable to store results. - If None, creates a new variable. - variable_name: Name for new variable if distance_var is None - - Returns: - MeshVariable with distance values (scalar, 1 component) - - Raises: - ValueError: If collection is empty or no faults are triangulated - ImportError: If pyvista not available + Parameters + ---------- + mesh : Mesh + The mesh to compute distances on. + distance_var : MeshVariable, optional + Existing MeshVariable to store results. + If None, creates a new variable. + variable_name : str, optional + Name for new variable if distance_var is None. + + Returns + ------- + MeshVariable + Scalar variable with distance values (1 component). + + Raises + ------ + ValueError + If collection is empty or no faults are triangulated. + ImportError + If pyvista is not available. """ pv = _require_pyvista() @@ -544,18 +553,27 @@ def transfer_normals( For each mesh point, finds the closest fault face (from any fault in the collection) and copies that face's normal vector. - Args: - mesh: The mesh to transfer normals to - coords: Optional coordinates to query. If None, uses mesh.X.coords - normal_var: Optional existing MeshVariable to store results. - If None, creates a new variable. - variable_name: Name for new variable if normal_var is None - - Returns: - MeshVariable with normal vectors (3 components) - - Raises: - ValueError: If collection is empty or no faults are triangulated + Parameters + ---------- + mesh : Mesh + The mesh to transfer normals to. + coords : ndarray, optional + Coordinates to query. If None, uses mesh.X.coords. + normal_var : MeshVariable, optional + Existing MeshVariable to store results. + If None, creates a new variable. + variable_name : str, optional + Name for new variable if normal_var is None. + + Returns + ------- + MeshVariable + Variable with normal vectors (3 components). + + Raises + ------ + ValueError + If collection is empty or no faults are triangulated. """ if len(self.faults) == 0: raise ValueError("Cannot transfer normals: no faults in collection") diff --git a/src/underworld3/meshing/geographic.py b/src/underworld3/meshing/geographic.py index 7d0fefde6..77836e94c 100644 --- a/src/underworld3/meshing/geographic.py +++ b/src/underworld3/meshing/geographic.py @@ -488,44 +488,29 @@ def RegionalGeographicBox( Examples -------- - # Create mesh for southeastern Australia, 0-400 km depth - mesh = uw.meshing.RegionalGeographicBox( - lon_range=(135, 140), - lat_range=(-35, -30), - depth_range=(0, 400), - ellipsoid='WGS84', - numElements=(20, 20, 10), - ) + Create mesh for southeastern Australia, 0-400 km depth:: + + mesh = uw.meshing.RegionalGeographicBox( + lon_range=(135, 140), + lat_range=(-35, -30), + depth_range=(0, 400), + ellipsoid='WGS84', + numElements=(20, 20, 10), + ) - # Access geographic coordinates - lon = mesh.geo.lon # Longitude array (degrees East) - lat = mesh.geo.lat # Latitude array (degrees North) - depth = mesh.geo.depth # Depth array (km below surface) - - # Use in equations - λ_lon, λ_lat, λ_d = mesh.geo[:] - T = 1600 - 0.5 * λ_d # Temperature decreasing with depth - - # Basis vectors for boundary conditions - v_surface = 0 * mesh.geo.unit_up # No vertical flow at surface - v_bottom = 10 * mesh.geo.unit_down # Downward flow at bottom - - # Mars example - mesh_mars = uw.meshing.RegionalGeographicBox( - lon_range=(0, 45), - lat_range=(-22.5, 22.5), - depth_range=(0, 200), - ellipsoid='Mars', - numElements=(15, 15, 8), - ) + Access geographic coordinates:: + + lon = mesh.geo.lon + lat = mesh.geo.lat + depth = mesh.geo.depth Notes ----- - - Uses geodetic latitude (GPS/map standard), not geocentric latitude - - Depth is measured from reference ellipsoid surface, not from center - - mesh.R provides spherical coordinates $(r, \\theta, \\phi)$ for backward compatibility - - mesh.geo provides geographic coordinates (lon, lat, depth) with ellipsoid geometry - - Right-handed coordinate system: WE × SN = down + Uses geodetic latitude (GPS/map standard), not geocentric latitude. + Depth is measured from reference ellipsoid surface, not from center. + ``mesh.R`` provides spherical coordinates for backward compatibility. + ``mesh.geo`` provides geographic coordinates with ellipsoid geometry. + Right-handed coordinate system: WE x SN = down. """ from underworld3.coordinates import ELLIPSOIDS, geographic_to_cartesian from underworld3.units import ( diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py index 1ee4db5f7..732705d86 100644 --- a/src/underworld3/meshing/spherical.py +++ b/src/underworld3/meshing/spherical.py @@ -289,56 +289,42 @@ def SphericalShellInternalBoundary( Generates a spherical shell with an internal boundary using Gmsh. The function creates a 3D mesh of a spherical shell defined by outer, internal, and inner radii. Mesh size, polynomial degree, and Gmsh verbosity can be customized. - Parameters: - ----------- + Parameters + ---------- radiusOuter : float, optional The outer radius of the spherical shell. Default is 1.0. - radiusInternal : float, optional - The radius of the internal boundary within the spherical shell. Default is 0.8. - + The radius of the internal boundary. Default is 0.8. radiusInner : float, optional The inner radius of the spherical shell. Default is 0.547. - cellSize : float, optional - The target size for the mesh elements. This controls the density of the mesh. Default is 0.1. - + The target size for the mesh elements. Default is 0.1. degree : int, optional - The polynomial degree of the finite elements used in the mesh. Default is 1. - + The polynomial degree of the finite elements. Default is 1. qdegree : int, optional - The quadrature degree for integration. Higher values may improve accuracy but increase computation time. Default is 2. - + The quadrature degree for integration. Default is 2. filename : str, optional - The name of the file where the mesh will be saved. If None, a default name is generated based on the radii and mesh size. Default is None. - + The name of the file where the mesh will be saved. Default is None. refinement : optional - Refinement level or method for the mesh. Used to increase the resolution of the mesh in certain regions. Default is None. - + Refinement level for the mesh. Default is None. gmsh_verbosity : int, optional - Controls the verbosity of Gmsh output. Set to 0 for minimal output, higher numbers for more detailed logs. Default is 0. - + Gmsh output verbosity (0=quiet). Default is 0. verbose : bool, optional - If True, the function prints additional information during execution. Default is False. + If True, print additional information. Default is False. - Returns: - -------- - None - The function generates and saves a mesh file according to the specified parameters. + Returns + ------- + Mesh + The generated spherical shell mesh with internal boundary. - Example: + Examples -------- - mesh = uw.meshing.SphericalShellInternalBoundary( - radiusOuter=2.0, - radiusInternal=1.5, - radiusInner=1.0, - cellSize=0.05, - degree=2, - qdegree=3, - filename="custom_spherical_shell.msh", - gmsh_verbosity=1, - verbose=True - ) + >>> mesh = uw.meshing.SphericalShellInternalBoundary( + ... radiusOuter=2.0, + ... radiusInternal=1.5, + ... radiusInner=1.0, + ... cellSize=0.05, + ... ) """ class boundaries(Enum): @@ -511,64 +497,48 @@ def SegmentofSphere( Generates a segment of a sphere using Gmsh. This function creates a 3D mesh of a spherical segment defined by outer and inner radii, and the extent in longitude and latitude. The mesh can be customized in terms of size, polynomial degree, and verbosity. - Parameters: - ----------- + Parameters + ---------- radiusOuter : float, optional The outer radius of the spherical segment. Default is 1.0. - radiusInner : float, optional The inner radius of the spherical segment. Default is 0.547. - longitudeExtent : float, optional - The angular extent of the segment in the longitudinal direction (in degrees). Default is 90.0. - + The angular extent in longitude (degrees). Default is 90.0. latitudeExtent : float, optional - The angular extent of the segment in the latitudinal direction (in degrees). Default is 90.0. - + The angular extent in latitude (degrees). Default is 90.0. cellSize : float, optional - The target size for the mesh elements. This controls the density of the mesh. Default is 0.1. - + The target size for the mesh elements. Default is 0.1. degree : int, optional - The polynomial degree of the finite elements used in the mesh. Default is 1. - + The polynomial degree of the finite elements. Default is 1. qdegree : int, optional - The quadrature degree for integration. Higher values may improve accuracy but increase computation time. Default is 2. - + The quadrature degree for integration. Default is 2. filename : str, optional - The name of the file where the mesh will be saved. If None, a default name is generated based on the parameters. Default is None. - + The name of the file where the mesh will be saved. Default is None. refinement : optional - Refinement level or method for the mesh. Used to increase the resolution of the mesh in certain regions. Default is None. - + Refinement level for the mesh. Default is None. gmsh_verbosity : int, optional - Controls the verbosity of Gmsh output. Set to 0 for minimal output, higher numbers for more detailed logs. Default is 0. - + Gmsh output verbosity (0=quiet). Default is 0. verbose : bool, optional - If True, the function prints additional information during execution. Default is False. - - centroid : Tuple[float, float, float], optional - The coordinates of the centroid (center) of the sphere segment. Default is (0.0, 0.0, 0.0). + If True, print additional information. Default is False. + centroid : tuple of float, optional + The coordinates of the centroid of the sphere segment. + Default is (0.0, 0.0, 0.0). - Returns: - -------- - None - The function generates and saves a mesh file according to the specified parameters. + Returns + ------- + Mesh + The generated spherical segment mesh. - Example: + Examples -------- - mesh = uw.meshing.SegmentofSphere( - radiusOuter=2.0, - radiusInner=1.0, - longitudeExtent=120.0, - latitudeExtent=60.0, - cellSize=0.05, - degree=2, - qdegree=3, - filename="custom_sphere_segment.msh", - centroid=(0.0, 0.0, 0.0), - gmsh_verbosity=1, - verbose=True - ) + >>> mesh = uw.meshing.SegmentofSphere( + ... radiusOuter=2.0, + ... radiusInner=1.0, + ... longitudeExtent=120.0, + ... latitudeExtent=60.0, + ... cellSize=0.05, + ... ) """ class boundaries(Enum): diff --git a/src/underworld3/meshing/surfaces.py b/src/underworld3/meshing/surfaces.py index e055b47ae..d1724fc56 100644 --- a/src/underworld3/meshing/surfaces.py +++ b/src/underworld3/meshing/surfaces.py @@ -359,15 +359,19 @@ def __init__( ): """Create a surface. - Args: - name: Identifier for this surface - mesh: Computational mesh (required for .sym access and distance field) - control_points: (N, 3) array of 3D points defining the surface. - If None, the surface is empty and must be loaded or - have points set later. - symbol: Short LaTeX-friendly symbol for math display (e.g., "F" for "fault"). - If None, defaults to first letter of name capitalized. - Used in expressions like d_F instead of {surf_fault_distance}. + Parameters + ---------- + name : str + Identifier for this surface. + mesh : Mesh, optional + Computational mesh (required for ``.sym`` access and distance field). + control_points : ndarray, optional + (N, 3) array of 3D points defining the surface. + If None, the surface is empty and must be loaded or + have points set later. + symbol : str, optional + Short LaTeX-friendly symbol for math display (e.g., ``"F"``). + If None, defaults to first letter of name capitalized. """ self.name = name self.mesh = mesh @@ -468,6 +472,9 @@ def _dimensionalise_coords(self, coords: np.ndarray) -> np.ndarray: if coords is None: return None + if self.mesh is None: + return coords + cs = getattr(self.mesh, "CoordinateSystem", None) if cs is not None and getattr(cs, "_scaled", False): coords = coords * cs._length_scale @@ -886,26 +893,31 @@ def influence_function( Uses the absolute value of the signed distance field, so the influence is symmetric on both sides of the surface. For asymmetric behavior, - access the signed distance directly via surface.distance.sym[0]. + access the signed distance directly via ``surface.distance.sym[0]``. - Args: - width: Characteristic width of the transition zone. Can include - Pint units (e.g., 500 * u.meter), which will be converted - using the model's scaling system. - value_near: Value at/near the surface (can be a scalar or sympy expression) - value_far: Value far from the surface (can be a scalar or sympy expression) - profile: Transition profile type: - - "step": Sharp transition at distance = width - - "linear": Linear ramp from 0 to width - - "gaussian": Smooth Gaussian decay - - "smoothstep": C1-continuous Hermite interpolation + Parameters + ---------- + width : float + Characteristic width of the transition zone. + value_near : float or sympy.Expr, optional + Value at/near the surface. Default is 1.0. + value_far : float or sympy.Expr, optional + Value far from the surface. Default is 0.0. + profile : str, optional + Transition profile type. One of ``"step"`` (sharp transition), + ``"linear"`` (linear ramp), ``"gaussian"`` (smooth decay), + or ``"smoothstep"`` (C1-continuous Hermite). Default is ``"step"``. - Returns: - sympy.Expr that can be used in Underworld expressions + Returns + ------- + sympy.Expr + Expression that can be used in Underworld equations. - Example: - >>> # Step function for fault zone viscosity - >>> eta = surface.influence_function( + Examples + -------- + Step function for fault zone viscosity: + + >>> eta = surface.influence_function( ... width=0.05, ... value_near=0.01, ... value_far=1.0, diff --git a/src/underworld3/model.py b/src/underworld3/model.py index f9f320773..7f9795fb8 100644 --- a/src/underworld3/model.py +++ b/src/underworld3/model.py @@ -120,11 +120,11 @@ def __init__(self, name: Optional[str] = None, **kwargs): """ Initialize a new Model instance. - Parameters: - ----------- + Parameters + ---------- name : str, optional Human-readable name for this model instance - **kwargs : dict + \*\*kwargs : dict Additional arguments for Pydantic BaseModel """ # Handle name generation before calling super().__init__ @@ -534,13 +534,13 @@ def define_parameter(self, name: str, ptype=None, **kwargs): NOTE: Parameter system not yet implemented. Use model.materials dict directly. - Parameters: - ----------- + Parameters + ---------- name : str Parameter path (e.g., 'material.viscosity', 'solver.tolerance') ptype : ParameterType, optional Parameter type for validation (not used yet) - **kwargs : dict + \*\*kwargs : dict Additional arguments """ # TODO: Implement when parameter system is ready @@ -581,8 +581,8 @@ def set_reference_quantities(self, verbose=False, nondimensional_scaling=True, * By default, this automatically enables non-dimensionalization for solvers, ensuring consistent behavior between user-facing units and solver internals. - Parameters: - ----------- + Parameters + ---------- verbose : bool, optional If True, print diagnostic information about dimensional analysis and scale derivation. Default: False. @@ -591,60 +591,24 @@ def set_reference_quantities(self, verbose=False, nondimensional_scaling=True, * When True (recommended), solver operations work in non-dimensional [0-1] space while user-facing values remain in physical units. Set to False for expert mode (dimensional units only, no scaling). - **Warning**: Disabling this may cause numerical conditioning issues - and inconsistencies in unit conversions. - **quantities : dict - Named reference quantities using Pint units or UWQuantity objects, e.g.: - - mantle_viscosity=1e21*uw.units.Pa*uw.units.s - - plate_velocity=5*uw.units.cm/uw.units.year - - domain_depth=3000*uw.units.km - OR using uw.quantity(): - - domain_depth=uw.quantity(2900, "km") - - Example: - -------- - >>> # Standard usage (recommended for most users) - >>> model.set_reference_quantities( - ... mantle_viscosity=1e21*uw.units.Pa*uw.units.s, - ... plate_velocity=5*uw.units.cm/uw.units.year, - ... domain_depth=3000*uw.units.km - ... ) - # ✓ Units system active with automatic non-dimensionalization + Disabling this may cause numerical conditioning issues. + \*\*quantities : dict + Named reference quantities using Pint units or UWQuantity objects, + e.g. ``domain_depth=uw.quantity(2900, "km")``. - >>> # Also accepts UWQuantity objects - >>> model.set_reference_quantities( - ... domain_depth=uw.quantity(2900, "km"), - ... mantle_viscosity=uw.quantity(1e21, "Pa*s") - ... ) - - >>> # Expert mode (not recommended - dimensional units without scaling) - >>> model.set_reference_quantities( - ... domain_depth=uw.quantity(2900, "km"), - ... nondimensional_scaling=False - ... ) - # ⚠ Expert mode: Units active WITHOUT non-dimensionalization - - >>> # With diagnostic output - >>> model.set_reference_quantities( - ... verbose=True, - ... domain_depth=uw.quantity(500, "m"), - ... mantle_temperature=uw.quantity(1300, "K") - ... ) - - Notes: + Raises ------ + RuntimeError + If called after a mesh has been created (units are locked) + + Notes + ----- This method creates a Pint-native registry with model-specific constants using the _constants pattern for optimal numerical conditioning. - The default behavior (nondimensional_scaling=True) ensures: - - User-facing values in physical units (km, Myr, Pa*s, etc.) - - Solver operations in well-conditioned non-dimensional [0-1] space - - Automatic conversions handled transparently - - Raises: - ------- - RuntimeError - If called after a mesh has been created (units are locked) + The default behavior (nondimensional_scaling=True) ensures user-facing + values in physical units, solver operations in well-conditioned + non-dimensional [0-1] space, and automatic conversions. """ # Check if units are locked self._check_units_locked() diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 9fdfd86ba..72048f8f6 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -6,15 +6,12 @@ representations of material history, composition, and other quantities that move with the flow. -Key Components --------------- -SwarmType : enum - PETSc swarm type specification (BASIC or PIC). -SwarmVariable : class - Variable storing values at particle locations with mesh-based proxy - for use in symbolic expressions. -IndexSwarmVariable : class - Integer-valued swarm variable for material indexing. +**SwarmType** -- PETSc swarm type specification (BASIC or PIC). + +**SwarmVariable** -- Variable storing values at particle locations with +mesh-based proxy for use in symbolic expressions. + +**IndexSwarmVariable** -- Integer-valued swarm variable for material indexing. The swarm module integrates with PETSc's DMSwarm for parallel particle management and provides automatic population, advection, and repopulation @@ -3047,14 +3044,17 @@ def migration_control(self, disable=False): If False (default), migration is deferred until context exit. If True, migration is completely disabled. - Usage: - # Defer migration until end (default) + Examples + -------- + Defer migration until end (default):: + with swarm.migration_control(): swarm.points[mask1] += delta1 swarm.points[mask2] *= scale # Migration happens HERE on exit - # Completely disable migration + Completely disable migration:: + with swarm.migration_control(disable=True): # Operations where migration should never happen # No migration on exit diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 4861fb5ec..b0a9bec8b 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -13,28 +13,23 @@ - ``bdf(order)`` — backward differentiation formula (returns Δψ, divide by Δt for rate) - ``adams_moulton_flux(order)`` — weighted flux for implicit integration -Classes -------- -Symbolic - Pure symbolic history — no mesh storage. Used for flux tracking - in SNES_Diffusion where the flux expression is a SymPy tree, not - a mesh variable. -Eulerian - Fixed-grid ∂φ/∂t with optional grid-based advection u·∇φ. When - ``V_fn`` is provided, ``update_pre_solve`` applies an explicit - advection correction so that ``bdf()`` approximates the full - material derivative D/Dt = ∂/∂t + u·∇. -SemiLagrangian - Characteristic-based D/Dt via departure points. Traces backward - along velocity field to sample upstream values. Unconditionally - stable for advection (no CFL constraint) but less accurate when - velocity is near zero. -Lagrangian - Full particle-following D/Dt. Creates and manages its own swarm. - The swarm is advected during ``update_post_solve``. -Lagrangian_Swarm - Specialized swarm-based Lagrangian using a user-provided swarm. - The swarm advection is the user's responsibility. +**Symbolic** -- Pure symbolic history, no mesh storage. Used for flux +tracking in SNES_Diffusion where the flux expression is a SymPy tree. + +**Eulerian** -- Fixed-grid time derivative with optional grid-based +advection. When ``V_fn`` is provided, ``update_pre_solve`` applies an +explicit advection correction so that ``bdf()`` approximates the full +material derivative. + +**SemiLagrangian** -- Characteristic-based D/Dt via departure points. +Traces backward along velocity field to sample upstream values. +Unconditionally stable for advection (no CFL constraint). + +**Lagrangian** -- Full particle-following D/Dt. Creates and manages +its own swarm, advected during ``update_post_solve``. + +**Lagrangian_Swarm** -- Specialized swarm-based Lagrangian using a +user-provided swarm. Swarm advection is the user's responsibility. Notes ----- diff --git a/src/underworld3/units.py b/src/underworld3/units.py index a5e83ba9d..023e73dc4 100644 --- a/src/underworld3/units.py +++ b/src/underworld3/units.py @@ -937,40 +937,42 @@ def dimensionalise(expression, target_dimensionality=None, model=None) -> Any: 1. **Auto mode**: Extract dimensionality from the expression itself (if preserved) 2. **Explicit mode**: Use provided target_dimensionality - Args: - expression: Non-dimensional value (UWQuantity, UnitAwareArray, or plain number) - with preserved dimensionality metadata - target_dimensionality: Optional dict specifying target dimensionality - (Pint format: e.g., {'[length]': 1, '[time]': -1} for velocity) - If None, uses dimensionality from expression - model: Model instance with reference quantities (uses default if None) - - Returns: - Dimensional quantity with appropriate units - - Raises: - ValueError: If no dimensionality information available - ValueError: If model has no reference quantities - - Examples: - >>> # Auto mode - dimensionality preserved from non_dimensionalise() - >>> velocity_qty = uw.quantity(5.0, "cm/year") - >>> nondim_vel = non_dimensionalise(velocity_qty, model) - >>> # nondim_vel remembers it was velocity - >>> dimensional_vel = dimensionalise(nondim_vel, model=model) - >>> # Result has appropriate units based on model scales - - >>> # Explicit mode - specify dimensionality - >>> plain_value = 2.5 # dimensionless number - >>> velocity_dimensionality = {'[length]': 1, '[time]': -1} - >>> velocity = dimensionalise(plain_value, velocity_dimensionality, model) - >>> # Result is 2.5 * (length_scale / time_scale) - - >>> # With arrays - >>> nondim_array = UnitAwareArray([1.0, 2.0, 3.0], - ... units="dimensionless", - ... dimensionality={'[length]': 1}) - >>> dimensional_array = dimensionalise(nondim_array, model=model) + Parameters + ---------- + expression : UWQuantity, UnitAwareArray, or number + Non-dimensional value with preserved dimensionality metadata. + target_dimensionality : dict, optional + Target dimensionality in Pint format, e.g. + ``{'[length]': 1, '[time]': -1}`` for velocity. + If None, uses dimensionality from the expression. + model : Model, optional + Model instance with reference quantities. Uses default if None. + + Returns + ------- + quantity + Dimensional quantity with appropriate units. + + Raises + ------ + ValueError + If no dimensionality information is available. + ValueError + If model has no reference quantities. + + Examples + -------- + Auto mode -- dimensionality preserved from ``non_dimensionalise()``: + + >>> velocity_qty = uw.quantity(5.0, "cm/year") + >>> nondim_vel = non_dimensionalise(velocity_qty, model) + >>> dimensional_vel = dimensionalise(nondim_vel, model=model) + + Explicit mode -- specify dimensionality: + + >>> plain_value = 2.5 + >>> velocity_dimensionality = {'[length]': 1, '[time]': -1} + >>> velocity = dimensionalise(plain_value, velocity_dimensionality, model) """ # Get model if not provided if model is None: diff --git a/src/underworld3/utilities/nd_array_callback.py b/src/underworld3/utilities/nd_array_callback.py index 3d153ada2..fc5ba488d 100644 --- a/src/underworld3/utilities/nd_array_callback.py +++ b/src/underworld3/utilities/nd_array_callback.py @@ -104,32 +104,17 @@ def add_delayed_callback(self, array, callback_func, change_info): class NDArray_With_Callback(np.ndarray): - """ - # NDArray_With_Callback - - A numpy ndarray subclass that triggers **callbacks** when array data is modified. - This class maintains full numpy array compatibility while providing reactive programming - capabilities for scientific computing applications. - - ## Mathematical Representation - - Given an array $\\mathbf{A} \\in \\mathbb{R}^{n \\times m}$, any modification operation - $\\mathcal{O}(\\mathbf{A}) \\rightarrow \\mathbf{A}'$ will trigger registered callbacks: + """A numpy ndarray subclass that triggers callbacks when array data is modified. - $$\\mathbf{A}' = \\mathcal{O}(\\mathbf{A}) \\implies \\text{callback}(\\mathbf{A}', \\text{change\\_info})$$ + This class maintains full numpy array compatibility while providing reactive + programming capabilities for scientific computing applications. - Where $\\mathcal{O}$ represents operations like assignment, in-place arithmetic, or array methods. + **Callback Function Signature**:: - ## Callback Interface - - ### Function Signature - ```python - def callback(array: NDArray_With_Callback, change_info: dict) -> None: - # Handle array modification - pass - ``` + def callback(array: NDArray_With_Callback, change_info: dict) -> None: + pass - ### change_info Dictionary + The ``change_info`` dictionary contains: - ``operation`` (str): Operation name ('setitem', 'iadd', 'fill', etc.) - ``indices`` (tuple/slice/None): Location of change (for setitem operations) @@ -138,109 +123,16 @@ def callback(array: NDArray_With_Callback, change_info: dict) -> None: - ``array_shape`` (tuple): Current shape of the array - ``array_dtype`` (np.dtype): Data type of the array - ## Usage Examples - - ### Basic Reactive Array - ```python - def on_change(array, change_info): - print(f"🔔 {change_info['operation']} at {change_info['indices']}") - print(f" New value: {change_info['new_value']}") - - # Create reactive array - arr = NDArray_With_Callback([1, 2, 3]) - arr.set_callback(on_change) - - # Modifications trigger callbacks - arr[0] = 99 # 🔔 setitem at 0, New value: 99 - arr += 10 # 🔔 iadd at None, New value: 10 - arr.fill(0) # 🔔 fill at None, New value: 0 - ``` - - ### Scientific Computing Integration - ```python - class Mesh: - def __init__(self, coordinates): - self._coords = coordinates - - @property - def data(self): - arr = NDArray_With_Callback(self._coords, owner=self) - arr.set_callback(self._on_coordinates_changed) - return arr - - def _on_coordinates_changed(self, array, info): - # Invalidate cached computations - self._jacobians = None - self._mesh_quality = None - # Trigger dependent updates - self._update_connectivity() - ``` - - ### Backend Abstraction Pattern - - The callback mechanism enables **backend-agnostic data storage**. Different storage - backends (PETSc, pyvista, HDF5, etc.) can be wrapped with the same interface by - providing different callbacks: - - ```python - # Example 1: PETSc backend (for MeshVariable) - def petsc_sync_callback(array, info): - '''Sync numpy array to PETSc Vec.''' - self._petsc_vec.setArray(array) - self._petsc_vec.assemble() - - mesh_var.data = NDArray_With_Callback( - self._petsc_vec.getArray(), - callback=petsc_sync_callback, - owner=self, - ) - - # Example 2: pyvista backend (for SurfaceVariable) - def pyvista_sync_callback(array, info): - '''Sync numpy array to pyvista point_data.''' - self.surface._pv_mesh.point_data[self.name] = np.asarray(array) - self._proxy_stale = True # Mark derived data as stale - - surface_var.data = NDArray_With_Callback( - self._pv_mesh.point_data[name], - callback=pyvista_sync_callback, - owner=self, - ) - ``` - - This pattern allows `UnitAwareArray` (which extends this class) to work with - any storage backend without modification - only the callback changes. - - ### Delayed Callback Context - ```python - # Batch multiple operations - with arr.delay_callback("batch update"): - arr[0] = 1 - arr[1] = 2 - arr[2] = 3 - # All callbacks fire here (with MPI synchronization) - - # Global coordination across arrays - with NDArray_With_Callback.delay_callbacks_global("mesh deformation"): - mesh.data += displacement - swarm.data += velocity * dt - # Synchronized callback execution across all arrays - ``` - - ## Advanced Features - - - **Multiple callbacks**: `add_callback()`, `remove_callback()`, `clear_callbacks()` - - **Enable/disable**: `enable_callbacks()`, `disable_callbacks()` - - **Delayed execution**: `delay_callback()`, `delay_callbacks_global()` + **Features**: + + - **Multiple callbacks**: ``add_callback()``, ``remove_callback()``, ``clear_callbacks()`` + - **Enable/disable**: ``enable_callbacks()``, ``disable_callbacks()`` + - **Delayed execution**: ``delay_callback()``, ``delay_callbacks_global()`` - **MPI synchronization**: Automatic barriers in parallel contexts - **Weak references**: Owner tracking without circular dependencies - - **Error resilience**: Callback exceptions don't break array operations - - **Global reductions**: MPI-aware `global_max()`, `global_min()`, `global_sum()`, etc. + - **Global reductions**: MPI-aware ``global_max()``, ``global_min()``, ``global_sum()``, etc. - ## Global Reduction Operations (MPI-aware) - - These methods perform reduction operations across all MPI ranks, essential for - parallel scientific computing where data is distributed across processes. + **Global Reduction Operations (MPI-aware)**: - ``global_max(axis=None)``: Maximum value across all MPI ranks - ``global_min(axis=None)``: Minimum value across all MPI ranks @@ -250,37 +142,9 @@ def pyvista_sync_callback(array, info): - ``global_norm(ord=2)``: 2-norm (Euclidean) across all ranks - ``global_rms()``: Root mean square across all ranks - ### Usage Example - ```python - # In parallel code, each rank has a portion of mesh coordinates - coords = mesh.X.coords # NDArray_With_Callback or subclass - - # Find global bounds (across all MPI ranks) - x_min = coords[:, 0].global_min() # True minimum across all ranks - x_max = coords[:, 0].global_max() # True maximum across all ranks - - # Compute global statistics - mean_coord = coords.global_mean() # True mean (not just local mean!) - total_size = coords.global_size() # Total elements across all ranks - - # Compute global norms - rms_value = coords.global_rms() # Root mean square - l2_norm = coords.global_norm() # Euclidean norm - ``` - - ### Important Notes - - - These methods use MPI collective operations (`allreduce`) - - **All ranks must call these methods** (they are collective operations) - - Subclasses like `UnitAwareArray` override these to preserve units - - ## Performance Notes - - - **Zero overhead** when callbacks disabled - - **Minimal impact** on array operations (< 5% typical) - - **Batch processing** via delayed contexts for optimal performance - - **Thread-safe** delayed callback management - - **Memory efficient** weak reference ownership tracking + These methods use MPI collective operations (``allreduce``). + All ranks must call these methods (they are collective operations). + Subclasses like ``UnitAwareArray`` override these to preserve units. """ def __new__(cls, input_array=None, owner=None, callback=None, disable_inplace_operators=False): @@ -296,8 +160,9 @@ def __new__(cls, input_array=None, owner=None, callback=None, disable_inplace_op callback : callable, optional Initial callback function to register disable_inplace_operators : bool, optional - If True, in-place operators (+=, -=, *=, /=, etc.) will raise RuntimeError - for parallel safety. Default is False for backward compatibility. + If True, in-place operators (``+=``, ``-=``, ``*=``, ``/=``, etc.) + will raise RuntimeError for parallel safety. + Default is False for backward compatibility. """ if input_array is None: input_array = [] diff --git a/src/underworld3/utilities/unit_aware_array.py b/src/underworld3/utilities/unit_aware_array.py index ef62b340d..bef1d283d 100644 --- a/src/underworld3/utilities/unit_aware_array.py +++ b/src/underworld3/utilities/unit_aware_array.py @@ -31,51 +31,15 @@ class UnitAwareArray(NDArray_With_Callback): - """ - A numpy ndarray subclass that combines callback functionality with unit awareness. - - This class extends NDArray_With_Callback to provide: - - Automatic unit tracking and propagation - - Unit compatibility checking for operations - - Integration with UW3 unit conversion system - - Preservation of all callback functionality - - Unit-preserving global reduction operations (MPI-aware) - - Mathematical Representation: - Given an array A with units [A], operations preserve dimensional consistency: - - - A [m] + B [m] -> C [m] (compatible units) - - A [m] + B [s] -> Error (incompatible units) - - A [m] * B [s] -> C [m*s] (unit multiplication) - - A [m] * 2 -> C [m] (scalar multiplication) - - Usage Examples: - ```python - # Create arrays with units - length = UnitAwareArray([1, 2, 3], units="m") - time = UnitAwareArray([0.1, 0.2, 0.3], units="s") + """A numpy ndarray subclass that combines callback functionality with unit awareness. - # Operations preserve units - velocity = length / time # Result has units m/s + Extends ``NDArray_With_Callback`` to provide automatic unit tracking, + compatibility checking, and integration with the UW3 unit conversion system. - # Unit checking prevents errors - total = length + time # Raises ValueError (incompatible units) + Operations preserve dimensional consistency: compatible units are added, + incompatible units raise errors, and multiplication combines units. - # Automatic conversion when possible - length_km = UnitAwareArray([1, 2, 3], units="km") - total_length = length + length_km # Converts km to m automatically - - # Callbacks still work - def on_change(array, info): - print(f"Array {array.units} changed: {info['operation']}") - length.set_callback(on_change) - length[0] = 5 # Triggers callback - ``` - - Global Reduction Operations (MPI-aware, unit-preserving): - - These methods extend the parent class's global reduction operations to preserve - units in the result. Essential for parallel simulations with physical quantities. + **Global Reduction Operations (MPI-aware, unit-preserving)**: - ``global_max()`` -> UWQuantity (same units as array) - ``global_min()`` -> UWQuantity (same units as array) @@ -87,53 +51,13 @@ def on_change(array, info): - ``global_var()`` -> UWQuantity (units squared) - ``global_std()`` -> UWQuantity (same units as array) - ```python - # Example: mesh coordinates with units - mesh.units = "kilometer" - coords = mesh.X.coords # UnitAwareArray with units="kilometer" - - # Global reductions preserve units - x_max = coords[:, 0].global_max() # Returns: 100.0 [kilometer] - x_min = coords[:, 0].global_min() # Returns: 0.0 [kilometer] - mean_coord = coords.global_mean() # Returns: UWQuantity in kilometers - - # Variance has squared units - variance = coords.global_var() # Returns: UWQuantity in (kilometer)**2 - std_dev = coords.global_std() # Returns: UWQuantity in kilometers - ``` - - Note: Tensor arrays (ndim > 2) raise NotImplementedError for global reductions. + Tensor arrays (ndim > 2) raise NotImplementedError for global reductions. Use component-wise operations or slice the array for tensors. - Backend Abstraction: - - UnitAwareArray inherits the callback mechanism from NDArray_With_Callback, - enabling it to work with **any storage backend** - only the callback changes. - This provides a consistent unit-aware interface regardless of underlying storage: - - ```python - # MeshVariable: PETSc backend with units - def petsc_sync(array, info): - self._petsc_vec.setArray(array) - - mesh_var.data = UnitAwareArray( - petsc_data, units="m/s", - callback=petsc_sync, owner=self, - ) - - # SurfaceVariable: pyvista backend with units - def pyvista_sync(array, info): - self._pv_mesh.point_data[name] = np.asarray(array) - self._proxy_stale = True - - surface_var.data = UnitAwareArray( - pv_data, units="Pa", - callback=pyvista_sync, owner=self, - ) - ``` - - The same UnitAwareArray class provides consistent `.units`, `.magnitude`, and - unit-aware arithmetic for both backends without any code changes. + Inherits the callback mechanism from ``NDArray_With_Callback``, enabling it + to work with any storage backend (PETSc, pyvista, etc.) by changing only + the callback. Provides consistent ``.units``, ``.magnitude``, and unit-aware + arithmetic regardless of underlying storage. """ def __new__( diff --git a/tests/test_0601_mesh_vector_calc.py b/tests/test_0601_mesh_vector_calc.py index 0e512e605..ecf0e4fb4 100644 --- a/tests/test_0601_mesh_vector_calc.py +++ b/tests/test_0601_mesh_vector_calc.py @@ -146,8 +146,6 @@ def test_mesh_variables_with_units(): # assert mesh.vector.to_vector(v1.sym) == v1.ijk -# This is to be lit up after we properly deprecate .ijk -@pytest.mark.xfail(raises=AttributeError) def test_no_ijk(): # This object should not exist - Attribute error v13.ijk diff --git a/tests/test_0610_constitutive_tensor_regression.py b/tests/test_0610_constitutive_tensor_regression.py index 5cfe1612b..9a7e64e5b 100644 --- a/tests/test_0610_constitutive_tensor_regression.py +++ b/tests/test_0610_constitutive_tensor_regression.py @@ -280,7 +280,6 @@ def test_viscous_flow_model_parameter_types(self, mesh_2d): except Exception as e: pytest.fail(f"ViscousFlowModel failed with parameter type {type(param)}: {e}") - @pytest.mark.xfail(reason="ViscoElasticPlasticFlowModel has copy() bug - known issue") def test_viscoelastic_plastic_model_tensors(self, mesh_2d): """Test ViscoElasticPlasticFlowModel tensor operations.""" u = uw.discretisation.MeshVariable("U_vep", mesh_2d, mesh_2d.dim, degree=2) diff --git a/tests/test_0850_units_propagation.py b/tests/test_0850_units_propagation.py index ae8500eed..dc6817ebd 100644 --- a/tests/test_0850_units_propagation.py +++ b/tests/test_0850_units_propagation.py @@ -198,9 +198,6 @@ def test_addition_requires_matching_units(self): units_str = str(units).replace(" ", "") # Remove spaces for comparison assert "meter" in units_str - @pytest.mark.xfail( - reason="Rayleigh number may not simplify to exact dimensionless due to Pint representation" - ) def test_rayleigh_number_is_dimensionless(self): """Test that Rayleigh number calculation yields dimensionless number.""" uw.reset_default_model() diff --git a/tests/test_quantities_simplified.py b/tests/test_quantities_simplified.py index 68b0e59c9..71d4c5ffa 100644 --- a/tests/test_quantities_simplified.py +++ b/tests/test_quantities_simplified.py @@ -7,11 +7,22 @@ 3. Unit conversions 4. .data property (non-dimensional values) 5. SymPy compatibility + +NOTE: These tests import from `underworld3.function.quantities_simplified` +which was planned but never implemented. The existing quantities module is +`underworld3.function.quantities`. This entire file is skipped until the +simplified API is implemented or these tests are rewritten against the +existing quantities module. """ import pytest import numpy as np +pytestmark = pytest.mark.skip( + reason="Tests import quantities_simplified module which does not exist. " + "Rewrite against underworld3.function.quantities or implement the module." +) + class TestSimplifiedQuantityBasics: """Test basic UWQuantity creation and properties."""