Profiling

conf/config.yaml

@@ -18,7 +18,7 @@ N: 32
 
 # Solver control
 tolerance: 1.0e-6
-max_iterations: 100000
+max_iterations: 500000
 
 # Experiment naming
 experiment_name: LDC-Dev

conf/experiment/validation/convergence/spectral-regu-1000.yaml

@@ -8,7 +8,7 @@ defaults:
 
 # MLflow experiment
 experiment_name: LDC-Convergence
-sweep_name: regu-Re1000
+sweep_name: regu-Re1000-b
 
 # Use Saad/polynomial corner treatment for regularized problem
 solver:
@@ -18,4 +18,4 @@ hydra:
   sweeper:
     params:
       Re: 1000
-      N: 48, 64, 80, 96, 112, 128
+      N: 26, 28, 32, 36, 40, 44, 48, 64, 80, 96, 112, 128

conf/experiment/validation/convergence/spectral-regu-400.yaml

@@ -18,4 +18,4 @@ hydra:
   sweeper:
     params:
       Re: 400
-      N: 24, 32, 40, 48, 56, 64
+      N: 12, 14, 16, 18, 24, 32, 40, 48, 56, 64

conf/experiment/validation/ghia/spectral.yaml

@@ -7,14 +7,14 @@ defaults:
 
 # MLflow experiment
 experiment_name: LDC-Validation
-sweep_name: ghia-Re${Re}
+sweep_name: ghia-HPC
 
 # Validation parameters
 
 hydra:
   sweeper:
     params:
       # Sweep over spectral solver types
-      solver: spectral/fsg #,spectral/fsg
-      Re: 100
-      N: 12, 14, 16, 18, 20
+      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

conf/solver/spectral/fsg.yaml

@@ -12,6 +12,6 @@ name: spectral_fsg
 # FSG Multigrid settings
 multigrid: fsg
 n_levels: 3
-coarse_tolerance_factor: 10.0
+coarse_tolerance_factor: 1.0
 prolongation_method: fft
 restriction_method: fft

conf/solver/spectral/sg.yaml

@@ -17,7 +17,7 @@ max_iterations: ${max_iterations}
 
 # Spectral-specific parameters
 basis_type: chebyshev    # "chebyshev" or "legendre"
-CFL: 2.5
+CFL: 2.0
 beta_squared: 5.0        # artificial compressibility parameter
 
 # Corner singularity treatment

scripts/hpc_submit.py

@@ -118,8 +118,8 @@ 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="6:00", help="Wall time (default: 1:00)")
-    parser.add_argument("--cores", "-n", type=int, default=4, help="Cores per job (default: 4)")
+    parser.add_argument("--time", "-W", default="0:30", 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")
     parser.add_argument("--test-index", type=int, help="Test: show command for specific index")

src/solvers/base.py

@@ -152,6 +152,13 @@ def downsample(data):
             palinstrophy=downsample(palinstrophy_history),
         )
 
+        # Compute vortex metrics (streamfunction-based)
+        try:
+            vortex_metrics = self.compute_vortex_metrics()
+        except Exception as e:
+            log.warning(f"Failed to compute vortex metrics: {e}")
+            vortex_metrics = {}
+
         # Update metrics with convergence info (use FINAL values, not downsampled)
         self.metrics = Metrics(
             iterations=final_iter_count,
@@ -170,6 +177,22 @@ def downsample(data):
             final_palinstrophy=palinstrophy_history[-1]
             if palinstrophy_history
             else 0.0,
+            # Vortex metrics
+            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_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),
+            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),
+            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),
+            psi_TL_x=vortex_metrics.get("psi_TL_x", 0.0),
+            psi_TL_y=vortex_metrics.get("psi_TL_y", 0.0),
         )
 
     def solve(self, tolerance: float = None, max_iter: int = None):
@@ -535,6 +558,205 @@ def compute_global_quantities(self) -> dict:
             "P": self._compute_palinstrophy(),
         }
 
