From bbf7e6404c92aa22de77cf2972ae614d122b666a Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 15:16:44 +0200 Subject: [PATCH 01/22] docs: Complete Phase 7 - documentation and examples - Update README with comprehensive API documentation for v0.1.6: - Backend availability API (Scalar, SIMD, OpenMP, CUDA) - Boundary conditions with all BC types and functions - Derived fields and statistics - CPU features detection (AVX2, NEON) - Error handling with exception classes - Migration guide from v0.1.0 - Update MIGRATION_PLAN.md to mark Phase 7 complete - Add new examples: - backend_detection.py - CPU features, solver/BC backend queries - boundary_conditions.py - BC types, backends, inlet/outlet - channel_flow.py - Parabolic inlet, Poiseuille flow comparison - derived_fields.py - Statistics, velocity magnitude, VTK output - error_handling.py - Exception classes, raise_for_status - lid_driven_cavity_advanced.py - Complete cavity simulation - solver_comparison.py - Backend performance benchmarking - vtk_output.py - Scalar/vector fields, time series, CSV output --- MIGRATION_PLAN.md | 120 +++++++--- README.md | 250 +++++++++++++++++++- examples/backend_detection.py | 198 ++++++++++++++++ examples/boundary_conditions.py | 193 ++++++++++++++++ examples/channel_flow.py | 250 ++++++++++++++++++++ examples/derived_fields.py | 227 ++++++++++++++++++ examples/error_handling.py | 228 ++++++++++++++++++ examples/lid_driven_cavity_advanced.py | 231 +++++++++++++++++++ examples/solver_comparison.py | 308 +++++++++++++++++++++++++ examples/vtk_output.py | 268 +++++++++++++++++++++ 10 files changed, 2236 insertions(+), 37 deletions(-) create mode 100644 examples/backend_detection.py create mode 100644 examples/boundary_conditions.py create mode 100644 examples/channel_flow.py create mode 100644 examples/derived_fields.py create mode 100644 examples/error_handling.py create mode 100644 examples/lid_driven_cavity_advanced.py create mode 100644 examples/solver_comparison.py create mode 100644 examples/vtk_output.py diff --git a/MIGRATION_PLAN.md b/MIGRATION_PLAN.md index 83e88da..0d9fcf6 100644 --- a/MIGRATION_PLAN.md +++ b/MIGRATION_PLAN.md @@ -4,9 +4,9 @@ This document outlines the required changes to update cfd-python bindings to wor ## Current State -- **cfd-python version:** 0.1.0 (outdated) +- **cfd-python version:** 0.2.0 (compatible with CFD v0.1.6) - **Target CFD library:** v0.1.6 -- **Status:** BROKEN - will not compile against current CFD library +- **Status:** ✅ COMPLETE - All migration phases finished ## What's New in v0.1.6 @@ -457,37 +457,44 @@ cpu_features_t cfd_get_cpu_features(void); **Actual effort:** < 0.5 days -### Phase 7: Documentation & Tests (Required) +### Phase 7: Documentation & Tests (Required) ✅ COMPLETED **Priority:** P1 - Required for release -**Tasks:** - -- [ ] **7.1 Update README** - - Installation instructions - - API changes documentation - - Migration guide for users - - Backend availability examples +**Status:** Completed on 2026-01-03 -- [ ] **7.2 Update Python docstrings** - - All new functions - - BC examples - - Error handling examples - - Backend detection examples - -- [ ] **7.3 Add comprehensive tests** - - Test all BC types - - Test error handling - - Test derived fields - - Test with different backends - - Test backend availability API +**Tasks:** -- [ ] **7.4 Update examples** - - BC usage examples - - Derived fields examples - - Backend detection examples +- [x] **7.1 Update README** + - Complete API reference for all v0.1.6 features + - Boundary conditions documentation with examples + - Backend availability API documentation + - Error handling with exception classes + - CPU features detection + - Migration guide from v0.1.0 + +- [x] **7.2 Update Python docstrings** + - Comprehensive module docstring in `__init__.py` + - All BC functions, types, and edges documented + - Error handling API documented + - Backend availability API documented + - CPU features API documented + +- [x] **7.3 Add comprehensive tests** + - `test_boundary_conditions.py` - All BC types and backends + - `test_backend_availability.py` - Backend constants and functions + - `test_derived_fields.py` - Statistics and velocity magnitude + - `test_errors.py` - Exception classes and raise_for_status + - `test_cpu_features.py` - SIMD detection and grid stretching + - `test_abi_compatibility.py` - NULL handling and stress tests + +- [x] **7.4 Update examples** + - `boundary_conditions.py` - BC types, backends, inlet/outlet + - `backend_detection.py` - CPU features, solver/BC backends + - `derived_fields.py` - Statistics, velocity magnitude, VTK output + - `error_handling.py` - Exception classes, raise_for_status -**Estimated effort:** 2 days +**Actual effort:** < 0.5 days --- @@ -569,9 +576,9 @@ find_library(CFD_LIBRARY cfd_library) # Unified library name | Phase 4: Error Handling | ~~1 day~~ ✅ < 0.5 days | 4 days | | Phase 5: Backend Availability (v0.1.6) | ✅ 0.5 days | 4.5 days | | Phase 6: CPU Features | ~~1 day~~ ✅ < 0.5 days | 5 days | -| Phase 7: Docs & Tests | 2 days | 7 days | +| Phase 7: Docs & Tests | ~~2 days~~ ✅ < 0.5 days | 5.5 days | -**Total estimated effort:** ~~9-10 days~~ ~7 days (5 days completed) +**Total estimated effort:** ~~9-10 days~~ ~5.5 days ✅ ALL PHASES COMPLETE --- @@ -694,3 +701,58 @@ While `_exceptions.py` has type hints, the C extension functions lack type stubs 2. **Standardize coordinate naming** across grid functions 3. **Consider IntEnum** for constant groups (SIMD, BC types, backends) 4. **Optional high-level wrappers** for BC operations (Phase 8 candidate) + +--- + +## Phase 7 (CFD ROADMAP) Status: Python Integration + +This section tracks progress on Phase 7 from the main CFD library ROADMAP.md, which covers Python integration. + +### 7.1 Python Bindings (P1) + +| Task | Status | Notes | +|------|--------|-------| +| Complete C extension module | ✅ Done | `cfd_python.c` with Stable ABI (Py_LIMITED_API 3.9+) | +| NumPy array integration | ⚠️ Partial | Functions accept/return Python lists; NumPy conversion at Python layer | +| Pythonic API design | ✅ Done | Snake_case naming, native Python types, exception hierarchy | +| Type hints and stubs | ❌ Pending | Phase 8 in cfd-python ROADMAP - `.pyi` stubs not yet created | +| Pre-built wheels (manylinux, macOS, Windows) | ✅ Done | CI builds wheels for all platforms; CUDA variant for Linux/Windows | + +**C Extension Module Features:** +- Grid creation (uniform and stretched) +- Flow field operations +- Solver registry and creation +- Simulation API (`run_simulation`, `run_simulation_step`) +- Boundary conditions (all types: Neumann, Dirichlet, no-slip, inlet, outlet) +- Derived fields and statistics +- Error handling with Python exceptions +- Backend availability API (Scalar, SIMD, OMP, CUDA) +- CPU features detection (AVX2, NEON) + +### 7.2 High-level Python API (P1) + +| Task | Status | Notes | +|------|--------|-------| +| Problem definition classes | ❌ Pending | Not yet implemented | +| Mesh generation helpers | ⚠️ Partial | `create_grid`, `create_grid_stretched` available | +| Post-processing utilities | ⚠️ Partial | `compute_velocity_magnitude`, `compute_flow_statistics` available | +| Jupyter notebook integration | ❌ Pending | Not yet implemented | + +**What's Available:** +- Low-level C bindings exposed directly to Python +- Exception hierarchy for error handling (`CFDError` and subclasses) +- `raise_for_status()` helper for checking C library return codes + +**What's Missing:** +- High-level `Simulation` class for problem setup +- `BoundaryConditions` builder/context manager pattern +- Integration with matplotlib/visualization in notebooks +- Problem templates (lid-driven cavity, channel flow, etc.) + +### Summary + +Phase 7.1 (Python Bindings) is **largely complete** - the C extension module is functional, builds wheels for all platforms, and provides a working Python API. The main gaps are: +1. Type stubs for IDE support (Phase 8 in cfd-python) +2. NumPy direct integration (currently uses list conversion) + +Phase 7.2 (High-level Python API) is **not started** - the current API exposes C functions directly without high-level Pythonic wrappers or classes. diff --git a/README.md b/README.md index d18612b..6bce276 100644 --- a/README.md +++ b/README.md @@ -4,11 +4,13 @@ Python bindings for high-performance CFD simulation library using CPython C-API ## Features -- **High Performance**: Direct bindings to optimized C library +- **High Performance**: Direct bindings to optimized C library with SIMD (AVX2/NEON) and OpenMP support - **Stable ABI**: Compatible across Python 3.9+ versions -- **Static Linking**: Self-contained wheels with no external dependencies +- **Multiple Backends**: Scalar, SIMD, OpenMP, and CUDA (GPU) backends with runtime detection +- **Boundary Conditions**: Full BC support (Neumann, Dirichlet, no-slip, inlet, outlet) - **Dynamic Solver Discovery**: New solvers automatically available - **Multiple Output Formats**: VTK and CSV export support +- **Error Handling**: Rich exception hierarchy with detailed error messages ## Installation @@ -162,6 +164,20 @@ Create a computational grid. **Returns:** Dictionary with `nx`, `ny`, `xmin`, `xmax`, `ymin`, `ymax`, `x_coords`, `y_coords` +#### `create_grid_stretched(nx, ny, xmin, xmax, ymin, ymax, beta)` + +Create a grid with non-uniform (stretched) spacing. + +**Parameters:** + +- `nx`, `ny`: Grid dimensions +- `xmin`, `xmax`, `ymin`, `ymax`: Domain bounds +- `beta`: Stretching parameter (higher = more clustering) + +**Returns:** Dictionary with grid info including `x_coords`, `y_coords` + +> **Note:** The stretched grid implementation has a known bug - see ROADMAP.md for details. + #### `get_default_solver_params()` Get default solver parameters. @@ -199,6 +215,204 @@ cfd_python.SOLVER_PROJECTION # 'projection' # ... more solvers as registered in C library ``` +### Backend Availability (v0.1.6+) + +Query and select compute backends at runtime: + +```python +import cfd_python + +# Check available backends +print(cfd_python.get_available_backends()) +# ['Scalar', 'SIMD', 'OpenMP'] # CUDA if GPU available + +# Check specific backend +if cfd_python.backend_is_available(cfd_python.BACKEND_SIMD): + print("SIMD backend available!") + +# Get backend name +print(cfd_python.backend_get_name(cfd_python.BACKEND_OMP)) # 'OpenMP' + +# List solvers for a backend +omp_solvers = cfd_python.list_solvers_by_backend(cfd_python.BACKEND_OMP) +print(f"OpenMP solvers: {omp_solvers}") +``` + +**Backend Constants:** + +- `BACKEND_SCALAR`: Basic scalar CPU implementation +- `BACKEND_SIMD`: SIMD-optimized (AVX2/NEON) +- `BACKEND_OMP`: OpenMP parallelized +- `BACKEND_CUDA`: CUDA GPU acceleration + +### Boundary Conditions + +Apply various boundary conditions to flow fields: + +```python +import cfd_python + +# Create velocity field (as flat lists) +nx, ny = 50, 50 +u = [0.0] * (nx * ny) +v = [0.0] * (nx * ny) + +# Apply uniform inlet on left edge +cfd_python.bc_apply_inlet_uniform( + u, v, nx, ny, + u_inlet=1.0, v_inlet=0.0, + edge=cfd_python.BC_EDGE_LEFT +) + +# Apply zero-gradient outlet on right edge +cfd_python.bc_apply_outlet_velocity(u, v, nx, ny, cfd_python.BC_EDGE_RIGHT) + +# Apply no-slip walls on top and bottom +cfd_python.bc_apply_noslip(u, v, nx, ny) + +# Check/set BC backend +print(f"BC Backend: {cfd_python.bc_get_backend_name()}") +if cfd_python.bc_backend_available(cfd_python.BC_BACKEND_OMP): + cfd_python.bc_set_backend(cfd_python.BC_BACKEND_OMP) +``` + +**BC Type Constants:** + +- `BC_TYPE_PERIODIC`: Periodic boundaries +- `BC_TYPE_NEUMANN`: Zero-gradient boundaries +- `BC_TYPE_DIRICHLET`: Fixed value boundaries +- `BC_TYPE_NOSLIP`: No-slip wall (zero velocity) +- `BC_TYPE_INLET`: Inlet velocity specification +- `BC_TYPE_OUTLET`: Outlet conditions + +**BC Edge Constants:** + +- `BC_EDGE_LEFT`, `BC_EDGE_RIGHT`, `BC_EDGE_BOTTOM`, `BC_EDGE_TOP` + +**BC Backend Constants:** + +- `BC_BACKEND_AUTO`: Auto-select best available +- `BC_BACKEND_SCALAR`: Single-threaded scalar +- `BC_BACKEND_OMP`: OpenMP parallel +- `BC_BACKEND_SIMD`: SIMD + OpenMP (AVX2/NEON) +- `BC_BACKEND_CUDA`: GPU acceleration + +**BC Functions:** + +- `bc_apply_scalar(field, nx, ny, bc_type)`: Apply BC to scalar field +- `bc_apply_velocity(u, v, nx, ny, bc_type)`: Apply BC to velocity fields +- `bc_apply_dirichlet(field, nx, ny, left, right, bottom, top)`: Fixed values +- `bc_apply_noslip(u, v, nx, ny)`: Zero velocity at walls +- `bc_apply_inlet_uniform(u, v, nx, ny, u_inlet, v_inlet, edge)`: Uniform inlet +- `bc_apply_inlet_parabolic(u, v, nx, ny, max_velocity, edge)`: Parabolic inlet +- `bc_apply_outlet_scalar(field, nx, ny, edge)`: Zero-gradient outlet +- `bc_apply_outlet_velocity(u, v, nx, ny, edge)`: Zero-gradient outlet + +### Derived Fields & Statistics + +Compute derived quantities from flow fields: + +```python +import cfd_python + +# Compute velocity magnitude +u = [1.0] * 100 +v = [0.5] * 100 +vel_mag = cfd_python.compute_velocity_magnitude(u, v, 10, 10) + +# Calculate field statistics +stats = cfd_python.calculate_field_stats(vel_mag) +print(f"Min: {stats['min']}, Max: {stats['max']}, Avg: {stats['avg']}") + +# Comprehensive flow statistics +p = [0.0] * 100 # pressure field +flow_stats = cfd_python.compute_flow_statistics(u, v, p, 10, 10) +print(f"Max velocity: {flow_stats['velocity_magnitude']['max']}") +``` + +### CPU Features Detection + +Detect SIMD capabilities at runtime: + +```python +import cfd_python + +# Check SIMD architecture +arch = cfd_python.get_simd_arch() +name = cfd_python.get_simd_name() # 'avx2', 'neon', or 'none' + +# Check specific capabilities +if cfd_python.has_avx2(): + print("AVX2 available!") +elif cfd_python.has_neon(): + print("ARM NEON available!") + +# General SIMD check +if cfd_python.has_simd(): + print(f"SIMD enabled: {name}") +``` + +**SIMD Constants:** + +- `SIMD_NONE`: No SIMD support +- `SIMD_AVX2`: x86-64 AVX2 +- `SIMD_NEON`: ARM NEON + +### Error Handling + +Handle errors with Python exceptions: + +```python +import cfd_python +from cfd_python import ( + CFDError, + CFDMemoryError, + CFDInvalidError, + CFDDivergedError, + raise_for_status, +) + +# Check status codes +status = cfd_python.get_last_status() +if status != cfd_python.CFD_SUCCESS: + error_msg = cfd_python.get_last_error() + print(f"Error: {error_msg}") + +# Use raise_for_status helper +try: + raise_for_status(status, context="simulation step") +except CFDDivergedError as e: + print(f"Solver diverged: {e}") +except CFDInvalidError as e: + print(f"Invalid parameter: {e}") +except CFDError as e: + print(f"CFD error: {e}") + +# Clear error state +cfd_python.clear_error() +``` + +**Error Constants:** + +- `CFD_SUCCESS`: Operation successful (0) +- `CFD_ERROR`: Generic error (-1) +- `CFD_ERROR_NOMEM`: Out of memory (-2) +- `CFD_ERROR_INVALID`: Invalid argument (-3) +- `CFD_ERROR_IO`: File I/O error (-4) +- `CFD_ERROR_UNSUPPORTED`: Operation not supported (-5) +- `CFD_ERROR_DIVERGED`: Solver diverged (-6) +- `CFD_ERROR_MAX_ITER`: Max iterations reached (-7) + +**Exception Classes:** + +- `CFDError`: Base exception class +- `CFDMemoryError(CFDError, MemoryError)`: Memory allocation failed +- `CFDInvalidError(CFDError, ValueError)`: Invalid argument +- `CFDIOError(CFDError, IOError)`: File I/O error +- `CFDUnsupportedError(CFDError, NotImplementedError)`: Unsupported operation +- `CFDDivergedError(CFDError)`: Solver diverged +- `CFDMaxIterError(CFDError)`: Max iterations reached + ### Output Functions #### `set_output_dir(directory)` @@ -231,14 +445,34 @@ Write simulation timeseries data to CSV file. ### Output Type Constants ```python -cfd_python.OUTPUT_PRESSURE # Pressure field (VTK) -cfd_python.OUTPUT_VELOCITY # Velocity field (VTK) -cfd_python.OUTPUT_FULL_FIELD # Complete flow field (VTK) -cfd_python.OUTPUT_CSV_TIMESERIES # Time series (CSV) -cfd_python.OUTPUT_CSV_CENTERLINE # Centerline profile (CSV) -cfd_python.OUTPUT_CSV_STATISTICS # Global statistics (CSV) +cfd_python.OUTPUT_VELOCITY_MAGNITUDE # Velocity magnitude scalar (VTK) +cfd_python.OUTPUT_VELOCITY # Velocity vector field (VTK) +cfd_python.OUTPUT_FULL_FIELD # Complete flow field (VTK) +cfd_python.OUTPUT_CSV_TIMESERIES # Time series (CSV) +cfd_python.OUTPUT_CSV_CENTERLINE # Centerline profile (CSV) +cfd_python.OUTPUT_CSV_STATISTICS # Global statistics (CSV) ``` +## Migration from v0.1.0 + +If you're upgrading from an older version, note these changes: + +### API Changes + +1. **Type names changed** (internal, transparent to Python users) +2. **Error handling improved**: Use `get_last_error()`, `get_last_status()`, `clear_error()` +3. **New solver backends**: SIMD, OpenMP, CUDA available via `backend_is_available()` + +### New Features in v0.1.6 + +- **Boundary conditions**: Full BC API with inlet/outlet support +- **Backend selection**: Query and select compute backends at runtime +- **Derived fields**: `compute_velocity_magnitude()`, `compute_flow_statistics()` +- **Error handling**: Python exception classes with `raise_for_status()` +- **CPU detection**: `has_avx2()`, `has_neon()`, `has_simd()` + +See [MIGRATION_PLAN.md](MIGRATION_PLAN.md) for detailed migration information. + ## Requirements - Python 3.9+ diff --git a/examples/backend_detection.py b/examples/backend_detection.py new file mode 100644 index 0000000..294e459 --- /dev/null +++ b/examples/backend_detection.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python3 +""" +Backend Detection Example + +This example demonstrates how to query and use different compute backends +in cfd_python. It covers solver backends, BC backends, and CPU feature detection. +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import cfd_python +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to build the package first:") + print(" pip install -e .") + sys.exit(1) + + +def main(): + print("CFD Python Backend Detection Example") + print("=" * 60) + print(f"Version: {cfd_python.__version__}") + + # ================================================================= + # 1. CPU Feature Detection + # ================================================================= + print("\n1. CPU Feature Detection") + print("-" * 60) + + # Get SIMD architecture + arch = cfd_python.get_simd_arch() + arch_name = cfd_python.get_simd_name() + + print(f" SIMD Architecture: {arch_name} (id={arch})") + + # Check specific capabilities + has_avx2 = cfd_python.has_avx2() + has_neon = cfd_python.has_neon() + has_any_simd = cfd_python.has_simd() + + print("\n Feature Support:") + print(f" AVX2: {'Yes' if has_avx2 else 'No'}") + print(f" NEON: {'Yes' if has_neon else 'No'}") + print(f" Any SIMD: {'Yes' if has_any_simd else 'No'}") + + # SIMD constant mapping + print("\n SIMD Constants:") + print(f" SIMD_NONE = {cfd_python.SIMD_NONE}") + print(f" SIMD_AVX2 = {cfd_python.SIMD_AVX2}") + print(f" SIMD_NEON = {cfd_python.SIMD_NEON}") + + # ================================================================= + # 2. Solver Backend Availability + # ================================================================= + print("\n2. Solver Backend Availability") + print("-" * 60) + + # Get all available backends + available_backends = cfd_python.get_available_backends() + print(f" Available backends: {available_backends}") + + # Check each backend + backends = [ + ("SCALAR", cfd_python.BACKEND_SCALAR), + ("SIMD", cfd_python.BACKEND_SIMD), + ("OMP", cfd_python.BACKEND_OMP), + ("CUDA", cfd_python.BACKEND_CUDA), + ] + + print("\n Backend details:") + for name, backend_id in backends: + available = cfd_python.backend_is_available(backend_id) + backend_name = cfd_python.backend_get_name(backend_id) + status = "Available" if available else "Not available" + print(f" {name} ({backend_name}): {status}") + + # ================================================================= + # 3. Solvers by Backend + # ================================================================= + print("\n3. Solvers by Backend") + print("-" * 60) + + for name, backend_id in backends: + if cfd_python.backend_is_available(backend_id): + solvers = cfd_python.list_solvers_by_backend(backend_id) + backend_name = cfd_python.backend_get_name(backend_id) + if solvers: + print(f"\n {backend_name.upper()} backend solvers:") + for solver in solvers: + print(f" - {solver}") + else: + print(f"\n {backend_name.upper()}: No solvers registered") + + # ================================================================= + # 4. All Available Solvers + # ================================================================= + print("\n4. All Available Solvers") + print("-" * 60) + + all_solvers = cfd_python.list_solvers() + print(f" Total solvers: {len(all_solvers)}") + for solver in all_solvers: + info = cfd_python.get_solver_info(solver) + print(f"\n {solver}:") + print(f" Description: {info.get('description', 'N/A')}") + print(f" Version: {info.get('version', 'N/A')}") + + # ================================================================= + # 5. Boundary Condition Backends + # ================================================================= + print("\n5. Boundary Condition Backends") + print("-" * 60) + + bc_backends = [ + ("AUTO", cfd_python.BC_BACKEND_AUTO), + ("SCALAR", cfd_python.BC_BACKEND_SCALAR), + ("OMP", cfd_python.BC_BACKEND_OMP), + ("SIMD", cfd_python.BC_BACKEND_SIMD), + ("CUDA", cfd_python.BC_BACKEND_CUDA), + ] + + current_bc = cfd_python.bc_get_backend() + current_bc_name = cfd_python.bc_get_backend_name() + print(f" Current BC backend: {current_bc_name}") + + print("\n BC backend availability:") + for name, backend_id in bc_backends: + available = cfd_python.bc_backend_available(backend_id) + status = "Available" if available else "Not available" + marker = " (current)" if backend_id == current_bc else "" + print(f" BC_BACKEND_{name}: {status}{marker}") + + # ================================================================= + # 6. Selecting Optimal Backend + # ================================================================= + print("\n6. Selecting Optimal Backend") + print("-" * 60) + + # Determine best available backend + if cfd_python.backend_is_available(cfd_python.BACKEND_CUDA): + best = "CUDA" + elif cfd_python.backend_is_available(cfd_python.BACKEND_SIMD): + best = "SIMD" + elif cfd_python.backend_is_available(cfd_python.BACKEND_OMP): + best = "OpenMP" + else: + best = "Scalar" + + print(f" Recommended solver backend: {best}") + + # Set BC backend based on availability + if cfd_python.bc_backend_available(cfd_python.BC_BACKEND_SIMD): + cfd_python.bc_set_backend(cfd_python.BC_BACKEND_SIMD) + print(f" Set BC backend to: {cfd_python.bc_get_backend_name()}") + elif cfd_python.bc_backend_available(cfd_python.BC_BACKEND_OMP): + cfd_python.bc_set_backend(cfd_python.BC_BACKEND_OMP) + print(f" Set BC backend to: {cfd_python.bc_get_backend_name()}") + + # ================================================================= + # 7. Performance Recommendations + # ================================================================= + print("\n7. Performance Recommendations") + print("-" * 60) + + recommendations = [] + + if has_avx2: + recommendations.append("Use SIMD backend for vectorized operations (4x speedup)") + if cfd_python.backend_is_available(cfd_python.BACKEND_OMP): + recommendations.append("Use OpenMP backend for multi-core parallelization") + if cfd_python.backend_is_available(cfd_python.BACKEND_CUDA): + recommendations.append("Use CUDA backend for GPU acceleration (10-100x speedup)") + + if not recommendations: + recommendations.append("Only scalar backend available - consider building with SIMD/OpenMP") + + print(" Based on your system configuration:") + for rec in recommendations: + print(f" - {rec}") + + # ================================================================= + # Summary + # ================================================================= + print("\n" + "=" * 60) + print("Summary") + print("-" * 60) + print(f" CPU SIMD: {arch_name}") + print(f" Solver backends: {', '.join(available_backends)}") + print(f" BC backend: {cfd_python.bc_get_backend_name()}") + print(f" Total solvers: {len(all_solvers)}") + + +if __name__ == "__main__": + main() diff --git a/examples/boundary_conditions.py b/examples/boundary_conditions.py new file mode 100644 index 0000000..6a1666f --- /dev/null +++ b/examples/boundary_conditions.py @@ -0,0 +1,193 @@ +#!/usr/bin/env python3 +""" +Boundary Conditions Example + +This example demonstrates how to use the boundary condition API in cfd_python. +It shows the various BC types, backend selection, and practical usage patterns. +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import cfd_python +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to build the package first:") + print(" pip install -e .") + sys.exit(1) + + +def main(): + print("CFD Python Boundary Conditions Example") + print("=" * 60) + + # Grid dimensions + nx, ny = 20, 20 + size = nx * ny + + # ================================================================= + # 1. BC Backend Information + # ================================================================= + print("\n1. Boundary Condition Backend Information") + print("-" * 60) + + # Show current backend + current_backend = cfd_python.bc_get_backend() + current_name = cfd_python.bc_get_backend_name() + print(f" Current BC backend: {current_name} (id={current_backend})") + + # Check available backends + backends = [ + ("AUTO", cfd_python.BC_BACKEND_AUTO), + ("SCALAR", cfd_python.BC_BACKEND_SCALAR), + ("OMP", cfd_python.BC_BACKEND_OMP), + ("SIMD", cfd_python.BC_BACKEND_SIMD), + ("CUDA", cfd_python.BC_BACKEND_CUDA), + ] + + print("\n Backend availability:") + for name, backend_id in backends: + available = cfd_python.bc_backend_available(backend_id) + status = "Available" if available else "Not available" + print(f" BC_BACKEND_{name}: {status}") + + # ================================================================= + # 2. Neumann (Zero-Gradient) Boundary Conditions + # ================================================================= + print("\n2. Neumann (Zero-Gradient) Boundary Conditions") + print("-" * 60) + + # Create a scalar field with a gradient + pressure = [float(i % nx) for i in range(size)] # Gradient in x-direction + print(f" Initial field: min={min(pressure):.2f}, max={max(pressure):.2f}") + + # Apply Neumann BCs (zero gradient at boundaries) + cfd_python.bc_apply_scalar(pressure, nx, ny, cfd_python.BC_TYPE_NEUMANN) + print(" Applied Neumann BC to scalar field") + + # ================================================================= + # 3. Dirichlet (Fixed Value) Boundary Conditions + # ================================================================= + print("\n3. Dirichlet (Fixed Value) Boundary Conditions") + print("-" * 60) + + # Create a temperature-like field + temperature = [0.0] * size + + # Set fixed values at boundaries + left_val = 100.0 # Hot wall on left + right_val = 0.0 # Cold wall on right + bottom_val = 50.0 + top_val = 50.0 + + cfd_python.bc_apply_dirichlet(temperature, nx, ny, left_val, right_val, bottom_val, top_val) + print(f" Applied Dirichlet BC: left={left_val}, right={right_val}") + print(f" Left edge value: {temperature[0]}") + print(f" Right edge value: {temperature[nx-1]}") + + # ================================================================= + # 4. No-Slip Wall Boundary Conditions + # ================================================================= + print("\n4. No-Slip Wall Boundary Conditions") + print("-" * 60) + + # Create velocity fields with some initial values + u = [1.0] * size # Initial u velocity + v = [0.5] * size # Initial v velocity + + print(f" Before no-slip: u[0]={u[0]}, v[0]={v[0]}") + + # Apply no-slip (zero velocity at all walls) + cfd_python.bc_apply_noslip(u, v, nx, ny) + + print(" Applied no-slip BC to all walls") + print(f" After no-slip: u[0]={u[0]}, v[0]={v[0]} (corner)") + print(f" Interior: u[{nx+1}]={u[nx+1]}, v[{nx+1}]={v[nx+1]}") + + # ================================================================= + # 5. Inlet Boundary Conditions + # ================================================================= + print("\n5. Inlet Boundary Conditions") + print("-" * 60) + + # Reset velocity fields + u = [0.0] * size + v = [0.0] * size + + # Uniform inlet on left edge + u_inlet = 1.0 + v_inlet = 0.0 + cfd_python.bc_apply_inlet_uniform(u, v, nx, ny, u_inlet, v_inlet, cfd_python.BC_EDGE_LEFT) + print(f" Applied uniform inlet (u={u_inlet}, v={v_inlet}) on left edge") + print(f" Left edge u velocity: {u[0]}") + + # Parabolic inlet on left edge (for channel flow) + u2 = [0.0] * size + v2 = [0.0] * size + max_velocity = 1.5 + + cfd_python.bc_apply_inlet_parabolic(u2, v2, nx, ny, max_velocity, cfd_python.BC_EDGE_LEFT) + print(f"\n Applied parabolic inlet (max={max_velocity}) on left edge") + # Check velocity at center vs edge of inlet + center_idx = (ny // 2) * nx + edge_idx = 0 + print(f" Edge velocity: u={u2[edge_idx]:.4f}") + print(f" Center velocity: u={u2[center_idx]:.4f}") + + # ================================================================= + # 6. Outlet Boundary Conditions + # ================================================================= + print("\n6. Outlet Boundary Conditions") + print("-" * 60) + + # Use existing velocity field and apply outlet on right edge + cfd_python.bc_apply_outlet_velocity(u, v, nx, ny, cfd_python.BC_EDGE_RIGHT) + print(" Applied zero-gradient outlet on right edge") + + # Scalar field outlet + scalar_field = [float(i) for i in range(size)] + cfd_python.bc_apply_outlet_scalar(scalar_field, nx, ny, cfd_python.BC_EDGE_RIGHT) + print(" Applied zero-gradient outlet for scalar field") + + # ================================================================= + # 7. Switching Backends + # ================================================================= + print("\n7. Switching BC Backends") + print("-" * 60) + + # Try to set OpenMP backend if available + if cfd_python.bc_backend_available(cfd_python.BC_BACKEND_OMP): + cfd_python.bc_set_backend(cfd_python.BC_BACKEND_OMP) + print(f" Switched to: {cfd_python.bc_get_backend_name()}") + + # Apply BC with new backend + test_field = [0.0] * size + cfd_python.bc_apply_scalar(test_field, nx, ny, cfd_python.BC_TYPE_NEUMANN) + print(" Applied Neumann BC with OpenMP backend") + else: + print(" OpenMP backend not available, skipping") + + # Switch back to auto + cfd_python.bc_set_backend(cfd_python.BC_BACKEND_AUTO) + print(f" Reset to: {cfd_python.bc_get_backend_name()}") + + # ================================================================= + # Summary + # ================================================================= + print("\n" + "=" * 60) + print("Boundary Condition Types Available:") + print(" - BC_TYPE_PERIODIC: Periodic boundaries") + print(" - BC_TYPE_NEUMANN: Zero-gradient (∂φ/∂n = 0)") + print(" - BC_TYPE_DIRICHLET: Fixed value") + print(" - BC_TYPE_NOSLIP: Zero velocity at walls") + print(" - BC_TYPE_INLET: Inlet velocity (uniform/parabolic)") + print(" - BC_TYPE_OUTLET: Outlet conditions (zero-gradient)") + print("\nEdge Constants:") + print(" - BC_EDGE_LEFT, BC_EDGE_RIGHT, BC_EDGE_BOTTOM, BC_EDGE_TOP") + + +if __name__ == "__main__": + main() diff --git a/examples/channel_flow.py b/examples/channel_flow.py new file mode 100644 index 0000000..ef30f18 --- /dev/null +++ b/examples/channel_flow.py @@ -0,0 +1,250 @@ +#!/usr/bin/env python3 +""" +Channel Flow Example + +This example demonstrates setting up a channel flow simulation with: +- Parabolic inlet velocity profile +- Zero-gradient outlet +- No-slip walls on top and bottom +- Periodic or developed flow analysis +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import cfd_python + from cfd_python import ( + BC_EDGE_LEFT, + BC_EDGE_RIGHT, + ) +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to build the package first:") + print(" pip install -e .") + sys.exit(1) + + +def create_parabolic_profile(ny, u_max): + """Create a parabolic velocity profile for channel flow. + + For Poiseuille flow: u(y) = u_max * (1 - (2y/H - 1)^2) + where H is the channel height. + + Args: + ny: Number of grid points in y-direction + u_max: Maximum centerline velocity + + Returns: + list: Velocity values at each y-position + """ + profile = [] + for j in range(ny): + # Normalized y-coordinate: 0 at bottom, 1 at top + y_norm = j / (ny - 1) + # Parabolic profile: max at center (y=0.5), zero at walls + u = u_max * 4.0 * y_norm * (1.0 - y_norm) + profile.append(u) + return profile + + +def setup_channel_flow(nx, ny, u_max=1.0): + """Set up initial conditions for channel flow. + + Args: + nx, ny: Grid dimensions + u_max: Maximum inlet velocity + + Returns: + tuple: (u, v, p) velocity and pressure fields + """ + size = nx * ny + + # Initialize fields + u = [0.0] * size + v = [0.0] * size + p = [0.0] * size + + # Apply inlet BC (parabolic profile on left edge) + cfd_python.bc_apply_inlet_parabolic(u, v, nx, ny, u_max, BC_EDGE_LEFT) + + # Apply outlet BC (zero-gradient on right edge) + cfd_python.bc_apply_outlet_velocity(u, v, nx, ny, BC_EDGE_RIGHT) + + # Apply no-slip walls on top and bottom + # Note: bc_apply_noslip applies to all boundaries, so we apply walls specifically + for i in range(nx): + # Bottom wall (j=0) + idx_bottom = i + u[idx_bottom] = 0.0 + v[idx_bottom] = 0.0 + + # Top wall (j=ny-1) + idx_top = (ny - 1) * nx + i + u[idx_top] = 0.0 + v[idx_top] = 0.0 + + return u, v, p + + +def compute_analytical_solution(nx, ny, u_max, dp_dx, mu, H): + """Compute analytical Poiseuille flow solution. + + For fully developed channel flow: + u(y) = (1/2mu) * (-dp/dx) * y * (H - y) + + With u_max = (1/8mu) * (-dp/dx) * H^2 + + Args: + nx, ny: Grid dimensions + u_max: Maximum centerline velocity + dp_dx: Pressure gradient + mu: Dynamic viscosity + H: Channel height + + Returns: + list: Analytical velocity field + """ + u_analytical = [0.0] * (nx * ny) + dy = H / (ny - 1) + + for j in range(ny): + y = j * dy + # Parabolic profile + u_y = u_max * 4.0 * (y / H) * (1.0 - y / H) + + for i in range(nx): + idx = j * nx + i + u_analytical[idx] = u_y + + return u_analytical + + +def main(): + print("CFD Python - Channel Flow Example") + print("=" * 60) + + # Parameters + nx, ny = 64, 32 # Longer in x for channel flow + xmin, xmax = 0.0, 4.0 # Channel length = 4 + ymin, ymax = 0.0, 1.0 # Channel height = 1 + u_max = 1.0 # Maximum inlet velocity + + print("\nChannel Configuration:") + print(f" Grid: {nx} x {ny}") + print(f" Length: {xmax - xmin}") + print(f" Height: {ymax - ymin}") + print(f" Max inlet velocity: {u_max}") + + # Show backend info + print("\nBackend Information:") + print(f" SIMD: {cfd_python.get_simd_name()}") + print(f" BC Backend: {cfd_python.bc_get_backend_name()}") + + # Set up flow + print("\nSetting up channel flow...") + u, v, p = setup_channel_flow(nx, ny, u_max) + + # Verify inlet profile + print("\nInlet velocity profile (left edge):") + print(f" {'j':>4} {'y':>8} {'u':>10}") + print(" " + "-" * 26) + for j in range(0, ny, max(1, ny // 8)): + y = ymin + j * (ymax - ymin) / (ny - 1) + idx = j * nx # Left edge + print(f" {j:4d} {y:8.4f} {u[idx]:10.6f}") + + # Run simulation + print("\nRunning simulation...") + steps = 100 + + result = cfd_python.run_simulation_with_params( + nx=nx, ny=ny, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, steps=steps, dt=0.001, cfl=0.2 + ) + + print(f" Completed {steps} steps") + print(f" Solver: {result['solver_name']}") + + # Compute statistics + vel_mag = result["velocity_magnitude"] + flow_stats = cfd_python.compute_flow_statistics(u, v, p, nx, ny) + + print("\nFlow Statistics:") + print(" Velocity magnitude:") + print(f" Min: {flow_stats['velocity_magnitude']['min']:.6f}") + print(f" Max: {flow_stats['velocity_magnitude']['max']:.6f}") + print(f" Avg: {flow_stats['velocity_magnitude']['avg']:.6f}") + + # Check centerline velocity + centerline_velocities = [] + j_center = ny // 2 + for i in range(nx): + idx = j_center * nx + i + centerline_velocities.append(u[idx]) + + print(f"\nCenterline velocity (j={j_center}):") + print(f" Inlet: {centerline_velocities[0]:.6f}") + print(f" Center: {centerline_velocities[nx//2]:.6f}") + print(f" Outlet: {centerline_velocities[-1]:.6f}") + + # Compare with analytical solution + print("\nAnalytical Comparison:") + u_analytical = compute_analytical_solution(nx, ny, u_max, 0, 1, ymax - ymin) + u_analytical_stats = cfd_python.calculate_field_stats(u_analytical) + print(f" Analytical max velocity: {u_analytical_stats['max']:.6f}") + print(f" Simulated max velocity: {flow_stats['velocity_magnitude']['max']:.6f}") + + # Write output + print("\nWriting VTK output...") + + # Velocity magnitude + cfd_python.write_vtk_scalar( + "channel_velocity_magnitude.vtk", + "velocity_magnitude", + vel_mag, + nx, + ny, + xmin, + xmax, + ymin, + ymax, + ) + + # Velocity vectors + cfd_python.write_vtk_vector( + "channel_velocity.vtk", "velocity", u, v, nx, ny, xmin, xmax, ymin, ymax + ) + + print(" channel_velocity_magnitude.vtk") + print(" channel_velocity.vtk") + + # Cross-section profiles + print("\n" + "=" * 60) + print("Cross-section velocity profiles (u vs y):") + print("=" * 60) + + x_positions = [0, nx // 4, nx // 2, 3 * nx // 4, nx - 1] + x_labels = ["Inlet", "x=L/4", "x=L/2", "x=3L/4", "Outlet"] + + print(f"\n{'y':>8}", end="") + for label in x_labels: + print(f" {label:>10}", end="") + print() + print("-" * (8 + 12 * len(x_labels))) + + for j in range(0, ny, max(1, ny // 10)): + y = ymin + j * (ymax - ymin) / (ny - 1) + print(f"{y:8.4f}", end="") + for i in x_positions: + idx = j * nx + i + print(f" {u[idx]:10.6f}", end="") + print() + + print("\n" + "=" * 60) + print("Channel flow example complete!") + + +if __name__ == "__main__": + main() diff --git a/examples/derived_fields.py b/examples/derived_fields.py new file mode 100644 index 0000000..65f8968 --- /dev/null +++ b/examples/derived_fields.py @@ -0,0 +1,227 @@ +#!/usr/bin/env python3 +""" +Derived Fields and Statistics Example + +This example demonstrates how to compute derived quantities and statistics +from flow fields using cfd_python. +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import numpy as np + + import cfd_python +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to build the package first:") + print(" pip install -e .") + sys.exit(1) + + +def main(): + print("CFD Python Derived Fields Example") + print("=" * 60) + + # ================================================================= + # 1. Run a Simulation to Get Flow Fields + # ================================================================= + print("\n1. Running Simulation") + print("-" * 60) + + nx, ny = 32, 32 + result = cfd_python.run_simulation_with_params( + nx=nx, ny=ny, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=50, dt=0.001, cfl=0.2 + ) + + print(f" Grid: {nx} x {ny}") + print(f" Steps: {result['steps']}") + print(f" Solver: {result['solver_name']}") + + # ================================================================= + # 2. Compute Velocity Magnitude + # ================================================================= + print("\n2. Computing Velocity Magnitude") + print("-" * 60) + + # Create sample velocity fields (since run_simulation returns magnitude) + # For demonstration, create synthetic fields + u = [0.0] * (nx * ny) + v = [0.0] * (nx * ny) + + # Create a vortex-like velocity field + cx, cy = nx // 2, ny // 2 # Center + for j in range(ny): + for i in range(nx): + idx = j * nx + i + dx = i - cx + dy = j - cy + r = (dx * dx + dy * dy) ** 0.5 + 0.1 + # Tangential velocity (vortex) + u[idx] = -dy / r * 0.5 + v[idx] = dx / r * 0.5 + + # Compute velocity magnitude using cfd_python + vel_mag = cfd_python.compute_velocity_magnitude(u, v, nx, ny) + + print(f" Velocity magnitude computed for {len(vel_mag)} points") + print(f" Min: {min(vel_mag):.6f}") + print(f" Max: {max(vel_mag):.6f}") + + # Compare with numpy (if available) + u_np = np.array(u) + v_np = np.array(v) + vel_mag_np = np.sqrt(u_np**2 + v_np**2) + print(f" Numpy verification - Max diff: {np.max(np.abs(np.array(vel_mag) - vel_mag_np)):.2e}") + + # ================================================================= + # 3. Calculate Field Statistics + # ================================================================= + print("\n3. Calculating Field Statistics") + print("-" * 60) + + # Statistics for velocity magnitude + stats = cfd_python.calculate_field_stats(vel_mag) + print(" Velocity magnitude statistics:") + print(f" Min: {stats['min']:.6f}") + print(f" Max: {stats['max']:.6f}") + print(f" Avg: {stats['avg']:.6f}") + print(f" Sum: {stats['sum']:.6f}") + + # Statistics for u-velocity + u_stats = cfd_python.calculate_field_stats(u) + print("\n U-velocity statistics:") + print(f" Min: {u_stats['min']:.6f}") + print(f" Max: {u_stats['max']:.6f}") + print(f" Avg: {u_stats['avg']:.6f}") + + # Statistics for v-velocity + v_stats = cfd_python.calculate_field_stats(v) + print("\n V-velocity statistics:") + print(f" Min: {v_stats['min']:.6f}") + print(f" Max: {v_stats['max']:.6f}") + print(f" Avg: {v_stats['avg']:.6f}") + + # ================================================================= + # 4. Compute Comprehensive Flow Statistics + # ================================================================= + print("\n4. Comprehensive Flow Statistics") + print("-" * 60) + + # Create a pressure field + p = [0.0] * (nx * ny) + for j in range(ny): + for i in range(nx): + idx = j * nx + i + dx = i - cx + dy = j - cy + r2 = dx * dx + dy * dy + 0.1 + # Pressure decreases toward center (vortex core) + p[idx] = 1.0 - 0.5 / r2 + + # Get all flow statistics at once + flow_stats = cfd_python.compute_flow_statistics(u, v, p, nx, ny) + + print(" Full flow statistics:") + for field_name in ["u", "v", "p", "velocity_magnitude"]: + if field_name in flow_stats: + field_stats = flow_stats[field_name] + print(f"\n {field_name}:") + print(f" Min: {field_stats['min']:.6f}") + print(f" Max: {field_stats['max']:.6f}") + print(f" Avg: {field_stats['avg']:.6f}") + + # ================================================================= + # 5. Practical Application: Convergence Monitoring + # ================================================================= + print("\n5. Convergence Monitoring Example") + print("-" * 60) + + print(" Simulating convergence monitoring...") + + # Track statistics over multiple steps + max_velocities = [] + for step in range(5): + # Run a few steps + result = cfd_python.run_simulation_with_params( + nx=nx, ny=ny, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=10, dt=0.001 + ) + vel_stats = cfd_python.calculate_field_stats(result["velocity_magnitude"]) + max_velocities.append(vel_stats["max"]) + print(f" Step {(step+1)*10}: max_vel = {vel_stats['max']:.6f}") + + # Check for convergence (simplified) + if len(max_velocities) >= 2: + change = abs(max_velocities[-1] - max_velocities[-2]) + print(f"\n Velocity change: {change:.2e}") + + # ================================================================= + # 6. VTK Output with Derived Fields + # ================================================================= + print("\n6. VTK Output with Derived Fields") + print("-" * 60) + + grid = cfd_python.create_grid(nx, ny, 0.0, 1.0, 0.0, 1.0) + + # Write velocity magnitude + cfd_python.write_vtk_scalar( + "velocity_magnitude.vtk", + "velocity_magnitude", + vel_mag, + nx, + ny, + grid["xmin"], + grid["xmax"], + grid["ymin"], + grid["ymax"], + ) + print(" Written: velocity_magnitude.vtk") + + # Write pressure + cfd_python.write_vtk_scalar( + "pressure.vtk", + "pressure", + p, + nx, + ny, + grid["xmin"], + grid["xmax"], + grid["ymin"], + grid["ymax"], + ) + print(" Written: pressure.vtk") + + # Write velocity vectors + cfd_python.write_vtk_vector( + "velocity_field.vtk", + "velocity", + u, + v, + nx, + ny, + grid["xmin"], + grid["xmax"], + grid["ymin"], + grid["ymax"], + ) + print(" Written: velocity_field.vtk") + + # ================================================================= + # Summary + # ================================================================= + print("\n" + "=" * 60) + print("Available Derived Field Functions:") + print(" - compute_velocity_magnitude(u, v, nx, ny): sqrt(u^2 + v^2)") + print(" - calculate_field_stats(data): min, max, avg, sum") + print(" - compute_flow_statistics(u, v, p, nx, ny): all components") + print("\nThese functions are useful for:") + print(" - Post-processing simulation results") + print(" - Monitoring convergence") + print(" - Computing derived quantities for visualization") + + +if __name__ == "__main__": + main() diff --git a/examples/error_handling.py b/examples/error_handling.py new file mode 100644 index 0000000..081525e --- /dev/null +++ b/examples/error_handling.py @@ -0,0 +1,228 @@ +#!/usr/bin/env python3 +""" +Error Handling Example + +This example demonstrates the error handling API in cfd_python, +including exception classes and the raise_for_status helper. +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import cfd_python + from cfd_python import ( + CFDDivergedError, + CFDError, + CFDInvalidError, + CFDIOError, + CFDMaxIterError, + CFDMemoryError, + CFDUnsupportedError, + raise_for_status, + ) +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to build the package first:") + print(" pip install -e .") + sys.exit(1) + + +def main(): + print("CFD Python Error Handling Example") + print("=" * 60) + + # ================================================================= + # 1. Error Status Constants + # ================================================================= + print("\n1. Error Status Constants") + print("-" * 60) + + error_codes = [ + ("CFD_SUCCESS", cfd_python.CFD_SUCCESS), + ("CFD_ERROR", cfd_python.CFD_ERROR), + ("CFD_ERROR_NOMEM", cfd_python.CFD_ERROR_NOMEM), + ("CFD_ERROR_INVALID", cfd_python.CFD_ERROR_INVALID), + ("CFD_ERROR_IO", cfd_python.CFD_ERROR_IO), + ("CFD_ERROR_UNSUPPORTED", cfd_python.CFD_ERROR_UNSUPPORTED), + ("CFD_ERROR_DIVERGED", cfd_python.CFD_ERROR_DIVERGED), + ("CFD_ERROR_MAX_ITER", cfd_python.CFD_ERROR_MAX_ITER), + ] + + print(" Status codes defined:") + for name, value in error_codes: + desc = cfd_python.get_error_string(value) + print(f" {name} = {value}: {desc}") + + # ================================================================= + # 2. Checking Error State + # ================================================================= + print("\n2. Checking Error State") + print("-" * 60) + + # Clear any existing error + cfd_python.clear_error() + print(" Cleared error state") + + # Check current status + status = cfd_python.get_last_status() + error_msg = cfd_python.get_last_error() + + print(f" Current status: {status}") + print(f" Error message: {error_msg if error_msg else '(none)'}") + + # Demonstrate error checking pattern + if status == cfd_python.CFD_SUCCESS: + print(" Status is OK - no errors") + else: + print(f" Error occurred: {cfd_python.get_error_string(status)}") + + # ================================================================= + # 3. Exception Classes + # ================================================================= + print("\n3. Exception Classes") + print("-" * 60) + + exceptions = [ + ("CFDError", CFDError, Exception), + ("CFDMemoryError", CFDMemoryError, MemoryError), + ("CFDInvalidError", CFDInvalidError, ValueError), + ("CFDIOError", CFDIOError, IOError), + ("CFDUnsupportedError", CFDUnsupportedError, NotImplementedError), + ("CFDDivergedError", CFDDivergedError, CFDError), + ("CFDMaxIterError", CFDMaxIterError, CFDError), + ] + + print(" Exception hierarchy:") + for name, exc_class, base in exceptions: + bases = ", ".join(b.__name__ for b in exc_class.__mro__[1:-1]) + print(f" {name} <- {bases}") + + # ================================================================= + # 4. Using raise_for_status + # ================================================================= + print("\n4. Using raise_for_status") + print("-" * 60) + + # Test with success + try: + raise_for_status(cfd_python.CFD_SUCCESS, context="test operation") + print(" CFD_SUCCESS: No exception raised (as expected)") + except CFDError as e: + print(f" Unexpected error: {e}") + + # Test with each error code + test_codes = [ + (cfd_python.CFD_ERROR, "generic operation"), + (cfd_python.CFD_ERROR_NOMEM, "allocation"), + (cfd_python.CFD_ERROR_INVALID, "parameter validation"), + (cfd_python.CFD_ERROR_IO, "file writing"), + (cfd_python.CFD_ERROR_UNSUPPORTED, "feature check"), + (cfd_python.CFD_ERROR_DIVERGED, "solver step"), + (cfd_python.CFD_ERROR_MAX_ITER, "convergence"), + ] + + print("\n Testing raise_for_status with error codes:") + for code, context in test_codes: + try: + raise_for_status(code, context=context) + except CFDMemoryError as e: + print(f" {code}: CFDMemoryError - {e}") + except CFDInvalidError as e: + print(f" {code}: CFDInvalidError - {e}") + except CFDIOError as e: + print(f" {code}: CFDIOError - {e}") + except CFDUnsupportedError as e: + print(f" {code}: CFDUnsupportedError - {e}") + except CFDDivergedError as e: + print(f" {code}: CFDDivergedError - {e}") + except CFDMaxIterError as e: + print(f" {code}: CFDMaxIterError - {e}") + except CFDError as e: + print(f" {code}: CFDError - {e}") + + # ================================================================= + # 5. Practical Error Handling Pattern + # ================================================================= + print("\n5. Practical Error Handling Pattern") + print("-" * 60) + + def safe_simulation(nx, ny, steps): + """Run simulation with proper error handling.""" + try: + # Clear any previous errors + cfd_python.clear_error() + + # Run simulation + result = cfd_python.run_simulation_with_params( + nx=nx, ny=ny, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=steps, dt=0.001 + ) + + # Check for errors + status = cfd_python.get_last_status() + if status != cfd_python.CFD_SUCCESS: + raise_for_status(status, context="simulation") + + return result + + except CFDDivergedError: + print(" Solver diverged - try smaller time step") + return None + except CFDMaxIterError: + print(" Max iterations reached - may need more steps") + return None + except CFDInvalidError as e: + print(f" Invalid parameters: {e}") + return None + except CFDError as e: + print(f" CFD error: {e}") + return None + + print(" Running safe_simulation(32, 32, 10)...") + result = safe_simulation(32, 32, 10) + if result: + print(f" Success! Got {len(result['velocity_magnitude'])} values") + + # ================================================================= + # 6. Handling Specific Errors + # ================================================================= + print("\n6. Handling Specific Errors") + print("-" * 60) + + # Use Python's standard exception handling + print(" CFD exceptions inherit from Python built-ins:") + + # CFDMemoryError can be caught as MemoryError + try: + raise CFDMemoryError(-2, "Out of memory during grid allocation") + except MemoryError as e: + print(f" Caught as MemoryError: {e}") + + # CFDInvalidError can be caught as ValueError + try: + raise CFDInvalidError(-3, "Invalid grid dimensions") + except ValueError as e: + print(f" Caught as ValueError: {e}") + + # CFDIOError can be caught as IOError + try: + raise CFDIOError(-4, "Cannot write to output file") + except OSError as e: + print(f" Caught as IOError: {e}") + + # ================================================================= + # Summary + # ================================================================= + print("\n" + "=" * 60) + print("Error Handling Best Practices:") + print(" 1. Use cfd_python.clear_error() before critical operations") + print(" 2. Check cfd_python.get_last_status() after operations") + print(" 3. Use raise_for_status() for automatic exception raising") + print(" 4. Catch specific exceptions (CFDDivergedError, etc.)") + print(" 5. CFD exceptions inherit from Python built-ins for compatibility") + + +if __name__ == "__main__": + main() diff --git a/examples/lid_driven_cavity_advanced.py b/examples/lid_driven_cavity_advanced.py new file mode 100644 index 0000000..317e906 --- /dev/null +++ b/examples/lid_driven_cavity_advanced.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python3 +""" +Lid-Driven Cavity Flow - Advanced Example + +This example demonstrates a complete lid-driven cavity simulation with: +- Proper boundary condition setup +- Time stepping with convergence monitoring +- Post-processing and visualization output +- Statistics tracking +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import cfd_python + from cfd_python import ( + CFDError, + raise_for_status, + ) +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to build the package first:") + print(" pip install -e .") + sys.exit(1) + + +def setup_lid_driven_cavity(nx, ny, lid_velocity=1.0): + """Set up initial conditions for lid-driven cavity. + + Args: + nx, ny: Grid dimensions + lid_velocity: Velocity of the moving lid (top wall) + + Returns: + tuple: (u, v, p) velocity and pressure fields + """ + size = nx * ny + + # Initialize fields to zero + u = [0.0] * size + v = [0.0] * size + p = [0.0] * size + + # Apply boundary conditions + # No-slip on all walls first + cfd_python.bc_apply_noslip(u, v, nx, ny) + + # Moving lid on top (override top boundary) + # Set u = lid_velocity on top row + for i in range(nx): + idx = (ny - 1) * nx + i + u[idx] = lid_velocity + + return u, v, p + + +def run_cavity_simulation(nx=32, ny=32, Re=100, steps=500, output_interval=100): + """Run lid-driven cavity simulation. + + Args: + nx, ny: Grid dimensions + Re: Reynolds number + steps: Total time steps + output_interval: Steps between VTK outputs + + Returns: + dict: Simulation results + """ + print("\nLid-Driven Cavity Simulation") + print(f"Grid: {nx}x{ny}, Re={Re}, Steps={steps}") + print("=" * 60) + + # Grid setup + xmin, xmax = 0.0, 1.0 + ymin, ymax = 0.0, 1.0 + dx = (xmax - xmin) / (nx - 1) + dy = (ymax - ymin) / (ny - 1) + + # Time step (CFL-based) + lid_velocity = 1.0 + dt = 0.25 * min(dx, dy) / lid_velocity # CFL ~ 0.25 + dt = min(dt, 0.5 * Re * dx * dx) # Stability for diffusion + + print(f"Domain: [{xmin}, {xmax}] x [{ymin}, {ymax}]") + print(f"dx={dx:.4f}, dy={dy:.4f}, dt={dt:.6f}") + + # Set up output directory + output_dir = "cavity_output" + os.makedirs(output_dir, exist_ok=True) + cfd_python.set_output_dir(output_dir) + + # Initialize fields + u, v, p = setup_lid_driven_cavity(nx, ny, lid_velocity) + + # Track convergence + max_velocities = [] + residuals = [] + + print("\nRunning simulation...") + print(f"{'Step':>6} {'Max U':>12} {'Max V':>12} {'Residual':>12}") + print("-" * 48) + + prev_u = u.copy() + + for step in range(steps): + # Run simulation step using the library + try: + result = cfd_python.run_simulation_with_params( + nx=nx, ny=ny, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, steps=1, dt=dt, cfl=0.25 + ) + + # Check status + status = cfd_python.get_last_status() + if status != cfd_python.CFD_SUCCESS: + raise_for_status(status, context=f"step {step}") + + except CFDError as e: + print(f"Simulation error at step {step}: {e}") + break + + # Get velocity magnitude from result + vel_mag = result["velocity_magnitude"] + + # Compute statistics + stats = cfd_python.calculate_field_stats(vel_mag) + max_vel = stats["max"] + max_velocities.append(max_vel) + + # Compute residual (change in velocity) + if step > 0: + residual = sum(abs(vel_mag[i] - prev_u[i]) for i in range(len(vel_mag))) / len(vel_mag) + residuals.append(residual) + else: + residual = float("inf") + + prev_u = vel_mag.copy() + + # Output progress + if step % (steps // 10) == 0 or step == steps - 1: + u_stats = cfd_python.calculate_field_stats(vel_mag) + print(f"{step:6d} {u_stats['max']:12.6f} {u_stats['avg']:12.6f} {residual:12.2e}") + + # Write VTK output + if step % output_interval == 0 or step == steps - 1: + filename = f"cavity_step_{step:05d}.vtk" + cfd_python.write_vtk_scalar( + filename, "velocity_magnitude", vel_mag, nx, ny, xmin, xmax, ymin, ymax + ) + + print("\nSimulation complete!") + print(f"Output files written to: {output_dir}/") + + # Final statistics + final_stats = cfd_python.calculate_field_stats(vel_mag) + + return { + "nx": nx, + "ny": ny, + "Re": Re, + "steps": steps, + "final_max_velocity": final_stats["max"], + "final_avg_velocity": final_stats["avg"], + "max_velocities": max_velocities, + "residuals": residuals, + "output_dir": output_dir, + } + + +def analyze_results(results): + """Analyze and display simulation results.""" + print("\n" + "=" * 60) + print("Results Analysis") + print("=" * 60) + + print(f"\nGrid: {results['nx']} x {results['ny']}") + print(f"Reynolds number: {results['Re']}") + print(f"Time steps: {results['steps']}") + + print("\nFinal Statistics:") + print(f" Max velocity magnitude: {results['final_max_velocity']:.6f}") + print(f" Avg velocity magnitude: {results['final_avg_velocity']:.6f}") + + if results["residuals"]: + final_residual = results["residuals"][-1] + print(f" Final residual: {final_residual:.2e}") + + # Check convergence + if final_residual < 1e-6: + print(" Status: CONVERGED") + elif final_residual < 1e-4: + print(" Status: PARTIALLY CONVERGED") + else: + print(" Status: NOT CONVERGED (may need more steps)") + + print(f"\nOutput directory: {results['output_dir']}") + + +def main(): + print("CFD Python - Lid-Driven Cavity Example") + print("=" * 60) + + # Show system info + print("\nSystem Information:") + print(f" SIMD: {cfd_python.get_simd_name()}") + print(f" BC Backend: {cfd_python.bc_get_backend_name()}") + print(f" Available backends: {cfd_python.get_available_backends()}") + + # Run simulation + results = run_cavity_simulation(nx=32, ny=32, Re=100, steps=200, output_interval=50) + + # Analyze results + analyze_results(results) + + # Run parameter study (optional) + print("\n" + "=" * 60) + print("Parameter Study: Effect of Reynolds Number") + print("=" * 60) + + for Re in [50, 100, 200]: + result = cfd_python.run_simulation_with_params( + nx=32, ny=32, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=100, dt=0.0001 + ) + stats = cfd_python.calculate_field_stats(result["velocity_magnitude"]) + print(f" Re={Re:3d}: max_vel={stats['max']:.4f}, avg_vel={stats['avg']:.4f}") + + +if __name__ == "__main__": + main() diff --git a/examples/solver_comparison.py b/examples/solver_comparison.py new file mode 100644 index 0000000..d243d5a --- /dev/null +++ b/examples/solver_comparison.py @@ -0,0 +1,308 @@ +#!/usr/bin/env python3 +""" +Solver Comparison Example + +This example compares different solver backends available in cfd_python: +- Scalar (baseline, pure C implementation) +- SIMD (vectorized using AVX2/NEON) +- OpenMP (multi-threaded) +- CUDA (GPU-accelerated, if available) + +The comparison includes: +- Performance timing for each backend +- Result validation (all should produce identical results) +- Scaling analysis with problem size +""" + +import os +import sys +import time + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import cfd_python +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to build the package first:") + print(" pip install -e .") + sys.exit(1) + + +def benchmark_simulation(nx, ny, steps, label=""): + """Run a simulation and measure execution time. + + Args: + nx, ny: Grid dimensions + steps: Number of time steps + label: Description for output + + Returns: + tuple: (elapsed_time, result_dict) + """ + start = time.perf_counter() + + result = cfd_python.run_simulation_with_params( + nx=nx, ny=ny, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=steps, dt=0.001, cfl=0.2 + ) + + elapsed = time.perf_counter() - start + + return elapsed, result + + +def compute_checksum(field): + """Compute a simple checksum for validation.""" + return sum(field) + + +def main(): + print("CFD Python - Solver Comparison Example") + print("=" * 60) + + # ================================================================= + # 1. System Information + # ================================================================= + print("\n1. System Information") + print("-" * 60) + + # CPU features + print(f" SIMD Architecture: {cfd_python.get_simd_name()}") + print(f" Has AVX2: {cfd_python.has_avx2()}") + print(f" Has NEON: {cfd_python.has_neon()}") + print(f" Has Any SIMD: {cfd_python.has_simd()}") + + # Available backends + available = cfd_python.get_available_backends() + print(f"\n Available solver backends: {available}") + + # BC backend + print(f" BC backend: {cfd_python.bc_get_backend_name()}") + + # ================================================================= + # 2. Backend Details + # ================================================================= + print("\n2. Backend Details") + print("-" * 60) + + backends = [ + ("Scalar", cfd_python.BACKEND_SCALAR), + ("SIMD", cfd_python.BACKEND_SIMD), + ("OpenMP", cfd_python.BACKEND_OMP), + ("CUDA", cfd_python.BACKEND_CUDA), + ] + + print(f" {'Backend':<12} {'Available':<12} {'Solvers'}") + print(" " + "-" * 50) + + for name, backend_id in backends: + is_available = cfd_python.backend_is_available(backend_id) + if is_available: + solvers = cfd_python.list_solvers_by_backend(backend_id) + solver_count = len(solvers) if solvers else 0 + print(f" {name:<12} {'Yes':<12} {solver_count} registered") + else: + print(f" {name:<12} {'No':<12} -") + + # ================================================================= + # 3. Solver Inventory + # ================================================================= + print("\n3. Registered Solvers") + print("-" * 60) + + all_solvers = cfd_python.list_solvers() + print(f" Total registered solvers: {len(all_solvers)}") + + for solver_name in all_solvers: + info = cfd_python.get_solver_info(solver_name) + desc = info.get("description", "No description") + version = info.get("version", "N/A") + print(f"\n {solver_name}:") + print(f" Description: {desc}") + print(f" Version: {version}") + + # ================================================================= + # 4. Performance Comparison + # ================================================================= + print("\n4. Performance Comparison") + print("-" * 60) + + # Test configurations + test_configs = [ + (32, 32, 100, "Small grid"), + (64, 64, 100, "Medium grid"), + (128, 128, 50, "Large grid"), + ] + + print(f"\n {'Grid':<15} {'Steps':<8} {'Time (s)':<12} {'Cells/s':<15}") + print(" " + "-" * 55) + + results = [] + for nx, ny, steps, label in test_configs: + elapsed, result = benchmark_simulation(nx, ny, steps, label) + cells = nx * ny * steps + cells_per_sec = cells / elapsed if elapsed > 0 else 0 + + results.append( + { + "config": f"{nx}x{ny}", + "label": label, + "time": elapsed, + "cells_per_sec": cells_per_sec, + "checksum": compute_checksum(result["velocity_magnitude"]), + } + ) + + print(f" {nx}x{ny:<10} {steps:<8} {elapsed:<12.4f} {cells_per_sec:<15.0f}") + + # ================================================================= + # 5. BC Backend Comparison + # ================================================================= + print("\n5. Boundary Condition Backend Comparison") + print("-" * 60) + + # Test BC application performance + nx, ny = 128, 128 + size = nx * ny + iterations = 1000 + + bc_backends = [ + ("Scalar", cfd_python.BC_BACKEND_SCALAR), + ("SIMD", cfd_python.BC_BACKEND_SIMD), + ("OpenMP", cfd_python.BC_BACKEND_OMP), + ] + + print(f"\n Testing BC application ({iterations} iterations on {nx}x{ny} grid)") + print(f"\n {'Backend':<12} {'Available':<12} {'Time (ms)':<12} {'Speedup'}") + print(" " + "-" * 50) + + baseline_time = None + for name, backend_id in bc_backends: + if cfd_python.bc_backend_available(backend_id): + cfd_python.bc_set_backend(backend_id) + + # Create test fields + u = [1.0] * size + v = [0.5] * size + + # Time BC applications + start = time.perf_counter() + for _ in range(iterations): + cfd_python.bc_apply_noslip(u, v, nx, ny) + elapsed = time.perf_counter() - start + + elapsed_ms = elapsed * 1000 + + if baseline_time is None: + baseline_time = elapsed + speedup = 1.0 + else: + speedup = baseline_time / elapsed if elapsed > 0 else 0 + + print(f" {name:<12} {'Yes':<12} {elapsed_ms:<12.2f} {speedup:.2f}x") + else: + print(f" {name:<12} {'No':<12} -") + + # Reset to auto + cfd_python.bc_set_backend(cfd_python.BC_BACKEND_AUTO) + + # ================================================================= + # 6. Scaling Analysis + # ================================================================= + print("\n6. Scaling Analysis") + print("-" * 60) + + grid_sizes = [(16, 16), (32, 32), (64, 64), (128, 128)] + steps = 50 + + print(f"\n Grid size scaling ({steps} steps each)") + print(f"\n {'Grid':<12} {'Cells':<10} {'Time (s)':<10} {'Efficiency'}") + print(" " + "-" * 45) + + base_efficiency = None + for nx, ny in grid_sizes: + elapsed, _ = benchmark_simulation(nx, ny, steps) + cells = nx * ny + time_per_cell = elapsed / (cells * steps) if cells > 0 else 0 + + if base_efficiency is None: + base_efficiency = time_per_cell + efficiency = 100.0 + else: + efficiency = (base_efficiency / time_per_cell) * 100 if time_per_cell > 0 else 0 + + print(f" {nx}x{ny:<8} {cells:<10} {elapsed:<10.4f} {efficiency:.1f}%") + + # ================================================================= + # 7. Result Validation + # ================================================================= + print("\n7. Result Validation") + print("-" * 60) + + # Run same simulation with different grid sizes and check consistency + print("\n Verifying result consistency...") + + checksums = [] + for config in results: + checksums.append(config["checksum"]) + + print("\n Checksums from performance tests:") + for i, result in enumerate(results): + print(f" {result['config']}: {result['checksum']:.6f}") + + # ================================================================= + # 8. Recommendations + # ================================================================= + print("\n8. Performance Recommendations") + print("-" * 60) + + recommendations = [] + + # Check SIMD + if cfd_python.has_avx2(): + recommendations.append("AVX2 detected - SIMD backend recommended for CPU-intensive work") + elif cfd_python.has_neon(): + recommendations.append("NEON detected - SIMD backend recommended for CPU-intensive work") + + # Check OpenMP + if cfd_python.backend_is_available(cfd_python.BACKEND_OMP): + recommendations.append("OpenMP available - use for large grids to leverage all CPU cores") + + # Check CUDA + if cfd_python.backend_is_available(cfd_python.BACKEND_CUDA): + recommendations.append( + "CUDA available - use GPU backend for very large simulations (10-100x speedup)" + ) + else: + recommendations.append("CUDA not available - consider GPU build for maximum performance") + + # BC backend + if cfd_python.bc_backend_available(cfd_python.BC_BACKEND_SIMD): + recommendations.append("SIMD BC backend available - auto selection will use it") + + if recommendations: + print("\n Based on your system:") + for rec in recommendations: + print(f" - {rec}") + else: + recommendations.append("Only scalar backends available") + + # ================================================================= + # Summary + # ================================================================= + print("\n" + "=" * 60) + print("Summary") + print("-" * 60) + print(f" SIMD: {cfd_python.get_simd_name()}") + print(f" Solver backends: {', '.join(available)}") + print(f" BC backend: {cfd_python.bc_get_backend_name()}") + print(f" Total solvers: {len(all_solvers)}") + + # Best performance from tests + if results: + best = max(results, key=lambda x: x["cells_per_sec"]) + print(f" Best throughput: {best['cells_per_sec']:.0f} cells/s ({best['config']})") + + +if __name__ == "__main__": + main() diff --git a/examples/vtk_output.py b/examples/vtk_output.py new file mode 100644 index 0000000..dba1e29 --- /dev/null +++ b/examples/vtk_output.py @@ -0,0 +1,268 @@ +#!/usr/bin/env python3 +""" +VTK Output Example + +This example demonstrates the VTK output capabilities of cfd_python for +visualizing simulation results in ParaView or other VTK-compatible tools. + +VTK (Visualization Toolkit) format is widely used for scientific visualization. +This example shows how to: +- Write scalar fields (pressure, temperature, velocity magnitude) +- Write vector fields (velocity, vorticity) +- Create time series output for animations +- Configure output directories +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import cfd_python +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to build the package first:") + print(" pip install -e .") + sys.exit(1) + + +def create_sample_fields(nx, ny, xmin, xmax, ymin, ymax): + """Create sample fields for visualization. + + Creates synthetic data representing typical CFD quantities: + - Velocity field with a vortex pattern + - Pressure field with gradients + - Temperature field with boundary layer + + Args: + nx, ny: Grid dimensions + xmin, xmax, ymin, ymax: Domain bounds + + Returns: + dict: Dictionary of field arrays + """ + import math + + size = nx * ny + dx = (xmax - xmin) / (nx - 1) + dy = (ymax - ymin) / (ny - 1) + + # Initialize fields + u = [0.0] * size # x-velocity + v = [0.0] * size # y-velocity + p = [0.0] * size # pressure + T = [0.0] * size # temperature + + # Center of domain + xc = (xmax + xmin) / 2 + yc = (ymax + ymin) / 2 + + for j in range(ny): + y = ymin + j * dy + for i in range(nx): + x = xmin + i * dx + idx = j * nx + i + + # Create a vortex pattern for velocity + rx = x - xc + ry = y - yc + r = math.sqrt(rx * rx + ry * ry) + 1e-10 + + # Rankine vortex velocity profile + vortex_strength = 1.0 + core_radius = 0.3 + if r < core_radius: + factor = r / core_radius + else: + factor = core_radius / r + + u[idx] = -vortex_strength * factor * ry / r + v[idx] = vortex_strength * factor * rx / r + + # Pressure decreases toward vortex center + p[idx] = 1.0 - 0.5 * (vortex_strength * factor) ** 2 + + # Temperature with thermal boundary layer + # Hot on bottom, cold on top + T[idx] = 100.0 * (1.0 - y / ymax) + # Add some perturbation from velocity + T[idx] += 10.0 * math.exp(-(r**2) / (2 * core_radius**2)) + + return {"u": u, "v": v, "pressure": p, "temperature": T} + + +def main(): + print("CFD Python - VTK Output Example") + print("=" * 60) + + # ================================================================= + # 1. Setup + # ================================================================= + print("\n1. Setup") + print("-" * 60) + + # Grid parameters + nx, ny = 64, 64 + xmin, xmax = 0.0, 2.0 + ymin, ymax = 0.0, 2.0 + + print(f" Grid: {nx} x {ny}") + print(f" Domain: [{xmin}, {xmax}] x [{ymin}, {ymax}]") + + # Create output directory + output_dir = "vtk_output" + os.makedirs(output_dir, exist_ok=True) + cfd_python.set_output_dir(output_dir) + print(f" Output directory: {output_dir}/") + + # Create sample data + fields = create_sample_fields(nx, ny, xmin, xmax, ymin, ymax) + print(" Created sample vortex flow data") + + # ================================================================= + # 2. Writing Scalar Fields + # ================================================================= + print("\n2. Writing Scalar Fields") + print("-" * 60) + + # Write pressure field + cfd_python.write_vtk_scalar( + "pressure.vtk", "pressure", fields["pressure"], nx, ny, xmin, xmax, ymin, ymax + ) + print(" Written: pressure.vtk") + + # Write temperature field + cfd_python.write_vtk_scalar( + "temperature.vtk", "temperature", fields["temperature"], nx, ny, xmin, xmax, ymin, ymax + ) + print(" Written: temperature.vtk") + + # Compute and write velocity magnitude + vel_mag = cfd_python.compute_velocity_magnitude(fields["u"], fields["v"], nx, ny) + cfd_python.write_vtk_scalar( + "velocity_magnitude.vtk", "velocity_magnitude", vel_mag, nx, ny, xmin, xmax, ymin, ymax + ) + print(" Written: velocity_magnitude.vtk") + + # ================================================================= + # 3. Writing Vector Fields + # ================================================================= + print("\n3. Writing Vector Fields") + print("-" * 60) + + # Write velocity vector field + cfd_python.write_vtk_vector( + "velocity.vtk", "velocity", fields["u"], fields["v"], nx, ny, xmin, xmax, ymin, ymax + ) + print(" Written: velocity.vtk") + + # ================================================================= + # 4. Time Series Output + # ================================================================= + print("\n4. Time Series Output") + print("-" * 60) + print(" Creating animated vortex sequence...") + + import math + + # Create a series of snapshots showing vortex evolution + n_frames = 10 + for frame in range(n_frames): + # Rotate the vortex pattern + angle = frame * 2 * math.pi / n_frames + + # Create rotated velocity field + u_rot = [0.0] * (nx * ny) + v_rot = [0.0] * (nx * ny) + + cos_a = math.cos(angle) + sin_a = math.sin(angle) + + for idx in range(nx * ny): + u_orig = fields["u"][idx] + v_orig = fields["v"][idx] + u_rot[idx] = u_orig * cos_a - v_orig * sin_a + v_rot[idx] = u_orig * sin_a + v_orig * cos_a + + # Compute velocity magnitude for this frame + vel_mag_frame = cfd_python.compute_velocity_magnitude(u_rot, v_rot, nx, ny) + + # Write with frame number in filename + filename = f"vortex_{frame:04d}.vtk" + cfd_python.write_vtk_scalar( + filename, "velocity_magnitude", vel_mag_frame, nx, ny, xmin, xmax, ymin, ymax + ) + + print(f" Written: vortex_0000.vtk through vortex_{n_frames-1:04d}.vtk") + print(" (Load in ParaView as a time series)") + + # ================================================================= + # 5. CSV Output for Analysis + # ================================================================= + print("\n5. CSV Output for Analysis") + print("-" * 60) + + # Write velocity magnitude as CSV + cfd_python.write_csv_field("velocity_magnitude.csv", vel_mag, nx, ny) + print(" Written: velocity_magnitude.csv") + + # Write pressure as CSV + cfd_python.write_csv_field("pressure.csv", fields["pressure"], nx, ny) + print(" Written: pressure.csv") + + # ================================================================= + # 6. Field Statistics + # ================================================================= + print("\n6. Field Statistics") + print("-" * 60) + + # Compute statistics for each field + for name, field in [ + ("velocity_magnitude", vel_mag), + ("pressure", fields["pressure"]), + ("temperature", fields["temperature"]), + ]: + stats = cfd_python.calculate_field_stats(field) + print(f"\n {name}:") + print(f" Min: {stats['min']:.4f}") + print(f" Max: {stats['max']:.4f}") + print(f" Avg: {stats['avg']:.4f}") + + # Flow statistics (includes all velocity components) + print("\n Flow statistics:") + flow_stats = cfd_python.compute_flow_statistics( + fields["u"], fields["v"], fields["pressure"], nx, ny + ) + print(f" Velocity magnitude max: {flow_stats['velocity_magnitude']['max']:.4f}") + print(f" Pressure max: {flow_stats['pressure']['max']:.4f}") + + # ================================================================= + # Summary + # ================================================================= + print("\n" + "=" * 60) + print("VTK Output Summary") + print("-" * 60) + print(f" Output directory: {output_dir}/") + print("\n Files created:") + print(" Scalar fields:") + print(" - pressure.vtk") + print(" - temperature.vtk") + print(" - velocity_magnitude.vtk") + print(" Vector fields:") + print(" - velocity.vtk") + print(" Time series:") + print(f" - vortex_0000.vtk to vortex_{n_frames-1:04d}.vtk") + print(" CSV files:") + print(" - velocity_magnitude.csv") + print(" - pressure.csv") + print("\n To visualize:") + print(" 1. Open ParaView") + print(" 2. File -> Open -> Select VTK file") + print(" 3. Click 'Apply' in Properties panel") + print(" 4. For vectors: Filters -> Glyph") + print(" 5. For time series: Use playback controls") + + +if __name__ == "__main__": + main() From c3b6df927a1139d27334616f8821b318feabb709 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 15:58:55 +0200 Subject: [PATCH 02/22] refactor(examples): Clean up channel_flow.py - Remove unused create_parabolic_profile function - Simplify compute_analytical_solution by removing unused dp_dx and mu parameters - Update docstring to clarify its purpose as a reference implementation for validation --- examples/channel_flow.py | 40 ++++++++-------------------------------- 1 file changed, 8 insertions(+), 32 deletions(-) diff --git a/examples/channel_flow.py b/examples/channel_flow.py index ef30f18..fc0c5ff 100644 --- a/examples/channel_flow.py +++ b/examples/channel_flow.py @@ -27,29 +27,6 @@ sys.exit(1) -def create_parabolic_profile(ny, u_max): - """Create a parabolic velocity profile for channel flow. - - For Poiseuille flow: u(y) = u_max * (1 - (2y/H - 1)^2) - where H is the channel height. - - Args: - ny: Number of grid points in y-direction - u_max: Maximum centerline velocity - - Returns: - list: Velocity values at each y-position - """ - profile = [] - for j in range(ny): - # Normalized y-coordinate: 0 at bottom, 1 at top - y_norm = j / (ny - 1) - # Parabolic profile: max at center (y=0.5), zero at walls - u = u_max * 4.0 * y_norm * (1.0 - y_norm) - profile.append(u) - return profile - - def setup_channel_flow(nx, ny, u_max=1.0): """Set up initial conditions for channel flow. @@ -89,23 +66,22 @@ def setup_channel_flow(nx, ny, u_max=1.0): return u, v, p -def compute_analytical_solution(nx, ny, u_max, dp_dx, mu, H): - """Compute analytical Poiseuille flow solution. +def compute_analytical_solution(nx, ny, u_max, H): + """Compute analytical Poiseuille flow solution for validation. - For fully developed channel flow: - u(y) = (1/2mu) * (-dp/dx) * y * (H - y) + This is a reference implementation for comparing simulation results + against the known analytical solution for fully developed channel flow. - With u_max = (1/8mu) * (-dp/dx) * H^2 + For Poiseuille flow, the velocity profile is parabolic: + u(y) = u_max * 4 * (y/H) * (1 - y/H) Args: nx, ny: Grid dimensions u_max: Maximum centerline velocity - dp_dx: Pressure gradient - mu: Dynamic viscosity H: Channel height Returns: - list: Analytical velocity field + list: Analytical velocity field for comparison with simulation """ u_analytical = [0.0] * (nx * ny) dy = H / (ny - 1) @@ -191,7 +167,7 @@ def main(): # Compare with analytical solution print("\nAnalytical Comparison:") - u_analytical = compute_analytical_solution(nx, ny, u_max, 0, 1, ymax - ymin) + u_analytical = compute_analytical_solution(nx, ny, u_max, ymax - ymin) u_analytical_stats = cfd_python.calculate_field_stats(u_analytical) print(f" Analytical max velocity: {u_analytical_stats['max']:.6f}") print(f" Simulated max velocity: {flow_stats['velocity_magnitude']['max']:.6f}") From 3a3910086ce40a21886a5a6ec404755543a7d956 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 16:01:04 +0200 Subject: [PATCH 03/22] fix(examples): Prevent division by zero in lid_driven_cavity_advanced Use max(1, steps // 10) to avoid ZeroDivisionError when steps < 10. --- examples/lid_driven_cavity_advanced.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/lid_driven_cavity_advanced.py b/examples/lid_driven_cavity_advanced.py index 317e906..c69f252 100644 --- a/examples/lid_driven_cavity_advanced.py +++ b/examples/lid_driven_cavity_advanced.py @@ -139,7 +139,7 @@ def run_cavity_simulation(nx=32, ny=32, Re=100, steps=500, output_interval=100): prev_u = vel_mag.copy() # Output progress - if step % (steps // 10) == 0 or step == steps - 1: + if step % max(1, steps // 10) == 0 or step == steps - 1: u_stats = cfd_python.calculate_field_stats(vel_mag) print(f"{step:6d} {u_stats['max']:12.6f} {u_stats['avg']:12.6f} {residual:12.2e}") From dd7e1b4d7202121f810b7a3fecaca8961e735fef Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 16:03:42 +0200 Subject: [PATCH 04/22] style(examples): Use consistent import style for cfd_python Consolidate imports to use only 'import cfd_python' and access all symbols through the module namespace, avoiding mixed import styles. --- examples/channel_flow.py | 8 ++------ examples/lid_driven_cavity_advanced.py | 8 ++------ 2 files changed, 4 insertions(+), 12 deletions(-) diff --git a/examples/channel_flow.py b/examples/channel_flow.py index fc0c5ff..4a402a1 100644 --- a/examples/channel_flow.py +++ b/examples/channel_flow.py @@ -16,10 +16,6 @@ try: import cfd_python - from cfd_python import ( - BC_EDGE_LEFT, - BC_EDGE_RIGHT, - ) except ImportError as e: print(f"Import error: {e}") print("Make sure to build the package first:") @@ -45,10 +41,10 @@ def setup_channel_flow(nx, ny, u_max=1.0): p = [0.0] * size # Apply inlet BC (parabolic profile on left edge) - cfd_python.bc_apply_inlet_parabolic(u, v, nx, ny, u_max, BC_EDGE_LEFT) + cfd_python.bc_apply_inlet_parabolic(u, v, nx, ny, u_max, cfd_python.BC_EDGE_LEFT) # Apply outlet BC (zero-gradient on right edge) - cfd_python.bc_apply_outlet_velocity(u, v, nx, ny, BC_EDGE_RIGHT) + cfd_python.bc_apply_outlet_velocity(u, v, nx, ny, cfd_python.BC_EDGE_RIGHT) # Apply no-slip walls on top and bottom # Note: bc_apply_noslip applies to all boundaries, so we apply walls specifically diff --git a/examples/lid_driven_cavity_advanced.py b/examples/lid_driven_cavity_advanced.py index c69f252..3da17cf 100644 --- a/examples/lid_driven_cavity_advanced.py +++ b/examples/lid_driven_cavity_advanced.py @@ -16,10 +16,6 @@ try: import cfd_python - from cfd_python import ( - CFDError, - raise_for_status, - ) except ImportError as e: print(f"Import error: {e}") print("Make sure to build the package first:") @@ -115,9 +111,9 @@ def run_cavity_simulation(nx=32, ny=32, Re=100, steps=500, output_interval=100): # Check status status = cfd_python.get_last_status() if status != cfd_python.CFD_SUCCESS: - raise_for_status(status, context=f"step {step}") + cfd_python.raise_for_status(status, context=f"step {step}") - except CFDError as e: + except cfd_python.CFDError as e: print(f"Simulation error at step {step}: {e}") break From 19473b1de6d77d6b46c5a2a6a71c80274653d78c Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 16:56:09 +0200 Subject: [PATCH 05/22] docs: Update CHANGELOG with all migration phases and fix import style CHANGELOG updates: - Add Backend Availability API (Phase 5) - Add Derived Fields & Statistics (Phase 3) - Add Error Handling API (Phase 4) - Add CPU Features Detection (Phase 6) - Add CI/Build System (Phase 2.5) - Add Documentation & Examples (Phase 7) - Add new test suites error_handling.py: - Use consistent import style (only 'import cfd_python') - Access all symbols through module namespace --- CHANGELOG.md | 50 +++++++++++++++++++++++++++++++ examples/error_handling.py | 60 ++++++++++++++++---------------------- 2 files changed, 75 insertions(+), 35 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2d826aa..6468803 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,10 +8,60 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added + +#### Backend Availability API (Phase 5) + +- `BACKEND_SCALAR`, `BACKEND_SIMD`, `BACKEND_OMP`, `BACKEND_CUDA` constants +- `backend_is_available(backend)` - Check if backend is available at runtime +- `backend_get_name(backend)` - Get human-readable backend name +- `list_solvers_by_backend(backend)` - List solvers for specific backend +- `get_available_backends()` - List all available backend names + +#### Derived Fields & Statistics (Phase 3) + +- `calculate_field_stats(data)` - Compute min, max, avg, sum for a field +- `compute_velocity_magnitude(u, v, nx, ny)` - Compute velocity magnitude field +- `compute_flow_statistics(u, v, p, nx, ny)` - Comprehensive flow statistics + +#### Error Handling API (Phase 4) + +- Python exception hierarchy: `CFDError`, `CFDMemoryError`, `CFDInvalidError`, `CFDIOError`, `CFDUnsupportedError`, `CFDDivergedError`, `CFDMaxIterError` +- `raise_for_status(status_code, context)` - Raise appropriate exception for status codes +- Exception classes integrate with Python standard exceptions (e.g., `CFDMemoryError` inherits from `MemoryError`) + +#### CPU Features Detection (Phase 6) + +- `SIMD_NONE`, `SIMD_AVX2`, `SIMD_NEON` constants +- `get_simd_arch()` - Get SIMD architecture constant +- `get_simd_name()` - Get SIMD architecture name ("avx2", "neon", "none") +- `has_avx2()`, `has_neon()`, `has_simd()` - Check CPU SIMD capabilities +- `create_grid_stretched(nx, ny, xmin, xmax, ymin, ymax, beta)` - Stretched grid with hyperbolic cosine distribution + +#### CI/Build System (Phase 2.5) + - Dual-variant wheel builds supporting both CPU-only and CUDA-enabled configurations - Matrix build strategy in CI for separate CPU and CUDA wheel artifacts - Support for CFD library v0.1.6 modular backend libraries +#### Documentation & Examples (Phase 7) + +- Comprehensive README with full API reference +- New example scripts demonstrating cfd_python features: + - `lid_driven_cavity_advanced.py`: Complete cavity simulation with convergence monitoring + - `channel_flow.py`: Channel flow with parabolic inlet and analytical validation + - `vtk_output.py`: VTK output for ParaView visualization + - `solver_comparison.py`: Backend performance comparison and benchmarking + - `backend_detection.py`: CPU feature and backend availability detection + - `boundary_conditions.py`: Boundary condition API usage patterns + - `derived_fields.py`: Computing velocity magnitude and flow statistics + - `error_handling.py`: Error handling best practices +- New test suites: + - `test_backend_availability.py` - Backend constants and functions + - `test_derived_fields.py` - Statistics and velocity magnitude + - `test_errors.py` - Exception classes and raise_for_status + - `test_cpu_features.py` - SIMD detection and grid stretching + - `test_abi_compatibility.py` - NULL handling and stress tests + ### Changed - Updated build system to link modular CFD libraries (cfd_api, cfd_core, cfd_scalar, cfd_simd, cfd_omp, cfd_cuda) - Migrated to CUDA 12.0.0 from 12.6.2 for better stability and compatibility diff --git a/examples/error_handling.py b/examples/error_handling.py index 081525e..aa61f47 100644 --- a/examples/error_handling.py +++ b/examples/error_handling.py @@ -13,16 +13,6 @@ try: import cfd_python - from cfd_python import ( - CFDDivergedError, - CFDError, - CFDInvalidError, - CFDIOError, - CFDMaxIterError, - CFDMemoryError, - CFDUnsupportedError, - raise_for_status, - ) except ImportError as e: print(f"Import error: {e}") print("Make sure to build the package first:") @@ -86,13 +76,13 @@ def main(): print("-" * 60) exceptions = [ - ("CFDError", CFDError, Exception), - ("CFDMemoryError", CFDMemoryError, MemoryError), - ("CFDInvalidError", CFDInvalidError, ValueError), - ("CFDIOError", CFDIOError, IOError), - ("CFDUnsupportedError", CFDUnsupportedError, NotImplementedError), - ("CFDDivergedError", CFDDivergedError, CFDError), - ("CFDMaxIterError", CFDMaxIterError, CFDError), + ("CFDError", cfd_python.CFDError, Exception), + ("CFDMemoryError", cfd_python.CFDMemoryError, MemoryError), + ("CFDInvalidError", cfd_python.CFDInvalidError, ValueError), + ("CFDIOError", cfd_python.CFDIOError, IOError), + ("CFDUnsupportedError", cfd_python.CFDUnsupportedError, NotImplementedError), + ("CFDDivergedError", cfd_python.CFDDivergedError, cfd_python.CFDError), + ("CFDMaxIterError", cfd_python.CFDMaxIterError, cfd_python.CFDError), ] print(" Exception hierarchy:") @@ -108,9 +98,9 @@ def main(): # Test with success try: - raise_for_status(cfd_python.CFD_SUCCESS, context="test operation") + cfd_python.raise_for_status(cfd_python.CFD_SUCCESS, context="test operation") print(" CFD_SUCCESS: No exception raised (as expected)") - except CFDError as e: + except cfd_python.CFDError as e: print(f" Unexpected error: {e}") # Test with each error code @@ -127,20 +117,20 @@ def main(): print("\n Testing raise_for_status with error codes:") for code, context in test_codes: try: - raise_for_status(code, context=context) - except CFDMemoryError as e: + cfd_python.raise_for_status(code, context=context) + except cfd_python.CFDMemoryError as e: print(f" {code}: CFDMemoryError - {e}") - except CFDInvalidError as e: + except cfd_python.CFDInvalidError as e: print(f" {code}: CFDInvalidError - {e}") - except CFDIOError as e: + except cfd_python.CFDIOError as e: print(f" {code}: CFDIOError - {e}") - except CFDUnsupportedError as e: + except cfd_python.CFDUnsupportedError as e: print(f" {code}: CFDUnsupportedError - {e}") - except CFDDivergedError as e: + except cfd_python.CFDDivergedError as e: print(f" {code}: CFDDivergedError - {e}") - except CFDMaxIterError as e: + except cfd_python.CFDMaxIterError as e: print(f" {code}: CFDMaxIterError - {e}") - except CFDError as e: + except cfd_python.CFDError as e: print(f" {code}: CFDError - {e}") # ================================================================= @@ -163,20 +153,20 @@ def safe_simulation(nx, ny, steps): # Check for errors status = cfd_python.get_last_status() if status != cfd_python.CFD_SUCCESS: - raise_for_status(status, context="simulation") + cfd_python.raise_for_status(status, context="simulation") return result - except CFDDivergedError: + except cfd_python.CFDDivergedError: print(" Solver diverged - try smaller time step") return None - except CFDMaxIterError: + except cfd_python.CFDMaxIterError: print(" Max iterations reached - may need more steps") return None - except CFDInvalidError as e: + except cfd_python.CFDInvalidError as e: print(f" Invalid parameters: {e}") return None - except CFDError as e: + except cfd_python.CFDError as e: print(f" CFD error: {e}") return None @@ -196,19 +186,19 @@ def safe_simulation(nx, ny, steps): # CFDMemoryError can be caught as MemoryError try: - raise CFDMemoryError(-2, "Out of memory during grid allocation") + raise cfd_python.CFDMemoryError(-2, "Out of memory during grid allocation") except MemoryError as e: print(f" Caught as MemoryError: {e}") # CFDInvalidError can be caught as ValueError try: - raise CFDInvalidError(-3, "Invalid grid dimensions") + raise cfd_python.CFDInvalidError(-3, "Invalid grid dimensions") except ValueError as e: print(f" Caught as ValueError: {e}") # CFDIOError can be caught as IOError try: - raise CFDIOError(-4, "Cannot write to output file") + raise cfd_python.CFDIOError(-4, "Cannot write to output file") except OSError as e: print(f" Caught as IOError: {e}") From 7fb3ec4339fbe7b88331d84fb213f6cc606726dc Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 17:13:37 +0200 Subject: [PATCH 06/22] fix(examples): Correct comment/code mismatch in error_handling.py Update comment to reference OSError (the Python 3 name) instead of IOError, and update print statement to match. --- examples/error_handling.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/error_handling.py b/examples/error_handling.py index aa61f47..39bedd5 100644 --- a/examples/error_handling.py +++ b/examples/error_handling.py @@ -196,11 +196,11 @@ def safe_simulation(nx, ny, steps): except ValueError as e: print(f" Caught as ValueError: {e}") - # CFDIOError can be caught as IOError + # CFDIOError can be caught as OSError (IOError is an alias in Python 3) try: raise cfd_python.CFDIOError(-4, "Cannot write to output file") except OSError as e: - print(f" Caught as IOError: {e}") + print(f" Caught as OSError: {e}") # ================================================================= # Summary From bfde1a0d53c9a5d6a3db2e8608b4f4d250e4bba4 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 17:14:45 +0200 Subject: [PATCH 07/22] docs(examples): Clarify code behavior in channel_flow and error_handling channel_flow.py: - Add comment explaining that setup_channel_flow demonstrates BC API but run_simulation_with_params creates its own internal fields error_handling.py: - Add comment explaining __mro__[1:-1] slice behavior - Use _ for unused base variable in loop --- examples/channel_flow.py | 9 ++++++--- examples/error_handling.py | 5 +++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/examples/channel_flow.py b/examples/channel_flow.py index 4a402a1..e373b6d 100644 --- a/examples/channel_flow.py +++ b/examples/channel_flow.py @@ -115,11 +115,14 @@ def main(): print(f" SIMD: {cfd_python.get_simd_name()}") print(f" BC Backend: {cfd_python.bc_get_backend_name()}") - # Set up flow - print("\nSetting up channel flow...") + # Set up flow - demonstrates how to configure boundary conditions + # Note: setup_channel_flow shows BC API usage; run_simulation_with_params + # creates its own internal fields. The u, v, p here are used to verify + # the BC setup and for comparison with analytical solutions. + print("\nSetting up channel flow (BC demonstration)...") u, v, p = setup_channel_flow(nx, ny, u_max) - # Verify inlet profile + # Verify inlet profile (parabolic velocity distribution) print("\nInlet velocity profile (left edge):") print(f" {'j':>4} {'y':>8} {'u':>10}") print(" " + "-" * 26) diff --git a/examples/error_handling.py b/examples/error_handling.py index 39bedd5..9f16e5c 100644 --- a/examples/error_handling.py +++ b/examples/error_handling.py @@ -85,8 +85,9 @@ def main(): ("CFDMaxIterError", cfd_python.CFDMaxIterError, cfd_python.CFDError), ] - print(" Exception hierarchy:") - for name, exc_class, base in exceptions: + print(" Exception hierarchy (showing Method Resolution Order):") + for name, exc_class, _ in exceptions: + # __mro__[1:-1] gives all base classes except the class itself and 'object' bases = ", ".join(b.__name__ for b in exc_class.__mro__[1:-1]) print(f" {name} <- {bases}") From 1e68816c879ebda40581e9d7deebc9b745a5b1e5 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 17:23:45 +0200 Subject: [PATCH 08/22] refactor: Use explicit MRO filtering instead of slice Replace __mro__[1:-1] with explicit filtering using a generator expression that excludes the class itself and object. This is more readable and Pythonic. --- examples/error_handling.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/error_handling.py b/examples/error_handling.py index 9f16e5c..561cbff 100644 --- a/examples/error_handling.py +++ b/examples/error_handling.py @@ -87,8 +87,8 @@ def main(): print(" Exception hierarchy (showing Method Resolution Order):") for name, exc_class, _ in exceptions: - # __mro__[1:-1] gives all base classes except the class itself and 'object' - bases = ", ".join(b.__name__ for b in exc_class.__mro__[1:-1]) + # Filter MRO to show base classes, excluding the class itself and 'object' + bases = ", ".join(b.__name__ for b in exc_class.__mro__ if b not in (exc_class, object)) print(f" {name} <- {bases}") # ================================================================= From 51010824dcc1b0f831a86e56c156682530c57e8e Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 17:28:31 +0200 Subject: [PATCH 09/22] docs: Clarify Poiseuille flow formula in channel_flow example Add explanation of the factor of 4 in the parabolic velocity profile formula, clarifying that it ensures u_max occurs at the centerline and velocity is zero at the walls. --- examples/channel_flow.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/examples/channel_flow.py b/examples/channel_flow.py index e373b6d..2237ec0 100644 --- a/examples/channel_flow.py +++ b/examples/channel_flow.py @@ -71,9 +71,13 @@ def compute_analytical_solution(nx, ny, u_max, H): For Poiseuille flow, the velocity profile is parabolic: u(y) = u_max * 4 * (y/H) * (1 - y/H) + The factor of 4 ensures that u_max occurs at the centerline (y = H/2), + where the formula evaluates to u_max * 4 * 0.5 * 0.5 = u_max. + At the walls (y = 0 and y = H), velocity is zero (no-slip condition). + Args: nx, ny: Grid dimensions - u_max: Maximum centerline velocity + u_max: Maximum velocity at the channel centerline (y = H/2) H: Channel height Returns: From 5a61b7e3165205b582c9923e0dcf512885076d38 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 17:46:16 +0200 Subject: [PATCH 10/22] style: Use consistent import style in test_errors.py Remove mixed import styles (both 'import cfd_python' and 'from cfd_python import') and use only 'import cfd_python' with explicit module prefix for all symbols. --- tests/test_errors.py | 78 +++++++++++++++++++------------------------- 1 file changed, 34 insertions(+), 44 deletions(-) diff --git a/tests/test_errors.py b/tests/test_errors.py index e5e277a..4894d4f 100644 --- a/tests/test_errors.py +++ b/tests/test_errors.py @@ -5,16 +5,6 @@ import pytest import cfd_python -from cfd_python import ( - CFDDivergedError, - CFDError, - CFDInvalidError, - CFDIOError, - CFDMaxIterError, - CFDMemoryError, - CFDUnsupportedError, - raise_for_status, -) class TestErrorHandling: @@ -76,7 +66,7 @@ class TestExceptionClasses: def test_cfd_error_base(self): """Test CFDError base class""" - err = CFDError("test error", -1) + err = cfd_python.CFDError("test error", -1) assert err.message == "test error" assert err.status_code == -1 assert "test error" in str(err) @@ -84,47 +74,47 @@ def test_cfd_error_base(self): def test_cfd_error_default_status(self): """Test CFDError with default status code""" - err = CFDError("test error") + err = cfd_python.CFDError("test error") assert err.status_code == -1 def test_cfd_memory_error_inheritance(self): """Test CFDMemoryError inherits from both CFDError and MemoryError""" - err = CFDMemoryError("out of memory", -2) - assert isinstance(err, CFDError) + err = cfd_python.CFDMemoryError("out of memory", -2) + assert isinstance(err, cfd_python.CFDError) assert isinstance(err, MemoryError) assert err.status_code == -2 def test_cfd_invalid_error_inheritance(self): """Test CFDInvalidError inherits from both CFDError and ValueError""" - err = CFDInvalidError("invalid argument", -3) - assert isinstance(err, CFDError) + err = cfd_python.CFDInvalidError("invalid argument", -3) + assert isinstance(err, cfd_python.CFDError) assert isinstance(err, ValueError) assert err.status_code == -3 def test_cfd_io_error_inheritance(self): """Test CFDIOError inherits from both CFDError and IOError""" - err = CFDIOError("file not found", -4) - assert isinstance(err, CFDError) + err = cfd_python.CFDIOError("file not found", -4) + assert isinstance(err, cfd_python.CFDError) assert isinstance(err, IOError) assert err.status_code == -4 def test_cfd_unsupported_error_inheritance(self): """Test CFDUnsupportedError inherits from both CFDError and NotImplementedError""" - err = CFDUnsupportedError("not supported", -5) - assert isinstance(err, CFDError) + err = cfd_python.CFDUnsupportedError("not supported", -5) + assert isinstance(err, cfd_python.CFDError) assert isinstance(err, NotImplementedError) assert err.status_code == -5 def test_cfd_diverged_error(self): """Test CFDDivergedError""" - err = CFDDivergedError("solver diverged", -6) - assert isinstance(err, CFDError) + err = cfd_python.CFDDivergedError("solver diverged", -6) + assert isinstance(err, cfd_python.CFDError) assert err.status_code == -6 def test_cfd_max_iter_error(self): """Test CFDMaxIterError""" - err = CFDMaxIterError("max iterations reached", -7) - assert isinstance(err, CFDError) + err = cfd_python.CFDMaxIterError("max iterations reached", -7) + assert isinstance(err, cfd_python.CFDError) assert err.status_code == -7 @@ -133,61 +123,61 @@ class TestRaiseForStatus: def test_success_does_not_raise(self): """Test that success status codes do not raise""" - raise_for_status(0) # CFD_SUCCESS - raise_for_status(1) # Any positive value + cfd_python.raise_for_status(0) # CFD_SUCCESS + cfd_python.raise_for_status(1) # Any positive value def test_generic_error_raises_cfd_error(self): """Test that -1 raises CFDError""" - with pytest.raises(CFDError) as exc_info: - raise_for_status(-1) + with pytest.raises(cfd_python.CFDError) as exc_info: + cfd_python.raise_for_status(-1) assert exc_info.value.status_code == -1 def test_nomem_raises_cfd_memory_error(self): """Test that -2 raises CFDMemoryError""" - with pytest.raises(CFDMemoryError) as exc_info: - raise_for_status(-2) + with pytest.raises(cfd_python.CFDMemoryError) as exc_info: + cfd_python.raise_for_status(-2) assert exc_info.value.status_code == -2 def test_invalid_raises_cfd_invalid_error(self): """Test that -3 raises CFDInvalidError""" - with pytest.raises(CFDInvalidError) as exc_info: - raise_for_status(-3) + with pytest.raises(cfd_python.CFDInvalidError) as exc_info: + cfd_python.raise_for_status(-3) assert exc_info.value.status_code == -3 def test_io_raises_cfd_io_error(self): """Test that -4 raises CFDIOError""" - with pytest.raises(CFDIOError) as exc_info: - raise_for_status(-4) + with pytest.raises(cfd_python.CFDIOError) as exc_info: + cfd_python.raise_for_status(-4) assert exc_info.value.status_code == -4 def test_unsupported_raises_cfd_unsupported_error(self): """Test that -5 raises CFDUnsupportedError""" - with pytest.raises(CFDUnsupportedError) as exc_info: - raise_for_status(-5) + with pytest.raises(cfd_python.CFDUnsupportedError) as exc_info: + cfd_python.raise_for_status(-5) assert exc_info.value.status_code == -5 def test_diverged_raises_cfd_diverged_error(self): """Test that -6 raises CFDDivergedError""" - with pytest.raises(CFDDivergedError) as exc_info: - raise_for_status(-6) + with pytest.raises(cfd_python.CFDDivergedError) as exc_info: + cfd_python.raise_for_status(-6) assert exc_info.value.status_code == -6 def test_max_iter_raises_cfd_max_iter_error(self): """Test that -7 raises CFDMaxIterError""" - with pytest.raises(CFDMaxIterError) as exc_info: - raise_for_status(-7) + with pytest.raises(cfd_python.CFDMaxIterError) as exc_info: + cfd_python.raise_for_status(-7) assert exc_info.value.status_code == -7 def test_unknown_error_raises_cfd_error(self): """Test that unknown negative codes raise CFDError""" - with pytest.raises(CFDError) as exc_info: - raise_for_status(-99) + with pytest.raises(cfd_python.CFDError) as exc_info: + cfd_python.raise_for_status(-99) assert exc_info.value.status_code == -99 def test_context_included_in_message(self): """Test that context is included in error message""" - with pytest.raises(CFDError) as exc_info: - raise_for_status(-1, "during simulation") + with pytest.raises(cfd_python.CFDError) as exc_info: + cfd_python.raise_for_status(-1, "during simulation") assert "during simulation" in str(exc_info.value) From c3e0e25a919ac2ba327ff8dd6423123782e986a4 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 17:49:44 +0200 Subject: [PATCH 11/22] fix: Use simulation results for statistics in channel_flow example The example was computing flow_stats from the initial BC setup fields (u, v, p) instead of the actual simulation results. This was misleading because the "Flow Statistics" and "Simulated max velocity" values didn't reflect the simulation output. Changes: - Use calculate_field_stats(vel_mag) for simulation result statistics - Remove unused flow_stats computation - Update comment to clarify u, v are for BC demonstration only - Replace unused p with _ in tuple unpacking --- examples/channel_flow.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/channel_flow.py b/examples/channel_flow.py index 2237ec0..2b1212b 100644 --- a/examples/channel_flow.py +++ b/examples/channel_flow.py @@ -121,10 +121,10 @@ def main(): # Set up flow - demonstrates how to configure boundary conditions # Note: setup_channel_flow shows BC API usage; run_simulation_with_params - # creates its own internal fields. The u, v, p here are used to verify - # the BC setup and for comparison with analytical solutions. + # creates its own internal fields. The u, v here are used to verify + # the BC setup (inlet profile, cross-section profiles). print("\nSetting up channel flow (BC demonstration)...") - u, v, p = setup_channel_flow(nx, ny, u_max) + u, v, _ = setup_channel_flow(nx, ny, u_max) # Verify inlet profile (parabolic velocity distribution) print("\nInlet velocity profile (left edge):") @@ -146,15 +146,15 @@ def main(): print(f" Completed {steps} steps") print(f" Solver: {result['solver_name']}") - # Compute statistics + # Compute statistics from simulation result vel_mag = result["velocity_magnitude"] - flow_stats = cfd_python.compute_flow_statistics(u, v, p, nx, ny) + vel_mag_stats = cfd_python.calculate_field_stats(vel_mag) - print("\nFlow Statistics:") + print("\nSimulation Results:") print(" Velocity magnitude:") - print(f" Min: {flow_stats['velocity_magnitude']['min']:.6f}") - print(f" Max: {flow_stats['velocity_magnitude']['max']:.6f}") - print(f" Avg: {flow_stats['velocity_magnitude']['avg']:.6f}") + print(f" Min: {vel_mag_stats['min']:.6f}") + print(f" Max: {vel_mag_stats['max']:.6f}") + print(f" Avg: {vel_mag_stats['avg']:.6f}") # Check centerline velocity centerline_velocities = [] @@ -173,7 +173,7 @@ def main(): u_analytical = compute_analytical_solution(nx, ny, u_max, ymax - ymin) u_analytical_stats = cfd_python.calculate_field_stats(u_analytical) print(f" Analytical max velocity: {u_analytical_stats['max']:.6f}") - print(f" Simulated max velocity: {flow_stats['velocity_magnitude']['max']:.6f}") + print(f" Simulated max velocity: {vel_mag_stats['max']:.6f}") # Write output print("\nWriting VTK output...") From d2665eae00c9593fb5e54d2281df27a9c4b50a23 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 18:08:35 +0200 Subject: [PATCH 12/22] style: Use OSError consistently instead of IOError In Python 3, IOError is an alias for OSError. Updated all references to use OSError for consistency: - examples/error_handling.py: exceptions tuple - tests/test_errors.py: docstring and isinstance check --- examples/error_handling.py | 2 +- tests/test_errors.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/error_handling.py b/examples/error_handling.py index 561cbff..6fb93b3 100644 --- a/examples/error_handling.py +++ b/examples/error_handling.py @@ -79,7 +79,7 @@ def main(): ("CFDError", cfd_python.CFDError, Exception), ("CFDMemoryError", cfd_python.CFDMemoryError, MemoryError), ("CFDInvalidError", cfd_python.CFDInvalidError, ValueError), - ("CFDIOError", cfd_python.CFDIOError, IOError), + ("CFDIOError", cfd_python.CFDIOError, OSError), ("CFDUnsupportedError", cfd_python.CFDUnsupportedError, NotImplementedError), ("CFDDivergedError", cfd_python.CFDDivergedError, cfd_python.CFDError), ("CFDMaxIterError", cfd_python.CFDMaxIterError, cfd_python.CFDError), diff --git a/tests/test_errors.py b/tests/test_errors.py index 4894d4f..dcf53db 100644 --- a/tests/test_errors.py +++ b/tests/test_errors.py @@ -92,10 +92,10 @@ def test_cfd_invalid_error_inheritance(self): assert err.status_code == -3 def test_cfd_io_error_inheritance(self): - """Test CFDIOError inherits from both CFDError and IOError""" + """Test CFDIOError inherits from both CFDError and OSError""" err = cfd_python.CFDIOError("file not found", -4) assert isinstance(err, cfd_python.CFDError) - assert isinstance(err, IOError) + assert isinstance(err, OSError) assert err.status_code == -4 def test_cfd_unsupported_error_inheritance(self): From 607b5786edba0df11735b307a9c8f239fa166200 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 18:54:31 +0200 Subject: [PATCH 13/22] docs: Clarify synthetic fields in derived_fields.py example Make it explicit that the u, v fields created in section 2 are for demonstrating compute_velocity_magnitude() API usage, since run_simulation_with_params() only returns the velocity_magnitude. --- examples/derived_fields.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/derived_fields.py b/examples/derived_fields.py index 65f8968..02f239d 100644 --- a/examples/derived_fields.py +++ b/examples/derived_fields.py @@ -47,8 +47,9 @@ def main(): print("\n2. Computing Velocity Magnitude") print("-" * 60) - # Create sample velocity fields (since run_simulation returns magnitude) - # For demonstration, create synthetic fields + # Create synthetic velocity fields for demonstrating compute_velocity_magnitude(). + # Note: run_simulation_with_params() returns only velocity_magnitude, not the + # component fields. We create these synthetic fields to show the API usage. u = [0.0] * (nx * ny) v = [0.0] * (nx * ny) From 255dd5f2363b9fc1cab6d20a20a4a28d2ede11cf Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 19:05:22 +0200 Subject: [PATCH 14/22] docs: Fix misleading no-slip BC comment in channel_flow.py Clarify that bc_apply_noslip() is intentionally not used because it would overwrite the inlet/outlet conditions. The manual loop applies no-slip only to top and bottom walls. --- examples/channel_flow.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/channel_flow.py b/examples/channel_flow.py index 2b1212b..d0ae59e 100644 --- a/examples/channel_flow.py +++ b/examples/channel_flow.py @@ -46,8 +46,9 @@ def setup_channel_flow(nx, ny, u_max=1.0): # Apply outlet BC (zero-gradient on right edge) cfd_python.bc_apply_outlet_velocity(u, v, nx, ny, cfd_python.BC_EDGE_RIGHT) - # Apply no-slip walls on top and bottom - # Note: bc_apply_noslip applies to all boundaries, so we apply walls specifically + # Apply no-slip walls on top and bottom manually. + # We don't use bc_apply_noslip() (no-slip) here because it applies to all boundaries, + # which would overwrite the inlet/outlet conditions we just set. for i in range(nx): # Bottom wall (j=0) idx_bottom = i From a4e585863829bac14c87f7bee0c06afae2a833d7 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 19:38:21 +0200 Subject: [PATCH 15/22] fix: Correct exception constructor parameter order in error_handling.py The CFDError constructor signature is (message, status_code), not (status_code, message). Fixed the example to match the implementation and test file. --- examples/error_handling.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/error_handling.py b/examples/error_handling.py index 6fb93b3..20feb1f 100644 --- a/examples/error_handling.py +++ b/examples/error_handling.py @@ -187,19 +187,19 @@ def safe_simulation(nx, ny, steps): # CFDMemoryError can be caught as MemoryError try: - raise cfd_python.CFDMemoryError(-2, "Out of memory during grid allocation") + raise cfd_python.CFDMemoryError("Out of memory during grid allocation", -2) except MemoryError as e: print(f" Caught as MemoryError: {e}") # CFDInvalidError can be caught as ValueError try: - raise cfd_python.CFDInvalidError(-3, "Invalid grid dimensions") + raise cfd_python.CFDInvalidError("Invalid grid dimensions", -3) except ValueError as e: print(f" Caught as ValueError: {e}") # CFDIOError can be caught as OSError (IOError is an alias in Python 3) try: - raise cfd_python.CFDIOError(-4, "Cannot write to output file") + raise cfd_python.CFDIOError("Cannot write to output file", -4) except OSError as e: print(f" Caught as OSError: {e}") From 732e775b5e1f7529bd02fe59d125239d30e6ebf5 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 19:50:10 +0200 Subject: [PATCH 16/22] fix: Fix example scripts for cross-platform compatibility - boundary_conditions.py: Replace Unicode partial derivative symbol with ASCII 'dphi/dn' for Windows console compatibility - output_formats.py: Handle case where CSV file is written to CWD instead of temp directory - vtk_output.py: Replace missing write_csv_field() with inline Python, fix flow_stats key from 'pressure' to 'p' --- examples/boundary_conditions.py | 2 +- examples/output_formats.py | 5 ++++- examples/vtk_output.py | 16 ++++++++++++---- 3 files changed, 17 insertions(+), 6 deletions(-) diff --git a/examples/boundary_conditions.py b/examples/boundary_conditions.py index 6a1666f..8f0e633 100644 --- a/examples/boundary_conditions.py +++ b/examples/boundary_conditions.py @@ -180,7 +180,7 @@ def main(): print("\n" + "=" * 60) print("Boundary Condition Types Available:") print(" - BC_TYPE_PERIODIC: Periodic boundaries") - print(" - BC_TYPE_NEUMANN: Zero-gradient (∂φ/∂n = 0)") + print(" - BC_TYPE_NEUMANN: Zero-gradient (dphi/dn = 0)") print(" - BC_TYPE_DIRICHLET: Fixed value") print(" - BC_TYPE_NOSLIP: Zero velocity at walls") print(" - BC_TYPE_INLET: Inlet velocity (uniform/parabolic)") diff --git a/examples/output_formats.py b/examples/output_formats.py index f4c56f4..cf1f365 100644 --- a/examples/output_formats.py +++ b/examples/output_formats.py @@ -115,7 +115,10 @@ def main(): create_new=False, ) print(" Appended 5 timesteps") - print(f" Final size: {csv_file.stat().st_size} bytes") + if csv_file.exists(): + print(f" Final size: {csv_file.stat().st_size} bytes") + else: + print(" Note: CSV file created in current working directory") # Output Type Constants print("\n5. Output Type Constants") diff --git a/examples/vtk_output.py b/examples/vtk_output.py index dba1e29..bc9d202 100644 --- a/examples/vtk_output.py +++ b/examples/vtk_output.py @@ -203,12 +203,20 @@ def main(): print("\n5. CSV Output for Analysis") print("-" * 60) - # Write velocity magnitude as CSV - cfd_python.write_csv_field("velocity_magnitude.csv", vel_mag, nx, ny) + # Write velocity magnitude as CSV using Python + with open("velocity_magnitude.csv", "w") as f: + f.write("i,j,value\n") + for j in range(ny): + for i in range(nx): + f.write(f"{i},{j},{vel_mag[j * nx + i]}\n") print(" Written: velocity_magnitude.csv") # Write pressure as CSV - cfd_python.write_csv_field("pressure.csv", fields["pressure"], nx, ny) + with open("pressure.csv", "w") as f: + f.write("i,j,value\n") + for j in range(ny): + for i in range(nx): + f.write(f"{i},{j},{fields['pressure'][j * nx + i]}\n") print(" Written: pressure.csv") # ================================================================= @@ -235,7 +243,7 @@ def main(): fields["u"], fields["v"], fields["pressure"], nx, ny ) print(f" Velocity magnitude max: {flow_stats['velocity_magnitude']['max']:.4f}") - print(f" Pressure max: {flow_stats['pressure']['max']:.4f}") + print(f" Pressure max: {flow_stats['p']['max']:.4f}") # ================================================================= # Summary From 8f395802febe421125bccffd33749741eb8c9307 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 19:56:55 +0200 Subject: [PATCH 17/22] feat: Add matplotlib visualization example New example demonstrating CFD result visualization with matplotlib: - Velocity magnitude contour plots - Multi-panel vortex flow analysis (velocity, pressure, vectors, streamlines) - Centerline profile plots - Convergence history monitoring - Grid resolution comparison - 3D surface visualization --- examples/visualization_matplotlib.py | 346 +++++++++++++++++++++++++++ 1 file changed, 346 insertions(+) create mode 100644 examples/visualization_matplotlib.py diff --git a/examples/visualization_matplotlib.py b/examples/visualization_matplotlib.py new file mode 100644 index 0000000..edda0a5 --- /dev/null +++ b/examples/visualization_matplotlib.py @@ -0,0 +1,346 @@ +#!/usr/bin/env python3 +""" +Matplotlib Visualization Example + +This example demonstrates how to visualize CFD simulation results using +matplotlib. It creates various plots including contour plots, streamlines, +vector fields, and convergence history. + +Requirements: + pip install matplotlib numpy +""" + +import os +import sys + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..")) + +try: + import matplotlib.pyplot as plt + import numpy as np + + import cfd_python +except ImportError as e: + print(f"Import error: {e}") + print("Make sure to install dependencies:") + print(" pip install matplotlib numpy") + print(" pip install -e .") + sys.exit(1) + + +def create_vortex_field(nx, ny, xmin, xmax, ymin, ymax): + """Create a synthetic vortex velocity field for visualization.""" + dx = (xmax - xmin) / (nx - 1) + dy = (ymax - ymin) / (ny - 1) + xc, yc = (xmax + xmin) / 2, (ymax + ymin) / 2 + + u = [0.0] * (nx * ny) + v = [0.0] * (nx * ny) + p = [0.0] * (nx * ny) + + for j in range(ny): + y = ymin + j * dy + for i in range(nx): + x = xmin + i * dx + idx = j * nx + i + + rx, ry = x - xc, y - yc + r = np.sqrt(rx * rx + ry * ry) + 1e-10 + + # Rankine vortex + core_radius = 0.3 + strength = 1.0 + if r < core_radius: + factor = r / core_radius + else: + factor = core_radius / r + + u[idx] = -ry / r * strength * factor + v[idx] = rx / r * strength * factor + p[idx] = 1.0 - 0.5 * (strength * factor) ** 2 + + return u, v, p + + +def main(): + print("CFD Python - Matplotlib Visualization Example") + print("=" * 60) + + # ================================================================= + # 1. Run Simulation + # ================================================================= + print("\n1. Running Simulation") + print("-" * 60) + + nx, ny = 64, 64 + xmin, xmax = 0.0, 1.0 + ymin, ymax = 0.0, 1.0 + + result = cfd_python.run_simulation_with_params( + nx=nx, + ny=ny, + xmin=xmin, + xmax=xmax, + ymin=ymin, + ymax=ymax, + steps=100, + dt=0.001, + cfl=0.2, + ) + + print(f" Grid: {nx} x {ny}") + print(f" Solver: {result['solver_name']}") + + # Get velocity magnitude from simulation + vel_mag_sim = np.array(result["velocity_magnitude"]).reshape(ny, nx) + + # Create coordinate arrays + x = np.linspace(xmin, xmax, nx) + y = np.linspace(ymin, ymax, ny) + X, Y = np.meshgrid(x, y) + + # ================================================================= + # 2. Velocity Magnitude Contour Plot + # ================================================================= + print("\n2. Creating Velocity Magnitude Contour Plot") + print("-" * 60) + + fig, ax = plt.subplots(figsize=(8, 6)) + contour = ax.contourf(X, Y, vel_mag_sim, levels=20, cmap="viridis") + plt.colorbar(contour, ax=ax, label="Velocity Magnitude") + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.set_title("Velocity Magnitude Field") + ax.set_aspect("equal") + plt.savefig("velocity_magnitude_contour.png", dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: velocity_magnitude_contour.png") + + # ================================================================= + # 3. Synthetic Vortex Visualization + # ================================================================= + print("\n3. Creating Vortex Flow Visualization") + print("-" * 60) + + # Create vortex field for richer visualization + u, v, p = create_vortex_field(nx, ny, xmin, xmax, ymin, ymax) + U = np.array(u).reshape(ny, nx) + V = np.array(v).reshape(ny, nx) + P = np.array(p).reshape(ny, nx) + vel_mag = np.sqrt(U**2 + V**2) + + # Create multi-panel figure + fig, axes = plt.subplots(2, 2, figsize=(12, 10)) + + # Panel 1: Velocity magnitude with contours + ax1 = axes[0, 0] + cf1 = ax1.contourf(X, Y, vel_mag, levels=20, cmap="plasma") + ax1.contour(X, Y, vel_mag, levels=10, colors="white", linewidths=0.5, alpha=0.5) + plt.colorbar(cf1, ax=ax1, label="Velocity Magnitude") + ax1.set_xlabel("X") + ax1.set_ylabel("Y") + ax1.set_title("Velocity Magnitude") + ax1.set_aspect("equal") + + # Panel 2: Pressure field + ax2 = axes[0, 1] + cf2 = ax2.contourf(X, Y, P, levels=20, cmap="coolwarm") + plt.colorbar(cf2, ax=ax2, label="Pressure") + ax2.set_xlabel("X") + ax2.set_ylabel("Y") + ax2.set_title("Pressure Field") + ax2.set_aspect("equal") + + # Panel 3: Vector field (quiver plot) + ax3 = axes[1, 0] + skip = 4 # Skip every N points for clarity + ax3.quiver( + X[::skip, ::skip], + Y[::skip, ::skip], + U[::skip, ::skip], + V[::skip, ::skip], + vel_mag[::skip, ::skip], + cmap="viridis", + scale=15, + ) + ax3.set_xlabel("X") + ax3.set_ylabel("Y") + ax3.set_title("Velocity Vectors") + ax3.set_aspect("equal") + + # Panel 4: Streamlines + ax4 = axes[1, 1] + strm = ax4.streamplot(X, Y, U, V, color=vel_mag, cmap="viridis", density=2, linewidth=1) + plt.colorbar(strm.lines, ax=ax4, label="Velocity Magnitude") + ax4.set_xlabel("X") + ax4.set_ylabel("Y") + ax4.set_title("Streamlines") + ax4.set_aspect("equal") + + plt.tight_layout() + plt.savefig("vortex_visualization.png", dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: vortex_visualization.png") + + # ================================================================= + # 4. Centerline Profiles + # ================================================================= + print("\n4. Creating Centerline Profile Plots") + print("-" * 60) + + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + + # Horizontal centerline (y = 0.5) + j_center = ny // 2 + ax1 = axes[0] + ax1.plot(x, vel_mag[j_center, :], "b-", linewidth=2, label="Velocity Magnitude") + ax1.plot(x, U[j_center, :], "r--", linewidth=1.5, label="U velocity") + ax1.plot(x, V[j_center, :], "g--", linewidth=1.5, label="V velocity") + ax1.set_xlabel("X") + ax1.set_ylabel("Velocity") + ax1.set_title(f"Horizontal Centerline (y = {y[j_center]:.2f})") + ax1.legend() + ax1.grid(True, alpha=0.3) + + # Vertical centerline (x = 0.5) + i_center = nx // 2 + ax2 = axes[1] + ax2.plot(y, vel_mag[:, i_center], "b-", linewidth=2, label="Velocity Magnitude") + ax2.plot(y, U[:, i_center], "r--", linewidth=1.5, label="U velocity") + ax2.plot(y, V[:, i_center], "g--", linewidth=1.5, label="V velocity") + ax2.set_xlabel("Y") + ax2.set_ylabel("Velocity") + ax2.set_title(f"Vertical Centerline (x = {x[i_center]:.2f})") + ax2.legend() + ax2.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig("centerline_profiles.png", dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: centerline_profiles.png") + + # ================================================================= + # 5. Convergence History + # ================================================================= + print("\n5. Creating Convergence History Plot") + print("-" * 60) + + # Run multiple steps and track convergence + max_velocities = [] + residuals = [] + steps_list = [] + prev_vel = None + + for step in range(50): + result = cfd_python.run_simulation_with_params( + nx=32, ny=32, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=1, dt=0.001 + ) + vel = np.array(result["velocity_magnitude"]) + stats = cfd_python.calculate_field_stats(result["velocity_magnitude"]) + max_velocities.append(stats["max"]) + steps_list.append(step) + + if prev_vel is not None: + residual = np.mean(np.abs(vel - prev_vel)) + residuals.append(residual) + else: + residuals.append(np.nan) + prev_vel = vel.copy() + + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + + # Max velocity history + ax1 = axes[0] + ax1.plot(steps_list, max_velocities, "b-", linewidth=2) + ax1.set_xlabel("Step") + ax1.set_ylabel("Max Velocity") + ax1.set_title("Maximum Velocity History") + ax1.grid(True, alpha=0.3) + + # Residual history (skip first NaN, filter zeros for log scale) + ax2 = axes[1] + valid_residuals = [r if r > 0 else 1e-16 for r in residuals[1:]] + ax2.semilogy(steps_list[1:], valid_residuals, "r-", linewidth=2) + ax2.set_xlabel("Step") + ax2.set_ylabel("Residual (L1 norm)") + ax2.set_title("Convergence History") + ax2.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig("convergence_history.png", dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: convergence_history.png") + + # ================================================================= + # 6. Grid Comparison + # ================================================================= + print("\n6. Creating Grid Resolution Comparison") + print("-" * 60) + + fig, axes = plt.subplots(1, 3, figsize=(15, 4)) + grid_sizes = [(16, 16), (32, 32), (64, 64)] + + for idx, (gx, gy) in enumerate(grid_sizes): + result = cfd_python.run_simulation_with_params( + nx=gx, ny=gy, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=50, dt=0.001 + ) + vel = np.array(result["velocity_magnitude"]).reshape(gy, gx) + xx = np.linspace(0, 1, gx) + yy = np.linspace(0, 1, gy) + XX, YY = np.meshgrid(xx, yy) + + ax = axes[idx] + cf = ax.contourf(XX, YY, vel, levels=15, cmap="viridis") + plt.colorbar(cf, ax=ax, label="Velocity") + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.set_title(f"Grid: {gx}x{gy}") + ax.set_aspect("equal") + + plt.tight_layout() + plt.savefig("grid_comparison.png", dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: grid_comparison.png") + + # ================================================================= + # 7. 3D Surface Plot + # ================================================================= + print("\n7. Creating 3D Surface Plot") + print("-" * 60) + + fig = plt.figure(figsize=(10, 8)) + ax = fig.add_subplot(111, projection="3d") + + # Use the vortex velocity magnitude + surf = ax.plot_surface(X, Y, vel_mag, cmap="viridis", edgecolor="none", alpha=0.8) + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.set_zlabel("Velocity Magnitude") + ax.set_title("3D Velocity Magnitude Surface") + fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label="Velocity") + plt.savefig("velocity_surface_3d.png", dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: velocity_surface_3d.png") + + # ================================================================= + # Summary + # ================================================================= + print("\n" + "=" * 60) + print("Visualization Example Complete!") + print("=" * 60) + print("\nGenerated files:") + print(" - velocity_magnitude_contour.png: Basic contour plot") + print(" - vortex_visualization.png: Multi-panel vortex analysis") + print(" - centerline_profiles.png: 1D profile plots") + print(" - convergence_history.png: Convergence monitoring") + print(" - grid_comparison.png: Grid resolution study") + print(" - velocity_surface_3d.png: 3D surface visualization") + print("\nVisualization tips:") + print(" - Use contourf for filled contours, contour for lines") + print(" - Use quiver for vector fields (skip points for clarity)") + print(" - Use streamplot for streamlines") + print(" - Use semilogy for convergence plots") + + +if __name__ == "__main__": + main() From cb73ed6296723c46403519d16c2fecafab6bb4f1 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 20:23:19 +0200 Subject: [PATCH 18/22] refactor: Write all example outputs to examples/output/ directory - Update all examples to write VTK, CSV, and PNG files to output/ subdirectory - Add examples/output/ to .gitignore - Update examples/README.md with complete list of all 14 examples - Fix convergence plot to show meaningful data (velocity vs steps) Examples updated: - basic_example.py - channel_flow.py - derived_fields.py - lid_driven_cavity_advanced.py - visualization_matplotlib.py - visualization_numpy.py - vtk_output.py --- .gitignore | 4 ++ examples/README.md | 41 ++++++++---- examples/basic_example.py | 11 ++-- examples/channel_flow.py | 21 ++++-- examples/derived_fields.py | 16 +++-- examples/lid_driven_cavity_advanced.py | 8 +-- examples/visualization_matplotlib.py | 88 ++++++++++++++------------ examples/visualization_numpy.py | 8 ++- examples/vtk_output.py | 81 +++++++++++++++++------- 9 files changed, 181 insertions(+), 97 deletions(-) diff --git a/.gitignore b/.gitignore index 5a02d47..b669c7f 100644 --- a/.gitignore +++ b/.gitignore @@ -154,6 +154,10 @@ src/cfd_lib/ *.vtu *.pvd +# Example output files +examples/output/ +*.csv + # IDE .vscode/ .idea/ diff --git a/examples/README.md b/examples/README.md index 190244b..24de722 100644 --- a/examples/README.md +++ b/examples/README.md @@ -2,7 +2,7 @@ This directory contains example scripts demonstrating how to use the CFD Python library. -## Examples +## Examples Overview ### Basic Usage @@ -10,20 +10,37 @@ This directory contains example scripts demonstrating how to use the CFD Python |---------|-------------| | [basic_example.py](basic_example.py) | Fundamental usage of all main functions | | [solver_discovery.py](solver_discovery.py) | Discovering and using different solvers | +| [error_handling.py](error_handling.py) | Error handling API with exceptions | + +### Simulation Examples + +| Example | Description | +|---------|-------------| +| [lid_driven_cavity.py](lid_driven_cavity.py) | Classic CFD benchmark problem | +| [lid_driven_cavity_advanced.py](lid_driven_cavity_advanced.py) | Advanced cavity simulation with time stepping | +| [channel_flow.py](channel_flow.py) | Channel flow with inlet/outlet BCs | +| [parameter_study.py](parameter_study.py) | Running parameter studies and comparisons | + +### Boundary Conditions + +| Example | Description | +|---------|-------------| +| [boundary_conditions.py](boundary_conditions.py) | Full boundary conditions API demo | ### Output and Visualization | Example | Description | |---------|-------------| | [output_formats.py](output_formats.py) | Exporting to VTK and CSV formats | +| [vtk_output.py](vtk_output.py) | VTK output for ParaView visualization | | [visualization_numpy.py](visualization_numpy.py) | NumPy analysis of simulation results | +| [visualization_matplotlib.py](visualization_matplotlib.py) | Matplotlib plots (contours, vectors, streamlines) | -### Advanced Usage +### Derived Fields | Example | Description | |---------|-------------| -| [parameter_study.py](parameter_study.py) | Running parameter studies and comparisons | -| [lid_driven_cavity.py](lid_driven_cavity.py) | Classic CFD benchmark problem | +| [derived_fields.py](derived_fields.py) | Computing derived quantities and statistics | ## Running Examples @@ -32,9 +49,6 @@ Make sure the package is installed first: ```bash # Install in development mode pip install -e . - -# Or using dev_build.py -python dev_build.py develop ``` Then run any example: @@ -42,16 +56,17 @@ Then run any example: ```bash cd examples python basic_example.py -python solver_discovery.py -python output_formats.py +python visualization_matplotlib.py ``` ## Output Files -Examples may generate output files: +All examples write output files to the `output/` subdirectory. This directory is automatically created when running examples and is excluded from git. -- `.vtk` files - Open with ParaView or VisIt for visualization +Output file types: +- `.vtk` files - Open with ParaView or VisIt for 3D visualization - `.csv` files - Open with Excel, Python pandas, or any text editor +- `.png` files - Matplotlib visualizations ## Requirements @@ -59,5 +74,5 @@ Examples may generate output files: - `numpy` (included as dependency) Optional for visualization: -- `matplotlib` for plotting -- ParaView or VisIt for VTK visualization +- `matplotlib` for plotting (visualization_matplotlib.py) +- ParaView or VisIt for VTK file visualization diff --git a/examples/basic_example.py b/examples/basic_example.py index 47ab09c..35d8a28 100644 --- a/examples/basic_example.py +++ b/examples/basic_example.py @@ -27,6 +27,10 @@ def main(): print("=" * 50) print(f"Version: {cfd_python.__version__}") + # Create output directory + output_dir = os.path.join(os.path.dirname(__file__), "output") + os.makedirs(output_dir, exist_ok=True) + # Show available solvers print("\nAvailable Solvers:") for solver in cfd_python.list_solvers(): @@ -62,10 +66,9 @@ def main(): print("2. Simulation with VTK Output") print("-" * 50) - vel_mag_vtk = cfd_python.run_simulation( - nx=nx, ny=ny, steps=steps, output_file="basic_output.vtk" - ) - print(f" Output: basic_output.vtk ({len(vel_mag_vtk)} points)") + output_file = os.path.join(output_dir, "basic_output.vtk") + vel_mag_vtk = cfd_python.run_simulation(nx=nx, ny=ny, steps=steps, output_file=output_file) + print(f" Output: output/basic_output.vtk ({len(vel_mag_vtk)} points)") # Method 3: Using run_simulation_with_params print("\n" + "-" * 50) diff --git a/examples/channel_flow.py b/examples/channel_flow.py index d0ae59e..08c7ee7 100644 --- a/examples/channel_flow.py +++ b/examples/channel_flow.py @@ -103,6 +103,10 @@ def main(): print("CFD Python - Channel Flow Example") print("=" * 60) + # Create output directory + output_dir = os.path.join(os.path.dirname(__file__), "output") + os.makedirs(output_dir, exist_ok=True) + # Parameters nx, ny = 64, 32 # Longer in x for channel flow xmin, xmax = 0.0, 4.0 # Channel length = 4 @@ -181,7 +185,7 @@ def main(): # Velocity magnitude cfd_python.write_vtk_scalar( - "channel_velocity_magnitude.vtk", + os.path.join(output_dir, "channel_velocity_magnitude.vtk"), "velocity_magnitude", vel_mag, nx, @@ -194,11 +198,20 @@ def main(): # Velocity vectors cfd_python.write_vtk_vector( - "channel_velocity.vtk", "velocity", u, v, nx, ny, xmin, xmax, ymin, ymax + os.path.join(output_dir, "channel_velocity.vtk"), + "velocity", + u, + v, + nx, + ny, + xmin, + xmax, + ymin, + ymax, ) - print(" channel_velocity_magnitude.vtk") - print(" channel_velocity.vtk") + print(" output/channel_velocity_magnitude.vtk") + print(" output/channel_velocity.vtk") # Cross-section profiles print("\n" + "=" * 60) diff --git a/examples/derived_fields.py b/examples/derived_fields.py index 02f239d..99005a4 100644 --- a/examples/derived_fields.py +++ b/examples/derived_fields.py @@ -26,6 +26,10 @@ def main(): print("CFD Python Derived Fields Example") print("=" * 60) + # Create output directory + output_dir = os.path.join(os.path.dirname(__file__), "output") + os.makedirs(output_dir, exist_ok=True) + # ================================================================= # 1. Run a Simulation to Get Flow Fields # ================================================================= @@ -169,7 +173,7 @@ def main(): # Write velocity magnitude cfd_python.write_vtk_scalar( - "velocity_magnitude.vtk", + os.path.join(output_dir, "velocity_magnitude.vtk"), "velocity_magnitude", vel_mag, nx, @@ -179,11 +183,11 @@ def main(): grid["ymin"], grid["ymax"], ) - print(" Written: velocity_magnitude.vtk") + print(" Written: output/velocity_magnitude.vtk") # Write pressure cfd_python.write_vtk_scalar( - "pressure.vtk", + os.path.join(output_dir, "pressure.vtk"), "pressure", p, nx, @@ -193,11 +197,11 @@ def main(): grid["ymin"], grid["ymax"], ) - print(" Written: pressure.vtk") + print(" Written: output/pressure.vtk") # Write velocity vectors cfd_python.write_vtk_vector( - "velocity_field.vtk", + os.path.join(output_dir, "velocity_field.vtk"), "velocity", u, v, @@ -208,7 +212,7 @@ def main(): grid["ymin"], grid["ymax"], ) - print(" Written: velocity_field.vtk") + print(" Written: output/velocity_field.vtk") # ================================================================= # Summary diff --git a/examples/lid_driven_cavity_advanced.py b/examples/lid_driven_cavity_advanced.py index 3da17cf..5825140 100644 --- a/examples/lid_driven_cavity_advanced.py +++ b/examples/lid_driven_cavity_advanced.py @@ -84,7 +84,7 @@ def run_cavity_simulation(nx=32, ny=32, Re=100, steps=500, output_interval=100): print(f"dx={dx:.4f}, dy={dy:.4f}, dt={dt:.6f}") # Set up output directory - output_dir = "cavity_output" + output_dir = os.path.join(os.path.dirname(__file__), "output") os.makedirs(output_dir, exist_ok=True) cfd_python.set_output_dir(output_dir) @@ -141,13 +141,13 @@ def run_cavity_simulation(nx=32, ny=32, Re=100, steps=500, output_interval=100): # Write VTK output if step % output_interval == 0 or step == steps - 1: - filename = f"cavity_step_{step:05d}.vtk" + filename = os.path.join(output_dir, f"cavity_step_{step:05d}.vtk") cfd_python.write_vtk_scalar( filename, "velocity_magnitude", vel_mag, nx, ny, xmin, xmax, ymin, ymax ) print("\nSimulation complete!") - print(f"Output files written to: {output_dir}/") + print("Output files written to: output/") # Final statistics final_stats = cfd_python.calculate_field_stats(vel_mag) @@ -191,7 +191,7 @@ def analyze_results(results): else: print(" Status: NOT CONVERGED (may need more steps)") - print(f"\nOutput directory: {results['output_dir']}") + print("\nOutput directory: output/") def main(): diff --git a/examples/visualization_matplotlib.py b/examples/visualization_matplotlib.py index edda0a5..858f763 100644 --- a/examples/visualization_matplotlib.py +++ b/examples/visualization_matplotlib.py @@ -66,6 +66,11 @@ def main(): print("CFD Python - Matplotlib Visualization Example") print("=" * 60) + # Create output directory for generated images + output_dir = os.path.join(os.path.dirname(__file__), "output") + os.makedirs(output_dir, exist_ok=True) + print(f"\nOutput directory: {output_dir}") + # ================================================================= # 1. Run Simulation # ================================================================= @@ -112,9 +117,11 @@ def main(): ax.set_ylabel("Y") ax.set_title("Velocity Magnitude Field") ax.set_aspect("equal") - plt.savefig("velocity_magnitude_contour.png", dpi=150, bbox_inches="tight") + plt.savefig( + os.path.join(output_dir, "velocity_magnitude_contour.png"), dpi=150, bbox_inches="tight" + ) plt.close() - print(" Saved: velocity_magnitude_contour.png") + print(" Saved: output/velocity_magnitude_contour.png") # ================================================================= # 3. Synthetic Vortex Visualization @@ -178,9 +185,9 @@ def main(): ax4.set_aspect("equal") plt.tight_layout() - plt.savefig("vortex_visualization.png", dpi=150, bbox_inches="tight") + plt.savefig(os.path.join(output_dir, "vortex_visualization.png"), dpi=150, bbox_inches="tight") plt.close() - print(" Saved: vortex_visualization.png") + print(" Saved: output/vortex_visualization.png") # ================================================================= # 4. Centerline Profiles @@ -215,61 +222,62 @@ def main(): ax2.grid(True, alpha=0.3) plt.tight_layout() - plt.savefig("centerline_profiles.png", dpi=150, bbox_inches="tight") + plt.savefig(os.path.join(output_dir, "centerline_profiles.png"), dpi=150, bbox_inches="tight") plt.close() - print(" Saved: centerline_profiles.png") + print(" Saved: output/centerline_profiles.png") # ================================================================= - # 5. Convergence History + # 5. Convergence History (Grid Resolution Study) # ================================================================= - print("\n5. Creating Convergence History Plot") + print("\n5. Creating Convergence Analysis Plot") print("-" * 60) - # Run multiple steps and track convergence + # Since each run_simulation_with_params() call is independent, + # we show convergence by comparing results at different step counts + step_counts = [10, 20, 30, 50, 75, 100, 150, 200] max_velocities = [] - residuals = [] - steps_list = [] - prev_vel = None + avg_velocities = [] - for step in range(50): + for steps in step_counts: result = cfd_python.run_simulation_with_params( - nx=32, ny=32, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=1, dt=0.001 + nx=32, ny=32, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=steps, dt=0.001 ) - vel = np.array(result["velocity_magnitude"]) stats = cfd_python.calculate_field_stats(result["velocity_magnitude"]) max_velocities.append(stats["max"]) - steps_list.append(step) + avg_velocities.append(stats["avg"]) - if prev_vel is not None: - residual = np.mean(np.abs(vel - prev_vel)) - residuals.append(residual) - else: - residuals.append(np.nan) - prev_vel = vel.copy() + # Compute relative change between consecutive step counts + relative_changes = [] + for i in range(1, len(max_velocities)): + change = abs(max_velocities[i] - max_velocities[i - 1]) / max_velocities[i - 1] + relative_changes.append(change) fig, axes = plt.subplots(1, 2, figsize=(12, 5)) - # Max velocity history + # Velocity vs steps ax1 = axes[0] - ax1.plot(steps_list, max_velocities, "b-", linewidth=2) - ax1.set_xlabel("Step") - ax1.set_ylabel("Max Velocity") - ax1.set_title("Maximum Velocity History") + ax1.plot(step_counts, max_velocities, "b-o", linewidth=2, label="Max Velocity") + ax1.plot(step_counts, avg_velocities, "g--s", linewidth=2, label="Avg Velocity") + ax1.set_xlabel("Number of Steps") + ax1.set_ylabel("Velocity") + ax1.set_title("Velocity vs Simulation Steps") + ax1.legend() ax1.grid(True, alpha=0.3) - # Residual history (skip first NaN, filter zeros for log scale) + # Relative change (convergence indicator) ax2 = axes[1] - valid_residuals = [r if r > 0 else 1e-16 for r in residuals[1:]] - ax2.semilogy(steps_list[1:], valid_residuals, "r-", linewidth=2) - ax2.set_xlabel("Step") - ax2.set_ylabel("Residual (L1 norm)") - ax2.set_title("Convergence History") + ax2.semilogy(step_counts[1:], relative_changes, "r-o", linewidth=2) + ax2.set_xlabel("Number of Steps") + ax2.set_ylabel("Relative Change") + ax2.set_title("Convergence (Relative Change in Max Velocity)") ax2.grid(True, alpha=0.3) + ax2.axhline(y=0.01, color="k", linestyle="--", alpha=0.5, label="1% threshold") + ax2.legend() plt.tight_layout() - plt.savefig("convergence_history.png", dpi=150, bbox_inches="tight") + plt.savefig(os.path.join(output_dir, "convergence_history.png"), dpi=150, bbox_inches="tight") plt.close() - print(" Saved: convergence_history.png") + print(" Saved: output/convergence_history.png") # ================================================================= # 6. Grid Comparison @@ -298,9 +306,9 @@ def main(): ax.set_aspect("equal") plt.tight_layout() - plt.savefig("grid_comparison.png", dpi=150, bbox_inches="tight") + plt.savefig(os.path.join(output_dir, "grid_comparison.png"), dpi=150, bbox_inches="tight") plt.close() - print(" Saved: grid_comparison.png") + print(" Saved: output/grid_comparison.png") # ================================================================= # 7. 3D Surface Plot @@ -318,9 +326,9 @@ def main(): ax.set_zlabel("Velocity Magnitude") ax.set_title("3D Velocity Magnitude Surface") fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label="Velocity") - plt.savefig("velocity_surface_3d.png", dpi=150, bbox_inches="tight") + plt.savefig(os.path.join(output_dir, "velocity_surface_3d.png"), dpi=150, bbox_inches="tight") plt.close() - print(" Saved: velocity_surface_3d.png") + print(" Saved: output/velocity_surface_3d.png") # ================================================================= # Summary @@ -328,7 +336,7 @@ def main(): print("\n" + "=" * 60) print("Visualization Example Complete!") print("=" * 60) - print("\nGenerated files:") + print(f"\nGenerated files in {output_dir}:") print(" - velocity_magnitude_contour.png: Basic contour plot") print(" - vortex_visualization.png: Multi-panel vortex analysis") print(" - centerline_profiles.png: 1D profile plots") diff --git a/examples/visualization_numpy.py b/examples/visualization_numpy.py index 1e391ef..8071f49 100644 --- a/examples/visualization_numpy.py +++ b/examples/visualization_numpy.py @@ -23,6 +23,10 @@ def main(): print("CFD Python - NumPy Analysis Example") print("=" * 50) + # Create output directory + output_dir = os.path.join(os.path.dirname(__file__), "output") + os.makedirs(output_dir, exist_ok=True) + # Run simulation nx, ny = 30, 30 print(f"\n1. Running simulation ({nx}x{ny} grid, 50 steps)...") @@ -124,10 +128,10 @@ def main(): print("-" * 50) # Save as CSV for external tools - output_file = "velocity_data.csv" + output_file = os.path.join(output_dir, "velocity_data.csv") flat_data = np.column_stack([X.flatten(), Y.flatten(), vel_mag.flatten()]) np.savetxt(output_file, flat_data, delimiter=",", header="x,y,velocity_magnitude", comments="") - print(f" Saved to: {output_file}") + print(" Saved to: output/velocity_data.csv") print(f" Shape: {flat_data.shape}") print("\n" + "=" * 50) diff --git a/examples/vtk_output.py b/examples/vtk_output.py index bc9d202..6d3d38b 100644 --- a/examples/vtk_output.py +++ b/examples/vtk_output.py @@ -111,10 +111,10 @@ def main(): print(f" Domain: [{xmin}, {xmax}] x [{ymin}, {ymax}]") # Create output directory - output_dir = "vtk_output" + output_dir = os.path.join(os.path.dirname(__file__), "output") os.makedirs(output_dir, exist_ok=True) cfd_python.set_output_dir(output_dir) - print(f" Output directory: {output_dir}/") + print(" Output directory: output/") # Create sample data fields = create_sample_fields(nx, ny, xmin, xmax, ymin, ymax) @@ -128,22 +128,46 @@ def main(): # Write pressure field cfd_python.write_vtk_scalar( - "pressure.vtk", "pressure", fields["pressure"], nx, ny, xmin, xmax, ymin, ymax + os.path.join(output_dir, "pressure.vtk"), + "pressure", + fields["pressure"], + nx, + ny, + xmin, + xmax, + ymin, + ymax, ) - print(" Written: pressure.vtk") + print(" Written: output/pressure.vtk") # Write temperature field cfd_python.write_vtk_scalar( - "temperature.vtk", "temperature", fields["temperature"], nx, ny, xmin, xmax, ymin, ymax + os.path.join(output_dir, "temperature.vtk"), + "temperature", + fields["temperature"], + nx, + ny, + xmin, + xmax, + ymin, + ymax, ) - print(" Written: temperature.vtk") + print(" Written: output/temperature.vtk") # Compute and write velocity magnitude vel_mag = cfd_python.compute_velocity_magnitude(fields["u"], fields["v"], nx, ny) cfd_python.write_vtk_scalar( - "velocity_magnitude.vtk", "velocity_magnitude", vel_mag, nx, ny, xmin, xmax, ymin, ymax + os.path.join(output_dir, "velocity_magnitude.vtk"), + "velocity_magnitude", + vel_mag, + nx, + ny, + xmin, + xmax, + ymin, + ymax, ) - print(" Written: velocity_magnitude.vtk") + print(" Written: output/velocity_magnitude.vtk") # ================================================================= # 3. Writing Vector Fields @@ -153,9 +177,18 @@ def main(): # Write velocity vector field cfd_python.write_vtk_vector( - "velocity.vtk", "velocity", fields["u"], fields["v"], nx, ny, xmin, xmax, ymin, ymax + os.path.join(output_dir, "velocity.vtk"), + "velocity", + fields["u"], + fields["v"], + nx, + ny, + xmin, + xmax, + ymin, + ymax, ) - print(" Written: velocity.vtk") + print(" Written: output/velocity.vtk") # ================================================================= # 4. Time Series Output @@ -189,12 +222,12 @@ def main(): vel_mag_frame = cfd_python.compute_velocity_magnitude(u_rot, v_rot, nx, ny) # Write with frame number in filename - filename = f"vortex_{frame:04d}.vtk" + filename = os.path.join(output_dir, f"vortex_{frame:04d}.vtk") cfd_python.write_vtk_scalar( filename, "velocity_magnitude", vel_mag_frame, nx, ny, xmin, xmax, ymin, ymax ) - print(f" Written: vortex_0000.vtk through vortex_{n_frames-1:04d}.vtk") + print(f" Written: output/vortex_0000.vtk through output/vortex_{n_frames-1:04d}.vtk") print(" (Load in ParaView as a time series)") # ================================================================= @@ -204,20 +237,20 @@ def main(): print("-" * 60) # Write velocity magnitude as CSV using Python - with open("velocity_magnitude.csv", "w") as f: + with open(os.path.join(output_dir, "velocity_magnitude.csv"), "w") as f: f.write("i,j,value\n") for j in range(ny): for i in range(nx): f.write(f"{i},{j},{vel_mag[j * nx + i]}\n") - print(" Written: velocity_magnitude.csv") + print(" Written: output/velocity_magnitude.csv") # Write pressure as CSV - with open("pressure.csv", "w") as f: + with open(os.path.join(output_dir, "pressure.csv"), "w") as f: f.write("i,j,value\n") for j in range(ny): for i in range(nx): f.write(f"{i},{j},{fields['pressure'][j * nx + i]}\n") - print(" Written: pressure.csv") + print(" Written: output/pressure.csv") # ================================================================= # 6. Field Statistics @@ -251,19 +284,19 @@ def main(): print("\n" + "=" * 60) print("VTK Output Summary") print("-" * 60) - print(f" Output directory: {output_dir}/") + print(" Output directory: output/") print("\n Files created:") print(" Scalar fields:") - print(" - pressure.vtk") - print(" - temperature.vtk") - print(" - velocity_magnitude.vtk") + print(" - output/pressure.vtk") + print(" - output/temperature.vtk") + print(" - output/velocity_magnitude.vtk") print(" Vector fields:") - print(" - velocity.vtk") + print(" - output/velocity.vtk") print(" Time series:") - print(f" - vortex_0000.vtk to vortex_{n_frames-1:04d}.vtk") + print(f" - output/vortex_0000.vtk to output/vortex_{n_frames-1:04d}.vtk") print(" CSV files:") - print(" - velocity_magnitude.csv") - print(" - pressure.csv") + print(" - output/velocity_magnitude.csv") + print(" - output/pressure.csv") print("\n To visualize:") print(" 1. Open ParaView") print(" 2. File -> Open -> Select VTK file") From 8ec62d04b982c7465c6919a4263080663a0ee48e Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 20:34:32 +0200 Subject: [PATCH 19/22] docs: Add Rankine vortex explanation to visualization example - Add detailed docstring explaining Rankine vortex physics - Update 3D plot title and add explanatory comments - Print note about vortex characteristics when running - Update README to mention Rankine vortex in description The 3D surface plot shows velocity magnitude of a Rankine vortex: - "Volcano" shape with dip at center (zero velocity at axis) - Peak ring at core boundary (r = 0.3) - Decay outside core following 1/r --- examples/README.md | 2 +- examples/visualization_matplotlib.py | 40 +++++++++++++++++++++++++--- 2 files changed, 37 insertions(+), 5 deletions(-) diff --git a/examples/README.md b/examples/README.md index 24de722..7b6e30a 100644 --- a/examples/README.md +++ b/examples/README.md @@ -34,7 +34,7 @@ This directory contains example scripts demonstrating how to use the CFD Python | [output_formats.py](output_formats.py) | Exporting to VTK and CSV formats | | [vtk_output.py](vtk_output.py) | VTK output for ParaView visualization | | [visualization_numpy.py](visualization_numpy.py) | NumPy analysis of simulation results | -| [visualization_matplotlib.py](visualization_matplotlib.py) | Matplotlib plots (contours, vectors, streamlines) | +| [visualization_matplotlib.py](visualization_matplotlib.py) | Matplotlib plots (contours, vectors, streamlines, 3D Rankine vortex) | ### Derived Fields diff --git a/examples/visualization_matplotlib.py b/examples/visualization_matplotlib.py index 858f763..6374cf5 100644 --- a/examples/visualization_matplotlib.py +++ b/examples/visualization_matplotlib.py @@ -29,7 +29,24 @@ def create_vortex_field(nx, ny, xmin, xmax, ymin, ymax): - """Create a synthetic vortex velocity field for visualization.""" + """Create a synthetic Rankine vortex velocity field for visualization. + + The Rankine vortex is a classical fluid dynamics model combining: + - Solid body rotation inside the core (r < core_radius): velocity ~ r + - Irrotational (potential) flow outside (r > core_radius): velocity ~ 1/r + + This produces the characteristic "volcano" shape in 3D plots: + - Zero velocity at the center (fluid rotates around, not through the axis) + - Maximum velocity at the core boundary (r = core_radius) + - Decay toward outer edges following 1/r + + Args: + nx, ny: Grid dimensions + xmin, xmax, ymin, ymax: Domain bounds + + Returns: + tuple: (u, v, p) - x-velocity, y-velocity, and pressure fields + """ dx = (xmax - xmin) / (nx - 1) dy = (ymax - ymin) / (ny - 1) xc, yc = (xmax + xmin) / 2, (ymax + ymin) / 2 @@ -47,16 +64,20 @@ def create_vortex_field(nx, ny, xmin, xmax, ymin, ymax): rx, ry = x - xc, y - yc r = np.sqrt(rx * rx + ry * ry) + 1e-10 - # Rankine vortex + # Rankine vortex velocity profile core_radius = 0.3 strength = 1.0 if r < core_radius: + # Inside core: solid body rotation (v ~ r) factor = r / core_radius else: + # Outside core: irrotational flow (v ~ 1/r) factor = core_radius / r + # Tangential velocity components (perpendicular to radius) u[idx] = -ry / r * strength * factor v[idx] = rx / r * strength * factor + # Pressure from Bernoulli (lower in high-velocity regions) p[idx] = 1.0 - 0.5 * (strength * factor) ** 2 return u, v, p @@ -311,11 +332,21 @@ def main(): print(" Saved: output/grid_comparison.png") # ================================================================= - # 7. 3D Surface Plot + # 7. 3D Surface Plot (Rankine Vortex) # ================================================================= print("\n7. Creating 3D Surface Plot") print("-" * 60) + # This plot visualizes a Rankine vortex velocity field: + # - The "volcano" shape shows velocity magnitude vs position + # - Center dip: zero velocity at vortex core (fluid rotates around, not through) + # - Peak ring: maximum velocity at core radius boundary (r = 0.3) + # - Outer decay: velocity decreases as 1/r outside the core + # + # Rankine vortex combines: + # - Solid body rotation inside core (v ~ r) + # - Irrotational flow outside core (v ~ 1/r) + fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection="3d") @@ -324,11 +355,12 @@ def main(): ax.set_xlabel("X") ax.set_ylabel("Y") ax.set_zlabel("Velocity Magnitude") - ax.set_title("3D Velocity Magnitude Surface") + ax.set_title("3D Velocity Magnitude Surface (Rankine Vortex)") fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label="Velocity") plt.savefig(os.path.join(output_dir, "velocity_surface_3d.png"), dpi=150, bbox_inches="tight") plt.close() print(" Saved: output/velocity_surface_3d.png") + print(" Note: Shows Rankine vortex - peak velocity ring around core, decay at edges") # ================================================================= # Summary From 1a4d4345026c048f32046e107025bb67718a3c2d Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 20:42:03 +0200 Subject: [PATCH 20/22] docs: Add physics explanations to all example files Enhanced docstrings with educational content explaining: - basic_example.py: CFD fundamentals (grids, time stepping, CFL condition) - channel_flow.py: Poiseuille flow physics, boundary conditions, parabolic profile - derived_fields.py: Velocity magnitude, vorticity, kinetic energy, statistics - vtk_output.py: VTK format types, data types, ParaView visualization - lid_driven_cavity_advanced.py: Cavity physics, Reynolds number effects, convergence - visualization_numpy.py: NumPy array operations, gradients, meshgrid Each example now includes relevant physics background and numerical concepts to help users understand what the code is doing and why. --- examples/basic_example.py | 22 +++++++++++++++++ examples/channel_flow.py | 19 ++++++++++++++ examples/derived_fields.py | 25 +++++++++++++++++++ examples/lid_driven_cavity_advanced.py | 26 ++++++++++++++++++++ examples/visualization_numpy.py | 28 +++++++++++++++++++++ examples/vtk_output.py | 34 +++++++++++++++++++++----- 6 files changed, 148 insertions(+), 6 deletions(-) diff --git a/examples/basic_example.py b/examples/basic_example.py index 35d8a28..cb106d4 100644 --- a/examples/basic_example.py +++ b/examples/basic_example.py @@ -4,6 +4,15 @@ This example demonstrates the fundamental usage of the CFD Python bindings for running fluid dynamics simulations. + +CFD (Computational Fluid Dynamics) simulates fluid flow by solving the +Navier-Stokes equations on a discrete grid. Key concepts: + +- Grid: The computational domain is divided into nx × ny cells +- Time stepping: The simulation advances in discrete time steps (dt) +- CFL number: Courant-Friedrichs-Lewy condition ensures numerical stability + (CFL = u*dt/dx, typically kept < 1 for explicit methods) +- Velocity magnitude: sqrt(u² + v²) where u,v are velocity components """ import os @@ -48,6 +57,9 @@ def main(): print(f" Steps: {steps}") # Method 1: Simple simulation + # run_simulation() is the simplest API - just specify grid size and steps. + # It uses default boundary conditions (lid-driven cavity: moving top wall) + # and returns velocity magnitude at each grid point. print("\n" + "-" * 50) print("1. Simple Simulation (run_simulation)") print("-" * 50) @@ -62,6 +74,9 @@ def main(): print(f" Mean velocity: {np.mean(vel_array):.6f}") # Method 2: With VTK output + # VTK (Visualization Toolkit) format is widely used for scientific visualization. + # The output file can be opened in ParaView, VisIt, or other VTK-compatible tools + # to visualize the flow field in 2D/3D with contours, vectors, streamlines, etc. print("\n" + "-" * 50) print("2. Simulation with VTK Output") print("-" * 50) @@ -71,6 +86,10 @@ def main(): print(f" Output: output/basic_output.vtk ({len(vel_mag_vtk)} points)") # Method 3: Using run_simulation_with_params + # This API gives more control over simulation parameters: + # - dt: Time step size (smaller = more accurate but slower) + # - cfl: CFL number for stability (typically 0.1-0.5 for explicit methods) + # Returns a dict with results and metadata including solver info. print("\n" + "-" * 50) print("3. Custom Parameters (run_simulation_with_params)") print("-" * 50) @@ -86,6 +105,9 @@ def main(): print(f" Iterations: {stats.get('iterations', 'N/A')}") # Method 4: Grid and parameter inspection + # create_grid() generates the computational mesh with coordinates. + # get_default_solver_params() shows the default numerical parameters. + # These are useful for setting up custom simulations or post-processing. print("\n" + "-" * 50) print("4. Grid and Parameter Inspection") print("-" * 50) diff --git a/examples/channel_flow.py b/examples/channel_flow.py index 08c7ee7..0358a40 100644 --- a/examples/channel_flow.py +++ b/examples/channel_flow.py @@ -7,6 +7,25 @@ - Zero-gradient outlet - No-slip walls on top and bottom - Periodic or developed flow analysis + +Physics Background: +------------------- +Channel flow (Poiseuille flow) is fluid flow between two parallel plates driven +by a pressure gradient. Key characteristics: + +1. Parabolic velocity profile: Due to viscosity, fluid at the walls has zero + velocity (no-slip), while maximum velocity occurs at the centerline. + u(y) = u_max * 4 * (y/H) * (1 - y/H) + +2. Fully developed flow: After sufficient distance from the inlet, the velocity + profile becomes invariant along the flow direction (du/dx = 0). + +3. Boundary conditions: + - Inlet: Prescribed velocity (parabolic profile for developed flow) + - Outlet: Zero-gradient (Neumann BC) allows flow to exit freely + - Walls: No-slip condition (velocity = 0 at solid surfaces) + +This is a fundamental CFD validation case with known analytical solutions. """ import os diff --git a/examples/derived_fields.py b/examples/derived_fields.py index 99005a4..74e184c 100644 --- a/examples/derived_fields.py +++ b/examples/derived_fields.py @@ -4,6 +4,31 @@ This example demonstrates how to compute derived quantities and statistics from flow fields using cfd_python. + +Derived Fields in CFD: +---------------------- +Raw simulation outputs are typically velocity components (u, v) and pressure (p). +From these, we compute derived quantities for analysis: + +1. Velocity Magnitude: |V| = sqrt(u² + v²) + - Shows flow speed regardless of direction + - Useful for identifying high-speed regions, jets, wakes + +2. Field Statistics (min, max, avg, sum): + - Monitor simulation health (unbounded values indicate instability) + - Track convergence (values should stabilize over time) + - Validate against expected physical behavior + +3. Vorticity: ω = ∂v/∂x - ∂u/∂y + - Measures local rotation in the flow + - Identifies vortices and shear layers + +4. Kinetic Energy: KE = 0.5 * ρ * (u² + v²) + - Total energy in the flow + - Should be conserved or decay (not grow) for physical solutions + +These post-processing tools are essential for validating simulations and +extracting meaningful physical insights from the raw data. """ import os diff --git a/examples/lid_driven_cavity_advanced.py b/examples/lid_driven_cavity_advanced.py index 5825140..28330a2 100644 --- a/examples/lid_driven_cavity_advanced.py +++ b/examples/lid_driven_cavity_advanced.py @@ -7,6 +7,32 @@ - Time stepping with convergence monitoring - Post-processing and visualization output - Statistics tracking + +Physics Background: +------------------- +The lid-driven cavity is a fundamental CFD benchmark problem. A square (or +rectangular) cavity has: +- Three stationary walls (no-slip: u = v = 0) +- One moving wall (lid) with prescribed velocity (u = U_lid, v = 0) + +The moving lid drags fluid along, creating: +1. A primary vortex in the center rotating with the lid +2. Secondary corner vortices at low Reynolds numbers +3. Complex vortex structures at high Reynolds numbers + +Reynolds Number (Re): +Re = U_lid * L / ν where L is cavity size, ν is kinematic viscosity + +- Re < 100: Single primary vortex, steady flow +- Re ~ 1000: Secondary vortices appear in corners +- Re > 3000: Flow becomes unsteady, then turbulent + +Numerical Considerations: +- Time step limited by CFL condition (advection) and diffusion stability +- Grid resolution affects vortex capture and accuracy +- Convergence monitored by velocity change between steps + +This is the standard validation case for incompressible Navier-Stokes solvers. """ import os diff --git a/examples/visualization_numpy.py b/examples/visualization_numpy.py index 8071f49..f1ef22e 100644 --- a/examples/visualization_numpy.py +++ b/examples/visualization_numpy.py @@ -7,6 +7,34 @@ Note: This example uses only NumPy. For actual plotting, you would use matplotlib (see visualization_matplotlib.py). + +NumPy for CFD Post-Processing: +------------------------------ +NumPy is essential for CFD analysis because: + +1. Array Operations: + - Simulation returns flat lists; reshape to 2D grids for spatial analysis + - Vectorized operations are much faster than Python loops + - Example: vel_mag = np.sqrt(u**2 + v**2) # computes for all points at once + +2. Statistical Analysis: + - np.min/max/mean/std: Basic field statistics + - np.percentile: Distribution analysis (e.g., 99th percentile for outliers) + - Spatial statistics: centerlines, profiles, quadrant analysis + +3. Gradient Computation: + - np.gradient: Numerical derivatives for vorticity, strain rate + - Identifies regions of high shear, flow separation + +4. Data Export: + - np.savetxt: Export to CSV for external tools (Excel, MATLAB) + - np.save/load: Binary format for large datasets + +5. Meshgrid for Visualization: + - np.meshgrid: Create X, Y coordinate arrays for contour plots + - Pairs spatial coordinates with field values + +This example shows common analysis patterns used in CFD workflows. """ import os diff --git a/examples/vtk_output.py b/examples/vtk_output.py index 6d3d38b..869c53c 100644 --- a/examples/vtk_output.py +++ b/examples/vtk_output.py @@ -5,12 +5,34 @@ This example demonstrates the VTK output capabilities of cfd_python for visualizing simulation results in ParaView or other VTK-compatible tools. -VTK (Visualization Toolkit) format is widely used for scientific visualization. -This example shows how to: -- Write scalar fields (pressure, temperature, velocity magnitude) -- Write vector fields (velocity, vorticity) -- Create time series output for animations -- Configure output directories +VTK (Visualization Toolkit) Format: +----------------------------------- +VTK is the standard format for scientific visualization. Key concepts: + +1. File Types: + - .vtk (legacy format): Simple ASCII/binary, widely supported + - .vtu (XML format): More features, better for large datasets + - .pvd (ParaView Data): Collection file for time series + +2. Data Types: + - Scalar: Single value per point (pressure, temperature, speed) + - Vector: 2D/3D vector per point (velocity, force) + - Tensor: Matrix per point (stress, strain) + +3. Grid Types: + - Structured: Regular nx × ny grid (used here) + - Unstructured: Arbitrary connectivity (complex geometries) + - Rectilinear: Non-uniform spacing + +Visualization in ParaView: +- Contours: Show iso-lines/surfaces of scalar fields +- Glyphs: Arrows showing vector direction and magnitude +- Streamlines: Particle paths following velocity field +- Slice: Cut through 3D data + +Time Series: +Writing files with sequential names (field_0000.vtk, field_0001.vtk, ...) +allows ParaView to animate the simulation evolution. """ import os From 40db79c6aa2b565b39aa3d2554c26dad531278f6 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 20:59:53 +0200 Subject: [PATCH 21/22] fix: Correct lid_driven_cavity_advanced convergence tracking - Change from ineffective single-step loop to convergence study sampling at key step counts (1, 5, 10, 25, 50, 100, 150, 200) - Document API limitation: run_simulation_with_params() runs fresh simulation each call, doesn't preserve state between calls - Replace unused Re parameter study with meaningful grid resolution study that shows actual variation - Update column headers to reflect actual data (Max Vel, Avg Vel) --- examples/lid_driven_cavity_advanced.py | 80 ++++++++++++++++---------- 1 file changed, 50 insertions(+), 30 deletions(-) diff --git a/examples/lid_driven_cavity_advanced.py b/examples/lid_driven_cavity_advanced.py index 28330a2..7811feb 100644 --- a/examples/lid_driven_cavity_advanced.py +++ b/examples/lid_driven_cavity_advanced.py @@ -82,10 +82,17 @@ def setup_lid_driven_cavity(nx, ny, lid_velocity=1.0): def run_cavity_simulation(nx=32, ny=32, Re=100, steps=500, output_interval=100): """Run lid-driven cavity simulation. + Note: This example demonstrates the API usage pattern for monitoring convergence + and writing output. The current `run_simulation_with_params()` API runs a fresh + simulation each call and returns velocity magnitude - it doesn't expose the + internal u, v, p fields for true time-stepping continuation. Each iteration + here runs an independent simulation with increasing step counts to show how + the solution evolves. + Args: nx, ny: Grid dimensions Re: Reynolds number - steps: Total time steps + steps: Total time steps (used for progress monitoring) output_interval: Steps between VTK outputs Returns: @@ -117,30 +124,44 @@ def run_cavity_simulation(nx=32, ny=32, Re=100, steps=500, output_interval=100): # Initialize fields u, v, p = setup_lid_driven_cavity(nx, ny, lid_velocity) - # Track convergence - max_velocities = [] - residuals = [] + # Track convergence across different step counts + max_velocities = [] # Max velocity at each sample point + residuals = [] # Relative change between samples - print("\nRunning simulation...") - print(f"{'Step':>6} {'Max U':>12} {'Max V':>12} {'Residual':>12}") + # Run simulation and track convergence by comparing results at different step counts + # Note: The API runs a fresh simulation each call, so we sample at key intervals + # to show how the solution evolves with more time steps. + print("\nRunning convergence study...") + print(f"{'Steps':>6} {'Max Vel':>12} {'Avg Vel':>12} {'Change':>12}") print("-" * 48) - prev_u = u.copy() + prev_max = 0.0 + sample_steps = [1, 5, 10, 25, 50, 100, 150, 200] + sample_steps = [s for s in sample_steps if s <= steps] + if steps not in sample_steps: + sample_steps.append(steps) - for step in range(steps): - # Run simulation step using the library + for sim_steps in sample_steps: try: result = cfd_python.run_simulation_with_params( - nx=nx, ny=ny, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, steps=1, dt=dt, cfl=0.25 + nx=nx, + ny=ny, + xmin=xmin, + xmax=xmax, + ymin=ymin, + ymax=ymax, + steps=sim_steps, + dt=dt, + cfl=0.25, ) # Check status status = cfd_python.get_last_status() if status != cfd_python.CFD_SUCCESS: - cfd_python.raise_for_status(status, context=f"step {step}") + cfd_python.raise_for_status(status, context=f"steps={sim_steps}") except cfd_python.CFDError as e: - print(f"Simulation error at step {step}: {e}") + print(f"Simulation error at steps={sim_steps}: {e}") break # Get velocity magnitude from result @@ -151,23 +172,20 @@ def run_cavity_simulation(nx=32, ny=32, Re=100, steps=500, output_interval=100): max_vel = stats["max"] max_velocities.append(max_vel) - # Compute residual (change in velocity) - if step > 0: - residual = sum(abs(vel_mag[i] - prev_u[i]) for i in range(len(vel_mag))) / len(vel_mag) - residuals.append(residual) + # Compute change from previous sample + if prev_max > 0: + change = abs(max_vel - prev_max) / prev_max + residuals.append(change) else: - residual = float("inf") + change = float("inf") - prev_u = vel_mag.copy() + prev_max = max_vel - # Output progress - if step % max(1, steps // 10) == 0 or step == steps - 1: - u_stats = cfd_python.calculate_field_stats(vel_mag) - print(f"{step:6d} {u_stats['max']:12.6f} {u_stats['avg']:12.6f} {residual:12.2e}") + print(f"{sim_steps:6d} {stats['max']:12.6f} {stats['avg']:12.6f} {change:12.2e}") - # Write VTK output - if step % output_interval == 0 or step == steps - 1: - filename = os.path.join(output_dir, f"cavity_step_{step:05d}.vtk") + # Write VTK output at final step + if sim_steps == steps: + filename = os.path.join(output_dir, "cavity_final.vtk") cfd_python.write_vtk_scalar( filename, "velocity_magnitude", vel_mag, nx, ny, xmin, xmax, ymin, ymax ) @@ -236,17 +254,19 @@ def main(): # Analyze results analyze_results(results) - # Run parameter study (optional) + # Run parameter study: grid resolution effect + # Note: The current API doesn't expose viscosity/Reynolds number parameters, + # so we demonstrate grid resolution effects instead. print("\n" + "=" * 60) - print("Parameter Study: Effect of Reynolds Number") + print("Parameter Study: Effect of Grid Resolution") print("=" * 60) - for Re in [50, 100, 200]: + for grid_size in [16, 32, 48]: result = cfd_python.run_simulation_with_params( - nx=32, ny=32, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=100, dt=0.0001 + nx=grid_size, ny=grid_size, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, steps=100, dt=0.0001 ) stats = cfd_python.calculate_field_stats(result["velocity_magnitude"]) - print(f" Re={Re:3d}: max_vel={stats['max']:.4f}, avg_vel={stats['avg']:.4f}") + print(f" {grid_size}x{grid_size}: max_vel={stats['max']:.4f}, avg_vel={stats['avg']:.4f}") if __name__ == "__main__": From 05dfb46d865c846f99269ec393a11451d0f3f804 Mon Sep 17 00:00:00 2001 From: shaia Date: Sat, 3 Jan 2026 21:01:15 +0200 Subject: [PATCH 22/22] docs: Add missing examples to examples/README.md - Add backend_detection.py (SIMD/OpenMP/CUDA detection) - Add solver_comparison.py (comparing solver implementations) --- examples/README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/README.md b/examples/README.md index 7b6e30a..39ef59f 100644 --- a/examples/README.md +++ b/examples/README.md @@ -10,6 +10,8 @@ This directory contains example scripts demonstrating how to use the CFD Python |---------|-------------| | [basic_example.py](basic_example.py) | Fundamental usage of all main functions | | [solver_discovery.py](solver_discovery.py) | Discovering and using different solvers | +| [solver_comparison.py](solver_comparison.py) | Comparing different solver implementations | +| [backend_detection.py](backend_detection.py) | Detecting SIMD/OpenMP/CUDA backend availability | | [error_handling.py](error_handling.py) | Error handling API with exceptions | ### Simulation Examples