"""Entropy production (sigma = Jq . grad(1/T) >= 0).

Assertion-based CAS audit block.
Pillar: Thermodynamics | Chain: local equil -> entropy balance -> Fourier -> non-negativity
"""


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

    print("=== CAS AUDIT: F0016 — Entropy production ===\n")

    pass_count = 0
    fail_count = 0
    total_steps = 0

    print("Section A: Inputs defined.\n")

    x, y, z = symbols("x y z", real=True)
    T_field = Function("T_field")(x, y, z)
    Jqx = Function("Jqx")(x, y, z)
    Jqy = Function("Jqy")(x, y, z)
    Jqz = Function("Jqz")(x, y, z)
    kappa_th = symbols("kappa_th", positive=True)

    print("Section B: Local equilibrium, T > 0, kappa >= 0.\n")
    print("Section C: Lemmas declared.\n")
    print("Section D: Step log")
    print("---------------------------------------------")

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

    invT = 1 / T_field
    grad_invT_x = diff(invT, x)
    expected_grad_invT_x = -(1 / T_field**2) * gradT_x
    step1_residual = simplify(grad_invT_x - expected_grad_invT_x)

    total_steps += 1
    if simplify(step1_residual) == 0:
        print("  Step 1  PASS — grad(1/T) = -(1/T^2)*grad(T) [x-component]")
        pass_count += 1
    else:
        print(f"  Step 1  FAIL — residual: {step1_residual}")
        fail_count += 1

    grad_invT_y = diff(invT, y)
    grad_invT_z = diff(invT, z)
    step1b_res = simplify(grad_invT_y - (-(1 / T_field**2) * gradT_y))
    step1c_res = simplify(grad_invT_z - (-(1 / T_field**2) * gradT_z))

    total_steps += 1
    if simplify(step1b_res) == 0 and simplify(step1c_res) == 0:
        print("  Step 1b PASS — grad(1/T) y,z components confirmed")
        pass_count += 1
    else:
        print("  Step 1b FAIL")
        fail_count += 1

    # Step 2: Product rule
    div_JqOverT = diff(Jqx / T_field, x) + diff(Jqy / T_field, y) + diff(Jqz / T_field, z)
    div_Jq = diff(Jqx, x) + diff(Jqy, y) + diff(Jqz, z)
    Jq_dot_grad_invT = Jqx * grad_invT_x + Jqy * grad_invT_y + Jqz * grad_invT_z
    expected_div = (1 / T_field) * div_Jq + Jq_dot_grad_invT
    step2_residual = simplify(div_JqOverT - expected_div)

    total_steps += 1
    if simplify(step2_residual) == 0:
        print("  Step 2  PASS — div(Jq/T) = (1/T)*div(Jq) + Jq.grad(1/T)")
        pass_count += 1
    else:
        print(f"  Step 2  FAIL — residual: {step2_residual}")
        fail_count += 1

    # Step 3: Energy balance
    dudt_sym, divJq_sym = symbols("dudt_sym divJq_sym", real=True)
    T_pos = symbols("T_pos", positive=True)
    combined = (1 / T_pos) * ((-divJq_sym) + divJq_sym)
    step3_residual = simplify(combined)

    total_steps += 1
    if simplify(step3_residual) == 0:
        print("  Step 3  PASS — (1/T)*(du/dt + div(Jq)) = 0")
        pass_count += 1
    else:
        print(f"  Step 3  FAIL — residual: {step3_residual}")
        fail_count += 1

    # Step 4: Sigma isolation
    Jq_dot_gradT = Jqx * gradT_x + Jqy * gradT_y + Jqz * gradT_z
    sigma_form1 = Jq_dot_grad_invT
    sigma_form2 = -(1 / T_field**2) * Jq_dot_gradT
    step4_residual = simplify(sigma_form1 - sigma_form2)

    total_steps += 1
    if simplify(step4_residual) == 0:
        print("  Step 4  PASS — sigma = Jq.grad(1/T) = -(1/T^2)*Jq.grad(T)")
        pass_count += 1
    else:
        print(f"  Step 4  FAIL — residual: {step4_residual}")
        fail_count += 1

    # Step 5: Fourier substitution
    gradT_sq = symbols("gradT_sq", positive=True)
    T_val = symbols("T_val", positive=True)
    sigma_fourier = kappa_th * gradT_sq / T_val**2
    Jq_dot_grad_invT_fourier = (-kappa_th) * gradT_sq * (-(1 / T_val**2))
    step5_residual = simplify(Jq_dot_grad_invT_fourier - sigma_fourier)

    total_steps += 1
    if simplify(step5_residual) == 0:
        print("  Step 5  PASS — sigma = kappa*|grad T|^2/T^2 (Fourier)")
        pass_count += 1
    else:
        print(f"  Step 5  FAIL — residual: {step5_residual}")
        fail_count += 1

    # Step 6: Non-negativity
    total_steps += 1
    if simplify(sigma_fourier) >= 0:
        print("  Step 6  PASS — sigma >= 0 (second law)")
        pass_count += 1
    else:
        print("  Step 6  FAIL — Non-negativity check")
        fail_count += 1

    # Step 7: Equilibrium
    sigma_equil = sigma_fourier.subs(gradT_sq, 0)
    step7_residual = simplify(sigma_equil)

    total_steps += 1
    if simplify(step7_residual) == 0:
        print("  Step 7  PASS — Equilibrium: grad T = 0 => sigma = 0")
        pass_count += 1
    else:
        print(f"  Step 7  FAIL — Equilibrium residual: {step7_residual}")
        fail_count += 1

    # Step 8: Concrete 1D test
    kappa_val = 200
    dTdx_val = 10
    T_val_num = 300
    sigma_num = kappa_val * dTdx_val**2 / T_val_num**2
    sigma_expected = 2 / 9
    rel_error = abs(sigma_num - sigma_expected) / sigma_expected

    total_steps += 1
    if rel_error < 1e-12:
        print(f"  Step 8  PASS — Numerical: sigma = {sigma_num:.6f} W/(m^3 K)")
        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*K)] — PASS\n")

    # Self-test: wrong Fourier
    sigma_wrong = -kappa_th * gradT_sq / T_val**2
    wrong_residual = simplify(sigma_wrong - sigma_fourier)

    total_steps += 1
    if simplify(wrong_residual) != 0:
        print("  Self-test: Wrong Fourier sign gives sigma < 0 (detected)  PASS")
        pass_count += 1
    else:
        print("  Self-test: FAIL (wrong sign not detected)")
        fail_count += 1

    expected_wrong_res = -2 * kappa_th * gradT_sq / T_val**2
    wrong_quant = simplify(wrong_residual - expected_wrong_res)

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

    print("---------------------------------------------\n")
    print("=============================================")
    print("  F0016 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 F0016.")
    print(f"  ✓ F0016 — {pass_count}/{total_steps} PASS")


if __name__ == "__main__":
    run()