+    # =========================================================================
+    # Vortex Detection (Streamfunction-based, for comparison with Saad/Botella)
+    # =========================================================================
+
+    def _compute_streamfunction(self) -> tuple:
+        """Compute streamfunction ψ by solving ∇²ψ = -ω.
+
+        Uses Poisson solver with ψ=0 on all boundaries (closed cavity).
+
+        Returns
+        -------
+        psi_2d : np.ndarray
+            Streamfunction on 2D grid (ny, nx)
+        x_unique : np.ndarray
+            Unique x coordinates
+        y_unique : np.ndarray
+            Unique y coordinates
+        """
+        from scipy.sparse import diags
+        from scipy.sparse.linalg import spsolve
+
+        # Get vorticity
+        omega = self._compute_vorticity()
+
+        # Get grid info
+        shape = getattr(self, "shape_full", (self.params.nx, self.params.ny))
+        ny, nx = shape
+        dx, dy = self.dx_min, self.dy_min
+
+        # Reshape omega to 2D
+        omega_2d = omega.reshape(shape)
+
+        # Build Laplacian matrix for interior points with Dirichlet BC (ψ=0 on boundaries)
+        # Interior grid is (ny-2) x (nx-2)
+        n_interior = (ny - 2) * (nx - 2)
+
+        # Coefficients for 5-point stencil
+        cx = 1.0 / (dx * dx)
+        cy = 1.0 / (dy * dy)
+        cc = -2.0 * (cx + cy)
+
+        # Build sparse matrix
+        diag_main = np.full(n_interior, cc)
+        diag_x = np.full(n_interior - 1, cx)
+        diag_y = np.full(n_interior - (nx - 2), cy)
+
+        # Remove connections across row boundaries for x-direction
+        for i in range(1, ny - 2):
+            idx = i * (nx - 2) - 1
+            if idx < len(diag_x):
+                diag_x[idx] = 0.0
+
+        A = diags(
+            [diag_y, diag_x, diag_main, diag_x, diag_y],
+            [-(nx - 2), -1, 0, 1, (nx - 2)],
+            format="csr",
+        )
+
+        # RHS: -omega on interior (boundary ψ=0 contributions are zero)
+        rhs = -omega_2d[1:-1, 1:-1].ravel()
+
+        # Solve
+        psi_interior = spsolve(A, rhs)
+
+        # Reconstruct full ψ with zero boundaries
+        psi_2d = np.zeros((ny, nx))
+        psi_2d[1:-1, 1:-1] = psi_interior.reshape(ny - 2, nx - 2)
+
+        # Get coordinates
+        x_unique = np.sort(np.unique(self.fields.x))
+        y_unique = np.sort(np.unique(self.fields.y))
+
+        return psi_2d, x_unique, y_unique
+
+    def _find_primary_vortex(self) -> dict:
+        """Find the primary vortex (global minimum of streamfunction).
+
+        Returns
+        -------
+        dict
+            {'psi_min': float, 'x': float, 'y': float}
+        """
+        psi_2d, x_unique, y_unique = self._compute_streamfunction()
+
+        # Find global minimum
+        min_idx = np.unravel_index(np.argmin(psi_2d), psi_2d.shape)
+        psi_min = psi_2d[min_idx]
+        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)}
+
+    def _find_corner_vortices(self) -> dict:
+        """Find secondary corner vortices (BR, BL, TL).
+
+        Corner vortices have opposite sign to primary vortex:
+        - Primary vortex has ψ < 0 (clockwise rotation)
+        - Secondary vortices have ψ > 0 (counter-clockwise rotation)
+
+        Search regions:
+        - BR (bottom-right): x > 0.5, y < 0.5
+        - BL (bottom-left): x < 0.5, y < 0.5
+        - TL (top-left): x < 0.5, y > 0.5
+
+        Returns
+        -------
+        dict
+            {'BR': {'psi': float, 'x': float, 'y': float}, 'BL': {...}, 'TL': {...}}
+        """
+        psi_2d, x_unique, y_unique = self._compute_streamfunction()
+
+        # Create 2D coordinate arrays
+        X, Y = np.meshgrid(x_unique, y_unique)
+
+        results = {}
+
+        # Define search regions (corners)
+        regions = {
+            "BR": (X > 0.5) & (Y < 0.5),  # Bottom-right
+            "BL": (X < 0.5) & (Y < 0.5),  # Bottom-left
+            "TL": (X < 0.5) & (Y > 0.5),  # Top-left
+        }
+
+        for name, mask in regions.items():
+            # Secondary vortices have ψ > 0 (opposite sign to primary)
+            psi_region = np.where(mask, psi_2d, -np.inf)
+            max_idx = np.unravel_index(np.argmax(psi_region), psi_2d.shape)
+            psi_val = psi_2d[max_idx]
+
+            # Only report if we found a positive ψ (secondary vortex exists)
+            if psi_val > 0:
+                results[name] = {
+                    "psi": float(psi_val),
+                    "x": float(x_unique[max_idx[1]]),
+                    "y": float(y_unique[max_idx[0]]),
+                }
+            else:
+                results[name] = {"psi": 0.0, "x": 0.0, "y": 0.0}
+
+        return results
+
+    def _find_max_vorticity(self) -> dict:
+        """Find maximum vorticity and its location.
+
+        Returns
+        -------
+        dict
+            {'omega_max': float, 'x': float, 'y': float}
+        """
+        omega = self._compute_vorticity()
+
+        # Get grid shape
+        shape = getattr(self, "shape_full", (self.params.nx, self.params.ny))
+        omega_2d = omega.reshape(shape)
+
+        # Find maximum (by absolute value, but track actual sign)
+        max_abs_idx = np.unravel_index(np.argmax(np.abs(omega_2d)), omega_2d.shape)
+        omega_max = omega_2d[max_abs_idx]
+
+        # Get coordinates
+        x_unique = np.sort(np.unique(self.fields.x))
+        y_unique = np.sort(np.unique(self.fields.y))
+
+        return {
+            "omega_max": float(omega_max),
+            "x": float(x_unique[max_abs_idx[1]]),
+            "y": float(y_unique[max_abs_idx[0]]),
+        }
+
+    def compute_vortex_metrics(self) -> dict:
+        """Compute all vortex-related metrics for validation.
+
+        Returns
+        -------
+        dict
+            Dictionary with primary vortex, corner vortices, and max vorticity
+        """
+        primary = self._find_primary_vortex()
+        corners = self._find_corner_vortices()
+        max_omega = self._find_max_vorticity()
+
+        return {
+            "psi_min": primary["psi_min"],
+            "psi_min_x": primary["x"],
+            "psi_min_y": primary["y"],
+            "omega_max": max_omega["omega_max"],
+            "omega_max_x": max_omega["x"],
+            "omega_max_y": max_omega["y"],
+            "psi_BR": corners["BR"]["psi"],
+            "psi_BR_x": corners["BR"]["x"],
+            "psi_BR_y": corners["BR"]["y"],
+            "psi_BL": corners["BL"]["psi"],
+            "psi_BL_x": corners["BL"]["x"],
+            "psi_BL_y": corners["BL"]["y"],
+            "psi_TL": corners["TL"]["psi"],
+            "psi_TL_x": corners["TL"]["x"],
+            "psi_TL_y": corners["TL"]["y"],
+        }
+
     def save_vtk(self, filepath: Path):
         """Save solution to VTS file.
 

src/solvers/datastructures.py

@@ -71,6 +71,28 @@ class Metrics:
     final_enstrophy: float = 0.0
     final_palinstrophy: float = 0.0
 
+    # Primary vortex (global minimum of streamfunction)
+    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
+
+    # Maximum vorticity
+    omega_max: float = 0.0  # Maximum vorticity value
+    omega_max_x: float = 0.0  # x-coordinate of max vorticity
+    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
+    psi_BR: float = 0.0
+    psi_BR_x: float = 0.0
+    psi_BR_y: float = 0.0
+    psi_BL: float = 0.0
+    psi_BL_x: float = 0.0
+    psi_BL_y: float = 0.0
+    psi_TL: float = 0.0
+    psi_TL_x: float = 0.0
+    psi_TL_y: float = 0.0
+
     def to_mlflow(self) -> dict:
         """Convert to MLflow-compatible dict (bools as int, skip inf)."""
         return {