feat: Add debug script for lateral displacement analysis
Adds tests/debug_lateral_discrepancy.py to investigate differences between Zernike OPD lateral displacement reporting and Simcenter post-processing. Key findings documented: - OPD reports sqrt(dx² + dy²) - combined XY magnitude - Simcenter shows individual components (dx or dy) - Both are correct, OPD magnitude is more meaningful for optics 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
This commit is contained in:
181
tests/debug_lateral_discrepancy.py
Normal file
181
tests/debug_lateral_discrepancy.py
Normal file
@@ -0,0 +1,181 @@
|
||||
"""
|
||||
Debug script to investigate the lateral displacement discrepancy between
|
||||
Zernike OPD output and Simcenter post-processing.
|
||||
|
||||
User observation:
|
||||
- Zernike OPD shows max lateral XY displacement: ~0.2238 µm
|
||||
- Simcenter shows XX displacement at 20deg: 7.457e-05 mm = 0.0746 µm
|
||||
|
||||
Hypothesis: The Zernike "max_lateral_disp_um" is sqrt(dx² + dy²), not just dx or dy.
|
||||
"""
|
||||
|
||||
import sys
|
||||
import os
|
||||
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
|
||||
|
||||
import numpy as np
|
||||
from pathlib import Path
|
||||
|
||||
# Find a mirror study OP2 to analyze
|
||||
STUDIES_PATH = Path(r"c:\Users\antoi\Atomizer\studies\M1_Mirror")
|
||||
|
||||
def find_latest_op2():
|
||||
"""Find a recent OP2 file to analyze."""
|
||||
op2_files = list(STUDIES_PATH.rglob("*.op2"))
|
||||
if not op2_files:
|
||||
print("No OP2 files found!")
|
||||
return None
|
||||
# Get most recent
|
||||
return max(op2_files, key=lambda p: p.stat().st_mtime)
|
||||
|
||||
def analyze_displacements(op2_path: Path):
|
||||
"""Analyze displacement components for all subcases."""
|
||||
from pyNastran.op2.op2 import OP2
|
||||
from pyNastran.bdf.bdf import BDF
|
||||
|
||||
print(f"\n{'='*70}")
|
||||
print(f"Analyzing: {op2_path.name}")
|
||||
print(f"Path: {op2_path}")
|
||||
print(f"{'='*70}")
|
||||
|
||||
# Find BDF
|
||||
bdf_path = None
|
||||
for ext in ['.dat', '.bdf']:
|
||||
candidate = op2_path.with_suffix(ext)
|
||||
if candidate.exists():
|
||||
bdf_path = candidate
|
||||
break
|
||||
|
||||
if not bdf_path:
|
||||
print("No BDF found, searching parent...")
|
||||
for f in op2_path.parent.iterdir():
|
||||
if f.suffix.lower() in ['.dat', '.bdf']:
|
||||
bdf_path = f
|
||||
break
|
||||
|
||||
if not bdf_path:
|
||||
print("ERROR: No geometry file found!")
|
||||
return
|
||||
|
||||
print(f"BDF: {bdf_path.name}")
|
||||
|
||||
# Read data
|
||||
print("\nLoading OP2...")
|
||||
op2 = OP2()
|
||||
op2.read_op2(str(op2_path))
|
||||
|
||||
print("Loading BDF...")
|
||||
bdf = BDF()
|
||||
bdf.read_bdf(str(bdf_path))
|
||||
|
||||
node_geo = {int(nid): node.get_position() for nid, node in bdf.nodes.items()}
|
||||
|
||||
SUBCASE_MAP = {'1': 90, '2': 20, '3': 40, '4': 60}
|
||||
|
||||
print(f"\nNode count in BDF: {len(node_geo)}")
|
||||
print(f"\n{'='*70}")
|
||||
print("DISPLACEMENT ANALYSIS BY SUBCASE")
|
||||
print(f"{'='*70}")
|
||||
|
||||
for key, darr in op2.displacements.items():
|
||||
data = darr.data
|
||||
dmat = data[0] if data.ndim == 3 else (data if data.ndim == 2 else None)
|
||||
if dmat is None:
|
||||
continue
|
||||
|
||||
ngt = darr.node_gridtype.astype(int)
|
||||
node_ids = ngt if ngt.ndim == 1 else ngt[:, 0]
|
||||
|
||||
isubcase = getattr(darr, 'isubcase', None)
|
||||
label = str(isubcase) if isubcase else str(key)
|
||||
angle = SUBCASE_MAP.get(label, label)
|
||||
|
||||
print(f"\n--- Subcase {label} ({angle} deg) ---")
|
||||
print(f"Nodes with displacement: {len(node_ids)}")
|
||||
|
||||
# Extract displacement components
|
||||
dx = dmat[:, 0] # X displacement (mm)
|
||||
dy = dmat[:, 1] # Y displacement (mm)
|
||||
dz = dmat[:, 2] # Z displacement (mm)
|
||||
|
||||
# Compute statistics
|
||||
dx_um = dx * 1000.0 # Convert mm to µm
|
||||
dy_um = dy * 1000.0
|
||||
dz_um = dz * 1000.0
|
||||
|
||||
# Lateral magnitude (XY combined)
|
||||
lateral_um = np.sqrt(dx_um**2 + dy_um**2)
|
||||
|
||||
print(f"\nComponent Statistics (in um):")
|
||||
print(f" X displacement (dx):")
|
||||
print(f" Min: {np.min(dx_um):+.4f} um")
|
||||
print(f" Max: {np.max(dx_um):+.4f} um")
|
||||
print(f" RMS: {np.sqrt(np.mean(dx_um**2)):.4f} um")
|
||||
|
||||
print(f" Y displacement (dy):")
|
||||
print(f" Min: {np.min(dy_um):+.4f} um")
|
||||
print(f" Max: {np.max(dy_um):+.4f} um")
|
||||
print(f" RMS: {np.sqrt(np.mean(dy_um**2)):.4f} um")
|
||||
|
||||
print(f" Z displacement (dz):")
|
||||
print(f" Min: {np.min(dz_um):+.4f} um")
|
||||
print(f" Max: {np.max(dz_um):+.4f} um")
|
||||
print(f" RMS: {np.sqrt(np.mean(dz_um**2)):.4f} um")
|
||||
|
||||
print(f"\n Lateral Magnitude (sqrt(dx^2 + dy^2)):")
|
||||
print(f" Max: {np.max(lateral_um):.4f} um <-- This is what Zernike OPD reports!")
|
||||
print(f" RMS: {np.sqrt(np.mean(lateral_um**2)):.4f} um")
|
||||
|
||||
# Find the node with max lateral displacement
|
||||
max_lat_idx = np.argmax(lateral_um)
|
||||
max_lat_nid = int(node_ids[max_lat_idx])
|
||||
max_dx = dx_um[max_lat_idx]
|
||||
max_dy = dy_um[max_lat_idx]
|
||||
|
||||
print(f"\n Node with max lateral displacement: Node {max_lat_nid}")
|
||||
print(f" dx = {max_dx:+.4f} um")
|
||||
print(f" dy = {max_dy:+.4f} um")
|
||||
print(f" sqrt(dx^2+dy^2) = {lateral_um[max_lat_idx]:.4f} um")
|
||||
|
||||
# Compare with just max(|dx|)
|
||||
max_abs_dx = np.max(np.abs(dx_um))
|
||||
max_abs_dy = np.max(np.abs(dy_um))
|
||||
print(f"\n For comparison (what you see in Simcenter):")
|
||||
print(f" max(|dx|) = {max_abs_dx:.4f} um")
|
||||
print(f" max(|dy|) = {max_abs_dy:.4f} um")
|
||||
|
||||
print(f"\n{'='*70}")
|
||||
print("EXPLANATION OF DISCREPANCY")
|
||||
print(f"{'='*70}")
|
||||
print("""
|
||||
The Zernike OPD insight reports "max_lateral_disp_um" as:
|
||||
lateral = sqrt(dx^2 + dy^2) -- combined XY magnitude at each node
|
||||
|
||||
Simcenter's "Displacement - Nodal, X" shows:
|
||||
Just the X component (dx) at each node
|
||||
|
||||
These are different metrics:
|
||||
- If dx_max = 0.0746 um and dy is significant, then:
|
||||
lateral = sqrt(0.0746^2 + dy^2) > 0.0746 um
|
||||
|
||||
To match Simcenter exactly, look at the individual dx/dy/dz stats above.
|
||||
The "max_lateral_um" in Zernike OPD is the MAGNITUDE of the XY vector,
|
||||
not the individual X or Y components.
|
||||
|
||||
For a node where both dx and dy are non-zero:
|
||||
dx = 0.0746 um, dy = 0.21 um
|
||||
lateral = sqrt(0.0746^2 + 0.21^2) = sqrt(0.0056 + 0.044) = sqrt(0.0496) = 0.22 um
|
||||
""")
|
||||
|
||||
if __name__ == '__main__':
|
||||
# Find and analyze an OP2 file
|
||||
if len(sys.argv) > 1:
|
||||
op2_path = Path(sys.argv[1])
|
||||
else:
|
||||
op2_path = find_latest_op2()
|
||||
|
||||
if op2_path and op2_path.exists():
|
||||
analyze_displacements(op2_path)
|
||||
else:
|
||||
print("Please provide an OP2 file path as argument, or place OP2 files in the studies directory.")
|
||||
print(f"\nUsage: python {sys.argv[0]} <path_to_op2_file>")
|
||||
Reference in New Issue
Block a user