Issue with absorber column

[tsbdjo]

I'm not sure what's wrong here.

def absorber(feed_gas, lean_amine, component_list, n_comp, fluid, TPflash,
T_top, T_bottom, P_top, P_bottom, N_trays, EFFICIENCY, mixRule, debug=False,
max_iter=20, tolerance=1e-6):
vap_flow = np.ones(N_trays + 1) * (feed_gas.getTotalNumberOfMoles() / (N_trays+1))
vap_flow[-1] = feed_gas.getTotalNumberOfMoles()
liq_flow = np.ones(N_trays + 1) * (lean_amine.getTotalNumberOfMoles() / (N_trays+1))
liq_flow[0] = lean_amine.getTotalNumberOfMoles()

# Initial guess for compositions (assume each is the entering stream composition)
vap_comp = np.zeros((N_trays + 1, n_comp))
liq_comp = np.zeros((N_trays + 1, n_comp))

for iteration in range(max_iter):
    # Store previous for convergence check
    vap_flow_prev = vap_flow.copy()
    liq_flow_prev = liq_flow.copy()
    vap_comp_prev = vap_comp.copy()
    liq_comp_prev = liq_comp.copy()

    if debug:
        print(f"\n=== Absorber Iteration {iteration+1} ===\n")

    # Apply vapor and liquid inlet conditions
    vap_flow[-1] = feed_gas.getTotalNumberOfMoles()
    vap_comp[-1, :] = [feed_gas.getComponent(comp).getz() if feed_gas.getComponent(comp) is not None else 0.0 for comp in component_list]
    liq_flow[0] = lean_amine.getTotalNumberOfMoles()
    liq_comp[0, :] = [lean_amine.getComponent(comp).getz() if lean_amine.getComponent(comp) is not None else 0.0 for comp in component_list]

    # Sweep from bottom tray up (bottom-up is standard in equilibrium stage columns)
    for tray in reversed(range(N_trays)):
        T_tray = T_bottom + (tray * (T_top - T_bottom) / N_trays)
        P_tray = P_bottom + (tray * (P_top - P_bottom) / N_trays)
        n_vap_in = vap_flow[tray + 1]
        n_liq_in = liq_flow[tray]
        n_total = n_vap_in + n_liq_in

        if debug:
            print(f"\n--- Tray {tray + 1} ---")
            print(f"Temperature: {T_tray:.2f} °C | Pressure: {P_tray:.2f} bara")
            print(f"Vapor in: {n_vap_in:.5f} mol | Liquid in: {n_liq_in:.5f} mol")

        vap_comp = np.nan_to_num(vap_comp, nan=0.0, posinf=0.0, neginf=0.0)
        liq_comp = np.nan_to_num(liq_comp, nan=0.0, posinf=0.0, neginf=0.0)
        vap_comp = np.clip(vap_comp, 0.0, 1.0)
        liq_comp = np.clip(liq_comp, 0.0, 1.0)

        tray_fluid = fluid(feed_gas.getFluidName())

        safe_total = 0.0
        for j, comp in enumerate(component_list):
            mole_vap = max(vap_comp[tray + 1, j], 0.0)
            mole_liq = max(liq_comp[tray, j], 0.0)
            moles = n_vap_in * mole_vap + n_liq_in * mole_liq
            if np.isnan(moles) or moles < 1e-20:
                moles = 1e-20  # enforce trace but never negative or nan
            if moles > 0:
                tray_fluid.addComponent(comp, moles)
                safe_total += moles  # will help with optional normalization debug

        if safe_total == 0:
            raise ValueError(f"All component moles on tray {tray+1} dropped to zero!")

        tray_fluid.setTemperature(T_tray, "C")
        tray_fluid.setPressure(P_tray, "bara")
        tray_fluid.setMixingRule(mixRule)
        tray_fluid.setMultiPhaseCheck(True)
        tray_fluid.chemicalReactionInit()
        print("component_list:", component_list)
        print("vap_comp.shape:", vap_comp.shape)
        print("liq_comp.shape:", liq_comp.shape)
        if debug:
            print("Sum of tray components (moles):", safe_total)
        TPflash(tray_fluid)

        # Identify output phases
        gas_idx, aq_idx = None, None
        for idx in range(tray_fluid.getNumberOfPhases()):
            phase_name = str(tray_fluid.getPhase(idx).getPhaseTypeName()).lower()
            if "gas" in phase_name:
                gas_idx = idx
            elif "aqueous" in phase_name or "liquid" in phase_name:
                aq_idx = idx

        if gas_idx is not None:
            vap_phase = tray_fluid.getPhase(gas_idx)
            n_vap_out = vap_phase.getNumberOfMolesInPhase()
            vapor_out = np.array([
                vap_phase.getComponent(comp).getx() if vap_phase.hasComponent(comp) else 0.0
                for comp in component_list
            ])
            vap_flow[tray] = n_vap_out
            vap_comp[tray, :] = EFFICIENCY * vapor_out + (1 - EFFICIENCY) * vap_comp[tray, :]
        else:
            vap_flow[tray] = 0.0
            vapor_out = np.zeros(n_comp)
        
        if aq_idx is not None:
            liq_phase = tray_fluid.getPhase(aq_idx)
            n_liq_out = liq_phase.getNumberOfMolesInPhase()
            liquid_out = np.array([
                liq_phase.getComponent(comp).getx() if liq_phase.hasComponent(comp) else 0.0
                for comp in component_list
            ])
            liq_flow[tray] = n_liq_out
            liq_comp[tray, :] = EFFICIENCY * liquid_out + (1 - EFFICIENCY) * liq_comp[tray, :]

        else:
            n_liq_out = 0.0
            liquid_out = np.zeros(n_comp)

        if np.sum(vapor_out) > 0:
            vapor_out /= np.sum(vapor_out)
        if np.sum(liquid_out) > 0:
            liquid_out /= np.sum(liquid_out)

    # Check convergence
    flow_diff = np.max(np.abs(vap_flow - vap_flow_prev)) + np.max(np.abs(liq_flow - liq_flow_prev))
    comp_diff = np.max(np.abs(vap_comp - vap_comp_prev)) + np.max(np.abs(liq_comp - liq_comp_prev))
    if debug:
        print(f"\nAbsorber recycle convergence check: Δflow={flow_diff:.6e}, Δcomp={comp_diff:.6e}")
    if flow_diff < tolerance and comp_diff < tolerance:
        if debug:
            print("Absorber converged!\n")
        break

# Output result
sweet_gas = {'flow': vap_flow[0], 'composition': vap_comp[0, :].copy()}
rich_amine = {'flow': liq_flow[-1], 'composition': liq_comp[-1, :].copy()}
return sweet_gas, rich_amine