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/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/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/README.md b/examples/README.md index 190244b..39ef59f 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,39 @@ 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 + +| 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, 3D Rankine vortex) | -### 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 +51,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 +58,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 +76,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/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/basic_example.py b/examples/basic_example.py index 47ab09c..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 @@ -27,6 +36,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(): @@ -44,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) @@ -58,16 +74,22 @@ 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) - 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 + # 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) @@ -83,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/boundary_conditions.py b/examples/boundary_conditions.py new file mode 100644 index 0000000..8f0e633 --- /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 (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)") + 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..0358a40 --- /dev/null +++ b/examples/channel_flow.py @@ -0,0 +1,262 @@ +#!/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 + +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 +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 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, cfd_python.BC_EDGE_LEFT) + + # 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 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 + 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, H): + """Compute analytical Poiseuille flow solution for validation. + + This is a reference implementation for comparing simulation results + against the known analytical solution for fully developed channel flow. + + 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 velocity at the channel centerline (y = H/2) + H: Channel height + + Returns: + list: Analytical velocity field for comparison with simulation + """ + 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) + + # 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 + 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 - 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 here are used to verify + # the BC setup (inlet profile, cross-section profiles). + print("\nSetting up channel flow (BC demonstration)...") + u, v, _ = setup_channel_flow(nx, ny, u_max) + + # Verify inlet profile (parabolic velocity distribution) + 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 from simulation result + vel_mag = result["velocity_magnitude"] + vel_mag_stats = cfd_python.calculate_field_stats(vel_mag) + + print("\nSimulation Results:") + print(" Velocity magnitude:") + 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 = [] + 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, 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: {vel_mag_stats['max']:.6f}") + + # Write output + print("\nWriting VTK output...") + + # Velocity magnitude + cfd_python.write_vtk_scalar( + os.path.join(output_dir, "channel_velocity_magnitude.vtk"), + "velocity_magnitude", + vel_mag, + nx, + ny, + xmin, + xmax, + ymin, + ymax, + ) + + # Velocity vectors + cfd_python.write_vtk_vector( + os.path.join(output_dir, "channel_velocity.vtk"), + "velocity", + u, + v, + nx, + ny, + xmin, + xmax, + ymin, + ymax, + ) + + print(" output/channel_velocity_magnitude.vtk") + print(" output/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..74e184c --- /dev/null +++ b/examples/derived_fields.py @@ -0,0 +1,257 @@ +#!/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. + +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 +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) + + # 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 + # ================================================================= + 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 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) + + # 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( + os.path.join(output_dir, "velocity_magnitude.vtk"), + "velocity_magnitude", + vel_mag, + nx, + ny, + grid["xmin"], + grid["xmax"], + grid["ymin"], + grid["ymax"], + ) + print(" Written: output/velocity_magnitude.vtk") + + # Write pressure + cfd_python.write_vtk_scalar( + os.path.join(output_dir, "pressure.vtk"), + "pressure", + p, + nx, + ny, + grid["xmin"], + grid["xmax"], + grid["ymin"], + grid["ymax"], + ) + print(" Written: output/pressure.vtk") + + # Write velocity vectors + cfd_python.write_vtk_vector( + os.path.join(output_dir, "velocity_field.vtk"), + "velocity", + u, + v, + nx, + ny, + grid["xmin"], + grid["xmax"], + grid["ymin"], + grid["ymax"], + ) + print(" Written: output/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..20feb1f --- /dev/null +++ b/examples/error_handling.py @@ -0,0 +1,219 @@ +#!/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 +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", cfd_python.CFDError, Exception), + ("CFDMemoryError", cfd_python.CFDMemoryError, MemoryError), + ("CFDInvalidError", cfd_python.CFDInvalidError, ValueError), + ("CFDIOError", cfd_python.CFDIOError, OSError), + ("CFDUnsupportedError", cfd_python.CFDUnsupportedError, NotImplementedError), + ("CFDDivergedError", cfd_python.CFDDivergedError, cfd_python.CFDError), + ("CFDMaxIterError", cfd_python.CFDMaxIterError, cfd_python.CFDError), + ] + + print(" Exception hierarchy (showing Method Resolution Order):") + for name, exc_class, _ in exceptions: + # 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}") + + # ================================================================= + # 4. Using raise_for_status + # ================================================================= + print("\n4. Using raise_for_status") + print("-" * 60) + + # Test with success + try: + cfd_python.raise_for_status(cfd_python.CFD_SUCCESS, context="test operation") + print(" CFD_SUCCESS: No exception raised (as expected)") + except cfd_python.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: + cfd_python.raise_for_status(code, context=context) + except cfd_python.CFDMemoryError as e: + print(f" {code}: CFDMemoryError - {e}") + except cfd_python.CFDInvalidError as e: + print(f" {code}: CFDInvalidError - {e}") + except cfd_python.CFDIOError as e: + print(f" {code}: CFDIOError - {e}") + except cfd_python.CFDUnsupportedError as e: + print(f" {code}: CFDUnsupportedError - {e}") + except cfd_python.CFDDivergedError as e: + print(f" {code}: CFDDivergedError - {e}") + except cfd_python.CFDMaxIterError as e: + print(f" {code}: CFDMaxIterError - {e}") + except cfd_python.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: + cfd_python.raise_for_status(status, context="simulation") + + return result + + except cfd_python.CFDDivergedError: + print(" Solver diverged - try smaller time step") + return None + except cfd_python.CFDMaxIterError: + print(" Max iterations reached - may need more steps") + return None + except cfd_python.CFDInvalidError as e: + print(f" Invalid parameters: {e}") + return None + except cfd_python.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 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("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("Cannot write to output file", -4) + except OSError as e: + print(f" Caught as OSError: {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..7811feb --- /dev/null +++ b/examples/lid_driven_cavity_advanced.py @@ -0,0 +1,273 @@ +#!/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 + +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 +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 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. + + 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 (used for progress monitoring) + 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 = os.path.join(os.path.dirname(__file__), "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 across different step counts + max_velocities = [] # Max velocity at each sample point + residuals = [] # Relative change between samples + + # 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_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 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=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"steps={sim_steps}") + + except cfd_python.CFDError as e: + print(f"Simulation error at steps={sim_steps}: {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 change from previous sample + if prev_max > 0: + change = abs(max_vel - prev_max) / prev_max + residuals.append(change) + else: + change = float("inf") + + prev_max = max_vel + + print(f"{sim_steps:6d} {stats['max']:12.6f} {stats['avg']:12.6f} {change:12.2e}") + + # 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 + ) + + print("\nSimulation complete!") + print("Output files written to: output/") + + # 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("\nOutput directory: output/") + + +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: 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 Grid Resolution") + print("=" * 60) + + for grid_size in [16, 32, 48]: + result = cfd_python.run_simulation_with_params( + 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" {grid_size}x{grid_size}: max_vel={stats['max']:.4f}, avg_vel={stats['avg']:.4f}") + + +if __name__ == "__main__": + main() 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/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/visualization_matplotlib.py b/examples/visualization_matplotlib.py new file mode 100644 index 0000000..6374cf5 --- /dev/null +++ b/examples/visualization_matplotlib.py @@ -0,0 +1,386 @@ +#!/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 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 + + 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 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 + + +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 + # ================================================================= + 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( + os.path.join(output_dir, "velocity_magnitude_contour.png"), dpi=150, bbox_inches="tight" + ) + plt.close() + print(" Saved: output/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(os.path.join(output_dir, "vortex_visualization.png"), dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: output/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(os.path.join(output_dir, "centerline_profiles.png"), dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: output/centerline_profiles.png") + + # ================================================================= + # 5. Convergence History (Grid Resolution Study) + # ================================================================= + print("\n5. Creating Convergence Analysis Plot") + print("-" * 60) + + # 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 = [] + avg_velocities = [] + + 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=steps, dt=0.001 + ) + stats = cfd_python.calculate_field_stats(result["velocity_magnitude"]) + max_velocities.append(stats["max"]) + avg_velocities.append(stats["avg"]) + + # 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)) + + # Velocity vs steps + ax1 = axes[0] + 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) + + # Relative change (convergence indicator) + ax2 = axes[1] + 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(os.path.join(output_dir, "convergence_history.png"), dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: output/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(os.path.join(output_dir, "grid_comparison.png"), dpi=150, bbox_inches="tight") + plt.close() + print(" Saved: output/grid_comparison.png") + + # ================================================================= + # 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") + + # 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 (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 + # ================================================================= + print("\n" + "=" * 60) + print("Visualization Example Complete!") + print("=" * 60) + 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") + 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() diff --git a/examples/visualization_numpy.py b/examples/visualization_numpy.py index 1e391ef..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 @@ -23,6 +51,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 +156,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 new file mode 100644 index 0000000..869c53c --- /dev/null +++ b/examples/vtk_output.py @@ -0,0 +1,331 @@ +#!/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: +----------------------------------- +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 +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 = os.path.join(os.path.dirname(__file__), "output") + os.makedirs(output_dir, exist_ok=True) + cfd_python.set_output_dir(output_dir) + print(" Output directory: output/") + + # 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( + os.path.join(output_dir, "pressure.vtk"), + "pressure", + fields["pressure"], + nx, + ny, + xmin, + xmax, + ymin, + ymax, + ) + print(" Written: output/pressure.vtk") + + # Write temperature field + cfd_python.write_vtk_scalar( + os.path.join(output_dir, "temperature.vtk"), + "temperature", + fields["temperature"], + nx, + ny, + xmin, + xmax, + ymin, + ymax, + ) + 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( + os.path.join(output_dir, "velocity_magnitude.vtk"), + "velocity_magnitude", + vel_mag, + nx, + ny, + xmin, + xmax, + ymin, + ymax, + ) + print(" Written: output/velocity_magnitude.vtk") + + # ================================================================= + # 3. Writing Vector Fields + # ================================================================= + print("\n3. Writing Vector Fields") + print("-" * 60) + + # Write velocity vector field + cfd_python.write_vtk_vector( + os.path.join(output_dir, "velocity.vtk"), + "velocity", + fields["u"], + fields["v"], + nx, + ny, + xmin, + xmax, + ymin, + ymax, + ) + print(" Written: output/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 = 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: output/vortex_0000.vtk through output/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 using Python + 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: output/velocity_magnitude.csv") + + # Write pressure as CSV + 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: output/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['p']['max']:.4f}") + + # ================================================================= + # Summary + # ================================================================= + print("\n" + "=" * 60) + print("VTK Output Summary") + print("-" * 60) + print(" Output directory: output/") + print("\n Files created:") + print(" Scalar fields:") + print(" - output/pressure.vtk") + print(" - output/temperature.vtk") + print(" - output/velocity_magnitude.vtk") + print(" Vector fields:") + print(" - output/velocity.vtk") + print(" Time series:") + print(f" - output/vortex_0000.vtk to output/vortex_{n_frames-1:04d}.vtk") + print(" CSV files:") + print(" - output/velocity_magnitude.csv") + print(" - output/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() diff --git a/tests/test_errors.py b/tests/test_errors.py index e5e277a..dcf53db 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) - assert isinstance(err, 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, OSError) 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)