"""Reduced quasi-static invasion front coupled to ideal-gas pockets.
NOT Navier-Stokes, NOT resolved VOF, NOT experimentally validated casting CFD.
Finite-volume gas inventory is conserved except explicit vent/source terms.
"""
from dataclasses import dataclass,asdict
import heapq,math
import numpy as np
from scipy import ndimage
from scipy.optimize import brentq
from filling import pipe_velocity

@dataclass
class FrontSettings:
    gas_temperature_C: float = 100.0
    gas_viscosity: float = 2.2e-5
    permeability_m2: float = 1e-12
    atmospheric_pa: float = 101325.0
    surface_tension: float = .85
    contact_angle_deg: float = 120.0
    gas_generation_cm3_m2_s: float = 0.0
    max_time_s: float = 60.0
    max_step_s: float = .1

FRONT_LIMITS=[
 'KVASISTATICKA FRONTOVA APROXIMACE. Nejde o Navier-Stokes/VOF ani spolehlivou predpoved zaliti.',
 'Poloha cela: invazni fronta v sousednich voxelech, prednostne nizsi bunky. Zadna hybnost, proudove paprsky, vlny nebo turbulence.',
 'Idealni plyn v kazde propojene kapse je izotermicky a ma prostorove jednotny tlak; teplota plynu je zadana.',
 'Plyn unika Poiseuilleovym modelem kanalku s orificnim omezenim a Darcyho modelem pres formu. Parametry jsou odhady.',
 'Smaceni je zastoupeno statickym kapilarnim tlakem ekvivalentniho kruhoveho pruchodu; zadny dynamicky kontaktni uhel.',
 'Oxidace: pouze integral plochy rozhrani kov-plyn v case [m2 s]. Nejde o hmotnost/tloustku blan ani jejich zavleceni.',
 'Neobsahuje vodik rozpusteny v kovu, odpar kovu, chemickou kinetiku, tuhnuti proudiciho kovu ani defekty po ztuhnuti.',
 'Kapalna frakce bunek popisuje naplneni, nikoli termodynamicky podil kapaliny. Castecne bunky mohou nadhodnotit propojeni plynu.',
]

def gas_step(n,V,dt,C,T,pamb):
    """Implicit compressible Darcy/Poiseuille outflow C*(p²-pamb²).
    C has units mol/(s Pa²); permits ambient inflow for p<pamb.
    """
    R=8.314462618;b=V/(R*T)
    if C==0:return n,n*R*T/V
    a=dt*C;c=n+a*pamb*pamb
    p=2*c/(b+math.sqrt(b*b+4*a*c))
    return b*p,p

