#%%
# -*- coding: utf-8 -*-
"""
Automated Verification Script for Shallow unconfined flow with rainfall

@author: Jan Machacek

OVERVIEW:
This script provides an automated workflow to verify a numerical simulation 
against the analytical solution for a 1D uncoupled transient infiltration 
problem, based on the work of Srivastava & Yeh (1991).

WORKFLOW:
1.  Run Simulation: For each specified case, it executes the `numgeo` 
    numerical solver in a temporary directory.
2.  Process Data: It reads and assembles the simulation output (pore water 
    pressure at specific nodes over time) into a pandas DataFrame.
3.  Quantify Error: It calculates the Root Mean Squared Error (RMSE) between 
    the a reference solution and the simulation result to provide a quantitative measure of accuracy.
4.  Summarize Results: It generates a final summary table in Markdown format 
    comparing the RMSE across all test cases.

DEPENDENCIES:
- Python Libraries: shutil, subprocess, sys, os, numpy, pandas, matplotlib
- External Software: `numgeo` (must be in the system's PATH).
- LaTeX Distribution: Required for the plot's text rendering (`text.usetex`).

SETUP:
- `case_dirs`: A list of directories containing `numgeo` input files.
- `sol.txt`: A text file containing the roots (eigenvalues) for the 
  analytical solution.
"""
import shutil, subprocess, sys, os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib import cm

fontsize = 10
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
# Note: 'text.usetex' requires a LaTeX installation on your system.
plt.rcParams['text.usetex'] = True
plt.rcParams['font.size'] = fontsize
plt.rcParams['axes.spines.right'] = False
plt.rcParams['axes.spines.top'] = False

#
# Set up
#

# A list of the case directories to be processed.
case_dirs = ["p8", 
             "p6", 
             "u8p8-stab",
             "u6p6-stab",
             "u8p4",
             "u6p3",
             "u10p10-stap"]
files_left = ["PRINT-OUT_node_1489.dat",
              "PRINT-OUT_node_2153.dat",
              "PRINT-OUT_node_1489.dat",
              "PRINT-OUT_node_2153.dat",
              "PRINT-OUT_node_1489.dat", 
              "PRINT-OUT_node_2153.dat",
              "PRINT-OUT_node_2496.dat"]
files_right = ["PRINT-OUT_node_3956.dat",
               "PRINT-OUT_node_5627.dat",
               "PRINT-OUT_node_3956.dat",
               "PRINT-OUT_node_5627.dat",
               "PRINT-OUT_node_3956.dat",
               "PRINT-OUT_node_5627.dat",
               "PRINT-OUT_node_6343.dat"]

# reference values taken from p8 simulation
TARGET_LEFT = -0.810697E+01
TARGET_RIGHT = -0.595250E+01

# List of errors
ERR_LEFT = [] ; ERR_RIGHT = []

#
# Main Processing Loop
#

for src, lfile, rfile in zip(case_dirs,files_left,files_right):
    print(f"\n--- Processing Case: {src} ---")
    
    # Set up a temporary directory for the simulation run.

    dst = f"{src}-tmp"
    if os.path.exists(dst): shutil.rmtree(dst)
    shutil.copytree(src, dst)

    # 1. Run the numerical simulation.

    print(f"[{src}] Running simulation...")
    try:
        subprocess.run(
            ["numgeo", f"inp={src}-calculation", "out=sim", "ncpus=2"],
            cwd=dst, check=True, capture_output=True, text=True
        )
        print(f"[{src}] Simulation finished successfully.")
    except subprocess.CalledProcessError as e:
        print(f"[{src}] SIMULATION FAILED. Aborting.")
        print(f"Error: {e.stderr}")
        sys.exit(1)

    # 2. Read and assemble simulation data from output files.

    dat = np.genfromtxt(src+'-tmp/sim-print-out/'+lfile, skip_header=1)
    err = (1. - dat[1]/TARGET_LEFT) * 100.
    ERR_LEFT.append(err)

    dat = np.genfromtxt(src+'-tmp/sim-print-out/'+rfile, skip_header=1)
    err = (1. - dat[1]/TARGET_RIGHT) * 100.
    ERR_RIGHT.append(err)

#
# Final Report Generation
#

print(f"\n--- Summary ---")

labels = [r'$\Delta p^w_{left}$ (%)', r'$\Delta p^w_{right}$ (%)']
rows   = [ERR_LEFT, ERR_RIGHT]

# Format numbers: up to 3 decimals, strip trailing zeros
def fmt(x):
    s = f"{float(x):.3f}"
    return s.rstrip("0").rstrip(".")

# Build cells (first column = labels, then one column per case)
cols  = [""] + case_dirs
data  = [[labels[i]] + [fmt(v) for v in rows[i]] for i in range(2)]

# Column widths for neat spacing
widths = [max(len(cols[j]), max(len(r[j]) for r in data)) for j in range(len(cols))]

# Row joiner: left-align first col, right-align the numeric cols
def join(items):
    parts = []
    for j, cell in enumerate(items):
        w = widths[j]
        parts.append(cell.ljust(w) if j == 0 else cell.rjust(w))
    return "| " + " | ".join(parts) + " |"

header = join(cols)
sep    = "|" + "|".join("-" * (w + 2) for w in widths) + "|"
lines  = [header, sep] + [join(r) for r in data]

md = "\n".join(lines) + "\n"
print(md)