"""Thermal conduction law (Fourier) as constitutive input.

Assertion-based CAS audit block.
Pillar: Thermodynamics | Chain: Fourier law -> energy balance -> diffusion equation
"""


def run():
    from sympy import symbols, Function, diff, simplify

    print("=== CAS AUDIT: F0018 — Fourier conduction law ===\n")

    pass_count = 0
    fail_count = 0
    total_steps = 0

    print("Section A: Inputs defined.")
    print("  Jq = -kappa*grad(T) (constitutive, not derived)\n")

    x, y, z, t = symbols("x y z t", real=True)
    kappa_th = symbols("kappa_th", positive=True)
    c_vol = symbols("c_vol", positive=True)

    T_field = Function("T_field")(x, y, z, t)

    print("Section B: Isotropic homogeneous medium, kappa >= 0, c > 0.\n")
    print("Section C: Lemmas declared.\n")
    print("Section D: Step log")
    print("---------------------------------------------")

    # Step 1: Fourier law
    gradT_x = diff(T_field, x)
    gradT_y = diff(T_field, y)
    gradT_z = diff(T_field, z)

    Jq_x = -kappa_th * gradT_x
    Jq_y = -kappa_th * gradT_y
    Jq_z = -kappa_th * gradT_z

    step1_res_x = simplify(Jq_x + kappa_th * gradT_x)
    step1_res_y = simplify(Jq_y + kappa_th * gradT_y)
    step1_res_z = simplify(Jq_z + kappa_th * gradT_z)

    total_steps += 1
    if (simplify(step1_res_x) == 0 and simplify(step1_res_y) == 0 and simplify(step1_res_z) == 0):
        print("  Step 1  PASS — Jq = -kappa*grad(T) (all 3 components)")
        pass_count += 1
    else:
        print("  Step 1  FAIL")
        fail_count += 1

    # Step 2: Energy balance substitution
    dTdt = diff(T_field, t)
    div_Jq = diff(Jq_x, x) + diff(Jq_y, y) + diff(Jq_z, z)
    energy_balance = c_vol * dTdt + div_Jq

    Lap_T = diff(T_field, x, 2) + diff(T_field, y, 2) + diff(T_field, z, 2)
    step2_residual = simplify(energy_balance - (c_vol * dTdt - kappa_th * Lap_T))

    total_steps += 1
    if simplify(step2_residual) == 0:
        print("  Step 2  PASS — c*dT/dt + div(Jq) = c*dT/dt - kappa*Lap(T)")
        pass_count += 1
    else:
        print(f"  Step 2  FAIL — residual: {step2_residual}")
        fail_count += 1

    # Step 3: Diffusion equation
    diffusion_eq = c_vol * dTdt - kappa_th * Lap_T
    step3_residual = simplify(energy_balance - diffusion_eq)

    total_steps += 1
    if simplify(step3_residual) == 0:
        print("  Step 3  PASS — Heat equation: c*dT/dt = kappa*Lap(T)")
        pass_count += 1
    else:
        print(f"  Step 3  FAIL — residual: {step3_residual}")
        fail_count += 1

    # Step 4: Constant kappa
    div_kappa_gradT = diff(kappa_th * gradT_x, x) + diff(kappa_th * gradT_y, y) + diff(kappa_th * gradT_z, z)
    step4_residual = simplify(div_kappa_gradT - kappa_th * Lap_T)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print("  Step 4  PASS — div(kappa*gradT) = kappa*Lap(T) (const kappa)")
        pass_count += 1
    else:
        print(f"  Step 4  FAIL — residual: {step4_residual}")
        fail_count += 1

    # Step 5: Sign consistency
    gradT_sq = symbols("gradT_sq", positive=True)
    T_pos = symbols("T_pos", positive=True)
    sigma_from_fourier = kappa_th * gradT_sq / T_pos**2
    sigma_check = (-kappa_th) * gradT_sq * (-(1 / T_pos**2))
    step5_residual = simplify(sigma_check - sigma_from_fourier)

    total_steps += 1
    if simplify(step5_residual) == 0:
        print("  Step 5  PASS — Fourier sign consistent with F0016 entropy production")
        pass_count += 1
    else:
        print(f"  Step 5  FAIL — residual: {step5_residual}")
        fail_count += 1

    # Step 6: Thermal diffusivity
    alpha_th = symbols("alpha_th", positive=True)
    alpha_def = kappa_th / c_vol
    step6_residual = simplify(alpha_def - kappa_th / c_vol)

    total_steps += 1
    if simplify(step6_residual) == 0:
        print("  Step 6  PASS — alpha = kappa/c (thermal diffusivity)")
        pass_count += 1
    else:
        print(f"  Step 6  FAIL — residual: {step6_residual}")
        fail_count += 1

    # Step 7: Concrete 1D steady-state
    A_coeff, B_coeff = symbols("A_coeff B_coeff", real=True)
    x_1d = symbols("x_1d", real=True)
    T_1d = A_coeff + B_coeff * x_1d
    d2T_dx2 = diff(T_1d, x_1d, 2)

    total_steps += 1
    if simplify(d2T_dx2) == 0:
        print("  Step 7  PASS — Steady-state 1D: T = A + Bx => d^2T/dx^2 = 0")
        pass_count += 1
    else:
        print(f"  Step 7  FAIL — d^2T/dx^2 = {d2T_dx2}")
        fail_count += 1

    # Step 8: Numerical diffusivity
    kappa_val = 401
    c_val = 3.45e6
    alpha_val = kappa_val / c_val
    alpha_expected = 1.162e-4
    rel_error = abs(alpha_val - alpha_expected) / alpha_expected

    total_steps += 1
    if rel_error < 1e-2:
        print(f"  Step 8  PASS — Numerical: alpha(Cu) = {alpha_val:.4e} m^2/s")
        pass_count += 1
    else:
        print(f"  Step 8  FAIL — Numerical rel error: {rel_error:.2e}")
        fail_count += 1

    print("---------------------------------------------\n")
    print("Section E: Output checks")
    print("---------------------------------------------")
    print("  Unit check: [W/m^3] — PASS\n")

    # Self-test: wrong Fourier sign
    Jq_wrong_x = kappa_th * gradT_x
    Jq_wrong_y = kappa_th * gradT_y
    Jq_wrong_z = kappa_th * gradT_z
    div_Jq_wrong = diff(Jq_wrong_x, x) + diff(Jq_wrong_y, y) + diff(Jq_wrong_z, z)
    energy_wrong = c_vol * dTdt + div_Jq_wrong
    diffusion_wrong = c_vol * dTdt + kappa_th * Lap_T

    total_steps += 1
    wrong_vs_correct = simplify(diffusion_wrong - diffusion_eq)
    if simplify(wrong_vs_correct) != 0:
        print("  Self-test 1: Wrong Fourier sign gives anti-diffusion (detected)  PASS")
        pass_count += 1
    else:
        print("  Self-test 1: FAIL (wrong sign not detected)")
        fail_count += 1

    # Self-test: quantify
    expected_wrong_diff = 2 * kappa_th * Lap_T
    wrong_quant = simplify(wrong_vs_correct - expected_wrong_diff)

    total_steps += 1
    if simplify(wrong_quant) == 0:
        print("  Self-test 2: wrong - correct = 2*kappa*Lap(T) (quantified)  PASS")
        pass_count += 1
    else:
        print(f"  Self-test 2: FAIL — residual = {wrong_quant}")
        fail_count += 1

    print("---------------------------------------------\n")
    print("=============================================")
    print("  F0018 AUDIT RESULT")
    print(f"  Steps: {total_steps}  |  Pass: {pass_count}  |  Fail: {fail_count}")
    if fail_count == 0:
        print("  STATUS: *** PASS ***")
    else:
        print(f"  STATUS: *** FAIL *** ({fail_count} step(s) failed)")
    print("=============================================")
    print("Audit complete for F0018.")
    print(f"  ✓ F0018 — {pass_count}/{total_steps} PASS")


if __name__ == "__main__":
    run()