def run_front(mask,origin,s,m,c,fs,gs,progress=None,cancel=None,snapshot_plane=False):
    from core import snapshot_interval
    interval=snapshot_interval(s,mask.shape,gs.max_time_s,plane=snapshot_plane)
    def snapshot():
        return alpha[:,mask.shape[1]//2,:].T.copy() if snapshot_plane else alpha.copy()
    vals=asdict(gs)
    if not all(np.isfinite(v) for v in vals.values()):raise ValueError('Invalid gas inputs.')
    if min(gs.gas_viscosity,gs.atmospheric_pa,gs.max_time_s,gs.max_step_s)<=0 or gs.gas_temperature_C<=-273.15:raise ValueError('Invalid gas temperature, viscosity, time or pressure.')
    if min(gs.permeability_m2,gs.surface_tension,gs.gas_generation_cm3_m2_s)<0 or not 0<=gs.contact_angle_deg<=180:raise ValueError('Invalid permeability or wetting parameters.')
    if fs.viscosity<=0 or fs.liquid_density<=0 or fs.minor_loss<0 or fs.roughness_mm<0 or fs.backpressure_pa<0:raise ValueError('Invalid hydraulics.')
    R=8.314462618;Tgas=gs.gas_temperature_C+273.15;pa=float(gs.atmospheric_pa);rho=fs.liquid_density;mu=fs.viscosity
    dx=s.pitch_mm/1000;dv=dx**3;shape=mask.shape
    coords=np.argwhere(mask);world=origin+coords*s.pitch_mm
    start=tuple(coords[np.argmin(np.linalg.norm(world-np.array(c['gate']),axis=1))])
    vent_nodes=[tuple(coords[np.argmin(np.linalg.norm(world-np.asarray(v),axis=1))]) for v in c['vents']]
    top=max(world[:,2].max()+s.pitch_mm/2+s.head_mm,c.get('feeding_geometry',{}).get('top_z_mm',0))/1000
    alpha=np.zeros(shape);n=np.where(mask,pa*dv/(R*Tgas),0.);initial_moles=float(n.sum())
    arrival=np.full(shape,np.nan);maxp=np.zeros(shape);diam=np.maximum(dx,2*ndimage.distance_transform_edt(mask)*dx-dx)
    offsets=[(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(0,0,1),(0,0,-1)]
    def neighbors(u):
        for off in offsets:
            v=tuple(u[i]+off[i] for i in range(3))
            if all(0<=v[i]<shape[i] for i in range(3)):yield v
    # Count contact area of each cavity cell with the surrounding mold.
    faces=np.zeros(shape)
    for ax in range(3):
        faces+=(mask&~np.roll(mask,1,axis=ax));faces+=(mask&~np.roll(mask,-1,axis=ax))
    distpath=np.full(shape,np.inf);distpath[start]=dx
    bottleneck=np.zeros(shape);bottleneck[start]=min(s.gate_mm/1000,diam[start])
    heap=[(start[2],0.,start)];queued=np.zeros(shape,bool);queued[start]=True
    time=0.;steps=0;escaped=0.;generated=0.;exposure=0.;history=[];snaps=[(0.,snapshot())];next_snapshot=interval
    struct=ndimage.generate_binary_structure(3,1)
    global_pressure_max=pa;last_notice=-1
    # Each iteration partially fills one selected frontier voxel.
    while heap and time<gs.max_time_s and steps<60000:
        if cancel and cancel():raise InterruptedError('Front model stopped.')
        _,_,u=heapq.heappop(heap)
        if alpha[u]>=1-1e-7:continue
        gasmask=mask&(alpha<1-1e-7)
        labels,count=ndimage.label(gasmask,structure=struct)
        lab=int(labels[u])
        if lab==0:continue
        gasvol=np.where(gasmask,dv*(1-alpha),0.)
        vol=np.bincount(labels.ravel(),weights=gasvol.ravel(),minlength=count+1)
        inventory=np.bincount(labels.ravel(),weights=n.ravel(),minlength=count+1)
        # Conductance through mold approximated by thickness equal to margin.
        wallarea=np.bincount(labels.ravel(),weights=(faces*dx*dx*gasmask).ravel(),minlength=count+1)
        C=gs.permeability_m2*wallarea/(2*gs.gas_viscosity*max(s.margin_mm/1000,dx)*R*Tgas)
        Cporous=C.copy(); Cvent=np.zeros(count+1)
        ventarea=np.zeros(count+1)
        for v in vent_nodes:
            labv=int(labels[v])
            if labv:
                length=max(top-(origin[2]+v[2]*s.pitch_mm)/1000,dx)
                rad=s.vent_mm/2000
                Cvent[labv]+=math.pi*rad**4/(16*gs.gas_viscosity*length*R*Tgas)
                C[labv]=Cporous[labv]+Cvent[labv]
                ventarea[labv]+=math.pi*rad*rad
        # Actual metal/gas face area (fractional faces approximated by alpha difference).
        interface=np.zeros(shape)
        for ax in range(3):
            for direction in (-1,1):
                neighbor=np.roll(alpha,direction,axis=ax)
                interface+=np.maximum(neighbor-alpha,0)*gasmask*dx*dx
        interarea=np.bincount(labels.ravel(),weights=interface.ravel(),minlength=count+1)
        source=gs.gas_generation_cm3_m2_s*1e-6*pa/(R*293.15)*interarea
        d=float(bottleneck[u]);L=float(distpath[u])+s.head_mm/1000+.02
        z=(origin[2]+u[2]*s.pitch_mm)/1000
        capillary=max(0.,-4*gs.surface_tension*math.cos(math.radians(gs.contact_angle_deg))/d)
        pressure_available=rho*9.80665*(top-z)-capillary-fs.backpressure_pa
        freehead=max(pressure_available,0)/(rho*9.80665)
        vfree=pipe_velocity(freehead,rho,mu,d,L,fs.minor_loss,fs.roughness_mm/1000)
        re=rho*vfree*d/mu
        blend=float(np.clip((re-2300)/1700,0,1));ff=(-1.8*math.log10((fs.roughness_mm/1000/d/3.7)**1.11+6.9/max(re,2300)))**-2
        a=.5*rho*(fs.minor_loss+blend*ff*L/d);b=(1-blend)*32*mu*L/d**2
        area=math.pi*d*d/4
        def rate(p):
            dp=max(0.,pressure_available-(p-pa))
            if dp<=0:return 0.
            if a==0:return area*dp/max(b,1e-12)
            return area*2*dp/(b+math.sqrt(b*b+4*a*dp))
        # At most a quarter of a cell per normal step; long enough to permit venting if blocked.
        pold=inventory[lab]*R*Tgas/vol[lab]
        qold=rate(pold)
        dt=min(gs.max_step_s,gs.max_time_s-time,next_snapshot-time,max(1e-5,.25*dv/max(qold,1e-15)))
        remaining=dv*(1-alpha[u])
        # Do not numerically delete compressed gas in the last gas cell.
        maxfill=min(remaining,vol[lab]*(1-1e-9))
        def dq(p):return min(rate(p)*dt,maxfill)
        def outflow(k,p):
            dp=p-pa
            tube=Cvent[k]*(p*p-pa*pa)
            gasrho=max(p,pa)*.02897/(R*Tgas)
            orifice=.65*ventarea[k]*math.sqrt(2*abs(dp)/gasrho)*max(p,pa)/(R*Tgas)
            limited=math.copysign(min(abs(tube),orifice),dp)
            return Cporous[k]*(p*p-pa*pa)+limited
        def balance(p):return (vol[lab]-dq(p))*p/(R*Tgas)-inventory[lab]-dt*source[lab]+dt*outflow(lab,p)
        upper=max(pa+max(pressure_available,0)+1000,pold*2,pa*2)
        while balance(upper)<0:upper*=2
        pressure=brentq(balance,1.,upper,xtol=1e-5)
        fill=dq(pressure)
        oldtotal=float(n.sum());before=float(alpha[u]);alpha[u]+=fill/dv
        if alpha[u]>1-1e-7:alpha[u]=1.
        actual=(alpha[u]-before)*dv
        newV=vol[lab]-actual
        gasvol=np.where(gasmask,dv*(1-alpha),0.)
        # Pocket inventory update; selected pocket couples pressure and displaced volume.
        newinventory=np.zeros(count+1);pressures=np.full(count+1,pa,dtype=float)
        for k in range(1,count+1):
            if k==lab:
                newinventory[k]=max(0.,newV)*pressure/(R*Tgas);pressures[k]=pressure
            else:
                def residual(p):return vol[k]*p/(R*Tgas)-inventory[k]-dt*source[k]+dt*outflow(k,p)
                upperk=max(pa*2,inventory[k]*R*Tgas/vol[k]*2)
                while residual(upperk)<0:upperk*=2
                pressures[k]=brentq(residual,1.,upperk,xtol=1e-5)
                newinventory[k]=vol[k]*pressures[k]/(R*Tgas)
        distribution=np.divide(newinventory, np.maximum(vol-np.eye(1,count+1,lab).ravel()*actual,1e-30))
        n=gasvol*distribution[labels]
        gen=float(source[1:].sum()*dt);generated+=gen;escaped+=dt*sum(outflow(k,pressures[k]) for k in range(1,count+1))
        maxp=np.maximum(maxp,np.where(gasmask,pressures[labels]-pa,0))
        global_pressure_max=max(global_pressure_max,float(pressures.max()))
        exposure+=float(interarea[1:].sum()*dt)
        time+=dt;steps+=1
        if alpha[u]>=1-1e-7:
            arrival[u]=time
            for v in neighbors(u):
                if mask[v] and not queued[v]:
                    queued[v]=True;distpath[v]=distpath[u]+dx;bottleneck[v]=min(bottleneck[u],diam[v])
                    heapq.heappush(heap,(v[2],distpath[v],v))
        else:
            # A blocked branch is retried after other frontier cells.
            priority=u[2] if fill>1e-8*dv else shape[2]+time
            heapq.heappush(heap,(priority,distpath[u]+time*1e-9,u))
        frac=float(alpha.sum()/mask.sum())
        if time>=next_snapshot-1e-9:
            snaps.append((time,snapshot()));next_snapshot=time+interval
        if steps%50==0 or not heap:
            sealed=float(vol[(C==0)&(np.arange(count+1)>0)].sum()*1e6)
            history.append([time,frac,float(pressures.max()-pa),sealed,exposure,float(n.sum()),escaped,generated])
        notice=int(time*2)
        if progress and notice!=last_notice:
            progress(min(.99,time/gs.max_time_s),f'Fronta {time:.2f} s | naplneno {frac*100:.1f} % | plyn {pressures.max()-pa:.0f} Pa');last_notice=notice
        if frac>=1-1e-6:break
    history.append([time,float(alpha.sum()/mask.sum()),global_pressure_max-pa,0.,exposure,float(n.sum()),escaped,generated])
    if time>snaps[-1][0]+1e-9:snaps.append((time,snapshot()))
    error=float(n.sum())+escaped-generated-initial_moles
    return {'fill_fraction':alpha,'arrival_s':arrival,'max_gas_gauge_pa':maxp,'history':np.array(history),
      'snapshot_times':np.array([t for t,a in snaps]),'snapshots':np.array([a for t,a in snaps]),
      'snapshot_plane':snapshot_plane,'snapshot_interval_s':interval,
      'settings':asdict(gs),'hydraulics':asdict(fs),'thermal_inputs':asdict(s),'material':asdict(m),'channels':c,
      'geometry':{'shape':list(shape),'origin_mm':origin.tolist(),'pitch_mm':s.pitch_mm},'limitations':FRONT_LIMITS,
      'summary':{'time_s':time,'steps':steps,'filled_fraction_reduced_model':float(alpha.sum()/mask.sum()),
       'maximum_gas_gauge_pa':global_pressure_max-pa,'interface_exposure_m2_s':exposure,
       'gas_balance_error_mol':error,'initial_gas_mol':initial_moles,
       'termination':'model volume filled' if alpha.sum()/mask.sum()>=1-1e-6 else 'time/step limit or unreachable region',
       'warning':'Not a prediction of actual filling, oxides, porosity or agreement with reality.'}}