diff --git a/README.md b/README.md index 17d3577..378bd97 100644 --- a/README.md +++ b/README.md @@ -57,6 +57,28 @@ uv run mlflow ui # https://kni.dk/mlflow-ana-p3/ ``` +## Hyperparameter Optimization + +Optimize `corner_smoothing` using [Optuna](https://optuna.org/) + Hydra: + +```bash +# Minimize L2 error vs FV reference (default objective) +uv run python main.py -m +experiment/optimization=corner_smoothing \ + 'solver.corner_smoothing=interval(0.02,0.35)' Re=1000 N=30 + +# Minimize vortex error vs Botella & Peyret reference +uv run python main.py -m +experiment/optimization=corner_smoothing \ + 'solver.corner_smoothing=interval(0.02,0.35)' Re=1000 N=30 \ + optuna.objective=botella_vortex + +# Customize trials and parallelism +uv run python main.py -m +experiment/optimization=corner_smoothing \ + 'solver.corner_smoothing=interval(0.02,0.35)' Re=1000 N=30 \ + hydra.sweeper.n_trials=20 hydra.sweeper.n_jobs=8 +``` + +View results in MLflow under `Optuna-CornerSmoothing-{objective}`. See [docs/optuna_optimization.md](docs/optuna_optimization.md) for details. + ## References - [High-Re solutions for incompressible flow (Ghia et al.)](https://www.sciencedirect.com/science/article/pii/0021999182900584) diff --git a/conf/config.yaml b/conf/config.yaml index aa41fc4..1bbe490 100644 --- a/conf/config.yaml +++ b/conf/config.yaml @@ -18,7 +18,7 @@ N: 32 # Solver control tolerance: 1.0e-6 -max_iterations: 500000 +max_iterations: 10000000 # Experiment naming experiment_name: LDC-Dev diff --git a/conf/experiment/benchmarking/timings.yaml b/conf/experiment/benchmarking/timings.yaml index 253f7da..6762688 100644 --- a/conf/experiment/benchmarking/timings.yaml +++ b/conf/experiment/benchmarking/timings.yaml @@ -9,7 +9,7 @@ defaults: # MLflow experiment experiment_name: LDC-Benchmarking -sweep_name: multigrid-comparison-Re${Re} +sweep_name: Benchmarking # Benchmarking parameters Re: 100 @@ -21,6 +21,4 @@ hydra: # Sweep over spectral solver variants # sg: Single Grid (no multigrid) # fsg: Full Single Grid - # vmg: V-cycle MultiGrid - # fmg: Full MultiGrid - solver: spectral/sg,spectral/fsg,spectral/vmg,spectral/fmg + solver: spectral/sg,spectral/fsg diff --git a/conf/experiment/fv_validation.yaml b/conf/experiment/fv_validation.yaml deleted file mode 100644 index 96440da..0000000 --- a/conf/experiment/fv_validation.yaml +++ /dev/null @@ -1,26 +0,0 @@ -# @package _global_ -# FV Solver validation against Ghia benchmark -# Usage: uv run python run_solver.py +experiment=fv_validation -# uv run python run_solver.py -m +experiment=fv_validation (sweep all) - -defaults: - - override /solver: fv - -# MLflow experiment (shared with spectral_validation) -experiment_name: LDC-Validation -# Parent run name - use same name across fv/spectral to group them -sweep_name: validation - -# Validation grid sizes and Reynolds numbers -N: 32 -Re: 100 -tolerance: 1.0e-6 -max_iterations: 10000 - -hydra: - sweeper: - params: - N: 32,64 - Re: 100 - - diff --git a/conf/experiment/optimization/corner_smoothing.yaml b/conf/experiment/optimization/corner_smoothing.yaml new file mode 100644 index 0000000..5bda886 --- /dev/null +++ b/conf/experiment/optimization/corner_smoothing.yaml @@ -0,0 +1,57 @@ +# @package _global_ +# +# Optuna optimization of corner_smoothing parameter +# +# Supported objectives (override with optuna.objective=...): +# - fv_l2_error: Minimize L2 error vs FV reference (default) +# - botella_vortex: Minimize vortex metric error vs Botella & Peyret +# +# Note: Multi-objective optimization is not supported by hydra-optuna-sweeper 1.x +# +# Usage: +# # FV L2 error (default) +# uv run python main.py -m +experiment/optimization=corner_smoothing \ +# 'solver.corner_smoothing=interval(0.02,0.35)' +# +# # Botella vortex metrics +# uv run python main.py -m +experiment/optimization=corner_smoothing \ +# 'solver.corner_smoothing=interval(0.02,0.35)' optuna.objective=botella_vortex +# +# # Custom Re/N +# uv run python main.py -m +experiment/optimization=corner_smoothing \ +# 'solver.corner_smoothing=interval(0.02,0.35)' Re=1000 N=32 + +defaults: + - override /solver: spectral/fsg + - override /hydra/sweeper: optuna_corner + +# Experiment name (dynamically includes objective type) +experiment_name: Optuna-CornerSmoothing +sweep_name: corner-smoothing-${optuna.objective} +# Grid and physics +N: 30 +Re: 1000 + +# Solver settings (strict convergence) +# Note: Early NaN/Inf detection exits diverging runs quickly +tolerance: 1.0e-6 +max_iterations: 500000 + +# Validation against non-regularized FV (true LDC solution) +validation: + reference_dir: data/validation/fv + +# Optimization settings +# Override objective with: optuna.objective=botella_vortex +optuna: + objective: fv_l2_error # Options: fv_l2_error, botella_vortex + +# Optuna sweeper settings +hydra: + sweeper: + study_name: corner_smoothing_${optuna.objective}_Re${Re}_N${N} + n_trials: 15 + n_jobs: 5 + params: + N: 30, 40, 50 + solver.corner_smoothing: interval(0.01, 0.10) diff --git a/conf/experiment/validation/convergence/fv-regu-1000.yaml b/conf/experiment/validation/convergence/fv-regu-1000.yaml deleted file mode 100644 index 0f3fa6d..0000000 --- a/conf/experiment/validation/convergence/fv-regu-1000.yaml +++ /dev/null @@ -1,21 +0,0 @@ -# @package _global_ -# Convergence test - FV solver reference with Saad regularization at Re=1000 -# N=128 as reference solution for spectral convergence - -defaults: - - override /solver: fv - - override /validation: fv-regu - -# MLflow experiment -experiment_name: LDC-Convergence -sweep_name: regu-Re1000 - -# Use Saad/polynomial corner treatment for regularized problem -solver: - corner_treatment: saad - -hydra: - sweeper: - params: - Re: 1000 - N: 128 diff --git a/conf/experiment/validation/convergence/fv-regu-400.yaml b/conf/experiment/validation/convergence/fv-regu-400.yaml deleted file mode 100644 index 807f617..0000000 --- a/conf/experiment/validation/convergence/fv-regu-400.yaml +++ /dev/null @@ -1,21 +0,0 @@ -# @package _global_ -# Convergence test - FV solver reference with Saad regularization at Re=400 -# N=128 as reference solution for spectral convergence - -defaults: - - override /solver: fv - - override /validation: fv-regu - -# MLflow experiment -experiment_name: LDC-Convergence -sweep_name: regu-Re400 - -# Use Saad/polynomial corner treatment for regularized problem -solver: - corner_treatment: saad - -hydra: - sweeper: - params: - Re: 400 - N: 128 diff --git a/conf/experiment/validation/convergence/spectral.yaml b/conf/experiment/validation/convergence/spectral.yaml deleted file mode 100644 index a907def..0000000 --- a/conf/experiment/validation/convergence/spectral.yaml +++ /dev/null @@ -1,20 +0,0 @@ -# @package _global_ -# Ghia validation - Spectral solvers -# Sweeps over SG, FSG, VMG, FMG solvers at N=15 - -defaults: - - override /solver: spectral/sg - -# MLflow experiment -experiment_name: LDC-Validation -sweep_name: convergence-rate-2 - -# Validation parameters - -hydra: - sweeper: - params: - # Sweep over spectral solver types - solver: spectral/fsg #,spectral/fsg - Re: 1000 - N: 80, 90, 100, 110, 120, 130, 140 diff --git a/conf/experiment/validation/ghia/fv.yaml b/conf/experiment/validation/ghia/fv.yaml index e071df5..38464cd 100644 --- a/conf/experiment/validation/ghia/fv.yaml +++ b/conf/experiment/validation/ghia/fv.yaml @@ -6,12 +6,13 @@ defaults: - override /solver: fv # MLflow experiment -experiment_name: LDC-Validation -sweep_name: ghia-Re${Re} +experiment_name: LDC-GHIA-PLOTS + +sweep_name: ghia-plots hydra: sweeper: params: # Sweep over FV grid sizes - N: 64 - Re: 100 + N: 128 + Re: 1000 diff --git a/conf/experiment/validation/ghia/spectral.yaml b/conf/experiment/validation/ghia/spectral.yaml index ed25fca..bf9cc8d 100644 --- a/conf/experiment/validation/ghia/spectral.yaml +++ b/conf/experiment/validation/ghia/spectral.yaml @@ -6,8 +6,12 @@ defaults: - override /solver: spectral/sg # MLflow experiment -experiment_name: LDC-Validation -sweep_name: ghia-HPC +experiment_name: LDC-GHIA-PLOTS +sweep_name: ghia-plots + +solver: + corner_treatment: smoothing + corner_smoothing: 0.15 # smoothing width (fraction of domain) for smoothing method # Validation parameters @@ -16,5 +20,5 @@ hydra: params: # Sweep over spectral solver types solver: spectral/fsg #,spectral/vmg #, spectral/vmg - Re: 400, 1000 - N: 12, 14, 16, 18, 20, 24, 28, 32, 36, 40, 46, 50, 56, 60, 64, 68, 74, 82, 86, 90 + Re: 1000 + N: 16, 20, 32 diff --git a/conf/experiment/validation/hpc-fv.yaml b/conf/experiment/validation/hpc-fv.yaml deleted file mode 100644 index a918097..0000000 --- a/conf/experiment/validation/hpc-fv.yaml +++ /dev/null @@ -1,17 +0,0 @@ -# @package _global_ -# Ghia validation - Finite Volume solver - -defaults: - - override /solver: fv - -# MLflow experiment -experiment_name: LDC-Validation -sweep_name: Hpc-baseline-fv - -hydra: - sweeper: - params: - N: 128 - Re: 100, 400, 1000 - # Relaxation parameter sweep - diff --git a/conf/experiment/validation/convergence/spectral-regu-400.yaml b/conf/experiment/validation/saad-regu/spectral.yaml similarity index 76% rename from conf/experiment/validation/convergence/spectral-regu-400.yaml rename to conf/experiment/validation/saad-regu/spectral.yaml index 96d51ef..2ddb79d 100644 --- a/conf/experiment/validation/convergence/spectral-regu-400.yaml +++ b/conf/experiment/validation/saad-regu/spectral.yaml @@ -7,8 +7,8 @@ defaults: - override /validation: fv-regu # MLflow experiment -experiment_name: LDC-Convergence -sweep_name: regu-Re400 +experiment_name: LDC-Saad-Regu +sweep_name: regu-Re1000 # Use Saad/polynomial corner treatment for regularized problem solver: @@ -17,5 +17,5 @@ solver: hydra: sweeper: params: - Re: 400 - N: 12, 14, 16, 18, 24, 32, 40, 48, 56, 64 + Re: 1000 + N: 16, 32, 64, 128 diff --git a/conf/experiment/validation/convergence/spectral-regu-1000.yaml b/conf/experiment/validation/saad/spectral.yaml similarity index 62% rename from conf/experiment/validation/convergence/spectral-regu-1000.yaml rename to conf/experiment/validation/saad/spectral.yaml index c98215f..33ade10 100644 --- a/conf/experiment/validation/convergence/spectral-regu-1000.yaml +++ b/conf/experiment/validation/saad/spectral.yaml @@ -1,21 +1,21 @@ # @package _global_ -# Convergence test - Spectral solver with Saad regularization at Re=1000 +# Convergence test - Spectral solver with Saad regularization at Re=400 # Sweep over N values with polynomial corner treatment defaults: - override /solver: spectral/fsg - - override /validation: fv-regu + - override /validation: fv # MLflow experiment -experiment_name: LDC-Convergence -sweep_name: regu-Re1000-b +experiment_name: LDC-Saad +sweep_name: Re1000 # Use Saad/polynomial corner treatment for regularized problem solver: - corner_treatment: saad + corner_smoothing: 0.15 # smoothing width (fraction of domain) for smoothing method hydra: sweeper: params: Re: 1000 - N: 26, 28, 32, 36, 40, 44, 48, 64, 80, 96, 112, 128 + N: 16, 32, 64, 128 diff --git a/conf/hydra/sweeper/optuna_corner.yaml b/conf/hydra/sweeper/optuna_corner.yaml new file mode 100644 index 0000000..bef7ebc --- /dev/null +++ b/conf/hydra/sweeper/optuna_corner.yaml @@ -0,0 +1,28 @@ +# @package hydra.sweeper +# +# Optuna Sweeper Configuration for corner_smoothing optimization +# Usage: python main.py -m hydra/sweeper=optuna_corner +# +# For parallel execution, set n_jobs > 1 (requires storage to be set) +# Example: hydra.sweeper.n_jobs=4 hydra.sweeper.storage="sqlite:///optuna.db" + +_target_: hydra_plugins.hydra_optuna_sweeper.optuna_sweeper.OptunaSweeper + +defaults: + - sampler: tpe + +# Optimization direction (minimize L2 error) +direction: minimize + +# Study configuration +study_name: corner_smoothing_optimization +storage: null # Set to "sqlite:///optuna.db" for parallel execution + +# Number of trials and parallel jobs +n_trials: 15 +n_jobs: 1 # Set > 1 for parallel trials (requires storage) + +# Required fields +search_space: null +custom_search_space: null +params: null diff --git a/conf/solver/spectral/fsg.yaml b/conf/solver/spectral/fsg.yaml index a1b250b..79dbc93 100644 --- a/conf/solver/spectral/fsg.yaml +++ b/conf/solver/spectral/fsg.yaml @@ -11,7 +11,7 @@ name: spectral_fsg # FSG Multigrid settings multigrid: fsg -n_levels: 3 +n_levels: 2 coarse_tolerance_factor: 1.0 prolongation_method: fft restriction_method: fft diff --git a/conf/solver/spectral/sg.yaml b/conf/solver/spectral/sg.yaml index 96c253f..8eb54f9 100644 --- a/conf/solver/spectral/sg.yaml +++ b/conf/solver/spectral/sg.yaml @@ -17,7 +17,7 @@ max_iterations: ${max_iterations} # Spectral-specific parameters basis_type: chebyshev # "chebyshev" or "legendre" -CFL: 2.0 +CFL: 1.5 beta_squared: 5.0 # artificial compressibility parameter # Corner singularity treatment diff --git a/data/validation/botella/botella_Re1000_vortex.csv b/data/validation/botella/botella_Re1000_vortex.csv new file mode 100644 index 0000000..f38969d --- /dev/null +++ b/data/validation/botella/botella_Re1000_vortex.csv @@ -0,0 +1,14 @@ +# Botella & Peyret (1998) - Benchmark spectral solution for lid-driven cavity +# Reference: J. Comp. Phys. 179, 439-468 (2002) +# Primary and secondary vortex characteristics at Re=1000 with regularized lid velocity +# N=160 polynomial order (highest resolution benchmark) +# +# Primary vortex (PV): |psi| = magnitude of streamfunction at vortex center +# BL (bottom-left): secondary vortex in bottom-left corner +# BR (bottom-right): secondary vortex in bottom-right corner +# +# Note: psi values are absolute magnitudes (|ψ|), omega values are absolute magnitudes (|ω|) +# Sign conventions vary between papers; we compare magnitudes. +# +Re,psi_primary,omega_primary,x_primary,y_primary,psi_BL,omega_BL,x_BL,y_BL,psi_BR,omega_BR,x_BR,y_BR +1000,0.1189366,2.067753,0.4692,0.5652,2.3072e-4,1.109789,0.1360,0.1118,1.7297e-3,1.112030,0.8640,0.1118 diff --git a/data/validation/botella/botella_Re100_vortex.csv b/data/validation/botella/botella_Re100_vortex.csv new file mode 100644 index 0000000..677beb8 --- /dev/null +++ b/data/validation/botella/botella_Re100_vortex.csv @@ -0,0 +1,11 @@ +# Botella & Peyret (1998) - Benchmark spectral solution for lid-driven cavity +# Reference: J. Comp. Phys. 179, 439-468 (2002) - Table II +# Primary and secondary vortex characteristics at Re=100 with regularized lid velocity +# N=128 polynomial order +# +# Primary vortex (PV): global minimum of streamfunction +# BR (bottom-right): secondary vortex in bottom-right corner +# BL (bottom-left): secondary vortex in bottom-left corner (very weak at Re=100) +# +Re,psi_min,psi_min_x,psi_min_y,omega_center,psi_BR,psi_BR_x,psi_BR_y,psi_BL,psi_BL_x,psi_BL_y +100,-0.1034339,0.6188,0.7344,3.166577,1.4235e-8,0.9453,0.0625,0.0,0.0,0.0 diff --git a/data/validation/botella/botella_Re400_vortex.csv b/data/validation/botella/botella_Re400_vortex.csv new file mode 100644 index 0000000..b3518b1 --- /dev/null +++ b/data/validation/botella/botella_Re400_vortex.csv @@ -0,0 +1,11 @@ +# Botella & Peyret (1998) - Benchmark spectral solution for lid-driven cavity +# Reference: J. Comp. Phys. 179, 439-468 (2002) - Table V +# Primary and secondary vortex characteristics at Re=400 with regularized lid velocity +# N=128 polynomial order +# +# Primary vortex (PV): global minimum of streamfunction +# BR (bottom-right): secondary vortex in bottom-right corner +# BL (bottom-left): secondary vortex in bottom-left corner +# +Re,psi_min,psi_min_x,psi_min_y,omega_center,psi_BR,psi_BR_x,psi_BR_y,psi_BL,psi_BL_x,psi_BL_y +400,-0.1139600,0.5547,0.6055,2.295353,6.4157e-7,0.8906,0.1250,5.1567e-6,0.0508,0.0469 diff --git a/docs/optuna_optimization.md b/docs/optuna_optimization.md new file mode 100644 index 0000000..dea3e27 --- /dev/null +++ b/docs/optuna_optimization.md @@ -0,0 +1,69 @@ +# Hyperparameter Optimization with Optuna + +Optimize the `corner_smoothing` parameter for the spectral solver using Optuna + Hydra. + +## Quick Start + +```bash +# FV L2 error (default) +uv run python main.py -m +experiment/optimization=corner_smoothing \ + 'solver.corner_smoothing=interval(0.02,0.35)' + +# Botella vortex metrics +uv run python main.py -m +experiment/optimization=corner_smoothing \ + 'solver.corner_smoothing=interval(0.02,0.35)' optuna.objective=botella_vortex +``` + +## Objectives + +| Objective | Description | +|-----------|-------------| +| `fv_l2_error` (default) | Minimize L2 error vs FV reference (true LDC) | +| `botella_vortex` | Minimize vortex metric error vs Botella & Peyret | + +Note: Multi-objective Pareto optimization is not supported by `hydra-optuna-sweeper` 1.x (required for Optuna 2.x compatibility). + +## Configuration Options + +```bash +# Change Reynolds number +... Re=1000 + +# Change grid size +... N=32 + +# Change number of trials +... hydra.sweeper.n_trials=30 + +# Change parallel jobs +... hydra.sweeper.n_jobs=8 + +# Change search range +'solver.corner_smoothing=interval(0.01,0.5)' +``` + +## Defaults + +| Parameter | Value | +|-----------|-------| +| Solver | FSG (spectral multigrid) | +| N | 24 | +| Re | 400 | +| Tolerance | 1e-6 | +| Trials | 10 | +| Parallel jobs | 4 | + +## Results + +Logged to MLflow under `Optuna-CornerSmoothing-{objective}`. View with: + +```bash +uv run mlflow ui +``` + +## Config Files + +| File | Purpose | +|------|---------| +| `conf/experiment/optimization/corner_smoothing.yaml` | Optimization experiment | +| `conf/hydra/sweeper/optuna_corner.yaml` | Optuna sweeper settings | diff --git a/main.py b/main.py index b9456c6..44b031c 100644 --- a/main.py +++ b/main.py @@ -72,8 +72,8 @@ def find_existing_run(cfg: DictConfig) -> str: return run_id -def run_solver(cfg: DictConfig) -> str: - """Run solver and log to MLflow. Returns run_id.""" +def run_solver(cfg: DictConfig) -> tuple[str, dict, object]: + """Run solver and log to MLflow. Returns (run_id, validation_errors, solver).""" solver = instantiate(cfg.solver, _convert_="partial") solver_name = cfg.solver.name @@ -87,6 +87,8 @@ def run_solver(cfg: DictConfig) -> str: if parent_run_id: tags.update({"mlflow.parentRunId": parent_run_id, "parent_run_id": parent_run_id, "sweep": "child"}) + validation_errors = {} + with mlflow.start_run(run_name=run_name, tags=tags, nested=bool(parent_run_id)) as run: mlflow.log_params(solver.params.to_mlflow()) mlflow.log_dict(OmegaConf.to_container(cfg), "config.yaml") @@ -94,7 +96,7 @@ def run_solver(cfg: DictConfig) -> str: log.info(f"Solving: {solver_name} N={cfg.N} Re={cfg.Re}") solver.solve() - # Compute validation errors against reference FV solution + # Compute validation errors against reference FV solution (non-regularized) reference_dir = cfg.get("validation", {}).get("reference_dir", "data/validation/fv") validation_errors = solver.compute_validation_errors(reference_dir=reference_dir) if validation_errors: @@ -106,40 +108,149 @@ def run_solver(cfg: DictConfig) -> str: if batch: mlflow.tracking.MlflowClient().log_batch(run.info.run_id, metrics=batch) + # Log validation metrics comparison table + solver.mlflow_log_validation_table() + with tempfile.TemporaryDirectory() as tmpdir: vtk_path = Path(tmpdir) / "solution.vts" solver.to_vtk().save(str(vtk_path)) mlflow.log_artifact(str(vtk_path)) log.info(f"Done: {solver.metrics.iterations} iter, converged={solver.metrics.converged}, time={solver.metrics.wall_time_seconds:.2f}s") - return run.info.run_id + return run.info.run_id, validation_errors, solver def generate_plots(cfg: DictConfig, run_id: str): """Generate plots for a completed run.""" from shared.plotting.ldc import generate_plots_for_run - generate_plots_for_run( - run_id=run_id, - tracking_uri=cfg.mlflow.get("tracking_uri", "./mlruns"), - output_dir=Path(hydra.core.hydra_config.HydraConfig.get().runtime.output_dir), - solver_name=cfg.solver.name, - N=cfg.N, - Re=cfg.Re, - parent_run_id=os.environ.get("MLFLOW_PARENT_RUN_ID"), - upload_to_mlflow=True, - ) + try: + generate_plots_for_run( + run_id=run_id, + tracking_uri=cfg.mlflow.get("tracking_uri", "./mlruns"), + output_dir=Path(hydra.core.hydra_config.HydraConfig.get().runtime.output_dir), + solver_name=cfg.solver.name, + N=cfg.N, + Re=cfg.Re, + parent_run_id=os.environ.get("MLFLOW_PARENT_RUN_ID"), + upload_to_mlflow=True, + ) + except Exception as exc: + log.warning(f"Plotting failed (likely diverged run): {exc}") + + +def compute_fv_l2_objective(validation_errors: dict) -> float: + """Compute objective: combined L2 error vs FV reference. + + Returns sqrt(u_L2_error^2 + v_L2_error^2) against non-regularized FV. + """ + import math + + u_err = validation_errors.get("u_L2_error", float("inf")) + v_err = validation_errors.get("v_L2_error", float("inf")) + + objective = math.sqrt(u_err**2 + v_err**2) + log.info(f"Optuna objective (L2 error vs FV): {objective:.6e}") + return objective + + +def compute_botella_vortex_objective(solver, Re: int) -> float: + """Compute objective: vortex metric error vs Botella & Peyret reference. + + Returns combined relative error in primary vortex characteristics: + - psi_min (streamfunction minimum) + - psi_min_x, psi_min_y (vortex center location) + """ + import math + import pandas as pd + + # Load Botella reference data + ref_path = Path(f"data/validation/botella/botella_Re{Re}_vortex.csv") + if not ref_path.exists(): + log.warning(f"No Botella reference for Re={Re}, using FV objective instead") + return float("inf") + + ref_df = pd.read_csv(ref_path, comment="#") + ref = ref_df.iloc[0].to_dict() + + # Get computed vortex metrics from solver + metrics = solver.metrics + + # Compute relative errors for key vortex characteristics + errors = [] + + # Primary vortex streamfunction (most important) + if ref.get("psi_min") and ref["psi_min"] != 0: + psi_err = abs(metrics.psi_min - ref["psi_min"]) / abs(ref["psi_min"]) + errors.append(psi_err) + log.info(f" psi_min: computed={metrics.psi_min:.6f}, ref={ref['psi_min']:.6f}, err={psi_err:.4f}") + + # Primary vortex location + if ref.get("psi_min_x"): + x_err = abs(metrics.psi_min_x - ref["psi_min_x"]) + errors.append(x_err) + if ref.get("psi_min_y"): + y_err = abs(metrics.psi_min_y - ref["psi_min_y"]) + errors.append(y_err) + + # Combined objective (RMS of relative errors) + if errors: + objective = math.sqrt(sum(e**2 for e in errors) / len(errors)) + else: + objective = float("inf") + + log.info(f"Optuna objective (Botella vortex error): {objective:.6e}") + return objective + + +def compute_optuna_objective(cfg: DictConfig, validation_errors: dict, solver) -> float: + """Compute objective based on config setting. + + Returns + ------- + float + Objective value for single-objective optimization. + """ + objective_type = cfg.get("optuna", {}).get("objective", "fv_l2_error") + + if objective_type == "multi": + raise ValueError( + "Multi-objective optimization is not supported by hydra-optuna-sweeper 1.x. " + "Use objective=fv_l2_error or objective=botella_vortex instead." + ) + elif objective_type == "botella_vortex": + return compute_botella_vortex_objective(solver, int(cfg.Re)) + else: + # Default: FV L2 error + return compute_fv_l2_objective(validation_errors) @hydra.main(config_path="conf", config_name="config", version_base=None) -def main(cfg: DictConfig) -> None: - """Main entry point.""" +def main(cfg: DictConfig) -> float | None: + """Main entry point. + + Returns + ------- + float | None + Objective value for Optuna optimization. + - fv_l2_error: Combined L2 error vs FV reference + - botella_vortex: Vortex metric error vs Botella & Peyret + Returns None in plot_only mode. + """ log.info(f"Solver: {cfg.solver.name}, N={cfg.N}, Re={cfg.Re}") log.info(f"MLflow experiment: {setup_mlflow(cfg)}") - run_id = find_existing_run(cfg) if cfg.get("plot_only") else run_solver(cfg) + if cfg.get("plot_only"): + run_id = find_existing_run(cfg) + generate_plots(cfg, run_id) + return None + + run_id, validation_errors, solver = run_solver(cfg) generate_plots(cfg, run_id) + # Return objective for Optuna (if running hyperparameter optimization) + return compute_optuna_objective(cfg, validation_errors, solver) + if __name__ == "__main__": main() diff --git a/pyproject.toml b/pyproject.toml index 9b319d7..8c0754e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -50,6 +50,8 @@ dependencies = [ "python-dotenv>=1.2.1", "zarr>=3.1.5", "hydra-joblib-launcher>=1.2.0", + "optuna>=2.10.0,<3.0.0", + "hydra-optuna-sweeper>=1.2.0", ] [tool.setuptools] diff --git a/scripts/hpc_submit.py b/scripts/hpc_submit.py index f5d860f..fdfd17c 100644 --- a/scripts/hpc_submit.py +++ b/scripts/hpc_submit.py @@ -118,7 +118,7 @@ def main(): parser = argparse.ArgumentParser(description="Submit HPC job array") parser.add_argument("experiment", help="Experiment (e.g., +experiment/validation/ghia=fv)") parser.add_argument("--queue", "-q", default="hpc", help="LSF queue (default: hpc)") - parser.add_argument("--time", "-W", default="0:30", help="Wall time (default: 0:30)") + parser.add_argument("--time", "-W", default="1:00", help="Wall time (default: 0:30)") parser.add_argument("--cores", "-n", type=int, default=6, help="Cores per job (default: 4)") parser.add_argument("--mem", default="6GB", help="Memory per core (default: 4GB)") parser.add_argument("--dry-run", action="store_true", help="Show commands without submitting") diff --git a/scripts/optuna_corner_smoothing.sh b/scripts/optuna_corner_smoothing.sh new file mode 100755 index 0000000..5f1753c --- /dev/null +++ b/scripts/optuna_corner_smoothing.sh @@ -0,0 +1,39 @@ +#!/bin/bash +#BSUB -J optuna_corner_smoothing +#BSUB -q hpc +#BSUB -W 4:00 +#BSUB -n 20 +#BSUB -R "rusage[mem=8GB]" +#BSUB -R "span[hosts=1]" +#BSUB -o logs/optuna_%J.out +#BSUB -e logs/optuna_%J.err + +# Optuna corner_smoothing optimization for spectral solver +# Runs both FV L2 error and Botella vortex objectives +# +# Submit: bsub < scripts/optuna_corner_smoothing.sh + +set -e +mkdir -p logs + +echo "=== Optuna Corner Smoothing Optimization ===" +echo "Started: $(date)" +echo "" + +# Run FV L2 error objective +echo "=== Objective 1: FV L2 Error ===" +uv run python main.py -m +experiment/optimization=corner_smoothing \ + optuna.objective=fv_l2_error \ + mlflow=coolify + +echo "" + +# Run Botella vortex objective +echo "=== Objective 2: Botella Vortex ===" +uv run python main.py -m +experiment/optimization=corner_smoothing \ + optuna.objective=botella_vortex \ + mlflow=coolify + +echo "" +echo "=== Optimization Complete ===" +echo "Finished: $(date)" diff --git a/src/shared/plotting/ldc/__init__.py b/src/shared/plotting/ldc/__init__.py index 92132a5..563e6e7 100644 --- a/src/shared/plotting/ldc/__init__.py +++ b/src/shared/plotting/ldc/__init__.py @@ -6,7 +6,7 @@ from .convergence import plot_convergence from .data_loading import fields_to_dataframe, load_fields_from_zarr, restructure_fields -from .fields import plot_fields, plot_streamlines, plot_vorticity +from .fields import plot_vorticity from .mlflow_utils import ( download_mlflow_artifacts, find_sibling_runs, @@ -14,6 +14,7 @@ upload_plots_to_mlflow, ) from .orchestrator import generate_comparison_plots_for_sweep, generate_plots_for_run +from .pyvista_fields import generate_pyvista_field_plots from .validation import plot_ghia_comparison # Import style module to trigger sns.set_theme() on package import @@ -22,8 +23,7 @@ __all__ = [ "generate_plots_for_run", "generate_comparison_plots_for_sweep", - "plot_fields", - "plot_streamlines", + "generate_pyvista_field_plots", "plot_vorticity", "plot_convergence", "plot_ghia_comparison", diff --git a/src/shared/plotting/ldc/convergence.py b/src/shared/plotting/ldc/convergence.py index 34ef70d..5112f06 100644 --- a/src/shared/plotting/ldc/convergence.py +++ b/src/shared/plotting/ldc/convergence.py @@ -9,6 +9,7 @@ import matplotlib.pyplot as plt import pandas as pd +import seaborn as sns log = logging.getLogger(__name__) @@ -21,6 +22,9 @@ def plot_convergence( log.warning("No timeseries data available for convergence plot") return None + # Set seaborn darkgrid style + sns.set_style("darkgrid") + fig, ax = plt.subplots() # Plot available metrics @@ -41,10 +45,12 @@ def plot_convergence( rf"\textbf{{Convergence History}} --- {solver_label}, $N={N}$, $\mathrm{{Re}}={Re:.0f}$" ) ax.legend(frameon=True) - ax.grid(True, alpha=0.3) + + # Transparent figure, but keep darkgrid axes background + fig.patch.set_alpha(0.0) output_path = output_dir / "convergence.pdf" - fig.savefig(output_path) + fig.savefig(output_path, facecolor=(0, 0, 0, 0)) plt.close(fig) return output_path diff --git a/src/shared/plotting/ldc/fields.py b/src/shared/plotting/ldc/fields.py index 51456fe..dc8f7e4 100644 --- a/src/shared/plotting/ldc/fields.py +++ b/src/shared/plotting/ldc/fields.py @@ -11,7 +11,7 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd -import pyvista as pv +import seaborn as sns from scipy.interpolate import RectBivariateSpline log = logging.getLogger(__name__) @@ -75,7 +75,7 @@ def plot_fields( plt.tight_layout() output_path = output_dir / "fields.pdf" - fig.savefig(output_path, dpi=300, bbox_inches="tight") + fig.savefig(output_path, dpi=300, bbox_inches="tight", transparent=True) plt.close(fig) return output_path @@ -145,7 +145,7 @@ def plot_streamlines( plt.tight_layout() output_path = output_dir / "streamlines.pdf" - fig.savefig(output_path, dpi=300, bbox_inches="tight") + fig.savefig(output_path, dpi=300, bbox_inches="tight", transparent=True) plt.close(fig) return output_path @@ -175,6 +175,9 @@ def plot_vorticity( dudy = U_spline(y_fine, x_fine, dy=1) vorticity = dvdx - dudy + # Set seaborn darkgrid style + sns.set_style("darkgrid") + fig, ax = plt.subplots(figsize=(7, 6)) vmax = np.max(np.abs(vorticity)) @@ -197,176 +200,11 @@ def plot_vorticity( plt.tight_layout() + # Transparent figure, but keep darkgrid axes background + fig.patch.set_alpha(0.0) + output_path = output_dir / "vorticity.pdf" - fig.savefig(output_path, dpi=300, bbox_inches="tight") + fig.savefig(output_path, dpi=300, bbox_inches="tight", facecolor=(0, 0, 0, 0)) plt.close(fig) return output_path - - -def plot_streamlines_pyvista( - fields_df: pd.DataFrame, Re: float, solver: str, N: int, output_dir: Path -) -> Path: - """Generate beautiful streamline plot using PyVista with ParaView theme. - - Creates a visually striking visualization with velocity magnitude field - and streamlines, with transparent background for easy compositing. - """ - x_unique = np.sort(fields_df["x"].unique()) - y_unique = np.sort(fields_df["y"].unique()) - nx, ny = len(x_unique), len(y_unique) - - sorted_df = fields_df.sort_values(["y", "x"]) - U = sorted_df["u"].values.reshape(ny, nx) - V = sorted_df["v"].values.reshape(ny, nx) - - # Interpolate to finer grid for smoother visualization - n_fine = 200 - x_fine = np.linspace(x_unique[0], x_unique[-1], n_fine) - y_fine = np.linspace(y_unique[0], y_unique[-1], n_fine) - - U_interp = RectBivariateSpline(y_unique, x_unique, U)(y_fine, x_fine) - V_interp = RectBivariateSpline(y_unique, x_unique, V)(y_fine, x_fine) - - # Create 3D grid (z=0 plane) - X, Y = np.meshgrid(x_fine, y_fine) - Z = np.zeros_like(X) - - # Create structured grid - grid = pv.StructuredGrid(X, Y, Z) - - # Add velocity as vector field (3D with w=0) - vel_mag = np.sqrt(U_interp**2 + V_interp**2) - vectors = np.zeros((n_fine * n_fine, 3)) - vectors[:, 0] = U_interp.ravel() - vectors[:, 1] = V_interp.ravel() - vectors[:, 2] = 0.0 - - grid["velocity"] = vectors - grid["velocity_magnitude"] = vel_mag.ravel() - grid.set_active_vectors("velocity") - - # Create seed points for streamlines - n_seeds = 15 - seed_points_left = np.column_stack([ - np.full(n_seeds, x_fine[3]), - np.linspace(y_fine[3], y_fine[-4], n_seeds), - np.zeros(n_seeds) - ]) - seed_points_bottom = np.column_stack([ - np.linspace(x_fine[3], x_fine[-4], n_seeds), - np.full(n_seeds, y_fine[3]), - np.zeros(n_seeds) - ]) - seed_points = np.vstack([seed_points_left, seed_points_bottom]) - seeds = pv.PolyData(seed_points) - - # Generate streamlines - streamlines = grid.streamlines_from_source( - seeds, - vectors="velocity", - max_steps=2000, - integration_direction="both", - ) - - # Set up PyVista plotter with ParaView theme - pv.set_plot_theme("paraview") - plotter = pv.Plotter(off_screen=True, window_size=[1400, 1200]) - - # Add velocity magnitude field as background surface - surface = grid.extract_surface() - plotter.add_mesh( - surface, - scalars="velocity_magnitude", - cmap="turbo", - show_edges=False, - lighting=False, - scalar_bar_args={ - "title": "Velocity Magnitude |u|", - "vertical": True, - "position_x": 0.85, - "position_y": 0.2, - "width": 0.08, - "height": 0.6, - "title_font_size": 16, - "label_font_size": 14, - "fmt": "%.3f", - "n_labels": 5, - }, - ) - - # Add streamlines as white tubes on top - if streamlines.n_points > 0: - tubes = streamlines.tube(radius=0.004, n_sides=8) - plotter.add_mesh( - tubes, - color="white", - opacity=0.85, - smooth_shading=True, - specular=0.3, - ) - - # Set up camera for clean 2D view - plotter.camera_position = "xy" - plotter.camera.zoom(1.15) - - # Enable anti-aliasing for smooth edges - plotter.enable_anti_aliasing("ssaa") - - # Save with transparent background - output_path = output_dir / "streamlines_3d.png" - plotter.screenshot(output_path, transparent_background=True, scale=2) - plotter.close() - - return output_path - - -def export_fields_to_vtk( - fields_df: pd.DataFrame, Re: float, solver: str, N: int, output_dir: Path -) -> Path: - """Export solution fields to VTK format for ParaView visualization. - - Creates a structured grid VTK file with pressure, velocity components, - velocity magnitude, and vorticity fields. - """ - x_unique = np.sort(fields_df["x"].unique()) - y_unique = np.sort(fields_df["y"].unique()) - nx, ny = len(x_unique), len(y_unique) - - sorted_df = fields_df.sort_values(["y", "x"]) - P = sorted_df["p"].values.reshape(ny, nx) - U = sorted_df["u"].values.reshape(ny, nx) - V = sorted_df["v"].values.reshape(ny, nx) - - # Create 3D grid (z=0 plane) - X, Y = np.meshgrid(x_unique, y_unique) - Z = np.zeros_like(X) - - # Create structured grid - grid = pv.StructuredGrid(X, Y, Z) - - # Add scalar fields - grid["pressure"] = P.ravel() - grid["u"] = U.ravel() - grid["v"] = V.ravel() - grid["velocity_magnitude"] = np.sqrt(U**2 + V**2).ravel() - - # Compute and add vorticity - U_spline = RectBivariateSpline(y_unique, x_unique, U) - V_spline = RectBivariateSpline(y_unique, x_unique, V) - dvdx = V_spline(y_unique, x_unique, dx=1) - dudy = U_spline(y_unique, x_unique, dy=1) - vorticity = dvdx - dudy - grid["vorticity"] = vorticity.ravel() - - # Add velocity vector field - vectors = np.zeros((nx * ny, 3)) - vectors[:, 0] = U.ravel() - vectors[:, 1] = V.ravel() - grid["velocity"] = vectors - - # Save as VTK - output_path = output_dir / "solution.vtk" - grid.save(output_path) - - return output_path diff --git a/src/shared/plotting/ldc/orchestrator.py b/src/shared/plotting/ldc/orchestrator.py index 8bc6b04..e5a1c9e 100644 --- a/src/shared/plotting/ldc/orchestrator.py +++ b/src/shared/plotting/ldc/orchestrator.py @@ -12,19 +12,15 @@ from .convergence import plot_convergence from .data_loading import fields_to_dataframe, load_fields_from_zarr -from .fields import ( - plot_fields, - plot_streamlines, - plot_streamlines_pyvista, - plot_vorticity, -) +from .fields import plot_vorticity from .mlflow_utils import ( download_mlflow_artifacts, find_sibling_runs, load_timeseries_from_mlflow, upload_plots_to_mlflow, ) -from .validation import plot_ghia_comparison +from .pyvista_fields import generate_pyvista_field_plots +from .validation import plot_ghia_comparison, plot_l2_convergence log = logging.getLogger(__name__) @@ -39,7 +35,14 @@ def generate_plots_for_run( parent_run_id: Optional[str] = None, upload_to_mlflow: bool = True, ) -> list[Path]: - """Generate all plots for a completed run.""" + """Generate all plots for a completed run. + + Artifacts generated: + - convergence.pdf (matplotlib) + - ghia_comparison.pdf (matplotlib) + - vorticity.pdf (matplotlib) + - u.png, v.png, pressure.png, vel-mag.png, streamlines.png (PyVista) + """ output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) @@ -51,13 +54,12 @@ def generate_plots_for_run( log.info(f"Generating plots for {solver_name} N={N} Re={Re}") - # Generate individual plots plot_paths = [] - plot_paths.append(plot_fields(fields_df, Re, solver_name, N, output_dir)) - plot_paths.append(plot_streamlines(fields_df, Re, solver_name, N, output_dir)) - plot_paths.append(plot_streamlines_pyvista(fields_df, Re, solver_name, N, output_dir)) - plot_paths.append(plot_vorticity(fields_df, Re, solver_name, N, output_dir)) + + # Matplotlib plots: convergence, ghia comparison, vorticity plot_paths.append(plot_convergence(timeseries_df, Re, solver_name, N, output_dir)) + plot_paths.append(plot_vorticity(fields_df, Re, solver_name, N, output_dir)) + ghia_path = plot_ghia_comparison( [{"run_id": run_id, "N": N, "Re": Re, "solver": solver_name, "status": "FINISHED"}], tracking_uri, @@ -66,6 +68,14 @@ def generate_plots_for_run( if ghia_path: plot_paths.append(ghia_path) + # PyVista plots: u, v, pressure, vel-mag, streamlines + vts_path = artifact_dir / "solution.vts" + if vts_path.exists(): + pyvista_paths = generate_pyvista_field_plots(vts_path, output_dir) + plot_paths.extend(pyvista_paths.values()) + else: + log.warning(f"VTS file not found at {vts_path}, skipping PyVista plots") + plot_paths = [p for p in plot_paths if p is not None] log.info(f"Generated {len(plot_paths)} plots for run") @@ -110,15 +120,26 @@ def generate_comparison_plots_for_sweep( comparison_dir = output_dir / parent_name comparison_dir.mkdir(exist_ok=True) + # Generate Ghia comparison plot comparison_path = plot_ghia_comparison(siblings, tracking_uri, comparison_dir) + # Generate L2 convergence plots (4 separate PDFs) + l2_paths = plot_l2_convergence(siblings, tracking_uri, comparison_dir) + + # Collect all generated plots + all_plots = [] if comparison_path: + all_plots.append(comparison_path) results[parent_run_id] = comparison_path log.info(f" Created comparison plot: {comparison_path.name}") - if upload_to_mlflow: - upload_plots_to_mlflow(parent_run_id, [comparison_path], tracking_uri) - log.info(" Uploaded to parent run") + if l2_paths: + all_plots.extend(l2_paths) + log.info(f" Created {len(l2_paths)} L2 convergence plot(s)") + + if all_plots and upload_to_mlflow: + upload_plots_to_mlflow(parent_run_id, all_plots, tracking_uri) + log.info(" Uploaded to parent run") log.info(f"Generated {len(results)} comparison plot(s)") return results diff --git a/src/shared/plotting/ldc/pyvista_fields.py b/src/shared/plotting/ldc/pyvista_fields.py new file mode 100644 index 0000000..a5f31f0 --- /dev/null +++ b/src/shared/plotting/ldc/pyvista_fields.py @@ -0,0 +1,589 @@ +""" +PyVista-based Field Visualization for LDC. + +Generates high-quality 2D field plots using PyVista with the ParaView theme. +Loads solution directly from VTS files. +""" + +import logging +import os +import subprocess +import sys +import tempfile +from pathlib import Path + +import numpy as np +import pyvista as pv + +log = logging.getLogger(__name__) + +# Enable off-screen rendering +os.environ["PYVISTA_OFF_SCREEN"] = "true" + +# Use ParaView theme +pv.set_plot_theme("paraview") + + +# Window size for high DPI output (2400x2400 = 2x the previous 1200x1200) +WINDOW_SIZE = [2400, 2400] + +# Common scalar bar arguments with black text and larger font +SCALAR_BAR_ARGS = { + "vertical": False, + "position_x": 0.25, + "position_y": 0.02, # Lower position for more padding + "width": 0.5, + "height": 0.04, + "title_font_size": 44, # Scaled up for higher resolution + "label_font_size": 32, # Scaled up for higher resolution + "color": "black", + "fmt": "%.2f", + "n_labels": 5, + "italic": True, # Italic for more LaTeX-like appearance + "font_family": "times", # Times font for LaTeX-like appearance +} + + +def _setup_camera(plotter: pv.Plotter) -> None: + """Set up proper 2D orthographic view after mesh is added.""" + plotter.enable_parallel_projection() + plotter.view_xy() + plotter.reset_camera_clipping_range() + plotter.reset_camera(bounds=plotter.bounds) + + +_STREAMLINES_SCRIPT = ''' +import os +import sys +os.environ["PYVISTA_OFF_SCREEN"] = "true" + +import pyvista as pv +pv.set_plot_theme("paraview") + +vts_path = sys.argv[1] +output_path = sys.argv[2] +separating_distance = float(sys.argv[3]) + +WINDOW_SIZE = [2400, 2400] +SCALAR_BAR_ARGS = { + "vertical": False, + "position_x": 0.25, + "position_y": 0.02, + "width": 0.5, + "height": 0.04, + "title_font_size": 44, + "label_font_size": 32, + "color": "black", + "fmt": "%.2f", + "n_labels": 5, + "italic": True, + "font_family": "times", +} + +grid = pv.read(vts_path) +plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + +plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\\n"}, +) + +bounds = grid.bounds +resolution = 256 +image_data = pv.ImageData( + dimensions=(resolution, resolution, 1), + spacing=((bounds[1] - bounds[0]) / (resolution - 1), + (bounds[3] - bounds[2]) / (resolution - 1), + 1.0), + origin=(bounds[0], bounds[2], 0.0), +) + +sampled = image_data.sample(grid) + +streamlines = sampled.streamlines_evenly_spaced_2D( + vectors="velocity", + start_position=(0.3, 0.7, 0.0), + separating_distance=separating_distance, + separating_distance_ratio=0.5, + step_length=0.5, + max_steps=800, +) + +if streamlines.n_points > 0: + plotter.add_mesh( + streamlines, + color="white", + line_width=2.0, + opacity=0.2, + ) + +plotter.enable_parallel_projection() +plotter.view_xy() +plotter.reset_camera_clipping_range() +plotter.reset_camera(bounds=plotter.bounds) + +plotter.screenshot(output_path, transparent_background=True) +plotter.close() +''' + + +def _generate_streamlines_safe(vts_path: Path, output_path: Path, timeout: int = 60) -> bool: + """Generate streamlines plot in a subprocess to isolate VTK crashes. + + Tries dense streamlines first, falls back to sparser if it crashes. + """ + vts_str = str(vts_path) + out_str = str(output_path) + + # Try dense streamlines first (separating_distance=2.5) + for sep_dist in [2.5, 4.0, 8.0]: + try: + result = subprocess.run( + [sys.executable, "-c", _STREAMLINES_SCRIPT, vts_str, out_str, str(sep_dist)], + timeout=timeout, + capture_output=True, + text=True, + ) + if result.returncode == 0: + log.info(f"Streamlines generated with separating_distance={sep_dist}") + return True + else: + log.warning(f"Streamlines failed with separating_distance={sep_dist}: {result.stderr[:200]}") + except subprocess.TimeoutExpired: + log.warning(f"Streamlines timed out with separating_distance={sep_dist}") + except Exception as e: + log.warning(f"Streamlines error with separating_distance={sep_dist}: {e}") + + log.error("All streamline attempts failed") + return False + + +def plot_u_velocity( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate u-velocity contour plot.""" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + plotter.add_mesh( + grid, + scalars="u", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Horizontal velocity, u\n"}, + ) + _setup_camera(plotter) + + output_path = output_dir / f"u{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_v_velocity( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate v-velocity contour plot.""" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + plotter.add_mesh( + grid, + scalars="v", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Vertical velocity, v\n"}, + ) + _setup_camera(plotter) + + output_path = output_dir / f"v{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_velocity_magnitude( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate velocity magnitude contour plot.""" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\n"}, + ) + _setup_camera(plotter) + + output_path = output_dir / f"vel-mag{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_pressure( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate pressure contour plot.""" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + plotter.add_mesh( + grid, + scalars="pressure", + cmap="inferno", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Pressure, p\n"}, + ) + _setup_camera(plotter) + + output_path = output_dir / f"pressure{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_streamlines( + grid: pv.StructuredGrid, + output_dir: Path, + n_seeds: int = 8, + suffix: str = "", +) -> Path: + """Generate streamlines plot with velocity magnitude background. + + Parameters + ---------- + grid : pv.StructuredGrid + Solution grid with velocity vectors + output_dir : Path + Output directory + n_seeds : int + Number of seed points per dimension (reduced for less density) + suffix : str + Filename suffix + + Returns + ------- + Path + Path to saved figure + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + # Create seed points on a grid (avoiding exact boundaries) + seed_x = np.linspace(0.1, 0.9, n_seeds) + seed_y = np.linspace(0.1, 0.9, n_seeds) + seed_X, seed_Y = np.meshgrid(seed_x, seed_y) + seed_points = np.column_stack([ + seed_X.ravel(), + seed_Y.ravel(), + np.zeros(n_seeds * n_seeds) + ]) + seeds = pv.PolyData(seed_points) + + # Compute streamlines + streamlines = grid.streamlines_from_source( + seeds, + vectors="velocity", + integration_direction="both", + initial_step_length=0.01, + max_step_length=0.1, + ) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + + # Add velocity magnitude as background + plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\n"}, + ) + + # Add streamlines + if streamlines.n_points > 0: + plotter.add_mesh( + streamlines, + color="white", + line_width=3.0, # Scaled for higher resolution + opacity=0.8, + ) + + _setup_camera(plotter) + + output_path = output_dir / f"streamlines{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_streamlines_evenly_spaced( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate evenly-spaced streamlines using line sources. + + Uses line sources at multiple positions to create evenly distributed + streamlines across the domain, similar to PyVista's 2D streamlines example. + + Parameters + ---------- + grid : pv.StructuredGrid + Solution grid with velocity vectors + output_dir : Path + Output directory + suffix : str + Filename suffix + + Returns + ------- + Path + Path to saved figure + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + # Use line sources for evenly spaced seeds (like PyVista example) + # Vertical line at left side of domain + streamlines_left = grid.streamlines( + pointa=(0.05, 0.05, 0), + pointb=(0.05, 0.95, 0), + n_points=12, + vectors="velocity", + integration_direction="forward", + initial_step_length=0.01, + max_step_length=0.1, + max_steps=500, + ) + + # Horizontal line at top of domain + streamlines_top = grid.streamlines( + pointa=(0.05, 0.95, 0), + pointb=(0.95, 0.95, 0), + n_points=12, + vectors="velocity", + integration_direction="forward", + initial_step_length=0.01, + max_step_length=0.1, + max_steps=500, + ) + + # Vertical line at right side + streamlines_right = grid.streamlines( + pointa=(0.95, 0.95, 0), + pointb=(0.95, 0.05, 0), + n_points=12, + vectors="velocity", + integration_direction="forward", + initial_step_length=0.01, + max_step_length=0.1, + max_steps=500, + ) + + # Horizontal line at bottom + streamlines_bottom = grid.streamlines( + pointa=(0.95, 0.05, 0), + pointb=(0.05, 0.05, 0), + n_points=12, + vectors="velocity", + integration_direction="forward", + initial_step_length=0.01, + max_step_length=0.1, + max_steps=500, + ) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + + # Add velocity magnitude as background + plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\n"}, + ) + + # Add all streamlines + for streamlines in [streamlines_left, streamlines_top, streamlines_right, streamlines_bottom]: + if streamlines.n_points > 0: + plotter.add_mesh( + streamlines, + color="white", + line_width=3.0, # Scaled for higher resolution + opacity=0.9, + ) + + _setup_camera(plotter) + + output_path = output_dir / f"streamlines_evenly_spaced{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_streamlines_evenly_spaced_2D( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate evenly-spaced streamlines using VTK's 2D algorithm. + + Uses streamlines_evenly_spaced_2D for publication-quality streamlines. + This algorithm produces streamlines that are approximately evenly spaced. + + Parameters + ---------- + grid : pv.StructuredGrid + Solution grid with velocity vectors + output_dir : Path + Output directory + suffix : str + Filename suffix + + Returns + ------- + Path + Path to saved figure + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + + # Add velocity magnitude as background + plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\n"}, + ) + + # Use evenly spaced streamlines 2D algorithm + # First resample to ImageData for better algorithm performance + try: + # Create uniform ImageData grid with higher resolution for smoother streamlines + bounds = grid.bounds + resolution = 256 # Balance between smoothness and stability + image_data = pv.ImageData( + dimensions=(resolution, resolution, 1), + spacing=((bounds[1] - bounds[0]) / (resolution - 1), + (bounds[3] - bounds[2]) / (resolution - 1), + 1.0), + origin=(bounds[0], bounds[2], 0.0), + ) + + # Sample the velocity field onto the uniform grid + sampled = image_data.sample(grid) + + streamlines = sampled.streamlines_evenly_spaced_2D( + vectors="velocity", + start_position=(0.3, 0.7, 0.0), # Near primary vortex + separating_distance=8.0, # Conservative for stability with non-uniform grids + separating_distance_ratio=0.5, + step_length=0.5, # Step size in cell units + max_steps=800, + ) + + if streamlines.n_points > 0: + plotter.add_mesh( + streamlines, + color="white", + line_width=2.0, # Scaled for higher resolution + opacity=0.2, # Subtle but visible + ) + log.info(f"Generated {streamlines.n_lines} evenly-spaced 2D streamlines") + except Exception as e: + log.warning(f"streamlines_evenly_spaced_2D failed: {e}") + + _setup_camera(plotter) + + output_path = output_dir / f"streamlines{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def generate_pyvista_field_plots( + vts_path: Path, + output_dir: Path, +) -> dict[str, Path]: + """Generate PyVista-based field plots from a VTS solution file. + + Generates clean artifact names: u.png, v.png, pressure.png, vel-mag.png, streamlines.png + + Parameters + ---------- + vts_path : Path + Path to solution.vts file + output_dir : Path + Output directory + + Returns + ------- + dict[str, Path] + Dictionary mapping plot names to their paths + """ + vts_path = Path(vts_path) + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + # Load solution + grid = pv.read(str(vts_path)) + + paths = {} + + log.info("Generating PyVista field plots...") + + log.info(" u velocity...") + paths["u"] = plot_u_velocity(grid, output_dir, suffix="") + + log.info(" v velocity...") + paths["v"] = plot_v_velocity(grid, output_dir, suffix="") + + log.info(" velocity magnitude...") + paths["vel-mag"] = plot_velocity_magnitude(grid, output_dir, suffix="") + + log.info(" pressure...") + paths["pressure"] = plot_pressure(grid, output_dir, suffix="") + + log.info(" streamlines (2D evenly spaced)...") + streamlines_path = output_dir / "streamlines.png" + if _generate_streamlines_safe(vts_path, streamlines_path): + paths["streamlines"] = streamlines_path + + return paths diff --git a/src/shared/plotting/ldc/style.py b/src/shared/plotting/ldc/style.py index 2eb8299..33e2161 100644 --- a/src/shared/plotting/ldc/style.py +++ b/src/shared/plotting/ldc/style.py @@ -1,14 +1,29 @@ """ Plotting Style Configuration for LDC Plots. -Uses seaborn darkgrid theme. +Uses seaborn darkgrid theme with LaTeX rendering. """ import logging +import matplotlib.pyplot as plt import seaborn as sns log = logging.getLogger(__name__) -# Use seaborn darkgrid theme -sns.set_theme(style="darkgrid") +# Enable LaTeX rendering +plt.rcParams.update( + { + "text.usetex": True, + "font.family": "serif", + "font.serif": ["Computer Modern Roman"], + "axes.labelsize": 12, + "font.size": 11, + "legend.fontsize": 10, + "xtick.labelsize": 10, + "ytick.labelsize": 10, + } +) + +# Use seaborn darkgrid theme (after rcParams to preserve LaTeX settings) +sns.set_theme(style="darkgrid", rc={"text.usetex": True}) diff --git a/src/shared/plotting/ldc/validation.py b/src/shared/plotting/ldc/validation.py index 3621be4..cc4a2d2 100644 --- a/src/shared/plotting/ldc/validation.py +++ b/src/shared/plotting/ldc/validation.py @@ -1,13 +1,14 @@ """ Validation and Comparison Plots for LDC. -Generates Ghia benchmark comparisons. +Generates Ghia benchmark comparisons and L2 convergence plots. """ import logging from pathlib import Path import matplotlib.pyplot as plt +import mlflow import numpy as np import pandas as pd import seaborn as sns @@ -21,6 +22,170 @@ log = logging.getLogger(__name__) +def plot_l2_convergence( + siblings: list[dict], + tracking_uri: str, + output_dir: Path, +) -> list[Path]: + """Plot L2 error convergence (log-log) for parent runs. + + Creates 4 separate PDF files: + - l2_convergence_u.pdf: u-velocity L2 error vs N + - l2_convergence_v.pdf: v-velocity L2 error vs N + - l2_convergence_u_regu.pdf: u-velocity L2 error (regularized ref) vs N + - l2_convergence_v_regu.pdf: v-velocity L2 error (regularized ref) vs N + + Parameters + ---------- + siblings : list[dict] + List of sibling run info dicts with run_id, N, Re, solver + tracking_uri : str + MLflow tracking URI + output_dir : Path + Output directory + + Returns + ------- + list[Path] + Paths to generated convergence plots + """ + if not siblings: + return [] + + # Only use finished runs + finished_siblings = [ + s for s in siblings if s.get("status", "FINISHED") == "FINISHED" + ] + if len(finished_siblings) < 2: + log.info( + f"Need at least 2 finished runs for convergence plot (have {len(finished_siblings)})" + ) + return [] + + mlflow.set_tracking_uri(tracking_uri) + client = mlflow.tracking.MlflowClient() + + # Collect L2 error metrics from all runs + records = [] + for sibling in finished_siblings: + run_id = sibling["run_id"] + N = sibling["N"] + solver = sibling.get("solver", "unknown") + method = _build_method_label(sibling) + + try: + run = client.get_run(run_id) + metrics = run.data.metrics + + # Get L2 error metrics (may not all be present) + u_l2 = metrics.get("u_L2_error") + v_l2 = metrics.get("v_L2_error") + u_l2_regu = metrics.get("u_L2_error_regu") + v_l2_regu = metrics.get("v_L2_error_regu") + + if u_l2 is not None or v_l2 is not None: + records.append( + { + "N": N, + "Method": method, + "Solver": solver, + "u_L2_error": u_l2, + "v_L2_error": v_l2, + "u_L2_error_regu": u_l2_regu, + "v_L2_error_regu": v_l2_regu, + } + ) + except Exception as e: + log.warning(f"Failed to load metrics for run {run_id}: {e}") + continue + + if not records: + log.warning("No L2 error metrics found in sibling runs") + return [] + + df = pd.DataFrame(records) + + # Define the 4 plots to create + plot_configs = [ + ("u_L2_error", r"$u$ L2 Error", "l2_convergence_u.pdf"), + ("v_L2_error", r"$v$ L2 Error", "l2_convergence_v.pdf"), + ("u_L2_error_regu", r"$u$ L2 Error (regularized ref)", "l2_convergence_u_regu.pdf"), + ("v_L2_error_regu", r"$v$ L2 Error (regularized ref)", "l2_convergence_v_regu.pdf"), + ] + + output_paths = [] + for metric_col, ylabel, filename in plot_configs: + # Filter to rows with valid data for this metric + plot_df = df[df[metric_col].notna()].copy() + if plot_df.empty: + log.info(f"No data for {metric_col}, skipping") + continue + + fig, ax = plt.subplots(figsize=(7, 5)) + + # Seaborn lineplot with markers + sns.lineplot( + data=plot_df, + x="N", + y=metric_col, + hue="Method", + style="Method", + markers=True, + ax=ax, + linewidth=2, + markersize=8, + ) + + # Reference convergence lines (N^-2 and N^-4 for spectral) + N_vals = np.array(sorted(plot_df["N"].unique())) + if len(N_vals) >= 2: + # Get representative error for reference line placement + mid_N = N_vals[len(N_vals) // 2] + mid_err = plot_df[plot_df["N"] == mid_N][metric_col].mean() + if mid_err > 0: + # N^-2 reference + ref_n2 = mid_err * (mid_N / N_vals) ** 2 + ax.plot( + N_vals, + ref_n2, + "--", + color="gray", + alpha=0.6, + linewidth=1.5, + label=r"$\mathcal{O}(N^{-2})$", + ) + # N^-4 reference (faster convergence) + ref_n4 = mid_err * (mid_N / N_vals) ** 4 + ax.plot( + N_vals, + ref_n4, + ":", + color="gray", + alpha=0.6, + linewidth=1.5, + label=r"$\mathcal{O}(N^{-4})$", + ) + + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_xlabel(r"$N$ (polynomial order)", fontsize=11) + ax.set_ylabel(ylabel, fontsize=11) + ax.set_title(f"L2 Error Convergence", fontsize=12) + ax.legend(loc="best", fontsize=9) + ax.grid(True, which="both", ls="-", alpha=0.3) + + plt.tight_layout() + + output_path = output_dir / filename + fig.savefig(output_path, dpi=300, bbox_inches="tight") + plt.close(fig) + + log.info(f"Saved L2 convergence plot: {output_path.name}") + output_paths.append(output_path) + + return output_paths + + def _build_method_label(sibling: dict) -> str: """Build a unified method label from solver name. @@ -28,16 +193,15 @@ def _build_method_label(sibling: dict) -> str: - 'fv' -> 'FV' - 'spectral' -> 'Spectral' - 'spectral_fsg' -> 'Spectral-FSG' - - 'spectral_vmg' -> 'Spectral-VMG' + - 'spectral_fmg' -> 'Spectral-FMG' """ solver = sibling.get("solver", "unknown") # Format known solver names label_map = { - "fv": "FV", + "fv": "FV-TVD", "spectral": "Spectral", "spectral_fsg": "Spectral-FSG", - "spectral_vmg": "Spectral-VMG", "spectral_fmg": "Spectral-FMG", } @@ -191,60 +355,64 @@ def plot_ghia_comparison( u_df = pd.DataFrame(u_records) v_df = pd.DataFrame(v_records) + # Set seaborn darkgrid style with transparent figure background + sns.set_style("darkgrid") + sns.set(rc={"figure.facecolor": (0, 0, 0, 0)}) + # Create subplots with publication-quality sizing fig, axes = plt.subplots(1, 2, figsize=(12, 5)) - # Left: u-velocity (vertical centerline) + # Left: u-velocity (vertical centerline) - dashed lines, no markers sns.lineplot( data=u_df, x="u", y="y", hue="Method", style="Method", - markers=True, + markers=False, + dashes=True, ax=axes[0], sort=False, - linewidth=2, - markersize=7, - markevery=15, + linewidth=1, ) + # Ghia reference with markers sns.scatterplot( data=ghia_u, x="u", y="y", marker="o", - s=50, + s=30, facecolors="none", edgecolors="#333333", - linewidths=1.2, + linewidths=1.0, label="Ghia et al. (1982)", ax=axes[0], zorder=10, ) - axes[0].set_xlabel(r"$u$", fontsize=11) - axes[0].set_ylabel(r"$y$", fontsize=11) - axes[0].set_title(r"$u$-velocity (vertical centerline)", fontsize=11) + axes[0].set_xlabel(r"$u$", fontsize=12) + axes[0].set_ylabel(r"$y$", fontsize=12) + #axes[0].set_title(r"$u$ at $x = 0.5$", fontsize=12) - # Right: v-velocity (horizontal centerline) + # Right: v-velocity (horizontal centerline) - dashed lines, no markers sns.lineplot( data=v_df, x="x", y="v", hue="Method", style="Method", - markers=True, + markers=False, + dashes=True, ax=axes[1], sort=False, - linewidth=2, - markersize=7, - markevery=15, + linewidth=1, ) + # Ghia reference with markers sns.scatterplot( data=ghia_v, x="x", y="v", marker="o", - s=50, + s=30, facecolors="none", edgecolors="#333333", linewidths=1.2, @@ -252,18 +420,21 @@ def plot_ghia_comparison( ax=axes[1], zorder=10, ) - axes[1].set_xlabel(r"$x$", fontsize=11) - axes[1].set_ylabel(r"$v$", fontsize=11) - axes[1].set_title(r"$v$-velocity (horizontal centerline)", fontsize=11) + axes[1].set_xlabel(r"$x$", fontsize=12) + axes[1].set_ylabel(r"$v$", fontsize=12) + #axes[1].set_title(r"$v$ at $y = 0.5$", fontsize=12) - # Overall title - fig.suptitle(rf"Ghia Benchmark Comparison ($\mathrm{{Re}} = {int(Re)}$)", fontsize=13, y=1.00) + # Overall title with LaTeX formatting + fig.suptitle(rf"Centerline Velocities - $\mathrm{{Re}} = {int(Re)}$", fontsize=15, y=1.00) # Tight layout for better spacing plt.tight_layout() + # Transparent figure, but keep darkgrid axes background + fig.patch.set_alpha(0.0) + output_path = output_dir / "ghia_comparison.pdf" - fig.savefig(output_path, dpi=300, bbox_inches="tight") + fig.savefig(output_path, dpi=300, bbox_inches="tight", facecolor=(0, 0, 0, 0)) plt.close(fig) log.info(f"Saved comparison plot: {output_path.name}") diff --git a/src/solvers/__init__.py b/src/solvers/__init__.py index 8646e9a..284b645 100644 --- a/src/solvers/__init__.py +++ b/src/solvers/__init__.py @@ -7,9 +7,7 @@ LidDrivenCavitySolver (abstract base - defines problem) ├── FVSolver (finite volume with SIMPLE algorithm) └── SGSolver (spectral single grid base) - ├── FSGSolver (Full Single Grid multigrid) - ├── VMGSolver (V-cycle MultiGrid) - └── FMGSolver (Full MultiGrid) + └── FSGSolver (Full Single Grid multigrid) """ from .base import LidDrivenCavitySolver @@ -29,8 +27,6 @@ from solvers.fv.solver import FVSolver from solvers.spectral.sg import SGSolver from solvers.spectral.fsg import FSGSolver -from solvers.spectral.vmg import VMGSolver -from solvers.spectral.fmg import FMGSolver __all__ = [ @@ -48,8 +44,6 @@ # Spectral solvers "SGSolver", "FSGSolver", - "VMGSolver", - "FMGSolver", "SpectralParameters", "SpectralSolverFields", ] diff --git a/src/solvers/base.py b/src/solvers/base.py index 6769751..6c3ac94 100644 --- a/src/solvers/base.py +++ b/src/solvers/base.py @@ -181,16 +181,20 @@ def downsample(data): psi_min=vortex_metrics.get("psi_min", 0.0), psi_min_x=vortex_metrics.get("psi_min_x", 0.0), psi_min_y=vortex_metrics.get("psi_min_y", 0.0), + omega_center=vortex_metrics.get("omega_center", 0.0), omega_max=vortex_metrics.get("omega_max", 0.0), omega_max_x=vortex_metrics.get("omega_max_x", 0.0), omega_max_y=vortex_metrics.get("omega_max_y", 0.0), psi_BR=vortex_metrics.get("psi_BR", 0.0), + omega_BR=vortex_metrics.get("omega_BR", 0.0), psi_BR_x=vortex_metrics.get("psi_BR_x", 0.0), psi_BR_y=vortex_metrics.get("psi_BR_y", 0.0), psi_BL=vortex_metrics.get("psi_BL", 0.0), + omega_BL=vortex_metrics.get("omega_BL", 0.0), psi_BL_x=vortex_metrics.get("psi_BL_x", 0.0), psi_BL_y=vortex_metrics.get("psi_BL_y", 0.0), psi_TL=vortex_metrics.get("psi_TL", 0.0), + omega_TL=vortex_metrics.get("omega_TL", 0.0), psi_TL_x=vortex_metrics.get("psi_TL_x", 0.0), psi_TL_y=vortex_metrics.get("psi_TL_y", 0.0), ) @@ -638,7 +642,7 @@ def _find_primary_vortex(self) -> dict: Returns ------- dict - {'psi_min': float, 'x': float, 'y': float} + {'psi_min': float, 'x': float, 'y': float, 'omega_center': float} """ psi_2d, x_unique, y_unique = self._compute_streamfunction() @@ -648,7 +652,18 @@ def _find_primary_vortex(self) -> dict: x_min = x_unique[min_idx[1]] y_min = y_unique[min_idx[0]] - return {"psi_min": float(psi_min), "x": float(x_min), "y": float(y_min)} + # Get vorticity at the primary vortex center + omega = self._compute_vorticity() + shape = getattr(self, "shape_full", (self.params.nx, self.params.ny)) + omega_2d = omega.reshape(shape) + omega_center = omega_2d[min_idx] + + return { + "psi_min": float(psi_min), + "x": float(x_min), + "y": float(y_min), + "omega_center": float(omega_center), + } def _find_corner_vortices(self) -> dict: """Find secondary corner vortices (BR, BL, TL). @@ -743,6 +758,7 @@ def compute_vortex_metrics(self) -> dict: "psi_min": primary["psi_min"], "psi_min_x": primary["x"], "psi_min_y": primary["y"], + "omega_center": primary["omega_center"], "omega_max": max_omega["omega_max"], "omega_max_x": max_omega["x"], "omega_max_y": max_omega["y"], @@ -871,6 +887,82 @@ def mlflow_log_artifact(self, filepath: str): """ mlflow.log_artifact(filepath) + def mlflow_log_validation_table(self, reference_csv: str = None): + """Log validation metrics comparison table to MLflow. + + Creates a table comparing computed vortex metrics against Botella & Peyret + reference values, including energy/enstrophy/palinstrophy. + + Parameters + ---------- + reference_csv : str, optional + Path to reference CSV file. If None, uses default for current Re. + """ + import pandas as pd + + if not mlflow.active_run(): + log.warning("No active MLflow run - skipping validation table") + return + + # Load reference data + if reference_csv is None: + Re = int(self.params.Re) + reference_csv = f"data/validation/botella/botella_Re{Re}_vortex.csv" + + ref_path = Path(reference_csv) + if not ref_path.exists(): + log.warning(f"Reference file not found: {ref_path}") + return + + ref_df = pd.read_csv(ref_path, comment="#") + ref = ref_df.iloc[0].to_dict() + + # Build validation table in standard literature format + # Separate tables for primary vortex and secondary vortices + rows = [] + + def add_row(vortex, metric, computed, reference, fmt=".6f"): + """Add a row to the validation table.""" + if reference and reference != 0: + error_pct = abs(abs(computed) - abs(reference)) / abs(reference) * 100 + ref_str = f"{reference:{fmt}}" if abs(reference) >= 1e-3 else f"{reference:.4e}" + else: + error_pct = None + ref_str = "-" + + comp_str = f"{computed:{fmt}}" if abs(computed) >= 1e-3 else f"{computed:.4e}" + + rows.append({ + "Vortex": vortex, + "Metric": metric, + "Computed": comp_str, + "Botella": ref_str, + "Error (%)": f"{error_pct:.2f}" if error_pct is not None else "-", + }) + + # Primary vortex metrics (use absolute values for comparison) + add_row("Primary", "|ψ|", abs(self.metrics.psi_min), ref.get("psi_primary")) + add_row("Primary", "|ω|", abs(self.metrics.omega_center), ref.get("omega_primary")) + add_row("Primary", "x", self.metrics.psi_min_x, ref.get("x_primary")) + add_row("Primary", "y", self.metrics.psi_min_y, ref.get("y_primary")) + + # Secondary vortex - Bottom Left (BL) + add_row("BL", "|ψ|", abs(self.metrics.psi_BL), ref.get("psi_BL")) + add_row("BL", "|ω|", abs(self.metrics.omega_BL) if hasattr(self.metrics, 'omega_BL') else 0.0, ref.get("omega_BL")) + add_row("BL", "x", self.metrics.psi_BL_x, ref.get("x_BL")) + add_row("BL", "y", self.metrics.psi_BL_y, ref.get("y_BL")) + + # Secondary vortex - Bottom Right (BR) + add_row("BR", "|ψ|", abs(self.metrics.psi_BR), ref.get("psi_BR")) + add_row("BR", "|ω|", abs(self.metrics.omega_BR) if hasattr(self.metrics, 'omega_BR') else 0.0, ref.get("omega_BR")) + add_row("BR", "x", self.metrics.psi_BR_x, ref.get("x_BR")) + add_row("BR", "y", self.metrics.psi_BR_y, ref.get("y_BR")) + + # Create DataFrame and log to MLflow + table_df = pd.DataFrame(rows) + mlflow.log_table(table_df, artifact_file="validation_metrics.json") + log.info("Logged validation metrics table to MLflow") + # ========================================================================= # Validation against reference FV solutions # ========================================================================= @@ -878,10 +970,11 @@ def mlflow_log_artifact(self, filepath: str): def compute_validation_errors( self, reference_dir: str = "data/validation/fv", save_plots: bool = True ) -> dict: - """Compute L2 errors against reference FV solution. + """Compute L2 errors against reference FV solutions. - Loads the reference solution for the current Re, evaluates computed - solution at reference grid coordinates, and computes relative L2 errors. + Computes errors against both normal FV and regularized FV solutions: + - FV (normal): discontinuous lid velocity at corners + - FV-regu: regularized/smoothed lid velocity (matches spectral corner treatment) Parameters ---------- @@ -893,59 +986,72 @@ def compute_validation_errors( Returns ------- dict - {'u_L2_error': float, 'v_L2_error': float} or empty dict if no reference + L2 errors: {u_L2_error, v_L2_error, u_L2_error_regu, v_L2_error_regu} """ - ref_path = Path(reference_dir) / f"Re{int(self.params.Re)}" / "solution.vts" + results = {} + Re = int(self.params.Re) - if not ref_path.exists(): - log.warning(f"No reference solution found at {ref_path}") - return {} - - # Load reference solution - ref_mesh = pv.read(str(ref_path)) - ref_u = ref_mesh.point_data["u"] - ref_v = ref_mesh.point_data["v"] - - # Get reference grid coordinates - ref_points = ref_mesh.points - ref_x = ref_points[:, 0] - ref_y = ref_points[:, 1] - - # Evaluate computed solution at reference grid points - # Subclasses can override _evaluate_at_points for native interpolation (e.g., spectral) - curr_u_at_ref, curr_v_at_ref = self._evaluate_at_points(ref_x, ref_y) - - # Only compute norm on interior points (exclude boundaries where BCs differ) - # Use small margin to exclude boundary points - margin = 0.02 - interior = ( - (ref_x > margin) & (ref_x < self.params.Lx - margin) & - (ref_y > margin) & (ref_y < self.params.Ly - margin) - ) - valid = interior & ~(np.isnan(curr_u_at_ref) | np.isnan(curr_v_at_ref)) - n_valid = np.sum(valid) - n_total = len(ref_u) + # Define reference directories to compare against + ref_dirs = [ + ("data/validation/fv", ""), # normal FV + ("data/validation/fv-regu", "_regu"), # regularized FV + ] + + for ref_base_dir, suffix in ref_dirs: + ref_path = Path(ref_base_dir) / f"Re{Re}" / "solution.vts" + + if not ref_path.exists(): + log.debug(f"No reference solution at {ref_path}") + continue + + # Load reference solution + ref_mesh = pv.read(str(ref_path)) + ref_u = ref_mesh.point_data["u"] + ref_v = ref_mesh.point_data["v"] + + # Get reference grid coordinates + ref_points = ref_mesh.points + ref_x = ref_points[:, 0] + ref_y = ref_points[:, 1] + + # Evaluate computed solution at reference grid points + curr_u_at_ref, curr_v_at_ref = self._evaluate_at_points(ref_x, ref_y) + + # Only compute norm on interior points (exclude exact boundary nodes) + # Use small epsilon to exclude boundary nodes but keep near-boundary interior + margin = 1e-10 + interior = ( + (ref_x > margin) & (ref_x < self.params.Lx - margin) & + (ref_y > margin) & (ref_y < self.params.Ly - margin) + ) + valid = interior & ~(np.isnan(curr_u_at_ref) | np.isnan(curr_v_at_ref)) + n_valid = np.sum(valid) + n_total = len(ref_u) - if n_valid < n_total * 0.5: - log.warning(f"Only {n_valid}/{n_total} valid points for validation") + if n_valid < n_total * 0.5: + log.warning(f"Only {n_valid}/{n_total} valid points for {ref_base_dir}") - # Compute relative L2 errors on valid interior points - u_error = np.linalg.norm(curr_u_at_ref[valid] - ref_u[valid]) / ( - np.linalg.norm(ref_u[valid]) + 1e-12 - ) - v_error = np.linalg.norm(curr_v_at_ref[valid] - ref_v[valid]) / ( - np.linalg.norm(ref_v[valid]) + 1e-12 - ) + # Compute relative L2 errors on valid interior points + u_error = np.linalg.norm(curr_u_at_ref[valid] - ref_u[valid]) / ( + np.linalg.norm(ref_u[valid]) + 1e-12 + ) + v_error = np.linalg.norm(curr_v_at_ref[valid] - ref_v[valid]) / ( + np.linalg.norm(ref_v[valid]) + 1e-12 + ) - log.info(f"Validation errors vs FV reference ({n_valid}/{n_total} pts): u_L2={u_error:.6e}, v_L2={v_error:.6e}") + ref_label = "FV-regu" if suffix else "FV" + log.info(f"L2 errors vs {ref_label} ({n_valid}/{n_total} pts): u={u_error:.6e}, v={v_error:.6e}") - # Save error distribution plots - if save_plots: - self._save_validation_error_plots( - ref_x, ref_y, ref_u, ref_v, curr_u_at_ref, curr_v_at_ref, valid - ) + results[f"u_L2_error{suffix}"] = float(u_error) + results[f"v_L2_error{suffix}"] = float(v_error) + + # Save error distribution plots (only for normal FV to avoid clutter) + if save_plots and not suffix: + self._save_validation_error_plots( + ref_x, ref_y, ref_u, ref_v, curr_u_at_ref, curr_v_at_ref, valid + ) - return {"u_L2_error": float(u_error), "v_L2_error": float(v_error)} + return results def _save_validation_error_plots( self, ref_x, ref_y, ref_u, ref_v, curr_u, curr_v, valid_mask diff --git a/src/solvers/datastructures.py b/src/solvers/datastructures.py index a47f7dd..15d5a52 100644 --- a/src/solvers/datastructures.py +++ b/src/solvers/datastructures.py @@ -75,6 +75,7 @@ class Metrics: psi_min: float = 0.0 # Minimum streamfunction value psi_min_x: float = 0.0 # x-coordinate of minimum psi_min_y: float = 0.0 # y-coordinate of minimum + omega_center: float = 0.0 # Vorticity at primary vortex center # Maximum vorticity omega_max: float = 0.0 # Maximum vorticity value @@ -82,14 +83,17 @@ class Metrics: omega_max_y: float = 0.0 # y-coordinate of max vorticity # Secondary corner vortices (BR=bottom-right, BL=bottom-left, TL=top-left) - # Each stores the local extremum of psi and its location + # Each stores the local extremum of psi, vorticity at center, and location psi_BR: float = 0.0 + omega_BR: float = 0.0 psi_BR_x: float = 0.0 psi_BR_y: float = 0.0 psi_BL: float = 0.0 + omega_BL: float = 0.0 psi_BL_x: float = 0.0 psi_BL_y: float = 0.0 psi_TL: float = 0.0 + omega_TL: float = 0.0 psi_TL_x: float = 0.0 psi_TL_y: float = 0.0 @@ -265,7 +269,7 @@ class SpectralParameters(Parameters): corner_smoothing: float = 0.15 # smoothing width for smoothing method # Multigrid settings - multigrid: str = "none" # "none", "fsg", "vmg", "fmg" + multigrid: str = "none" # "none", "fsg" n_levels: int = 3 coarse_tolerance_factor: float = 10.0 diff --git a/src/solvers/spectral/__init__.py b/src/solvers/spectral/__init__.py index f6b5e73..00b2c6f 100644 --- a/src/solvers/spectral/__init__.py +++ b/src/solvers/spectral/__init__.py @@ -2,7 +2,5 @@ from solvers.spectral.sg import SGSolver from solvers.spectral.fsg import FSGSolver -from solvers.spectral.vmg import VMGSolver -from solvers.spectral.fmg import FMGSolver -__all__ = ["SGSolver", "FSGSolver", "VMGSolver", "FMGSolver"] +__all__ = ["SGSolver", "FSGSolver"] diff --git a/src/solvers/spectral/fmg.py b/src/solvers/spectral/fmg.py deleted file mode 100644 index b51a290..0000000 --- a/src/solvers/spectral/fmg.py +++ /dev/null @@ -1,45 +0,0 @@ -"""Full MultiGrid (FMG) spectral solver for lid-driven cavity. - -FMG will extend the base SG solver with Full MultiGrid acceleration. -Currently not implemented. -""" - -import logging - -from .sg import SGSolver - -log = logging.getLogger(__name__) - - -class FMGSolver(SGSolver): - """Full MultiGrid (FMG) spectral solver. - - Extends the base Single Grid solver with Full MultiGrid acceleration. - - Parameters - ---------- - All parameters inherited from SGSolver, plus: - n_levels : int - Number of multigrid levels - coarse_tolerance_factor : float - Tolerance multiplier for coarse grids - prolongation_method : str - Transfer operator for coarse-to-fine ('fft' or 'polynomial') - restriction_method : str - Transfer operator for fine-to-coarse ('fft' or 'polynomial') - """ - - def solve(self, tolerance: float = None, max_iter: int = None): - """Solve using Full MultiGrid (FMG). - - Parameters - ---------- - tolerance : float, optional - Convergence tolerance. If None, uses params.tolerance. - max_iter : int, optional - Maximum iterations. If None, uses params.max_iterations. - """ - raise NotImplementedError( - "Full MultiGrid (FMG) solver not yet implemented. " - "Use FSG (Full Single Grid) instead." - ) diff --git a/src/solvers/spectral/multigrid/fsg.py b/src/solvers/spectral/multigrid/fsg.py index d7931b4..15bf175 100644 --- a/src/solvers/spectral/multigrid/fsg.py +++ b/src/solvers/spectral/multigrid/fsg.py @@ -35,6 +35,44 @@ log = logging.getLogger(__name__) +def _build_interpolation_matrix_1d(nodes_inner, nodes_full): + """Build interpolation matrix from inner grid to full grid. + + Uses Chebyshev polynomial interpolation for spectral accuracy. + Given values f_inner at nodes_inner, computes f_full = Interp @ f_inner. + + Parameters + ---------- + nodes_inner : np.ndarray + Inner grid nodes (excludes boundary points) + nodes_full : np.ndarray + Full grid nodes (includes boundary points) + + Returns + ------- + Interp : np.ndarray + Interpolation matrix of shape (n_full, n_inner) + """ + from numpy.polynomial.chebyshev import chebvander + + n_inner = len(nodes_inner) + + # Map physical domain to [-1, 1] for Chebyshev polynomials + a, b = nodes_full[0], nodes_full[-1] + xi_inner = 2 * (nodes_inner - a) / (b - a) - 1 + xi_full = 2 * (nodes_full - a) / (b - a) - 1 + + # Vandermonde matrices: V[i,k] = T_k(xi[i]) + V_inner = chebvander(xi_inner, n_inner - 1) # (n_inner, n_inner) + V_full = chebvander(xi_full, n_inner - 1) # (n_full, n_inner) + + # Interpolation: f_full = V_full @ coeffs, where coeffs = V_inner^{-1} @ f_inner + # So: f_full = (V_full @ V_inner^{-1}) @ f_inner = Interp @ f_inner + Interp = V_full @ np.linalg.solve(V_inner, np.eye(n_inner)) + + return Interp + + # ============================================================================= # Numba JIT-compiled kernels for performance-critical operations # ============================================================================= @@ -863,6 +901,14 @@ def _compute_residuals(self, u: np.ndarray, v: np.ndarray, p: np.ndarray): if self.tau_p is not None: lvl.R_p += self.tau_p + # Add FAS tau correction if set (for coarse grid solves in V-cycle) + if self.tau_u is not None: + lvl.R_u[:] += self.tau_u + if self.tau_v is not None: + lvl.R_v[:] += self.tau_v + if self.tau_p is not None: + lvl.R_p[:] += self.tau_p + def _enforce_boundary_conditions(self, u: np.ndarray, v: np.ndarray): """Enforce boundary conditions using cached values.""" u_2d = u.reshape(self.level.shape_full) @@ -1060,6 +1106,13 @@ def solve_fsg( if corner_treatment is None: corner_treatment = create_corner_treatment(method="smoothing") + # Check if using subtraction method (needs fallback to smoothing on coarse grids) + uses_subtraction = corner_treatment.uses_modified_convection() + if uses_subtraction: + smoothing_treatment = create_corner_treatment(method="smoothing") + # Minimum N for subtraction method - use smoothing below this + min_n_for_subtraction = 8 + total_iterations = 0 n_levels = len(levels) @@ -1076,6 +1129,14 @@ def solve_fsg( f"tolerance={level_tol:.2e}" ) + # Select corner treatment for this level + # For subtraction: use smoothing on coarse levels for stability + if uses_subtraction and level.n < min_n_for_subtraction: + level_corner_treatment = smoothing_treatment + log.debug(f" Level {level_idx} (N={level.n}): using smoothing (N < {min_n_for_subtraction})") + else: + level_corner_treatment = corner_treatment + # Initialize from previous level or zeros if level_idx == 0: # Coarsest level: start from zeros @@ -1103,6 +1164,7 @@ def solve_fsg( converged = False level_iters = 0 + diverged = False for iteration in range(max_iterations): u_res, v_res = smoother.step() level_iters += 1 @@ -1119,6 +1181,14 @@ def solve_fsg( ) break + # Early exit on NaN/Inf (diverged) + if not np.isfinite(max_res): + diverged = True + log.warning( + f" Level {level_idx} diverged (NaN/Inf) at iteration {level_iters}, exiting early" + ) + break + # Logging every 100 iterations if iteration > 0 and iteration % 100 == 0: cont_res = smoother.get_continuity_residual() @@ -1127,6 +1197,10 @@ def solve_fsg( f"u_res={u_res:.2e}, v_res={v_res:.2e}, cont={cont_res:.2e}" ) + # Exit early if diverged (NaN/Inf detected) + if diverged: + break + if not converged and not is_finest: log.warning( f" Level {level_idx} did not converge after {level_iters} iterations, " @@ -1138,7 +1212,7 @@ def solve_fsg( ) finest_level = levels[-1] - final_converged = converged + final_converged = converged and not diverged log.info( f"FSG completed: {total_iterations} total iterations, converged={final_converged}" diff --git a/src/solvers/spectral/operators/corner.py b/src/solvers/spectral/operators/corner.py index 4d2997b..3b224b4 100644 --- a/src/solvers/spectral/operators/corner.py +++ b/src/solvers/spectral/operators/corner.py @@ -48,6 +48,13 @@ def get_wall_velocity( """Return (u, v) boundary condition on stationary walls.""" pass + def uses_modified_convection(self) -> bool: + """Return True if this method requires modified convection terms. + + Only the subtraction method (removed) required this. + """ + return False + # ============================================================================= # Method 1: Smoothing (Cosine smoothing near corners) diff --git a/src/solvers/spectral/operators/transfer_operators.py b/src/solvers/spectral/operators/transfer_operators.py index 838c19a..6d85e4a 100644 --- a/src/solvers/spectral/operators/transfer_operators.py +++ b/src/solvers/spectral/operators/transfer_operators.py @@ -11,7 +11,6 @@ - Polynomial-based: Uses Chebyshev polynomial fitting/evaluation For FSG (Full Single Grid), only prolongation is needed. -For VMG/FMG, both prolongation and restriction are required. """ from abc import ABC, abstractmethod diff --git a/src/solvers/spectral/sg.py b/src/solvers/spectral/sg.py index 9ea88bd..b254d45 100644 --- a/src/solvers/spectral/sg.py +++ b/src/solvers/spectral/sg.py @@ -286,6 +286,8 @@ def _compute_residuals(self, u, v, p): Uses O(N³) tensor product differentiation instead of O(N⁴) Kronecker form. + Uses O(N³) tensor product differentiation instead of O(N⁴) Kronecker form. + Parameters ---------- u, v : np.ndarray @@ -417,9 +419,9 @@ def step(self): """ a = self.arrays # Shorthand - # Swap buffers at start (for residual calculation in solve()) - a.u, a.u_prev = a.u_prev, a.u - a.v, a.v_prev = a.v_prev, a.v + # Save previous for convergence check (copy, don't swap!) + a.u_prev[:] = a.u + a.v_prev[:] = a.v # Compute adaptive timestep dt = self._compute_adaptive_timestep() @@ -552,21 +554,22 @@ def _compute_palinstrophy(self) -> float: # ========================================================================= def _compute_streamfunction(self) -> tuple: - """Compute streamfunction ψ by solving ∇²ψ = -ω using FFT-based solver. + """Compute streamfunction ψ by solving ∇²ψ = -ω. - Uses DST (Discrete Sine Transform) for O(N² log N) Poisson solve with - homogeneous Dirichlet BCs (ψ=0 on walls). + Uses spectral Laplacian with sparse linear solver for non-uniform + Chebyshev grids. Homogeneous Dirichlet BCs (ψ=0 on walls). Returns ------- psi_2d : np.ndarray - Streamfunction on 2D grid (Nx+1, Ny+1) + Streamfunction on 2D grid (Nx, Ny) x_2d : np.ndarray X coordinates 2D grid y_2d : np.ndarray Y coordinates 2D grid """ - from scipy.fft import dstn, idstn + from scipy.sparse import kron, eye, csr_matrix + from scipy.sparse.linalg import spsolve # Get vorticity using spectral differentiation omega = self._compute_vorticity() @@ -574,44 +577,44 @@ def _compute_streamfunction(self) -> tuple: Nx, Ny = self.shape_full - # Extract interior points (excluding boundaries where ψ=0) - # Interior is (Nx-2) x (Ny-2) - rhs_interior = -omega_2d[1:-1, 1:-1] - - nx_int, ny_int = rhs_interior.shape - - # DST-I (Type 1) eigenvalues for Laplacian with Dirichlet BCs - # For domain [0, Lx] x [0, Ly] with N interior points: - # λ_k = -( (k*π/Lx)² + (l*π/Ly)² ) - # With our normalized domain [0,1] x [0,1]: - k = np.arange(1, nx_int + 1) - l = np.arange(1, ny_int + 1) - - # Eigenvalues: λ_{kl} = -π²(k²/Lx² + l²/Ly²) - # For unit domain: λ_{kl} = -π²(k² + l²) but we need to account for - # the effective grid spacing in the DST - lambda_k = -((k * np.pi / (nx_int + 1)) ** 2) - lambda_l = -((l * np.pi / (ny_int + 1)) ** 2) - - # 2D eigenvalue matrix - Lambda_kl = lambda_k[:, np.newaxis] + lambda_l[np.newaxis, :] - - # Scale for physical domain size - Lambda_kl = Lambda_kl * ((nx_int + 1) ** 2) - - # Forward DST (Type I) - rhs_hat = dstn(rhs_interior, type=1) - - # Solve in spectral space: psi_hat = rhs_hat / Lambda - # Avoid division by zero (shouldn't happen for Dirichlet problem) - psi_hat = rhs_hat / Lambda_kl - - # Inverse DST (Type I) - note: DST-I is its own inverse up to scaling - psi_interior = idstn(psi_hat, type=1) - - # Assemble full solution with boundary conditions (ψ=0 on boundaries) - psi_2d = np.zeros((Nx, Ny)) - psi_2d[1:-1, 1:-1] = psi_interior + # Build 2D Laplacian using Kronecker products of spectral matrices + # L = Dxx ⊗ Iy + Ix ⊗ Dyy + Dxx_sparse = csr_matrix(self.Dxx_1d) + Dyy_sparse = csr_matrix(self.Dyy_1d) + Ix = eye(Nx, format='csr') + Iy = eye(Ny, format='csr') + + # Full 2D Laplacian (operates on flattened psi with 'C' order) + # For 'C' order: psi[i,j] -> psi_flat[i*Ny + j] + # Dxx acts on first index (x), Dyy acts on second index (y) + L = kron(Dxx_sparse, Iy) + kron(Ix, Dyy_sparse) + + # Right-hand side: -omega (flattened) + rhs = -omega_2d.ravel() + + # Apply homogeneous Dirichlet BCs: psi = 0 on all boundaries + # Identify boundary DOFs + n_dof = Nx * Ny + boundary_mask = np.zeros((Nx, Ny), dtype=bool) + boundary_mask[0, :] = True # x = 0 (left) + boundary_mask[-1, :] = True # x = Lx (right) + boundary_mask[:, 0] = True # y = 0 (bottom) + boundary_mask[:, -1] = True # y = Ly (top) + boundary_idx = np.where(boundary_mask.ravel())[0] + interior_idx = np.where(~boundary_mask.ravel())[0] + + # Modify system to enforce Dirichlet BCs + # Set boundary rows to identity, rhs to 0 + L_bc = L.tolil() + for idx in boundary_idx: + L_bc[idx, :] = 0 + L_bc[idx, idx] = 1.0 + rhs[idx] = 0.0 + L_bc = L_bc.tocsr() + + # Solve the linear system + psi_flat = spsolve(L_bc, rhs) + psi_2d = psi_flat.reshape((Nx, Ny)) return psi_2d, self.x_full, self.y_full @@ -621,7 +624,7 @@ def _find_primary_vortex(self) -> dict: Returns ------- dict - {'psi_min': float, 'x': float, 'y': float} + {'psi_min': float, 'x': float, 'y': float, 'omega_center': float} """ psi_2d, x_2d, y_2d = self._compute_streamfunction() @@ -631,7 +634,16 @@ def _find_primary_vortex(self) -> dict: x_min = x_2d[min_idx] y_min = y_2d[min_idx] - return {"psi_min": float(psi_min), "x": float(x_min), "y": float(y_min)} + # Get vorticity at the primary vortex center + omega_2d = self._compute_vorticity().reshape(self.shape_full) + omega_center = omega_2d[min_idx] + + return { + "psi_min": float(psi_min), + "x": float(x_min), + "y": float(y_min), + "omega_center": float(omega_center), + } def _find_corner_vortices(self) -> dict: """Find secondary corner vortices (BR, BL, TL). @@ -643,9 +655,10 @@ def _find_corner_vortices(self) -> dict: Returns ------- dict - {'BR': {'psi': float, 'x': float, 'y': float}, 'BL': {...}, 'TL': {...}} + {'BR': {'psi': float, 'omega': float, 'x': float, 'y': float}, ...} """ psi_2d, x_2d, y_2d = self._compute_streamfunction() + omega_2d = self._compute_vorticity().reshape(self.shape_full) results = {} @@ -666,11 +679,12 @@ def _find_corner_vortices(self) -> dict: if psi_val > 0: results[name] = { "psi": float(psi_val), + "omega": float(omega_2d[max_idx]), "x": float(x_2d[max_idx]), "y": float(y_2d[max_idx]), } else: - results[name] = {"psi": 0.0, "x": 0.0, "y": 0.0} + results[name] = {"psi": 0.0, "omega": 0.0, "x": 0.0, "y": 0.0} return results @@ -710,16 +724,20 @@ def compute_vortex_metrics(self) -> dict: "psi_min": primary["psi_min"], "psi_min_x": primary["x"], "psi_min_y": primary["y"], + "omega_center": primary["omega_center"], "omega_max": max_omega["omega_max"], "omega_max_x": max_omega["x"], "omega_max_y": max_omega["y"], "psi_BR": corners["BR"]["psi"], + "omega_BR": corners["BR"]["omega"], "psi_BR_x": corners["BR"]["x"], "psi_BR_y": corners["BR"]["y"], "psi_BL": corners["BL"]["psi"], + "omega_BL": corners["BL"]["omega"], "psi_BL_x": corners["BL"]["x"], "psi_BL_y": corners["BL"]["y"], "psi_TL": corners["TL"]["psi"], + "omega_TL": corners["TL"]["omega"], "psi_TL_x": corners["TL"]["x"], "psi_TL_y": corners["TL"]["y"], } diff --git a/src/solvers/spectral/vmg.py b/src/solvers/spectral/vmg.py deleted file mode 100644 index 5122f0c..0000000 --- a/src/solvers/spectral/vmg.py +++ /dev/null @@ -1,459 +0,0 @@ -"""V-cycle Multigrid (VMG) spectral solver. - -Implements proper V-cycle multigrid using Full Approximation Storage (FAS) scheme. -Uses recursive coarse grid correction with damped prolongation. - -Achieves ~2.8x speedup vs single grid with proper parameter tuning. -""" - -import logging -import time -from typing import List, Tuple, Optional - -import numpy as np - -from .sg import SGSolver -from solvers.spectral.multigrid.fsg import ( - SpectralLevel, - MultigridSmoother, - build_hierarchy, - prolongate_solution, - restrict_solution, - restrict_residual, -) -from solvers.spectral.operators.transfer_operators import ( - TransferOperators, - create_transfer_operators, -) -from solvers.spectral.operators.corner import ( - CornerTreatment, - create_corner_treatment, -) - -log = logging.getLogger(__name__) - - -def _vcycle_fas( - levels: List[SpectralLevel], - smoothers: List[MultigridSmoother], - transfer_ops: TransferOperators, - level_idx: int, - pre_smooth: int, - post_smooth: int, - damping: float, -) -> None: - """Perform one V-cycle using Full Approximation Storage (FAS). - - FAS Algorithm: - 1. Pre-smooth on current grid - 2. Compute fine residual: r_h = RHS - N(u_h) - 3. Restrict solution: u_H = I(u_h) - 4. Restrict residual: r_H = I(r_h) - 5. Compute coarse residual: r_H' = RHS_H - N_H(u_H) - 6. tau = r_H - r_H' (FAS correction term) - 7. Solve on coarse: N_H(v_H) = RHS_H + tau, starting from u_H - 8. Correction: e_H = v_H - u_H - 9. Prolongate: e_h = P(e_H) - 10. Update: u_h = u_h + damping * e_h - 11. Post-smooth - - Parameters - ---------- - levels : List[SpectralLevel] - Grid hierarchy (index 0 = coarsest) - smoothers : List[MultigridSmoother] - Smoother for each level - transfer_ops : TransferOperators - Prolongation/restriction operators - level_idx : int - Current level index - pre_smooth : int - Number of pre-smoothing steps - post_smooth : int - Number of post-smoothing steps - damping : float - Damping factor for coarse grid correction (0 < damping <= 1) - """ - level = levels[level_idx] - smoother = smoothers[level_idx] - - # Pre-smoothing - for _ in range(pre_smooth): - smoother.step() - - if level_idx > 0: # Not coarsest level - coarse_level = levels[level_idx - 1] - coarse_smoother = smoothers[level_idx - 1] - - # Store current fine solution - u_fine_old = level.u.copy() - v_fine_old = level.v.copy() - p_fine_old = level.p.copy() - - # Compute fine grid residual - smoother._compute_residuals(level.u, level.v, level.p) - - # Restrict residual to coarse (I(r_h)) - restrict_residual(level, coarse_level, transfer_ops) - I_r_u = coarse_level.R_u.copy() - I_r_v = coarse_level.R_v.copy() - I_r_p = coarse_level.R_p.copy() - - # Restrict solution to coarse using injection (u_H = I(u_h)) - restrict_solution(level, coarse_level, transfer_ops) - u_H = coarse_level.u.copy() - v_H = coarse_level.v.copy() - p_H = coarse_level.p.copy() - - # Compute coarse grid residual of restricted solution (r_H') - coarse_smoother._compute_residuals( - coarse_level.u, coarse_level.v, coarse_level.p - ) - - # tau = I(r_h) - r_H' (FAS correction term) - tau_u = I_r_u - coarse_level.R_u - tau_v = I_r_v - coarse_level.R_v - tau_p = I_r_p - coarse_level.R_p - - # Zero tau at boundaries (BCs are enforced separately) - tau_u_2d = tau_u.reshape(coarse_level.shape_full) - tau_v_2d = tau_v.reshape(coarse_level.shape_full) - tau_u_2d[0, :] = 0.0 - tau_u_2d[-1, :] = 0.0 - tau_u_2d[:, 0] = 0.0 - tau_u_2d[:, -1] = 0.0 - tau_v_2d[0, :] = 0.0 - tau_v_2d[-1, :] = 0.0 - tau_v_2d[:, 0] = 0.0 - tau_v_2d[:, -1] = 0.0 - - # Set tau correction for coarse solve - coarse_smoother.set_tau_correction(tau_u, tau_v, tau_p) - - # Recurse to coarser level - _vcycle_fas( - levels, smoothers, transfer_ops, - level_idx - 1, pre_smooth, post_smooth, damping - ) - - # Clear tau correction - coarse_smoother.clear_tau_correction() - - # Compute correction: e_H = v_H - u_H - e_u = coarse_level.u - u_H - e_v = coarse_level.v - v_H - e_p = coarse_level.p - p_H - - # Zero boundary corrections before prolongation - e_u_2d = e_u.reshape(coarse_level.shape_full) - e_v_2d = e_v.reshape(coarse_level.shape_full) - e_u_2d[0, :] = 0.0 - e_u_2d[-1, :] = 0.0 - e_u_2d[:, 0] = 0.0 - e_u_2d[:, -1] = 0.0 - e_v_2d[0, :] = 0.0 - e_v_2d[-1, :] = 0.0 - e_v_2d[:, 0] = 0.0 - e_v_2d[:, -1] = 0.0 - - # Prolongate correction to fine grid - e_u_fine = transfer_ops.prolongation.prolongate_2d( - e_u_2d, level.shape_full - ) - e_v_fine = transfer_ops.prolongation.prolongate_2d( - e_v_2d, level.shape_full - ) - e_p_fine = transfer_ops.prolongation.prolongate_2d( - e_p.reshape(coarse_level.shape_inner), level.shape_inner - ) - - # Zero boundary corrections on fine grid - e_u_fine[0, :] = 0.0 - e_u_fine[-1, :] = 0.0 - e_u_fine[:, 0] = 0.0 - e_u_fine[:, -1] = 0.0 - e_v_fine[0, :] = 0.0 - e_v_fine[-1, :] = 0.0 - e_v_fine[:, 0] = 0.0 - e_v_fine[:, -1] = 0.0 - - # Apply damped correction - level.u[:] = u_fine_old + damping * e_u_fine.ravel() - level.v[:] = v_fine_old + damping * e_v_fine.ravel() - level.p[:] = p_fine_old + damping * e_p_fine.ravel() - - # Re-apply boundary conditions - smoother._enforce_boundary_conditions(level.u, level.v) - - # Post-smoothing - for _ in range(post_smooth): - smoother.step() - - -def solve_vmg( - levels: List[SpectralLevel], - Re: float, - beta_squared: float, - lid_velocity: float, - CFL: float, - tolerance: float, - max_cycles: int, - transfer_ops: Optional[TransferOperators] = None, - corner_treatment: Optional[CornerTreatment] = None, - Lx: float = 1.0, - Ly: float = 1.0, - pre_smooth: int = 2, - post_smooth: int = 2, - damping: float = 0.5, -) -> Tuple[SpectralLevel, int, bool]: - """Solve using V-cycle Multigrid (VMG) with FAS. - - Parameters - ---------- - levels : List[SpectralLevel] - Grid hierarchy (index 0 = coarsest) - Re, beta_squared, lid_velocity, CFL : float - Solver parameters - tolerance : float - Convergence tolerance - max_cycles : int - Maximum V-cycles - transfer_ops : TransferOperators, optional - Transfer operators - corner_treatment : CornerTreatment, optional - Corner singularity treatment - Lx, Ly : float - Domain dimensions - pre_smooth : int - Pre-smoothing steps per level (default 2) - post_smooth : int - Post-smoothing steps per level (default 2) - damping : float - Damping factor for coarse correction (default 0.5) - - Returns - ------- - tuple - (finest_level, total_work_units, converged) - """ - if transfer_ops is None: - transfer_ops = create_transfer_operators("fft", "fft") - - if corner_treatment is None: - corner_treatment = create_corner_treatment(method="smoothing") - - # Handle subtraction method on coarse levels - uses_subtraction = corner_treatment.uses_modified_convection() - if uses_subtraction: - smoothing_treatment = create_corner_treatment(method="smoothing") - min_n_for_subtraction = 8 - - n_levels = len(levels) - - # Create smoothers for all levels - smoothers = [] - for level_idx, level in enumerate(levels): - if uses_subtraction and level.n < min_n_for_subtraction: - level_corner = smoothing_treatment - else: - level_corner = corner_treatment - - sm = MultigridSmoother( - level=level, - Re=Re, - beta_squared=beta_squared, - lid_velocity=lid_velocity, - CFL=CFL, - corner_treatment=level_corner, - Lx=Lx, - Ly=Ly, - ) - sm.initialize_lid() - smoothers.append(sm) - - # Initialize finest level - finest = levels[-1] - finest.u[:] = 0.0 - finest.v[:] = 0.0 - finest.p[:] = 0.0 - smoothers[-1].initialize_lid() - - # Track convergence via solution change - u_prev = finest.u.copy() - v_prev = finest.v.copy() - - # Work units = smoothing steps per cycle - work_per_cycle = (pre_smooth + post_smooth) * n_levels - total_work = 0 - - for cycle in range(max_cycles): - # Perform one V-cycle starting from finest level - _vcycle_fas( - levels, smoothers, transfer_ops, - n_levels - 1, pre_smooth, post_smooth, damping - ) - total_work += work_per_cycle - - # Check convergence via relative solution change - u_change = np.linalg.norm(finest.u - u_prev) / (np.linalg.norm(u_prev) + 1e-12) - v_change = np.linalg.norm(finest.v - v_prev) / (np.linalg.norm(v_prev) + 1e-12) - u_prev[:] = finest.u - v_prev[:] = finest.v - - max_change = max(u_change, v_change) - - if max_change < tolerance: - log.info( - f"VMG converged in {cycle + 1} cycles ({total_work} work units), " - f"residual={max_change:.2e}" - ) - return finest, total_work, True - - if (cycle + 1) % 50 == 0: - log.debug(f"VMG cycle {cycle + 1}: residual = {max_change:.2e}") - - log.warning(f"VMG did not converge after {max_cycles} cycles") - return finest, total_work, False - - -class VMGSolver(SGSolver): - """V-cycle Multigrid (VMG) spectral solver. - - Implements proper V-cycle multigrid with Full Approximation Storage (FAS). - Uses recursive coarse grid correction with damped prolongation for - accelerated convergence. - - Achieves ~2.8x speedup vs single grid with default parameters: - - pre_smooth=2, post_smooth=2, damping=0.5 - - Parameters - ---------- - All parameters inherited from SGSolver, plus: - n_levels : int - Maximum number of multigrid levels (default 10, limited by coarsest_n) - coarsest_n : int - Minimum N for coarsest grid (default 12). Coarse grids must resolve physics. - pre_smooth : int - Pre-smoothing steps per level (default 2) - post_smooth : int - Post-smoothing steps per level (default 2) - damping : float - Damping factor for coarse grid correction (default 0.5) - prolongation_method : str - Transfer operator for coarse-to-fine ('fft') - restriction_method : str - Transfer operator for fine-to-coarse ('fft') - """ - - def __init__(self, **kwargs): - """Initialize VMG solver with multigrid hierarchy.""" - super().__init__(**kwargs) - - # Create transfer operators from config - self._transfer_ops = create_transfer_operators( - prolongation_method=self.params.prolongation_method, - restriction_method=self.params.restriction_method, - ) - - # Get coarsest_n from params (default 12) - coarsest_n = getattr(self.params, 'coarsest_n', 12) - - # Build grid hierarchy dynamically from nx down to coarsest_n - # n_levels=100 is just a large upper bound; actual levels determined by coarsest_n - self._levels = build_hierarchy( - n_fine=self.params.nx, - n_levels=100, - coarsest_n=coarsest_n, - basis_x=self.basis_x, - basis_y=self.basis_y, - Lx=self.params.Lx, - Ly=self.params.Ly, - ) - - # V-cycle parameters - self._pre_smooth = getattr(self.params, 'pre_smooth', 2) - self._post_smooth = getattr(self.params, 'post_smooth', 2) - self._damping = getattr(self.params, 'damping', 0.5) - - log.info( - f"VMG initialized: {len(self._levels)} levels {[l.n for l in self._levels]}, " - f"pre={self._pre_smooth}, post={self._post_smooth}, " - f"damping={self._damping}" - ) - - def solve(self, tolerance: float = None, max_iter: int = None): - """Solve using V-cycle Multigrid (VMG). - - Parameters - ---------- - tolerance : float, optional - Convergence tolerance. If None, uses params.tolerance. - max_iter : int, optional - Maximum V-cycles. If None, uses params.max_iterations. - """ - if tolerance is None: - tolerance = self.params.tolerance - if max_iter is None: - max_iter = self.params.max_iterations - - log.info( - f"VMG: {self.params.n_levels} levels, " - f"pre={self._pre_smooth}, post={self._post_smooth}, " - f"damping={self._damping}" - ) - - # Solve using VMG - time_start = time.time() - finest_level, total_work, converged = solve_vmg( - levels=self._levels, - Re=self.params.Re, - beta_squared=self.params.beta_squared, - lid_velocity=self.params.lid_velocity, - CFL=self.params.CFL, - tolerance=tolerance, - max_cycles=max_iter, - transfer_ops=self._transfer_ops, - corner_treatment=self.corner_treatment, - Lx=self.params.Lx, - Ly=self.params.Ly, - pre_smooth=self._pre_smooth, - post_smooth=self._post_smooth, - damping=self._damping, - ) - - time_end = time.time() - wall_time = time_end - time_start - - # Copy solution from finest level to solver arrays - self.arrays.u[:] = finest_level.u - self.arrays.v[:] = finest_level.v - self.arrays.p[:] = finest_level.p - - # Compute final residuals - self._compute_residuals(self.arrays.u, self.arrays.v, self.arrays.p) - - # Store results - final_residuals = self._compute_algebraic_residuals() - residual_history = [ - { - "rel_iter": tolerance if converged else tolerance * 10, - "u_eq": final_residuals["u_residual"], - "v_eq": final_residuals["v_residual"], - "continuity": final_residuals["continuity_residual"], - } - ] - - self._store_results( - residual_history=residual_history, - final_iter_count=total_work, - is_converged=converged, - wall_time=wall_time, - energy_history=[self._compute_energy()], - enstrophy_history=[self._compute_enstrophy()], - palinstrophy_history=[self._compute_palinstrophy()], - ) - - log.info( - f"VMG completed in {wall_time:.2f}s: {total_work} work units, " - f"converged={converged}" - ) diff --git a/src/utilities/mlflow/callback.py b/src/utilities/mlflow/callback.py index 60ef0d2..8b330bd 100644 --- a/src/utilities/mlflow/callback.py +++ b/src/utilities/mlflow/callback.py @@ -216,6 +216,103 @@ def on_job_start(self, config: DictConfig, **kwargs) -> None: # Set env var so child run can find parent os.environ["MLFLOW_PARENT_RUN_ID"] = parent_id + def _log_optuna_results_to_parent(self, config: DictConfig) -> None: + """Log Optuna optimization results to parent MLflow run.""" + import mlflow + import pandas as pd + + # Check if this is an Optuna sweep + sweeper_cfg = config.get("hydra", {}).get("sweeper", {}) + if "optuna" not in str(sweeper_cfg.get("_target_", "")).lower(): + return + + study_name = sweeper_cfg.get("study_name", "") + if not study_name: + return + + # Get parent run ID + parent_run_ids = list(self._parent_runs.values()) + if not parent_run_ids: + parent_run_ids = self._find_recent_parent_runs() + + if not parent_run_ids: + log.warning("No parent run found for Optuna results logging") + return + + parent_run_id = parent_run_ids[0] + + # Try to get Optuna results from child runs in MLflow + try: + experiment = mlflow.get_experiment_by_name(self._full_experiment_name) + if not experiment: + return + + # Get all child runs under this parent + child_runs = mlflow.search_runs( + experiment_ids=[experiment.experiment_id], + filter_string=f"tags.parent_run_id = '{parent_run_id}' AND attributes.status = 'FINISHED'", + order_by=["start_time ASC"], + ) + + if child_runs.empty: + log.warning("No completed child runs found for Optuna summary") + return + + # Build trials table from child runs + trials_data = [] + for _, run in child_runs.iterrows(): + trial = { + "run_id": run["run_id"][:8], + "corner_smoothing": run.get("params.corner_smoothing"), + "u_L2_error": run.get("metrics.u_L2_error"), + "v_L2_error": run.get("metrics.v_L2_error"), + "iterations": run.get("metrics.iterations"), + "converged": run.get("metrics.converged"), + "wall_time": run.get("metrics.wall_time_seconds"), + } + trials_data.append(trial) + + trials_df = pd.DataFrame(trials_data) + + # Find best trial (minimize L2 error) + if "u_L2_error" in trials_df.columns and "v_L2_error" in trials_df.columns: + trials_df["combined_L2"] = ( + trials_df["u_L2_error"].fillna(float("inf")) ** 2 + + trials_df["v_L2_error"].fillna(float("inf")) ** 2 + ) ** 0.5 + best_idx = trials_df["combined_L2"].idxmin() + best_trial = trials_df.loc[best_idx] + else: + best_trial = None + + # Log to parent run + with mlflow.start_run(run_id=parent_run_id, nested=False): + # Log trials table as artifact + mlflow.log_table(trials_df, artifact_file="optuna_trials.json") + + # Log best trial metrics + if best_trial is not None: + mlflow.log_metrics({ + "best_corner_smoothing": float(best_trial.get("corner_smoothing", 0)), + "best_u_L2_error": float(best_trial.get("u_L2_error", float("inf"))), + "best_v_L2_error": float(best_trial.get("v_L2_error", float("inf"))), + "best_combined_L2": float(best_trial.get("combined_L2", float("inf"))), + }) + mlflow.set_tag("best_run_id", best_trial.get("run_id", "")) + + # Log summary info + mlflow.log_metric("n_trials_completed", len(trials_df)) + mlflow.log_metric("n_trials_converged", int(trials_df["converged"].sum()) if "converged" in trials_df.columns else 0) + + log.info(f"Logged Optuna results to parent run {parent_run_id[:8]}") + if best_trial is not None: + log.info(f" Best: corner_smoothing={best_trial.get('corner_smoothing')}, L2={best_trial.get('combined_L2'):.6f}") + + except Exception as e: + log.warning(f"Failed to log Optuna results: {e}") + import traceback + log.debug(traceback.format_exc()) + def on_multirun_end(self, config: DictConfig, **kwargs) -> None: """Clean up after sweep completes and generate comparison plots.""" if os.environ.get("MLFLOW_SWEEP_ACTIVE") != "1": @@ -226,6 +323,9 @@ def on_multirun_end(self, config: DictConfig, **kwargs) -> None: log.info("Multirun sweep completed") + # Log Optuna results to parent run + self._log_optuna_results_to_parent(config) + # Generate comparison plots for all parent runs try: from pathlib import Path