"""
Non-Linear Drawdown Diagnostic Test
Investigates why Period 3 pumping shows linear decline instead of accelerating drawdown
"""

import sys
import os

# Add current directory to Python path
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))

from gspy_interface import run_modflow_step, reset_state

def test_pumping_response_curve():
    """Test pumping response to identify linear vs non-linear behavior"""
    
    print("=" * 80)
    print("NON-LINEAR DRAWDOWN DIAGNOSTIC TEST")
    print("Investigating Period 3 pumping response curve")
    print("=" * 80)
    
    # Reset state for clean test
    reset_state()
    
    # Parameters
    hk = 10.0
    sy = 0.3
    
    print(f"Test Parameters: HK={hk}, Sy={sy}")
    print(f"Expected: Non-linear accelerating drawdown during pumping")
    
    # Simulate transition through periods
    results = []
    
    print(f"\n{'Day':<5} {'Period':<8} {'Well Rate':<12} {'Head (m)':<12} {'Δ Head':<12}")
    print("-" * 60)
    
    # Period 1: Day 1 (equilibrium)
    for day in [1]:
        result = run_modflow_step(1, 0.0, hk, sy)
        head = result[0]
        results.append((day, 1, 0.0, head))
        print(f"{day:<5} {'1':<8} {'0.0':<12} {head:<12.6f} {'--':<12}")
    
    # Period 2: Days 2-5 (gradient establishment) 
    for day in range(2, 6):
        result = run_modflow_step(2, 0.0, hk, sy)
        head = result[0]
        prev_head = results[-1][3]
        delta = head - prev_head
        results.append((day, 2, 0.0, head))
        print(f"{day:<5} {'2':<8} {'0.0':<12} {head:<12.6f} {delta:<12.6f}")
    
    # Period 3: Days 6-15 (pumping - should show accelerating decline)
    print("\n--- PERIOD 3 PUMPING ANALYSIS ---")
    print(f"{'Day':<5} {'Rate':<12} {'Head (m)':<12} {'Δ Head':<12} {'Δ²Head':<12}")
    print("-" * 65)
    
    for day in range(6, 16):
        result = run_modflow_step(3, -500.0, hk, sy)
        head = result[0]
        prev_head = results[-1][3]
        delta = head - prev_head
        
        # Calculate second derivative (acceleration in drawdown)
        if len(results) >= 2 and results[-1][1] == 3:  # Previous was also Period 3
            prev_delta = results[-1][3] - results[-2][3]
            delta2 = delta - prev_delta
        else:
            delta2 = 0.0
            
        results.append((day, 3, -500.0, head))
        print(f"{day:<5} {'-500.0':<12} {head:<12.6f} {delta:<12.6f} {delta2:<12.6f}")
    
    print(f"\n" + "=" * 80)
    print("DRAWDOWN ANALYSIS")
    print("=" * 80)
    
    # Analyze pumping period drawdown pattern
    pumping_results = [r for r in results if r[1] == 3]  # Period 3 only
    
    if len(pumping_results) >= 3:
        # Check if drawdown rate is increasing (non-linear)
        deltas = []
        for i in range(1, len(pumping_results)):
            delta = pumping_results[i][3] - pumping_results[i-1][3]
            deltas.append(abs(delta))
        
        print(f"Drawdown rates during pumping: {[f'{d:.6f}' for d in deltas]}")
        
        # Check for acceleration (increasing drawdown rate)
        accelerating = True
        for i in range(1, len(deltas)):
            if deltas[i] <= deltas[i-1]:
                accelerating = False
                break
        
        if accelerating and len(deltas) > 2:
            print(f"✓ ACCELERATING drawdown detected (non-linear)")
            ratio = deltas[-1] / deltas[0] if deltas[0] > 0 else float('inf')
            print(f"  Acceleration ratio: {ratio:.2f}x")
        else:
            print(f"✗ LINEAR or decelerating drawdown (unexpected)")
            print(f"  Expected: increasing drawdown rate over time")
    
    return results

def analyze_time_discretization_issue():
    """Analyze if daily time steps are too large for capturing transient effects"""
    
    print(f"\n" + "=" * 80)
    print("TIME DISCRETIZATION ANALYSIS")
    print("=" * 80)
    
    print("Potential causes of linear response:")
    print("1. Daily time steps may be too large for transient capture")
    print("2. Single time step per day (nstp=1) may miss intra-day dynamics")
    print("3. Storage coefficient may need adjustment")
    print("4. Pumping rate may be too small relative to aquifer properties")
    
    # Calculate time scales
    print(f"\nTime Scale Analysis:")
    
    # Aquifer properties
    hk = 10.0      # m/day
    sy = 0.3       # dimensionless  
    L = 1000.0     # cell size (m)
    
    # Characteristic time for diffusion
    t_char = sy * L**2 / hk  # days
    print(f"  Characteristic diffusion time: {t_char:.1f} days")
    print(f"  Daily time step: 1.0 days")
    print(f"  Ratio (Δt/t_char): {1.0/t_char:.3f}")
    
    if 1.0/t_char > 0.1:
        print(f"  ⚠ Time step may be too large for accurate transient capture")
        print(f"  ⚠ Consider multiple time steps per day (nstp > 1)")
    else:
        print(f"  ✓ Time step appears reasonable for transient effects")
    
    # Storage effects
    total_volume = (10 * 1000) * (10 * 1000) * 25  # L×W×thickness (m³)
    storable_volume = total_volume * sy
    daily_pumping = 500.0  # m³/day
    
    print(f"\nStorage Analysis:")
    print(f"  Total aquifer volume: {total_volume:,.0f} m³")
    print(f"  Storable volume (Sy×Volume): {storable_volume:,.0f} m³")
    print(f"  Daily pumping: {daily_pumping} m³/day")
    print(f"  Fraction pumped per day: {daily_pumping/storable_volume:.6f}")
    
    if daily_pumping/storable_volume < 0.001:
        print(f"  ⚠ Small pumping fraction may cause nearly-linear response")
    else:
        print(f"  ✓ Pumping fraction should cause observable non-linear effects")

def suggest_fixes():
    """Suggest potential fixes for linear response issue"""
    
    print(f"\n" + "=" * 80)
    print("SUGGESTED FIXES")
    print("=" * 80)
    
    print("To achieve non-linear drawdown response:")
    print()
    print("1. INCREASE TIME STEP RESOLUTION:")
    print("   - Modify DIS package: nstp=[10] instead of nstp=[1]")
    print("   - This creates 10 sub-steps per day for better transient capture")
    print()
    print("2. ADJUST AQUIFER THICKNESS:")
    print("   - Reduce thickness from 25m to 10m (more sensitive to pumping)")
    print("   - Or reduce specific yield to increase drawdown sensitivity")
    print()
    print("3. INCREASE PUMPING RATE:")
    print("   - Use -1000 m³/day instead of -500 m³/day")
    print("   - More aggressive pumping = more pronounced non-linear effects")
    print()
    print("4. VERIFY BOUNDARY CONDITIONS:")
    print("   - Ensure GHB conductance isn't too high (dampening response)")
    print("   - Current: 50000 m²/day - may be too high")
    print()
    print("Most likely fix: Increase nstp to capture intra-day transient dynamics")

if __name__ == "__main__":
    results = test_pumping_response_curve()
    analyze_time_discretization_issue()
    suggest_fixes()
    
    print(f"\n🎯 DIAGNOSTIC COMPLETE")
    print("Check results above to identify cause of linear response")