diff --git a/app/simulation/native_codegen/compiler.py b/app/simulation/native_codegen/compiler.py index 57205f4..a0878a0 100644 --- a/app/simulation/native_codegen/compiler.py +++ b/app/simulation/native_codegen/compiler.py @@ -260,7 +260,7 @@ def _compile_storage_anchored_program(network: SimulationNetwork) -> NativeProgr initial_volume = max(base, c.cvol0 / 100) if c.model_type == "amesim_pnch012" else base init.append(f"if (!native_gas_init({_number(c.p0)}, {_number(c.T0)}, {_number(initial_volume)}, &y[{si(c, 'm')}])) return 0;") lines.append(f"NativeGas gas_{gi};") - lines.append(f"if (!native_gas(y[{si(c, 'm')}], y[{si(c, 'U')}], {volume}, &gas_{gi})) return 0;") + lines.append(f"if (!native_gas_context(properties, y[{si(c, 'm')}], y[{si(c, 'U')}], {volume}, &gas_{gi})) return 0;") for field in ("m", "U"): put(c, field, f"y[{si(c, field)}]") for field in ("p", "T", "rho", "u", "h"): @@ -277,7 +277,7 @@ def _compile_storage_anchored_program(network: SimulationNetwork) -> NativeProgr continue a, b = chamber_at(c, "port_2"), chamber_at(c, "port_3") put(c, "xv", f"fmax(0.0, fmin(1.0, {get(c, 'res.signal')}))") - lines.append(f"if (!native_orifice({get(a, 'p')}, {get(b, 'p')}, {get(a, 'h')}, {get(b, 'h')}, {_number(c.effective_cq * c.maximum_area)}, {get(c, 'xv')}, &{get(c, 'port_2.m_flow')}, &{get(c, 'cm')}, &{get(c, 'gasvel')})) return 0;") + lines.append(f"if (!native_orifice_context(properties, {get(a, 'p')}, {get(b, 'p')}, {get(a, 'h')}, {get(b, 'h')}, {_number(c.effective_cq * c.maximum_area)}, {get(c, 'xv')}, &{get(c, 'port_2.m_flow')}, &{get(c, 'cm')}, &{get(c, 'gasvel')})) return 0;") assigned.update((key(c, "cm"), key(c, "gasvel"))) put(c, "port_3.m_flow", f"-{get(c, 'port_2.m_flow')}") put(c, "port_2.h_outflow", get(b, "h")) @@ -368,7 +368,10 @@ def _compile_storage_anchored_program(network: SimulationNetwork) -> NativeProgr f"const double model_atol[NSTATES] = {{{','.join(map(state_absolute_tolerance, state_keys))}}};", "const char *const model_output_keys[NOUTPUTS] = {" + ",".join(json.dumps(v.key, ensure_ascii=True) for v in variables) + "};", "int model_init(double *y) {", *init, "return 1; }", - "int model_eval(double t, const double *y, double *dy, double *w) {", "(void)t;", *lines, + "int model_eval(double t, const double *y, double *dy, double *w) {", "(void)t;", + f"NativePropertyState property_states[{min(256,max(16,4*len(components)))}];", + "NativePropertyCache property_cache, *properties=&property_cache;", + "native_properties_init(properties,property_states,sizeof(property_states)/sizeof(property_states[0]));", *lines, "for (int i=0;i={pb}?{hin(c,a,True)}:{hin(c,b,True)})' + T=f'native_temperature_ph_context(properties,{medium(c)},fmax(fmax({pa},{pb}),1),{pa}>={pb}?{hin(c,a,True)}:{hin(c,b,True)})' flow(c,a,pipe_flow(f'{medium(c)},{pa},{pb},{T},{num(c.diam)},{num(c.le)},{num(c.rr)},0')) else: opening = w(c,'xv') if kind!='amesim_pnor001' else '1.0' @@ -351,7 +351,7 @@ def compile_extended_program(network): area = c.effective_cq*(c.effective_area if kind=='amesim_pnor001' else c.maximum_area) inputs=f'{medium(c)},{pa},{pb},{hin(c,a)},{hin(c,b)},{num(area)},{opening}' outputs=(q(c,a),w(c,'cm'),w(c,'gasvel')) - code=f'if(!native_medium_orifice({inputs},&{outputs[0]},&{outputs[1]},&{outputs[2]})) return 0;' + code=f'if(!native_medium_orifice_context(properties,{inputs},&{outputs[0]},&{outputs[1]},&{outputs[2]})) return 0;' operations.append(Computation(f'flow:{c.name}.{a}',outputs,references(inputs),(code,))) flow_known.add(q(c,a));assigned.update((c.name+'.cm',c.name+'.gasvel')) flow(c,b,f'-{q(c,a)}') @@ -360,7 +360,7 @@ def compile_extended_program(network): for name in (['port_1'] if kind.endswith('1') else names): # Resistance uses the gas arriving from the upstream side, # including inflow into a PNL0001 storage volume. - T=f'({p(c,name)}>{gas}.p?native_temperature_ph({medium(c)},fmax({p(c,name)},1),{hin(c,name,True)}):{gas}.T)' + T=f'({p(c,name)}>{gas}.p?native_temperature_ph_context(properties,{medium(c)},fmax({p(c,name)},1),{hin(c,name,True)}):{gas}.T)' flow(c,name,pipe_flow(f'{medium(c)},{p(c,name)},{gas}.p,{T},{num(c.diam)},{num(c.le/(2 if kind.endswith("2") else 1))},{num(c.rr)},1')) for edge in network.connections: if edge.domain=='pneumatic': @@ -513,7 +513,7 @@ def compile_extended_program(network): if kind=='amesim_pnl0003': a,b=gases[c.name,1],gases[c.name,2] center=w(c,'dmctr') - put(c,'dmctr',f'native_pipe_flow({medium(c)},{a}.p,{b}.p,{a}.p>={b}.p?{a}.T:{b}.T,{num(c.diam)},{num(c.le)},{num(c.rr)},3)') + put(c,'dmctr',f'native_pipe_flow_context(properties,{medium(c)},{a}.p,{b}.p,{a}.p>={b}.p?{a}.T:{b}.T,{num(c.diam)},{num(c.le)},{num(c.rr)},3)') for half in halves: gas=gases[c.name,half];suffix=str(half) if half else '' names=['port_'+suffix] if half else pnames(c) @@ -536,7 +536,7 @@ def compile_extended_program(network): gas=gases[c.name,0] for name in pnames(c): flow=q(c,name);pp=f'({flow}>=0?fmax({p(c,name)},1):fmax({gas}.p,1))' - temp=f'({flow}>=0?fmax(native_temperature_ph({medium(c)},{pp},{hin(c,name,True)}),1):{gas}.T)' + temp=f'({flow}>=0?fmax(native_temperature_ph_context(properties,{medium(c)},{pp},{hin(c,name,True)}),1):{gas}.T)' diag.append((flow,pp,temp,c.le/2,0)) elif kind=='amesim_pnl0003': a,b=gases[c.name,1],gases[c.name,2];flow=w(c,'dmctr') @@ -545,13 +545,13 @@ def compile_extended_program(network): pa,pb=p(c,'port_1'),p(c,'port_2');pp=f'fmax(fmax({pa},{pb}),1)' if kind=='amesim_pnl0001': gas=gases[c.name,0] - temp=f'({pa}>{gas}.p?fmax(native_temperature_ph({medium(c)},{pp},{hin(c,"port_1",True)}),1):{gas}.T)' + temp=f'({pa}>{gas}.p?fmax(native_temperature_ph_context(properties,{medium(c)},{pp},{hin(c,"port_1",True)}),1):{gas}.T)' else: - temp=f'fmax(native_temperature_ph({medium(c)},{pp},{pa}>={pb}?{hin(c,"port_1",True)}:{hin(c,"port_2",True)}),1)' + temp=f'fmax(native_temperature_ph_context(properties,{medium(c)},{pp},{pa}>={pb}?{hin(c,"port_1",True)}:{hin(c,"port_2",True)}),1)' diag=[(q(c,'port_1'),pp,temp,c.le,0)] lines.append('{ double d[4],acc[4]={0};') for flow,pp,temp,length,diagnostic in diag: - lines.append(f'native_pipe_diagnostics({medium(c)},{flow},{pp},{temp},{num(c.diam)},{num(length)},{num(c.rr)},{diagnostic},d);') + lines.append(f'native_pipe_diagnostics_context(properties,{medium(c)},{flow},{pp},{temp},{num(c.diam)},{num(length)},{num(c.rr)},{diagnostic},d);') if len(diag)>1: lines.append('d[2]=fabs(d[2]);') lines.append('for(int i=0;i<4;i++) acc[i]+=d[i];') for i,field in enumerate(('re','cm','v','ff')): @@ -574,6 +574,9 @@ def compile_extended_program(network): *schedule_helpers, 'int model_init(double *y) {',*[f'y[{i}]={num(v)};' for i,v in enumerate(initial)],*gas_initializers,'return 1;}', 'int model_eval(double t,const double *y,double *dy,double *w) {', + f'NativePropertyState property_states[{min(256,max(16,4*gas_count+2*len(components)))}];', + 'NativePropertyCache property_cache, *properties=&property_cache;', + 'native_properties_init(properties,property_states,sizeof(property_states)/sizeof(property_states[0]));', *(['double projected[NSTATES];for(int i=0;ireal_helium==b->real_helium && a->R==b->R && a->cp==b->cp && + a->Tref==b->Tref && a->slope==b->slope && a->mu==b->mu && a->muT==b->muT && a->S==b->S; +} +void native_properties_init(NativePropertyCache *cache,NativePropertyState *states,size_t capacity) { + cache->states=states;cache->count=0;cache->capacity=states?capacity:0; +} +static NativePropertyState *property_new(NativePropertyCache *cache,const NativeMedium *m, + double p,double T,NativePropertyState *scratch) { + int valid=p>0 && T>0 && isfinite(p) && isfinite(T); + NativePropertyState *s=valid && cache && cache->countcapacity ? &cache->states[cache->count++] : scratch; + *s=(NativePropertyState){0};s->medium=*m;s->p=p;s->T=T;s->valid=valid?NATIVE_PROPERTY_PT:0; + return s; +} +static NativePropertyState *property_pt(NativePropertyCache *cache,const NativeMedium *m, + double p,double T,NativePropertyState *scratch) { + if(cache)for(size_t i=0;icount;i++) { + NativePropertyState *s=&cache->states[i]; + if((s->valid&NATIVE_PROPERTY_PT) && s->p==p && s->T==T && same_medium(&s->medium,m))return s; + } + return property_new(cache,m,p,T,scratch); +} +static double property_density(NativePropertyState *s) { + if(!(s->valid&NATIVE_PROPERTY_RHO)) { + s->rho=native_density(&s->medium,s->p,s->T); + if(s->rho>0 && isfinite(s->rho))s->valid|=NATIVE_PROPERTY_RHO; + } + return s->rho; +} +static double property_viscosity(NativePropertyState *s) { + if(!(s->valid&NATIVE_PROPERTY_MU)) { + s->mu=native_viscosity(&s->medium,s->T,0); + if(s->mu>0 && isfinite(s->mu))s->valid|=NATIVE_PROPERTY_MU; + } + return s->mu; +} +static void remember_gas(NativePropertyCache *cache,const NativeMedium *m,const NativeGas *gas) { + /* Below helium's critical temperature retain the existing vapor-root/PH + selection. A future phase-aware medium contract can carry that state. */ + if(!cache || (m->real_helium && gas->T<=TC))return; + NativePropertyState scratch,*s=property_pt(cache,m,gas->p,gas->T,&scratch); + if(!(s->valid&NATIVE_PROPERTY_PT) || !isfinite(gas->h))return; + if((s->valid&NATIVE_PROPERTY_H) && s->h!=gas->h) + s=property_new(cache,m,gas->p,gas->T,&scratch); + s->h=gas->h;s->valid|=NATIVE_PROPERTY_H; + if(gas->rho>0 && isfinite(gas->rho)){s->rho=gas->rho;s->valid|=NATIVE_PROPERTY_RHO;} +} +double native_temperature_ph_context(NativePropertyCache *cache,const NativeMedium *m,double p,double h) { + if(cache)for(size_t i=0;icount;i++) { + NativePropertyState *s=&cache->states[i]; + if((s->valid&NATIVE_PROPERTY_H) && s->p==p && s->h==h && same_medium(&s->medium,m))return s->T; + } + double T=native_temperature_ph(m,p,h); + if(cache && isfinite(h) && T>0 && isfinite(T) && p>0 && isfinite(p)) { + NativePropertyState scratch,*s=property_pt(cache,m,p,T,&scratch); + if((s->valid&NATIVE_PROPERTY_H) && s->h!=h)s=property_new(cache,m,p,T,&scratch); + s->h=h;s->valid|=NATIVE_PROPERTY_H; + } + return T; +} +double native_density_context(NativePropertyCache *cache,const NativeMedium *m,double p,double T) { + NativePropertyState scratch;return property_density(property_pt(cache,m,p,T,&scratch)); +} static double cube_root(double x) { return x==0 ? 0 : copysign(pow(fabs(x),1.0/3.0),x); } static void attraction(double T, double *a, double *da, double *dda) { @@ -68,39 +133,52 @@ static NativeGas gas_properties(double m,double U,double V) { double T=fmax(temperature_u(g.u),2.2); for(int i=0;i<16;i++) { double next=fmax(temperature_u(g.u-u_departure(T,g.rho)),2.2); int done=fabs(next-T)<=1e-10*fmax(T,1);T=next;if(done) break; } - g.T=T;g.p=pressure_rho(T,g.rho);g.h=h_ideal(T)+h_departure(g.p,T);return g; + g.T=T;g.p=pressure_rho(T,g.rho); + /* Same energy reference as u_ideal/h_ideal. Avoid a redundant cubic solve + in the single-root region; preserve vapor-root semantics below TC. */ + g.h=T>TC ? g.u+g.p/g.rho : h_ideal(T)+h_departure(g.p,T);return g; } static double temperature_ph(double p,double h) { double T=fmax(temperature_h(h),2.2); for(int i=0;i<16;i++) { double next=fmax(temperature_h(h-h_departure(p,T)),2.2); int done=fabs(next-T)<=1e-10*fmax(T,1);T=next;if(done) break; }return T; } -static void local_isentropic(double p,double T,double *factor,double *exponent) { - double rho=density(p,T),v=MOLAR_MASS/rho,a,da,dda;attraction(T,&a,&da,&dda); +static void local_isentropic(NativePropertyState *s) { + if(s->valid&NATIVE_PROPERTY_ISENTROPIC)return; + double p=s->p,T=s->T,rho=property_density(s),v=MOLAR_MASS/rho,a,da,dda;attraction(T,&a,&da,&dda); double d=v*(v+pr_b)+pr_b*(v-pr_b); double dpT=RU/(v-pr_b)-da/d; double dpR=(-RU*T/pow(v-pr_b,2)+a*2*(v+pr_b)/(d*d))*(-MOLAR_MASS/(rho*rho)); double cv=1.5*rg+T*dda*log_volume(rho); double cp=cv+T*dpT*dpT/(rho*rho*dpR),gamma=cp/cv; - *factor=p/(rho*dpR*gamma);*exponent=p*(gamma-1)/(gamma*T*dpT); + s->isentropic_factor=p/(rho*dpR*gamma);s->isentropic_exponent=p*(gamma-1)/(gamma*T*dpT); + if(isfinite(s->isentropic_factor) && isfinite(s->isentropic_exponent))s->valid|=NATIVE_PROPERTY_ISENTROPIC; } -static double isentropic(double p,double T,double pd) { - double f,e,fd,ed;local_isentropic(p,T,&f,&e);if(pd>=p) return f; - double Td=fmax(T*pow(fmax(pd/p,1e-12),e),2.2); - local_isentropic(fmax(pd,1),Td,&fd,&ed);return .5*(f+fd); +static double isentropic(NativePropertyCache *cache,NativePropertyState *up,double pd) { + local_isentropic(up);if(pd>=up->p)return up->isentropic_factor; + double Td=fmax(up->T*pow(fmax(pd/up->p,1e-12),up->isentropic_exponent),2.2); + NativePropertyState scratch,*down=property_pt(cache,&up->medium,fmax(pd,1),Td,&scratch); + local_isentropic(down);return .5*(up->isentropic_factor+down->isentropic_factor); } static double subsonic_cm(double r,double gamma,double rho,double T,double p) { return sqrt(fmax(2/(1-gamma)*rho*T/p*(pow(r,2*gamma)-pow(r,1+gamma)),0)); } -static void valve(double p,double pd,double T,double *cm,double *vel) { - p=fmax(p,1);pd=fmax(fmin(pd,p),0);T=fmax(T,1); - double g=fmax(1e-9,fmin(1-1e-9,isentropic(p,T,pd))),rho=fmax(density(p,T),1e-12); +static void state_valve(NativePropertyCache *cache,NativePropertyState *up,double pd,double *cm,double *vel) { + double p=up->p,T=up->T;pd=fmax(fmin(pd,p),0); + const NativeMedium *m=&up->medium; + double cp=m->cp+m->slope*(T-m->Tref); + double factor=m->real_helium?isentropic(cache,up,pd):(cp-m->R)/cp; + double g=fmax(1e-9,fmin(1-1e-9,factor)),rho=fmax(property_density(up),1e-12); double r=fmax(pd/p,0),critical=pow(2*g/(g+1),1/(1-g)),eff; if(r<=critical) { eff=critical;*cm=sqrt(2/(1+g)*rho*T/p)*pow(2*g/(g+1),g/(1-g));*vel=sqrt(2/(1+g)*p/rho); } else { eff=r;*cm=subsonic_cm(r,g,rho,T,p);*vel=sqrt(fmax(2/(1-g)*p/rho*(1-pow(r,1-g)),0)); } double ref=subsonic_cm(.9999,g,rho,T,p); if(*cm>0 && ref>0) { double smooth=tanh(fmax(12*fabs(*cm/ref)*log(eff)/log(.9999),0));*cm*=smooth;*vel*=smooth; } } +static void medium_valve(NativePropertyCache *cache,const NativeMedium *m,double p,double pd,double T,double *cm,double *vel) { + NativePropertyState scratch,*up=property_pt(cache,m,fmax(p,1),fmax(T,1),&scratch); + state_valve(cache,up,pd,cm,vel); +} int native_gas_init(double p, double T, double volume, double *mU) { if (!(p > 0 && T >= 2.2 && volume > 0)) return 0; @@ -115,14 +193,21 @@ int native_gas(double m, double U, double volume, NativeGas *gas) { *gas = gas_properties(m, U, volume); return gas->p > 0 && isfinite(gas->p) && isfinite(gas->h); } +int native_gas_context(NativePropertyCache *cache,double m,double U,double V,NativeGas *gas) { + int ok=native_gas(m,U,V,gas);if(ok)remember_gas(cache,&helium_medium,gas);return ok; +} int native_orifice(double p1, double p2, double h1, double h2, double cq_area, double opening, double *flow, double *cm, double *velocity) { + return native_orifice_context(NULL,p1,p2,h1,h2,cq_area,opening,flow,cm,velocity); +} +int native_orifice_context(NativePropertyCache *cache,double p1,double p2,double h1,double h2, + double cq_area,double opening,double *flow,double *cm,double *velocity) { int forward = p1 >= p2; double p = forward ? p1 : p2, pd = forward ? p2 : p1; - double T = temperature_ph(fmax(p, 1), forward ? h1 : h2); - valve(p, pd, T, cm, velocity); + double T = native_temperature_ph_context(cache,&helium_medium,fmax(p,1),forward?h1:h2); + medium_valve(cache,&helium_medium,p,pd,T,cm,velocity); double sign = forward ? 1 : -1; *velocity *= sign; *flow = opening == 0 || fabs(p1-p2) <= 1e-8 ? 0 : @@ -218,22 +303,18 @@ int native_medium_gas(const NativeMedium *medium, double m, double U, double V, g->h=medium->cp*g->T+.5*medium->slope*dt*dt; return g->T>0 && isfinite(g->p) && isfinite(g->h); } -static void medium_valve(const NativeMedium *m, double p, double pd, double T, double *cm, double *vel) { - if (m->real_helium) { valve(p,pd,T,cm,vel); return; } - p=fmax(p,1); pd=fmax(fmin(pd,p),0); T=fmax(T,1); - double cp=m->cp+m->slope*(T-m->Tref); - double g=fmax(1e-9,fmin(1-1e-9,(cp-m->R)/cp)),rho=fmax(native_density(m,p,T),1e-12); - double r=pd/p,critical=pow(2*g/(g+1),1/(1-g)),eff; - if(r<=critical) { eff=critical;*cm=sqrt(2/(1+g)*rho*T/p)*pow(2*g/(g+1),g/(1-g));*vel=sqrt(2/(1+g)*p/rho); } - else { eff=r;*cm=subsonic_cm(r,g,rho,T,p);*vel=sqrt(fmax(2/(1-g)*p/rho*(1-pow(r,1-g)),0)); } - double ref=subsonic_cm(.9999,g,rho,T,p); - if(*cm>0 && ref>0) { double smooth=tanh(fmax(12*fabs(*cm/ref)*log(eff)/log(.9999),0));*cm*=smooth;*vel*=smooth; } +int native_medium_gas_context(NativePropertyCache *cache,const NativeMedium *m,double mass,double U,double V,NativeGas *g) { + int ok=native_medium_gas(m,mass,U,V,g);if(ok)remember_gas(cache,m,g);return ok; } int native_medium_orifice(const NativeMedium *m, double p1, double p2, double h1, double h2, double area, double opening, double *q, double *cm, double *v) { + return native_medium_orifice_context(NULL,m,p1,p2,h1,h2,area,opening,q,cm,v); +} +int native_medium_orifice_context(NativePropertyCache *cache,const NativeMedium *m,double p1,double p2,double h1,double h2, + double area,double opening,double *q,double *cm,double *v) { double p=fmax(p1,p2),pd=fmin(p1,p2),sign=p1>=p2?1:-1; - double T=fmax(native_temperature_ph(m,fmax(p,1),p1>=p2?h1:h2),1); - medium_valve(m,p,pd,T,cm,v); + double T=fmax(native_temperature_ph_context(cache,m,fmax(p,1),p1>=p2?h1:h2),1); + medium_valve(cache,m,p,pd,T,cm,v); *q=fabs(p1-p2)<=1e-8?0:sign*area*opening*fmax(p,1)*(*cm)/sqrt(T); *v*=fabs(opening)<=1e-12?0:sign; return isfinite(*q) && isfinite(*cm) && isfinite(*v); @@ -253,38 +334,77 @@ static double pipe_friction_prepared(double re, double rr, double rough_limit) { static double pipe_friction(double re, double rr) { return pipe_friction_prepared(re,rr,pipe_rough_limit(rr)); } +/* Derivative of the existing friction blend, local to this scalar solve. */ +static double pipe_friction_derivative(double re,double rr,double rough,double *df) { + double lam=64/re,dl=-lam/re; + if(re<=89.96829989){*df=dl;return lam;} + double a=-1.8*log10(6.9/re),smooth=1/(a*a); + double ds=-2*smooth/a*1.8/(log(10.0)*re),turb=smooth,dt=ds; + if(rr>0){double r=re*rr,w=r*r/(r*r+180*180),dw=2*w*(1-w)/re; + turb+=w*(rough-smooth);dt=ds*(1-w)+dw*(rough-smooth);} + double z=pow((re-89.96829989)/2741.96700831,8.37293695),b=z/(1+z); + double db=8.37293695*b*(1-b)/(re-89.96829989); + *df=dl*(1-b)+b*dt+db*(turb-lam);return lam+b*(turb-lam); +} +double native_pipe_resistance(double K,double rr,double flow_per_re,NativePipeSolve *status) { + NativePipeSolve local={0,0,0,INFINITY};if(!status)status=&local; + *status=local; + if(!(K>=0 && rr>=0 && flow_per_re>0) || !isfinite(K) || !isfinite(rr) || !isfinite(flow_per_re))return NAN; + if(K/64<=89.96829989){status->converged=1;status->relative_residual=0;return K/64;} + double rough=pipe_rough_limit(rr),lo=0,hi=fmax(sqrt(K/.02),1),df; + for(int i=0;i<128 && hi*hi*pipe_friction_prepared(hi,rr,rough)iterations=i+1; + double f=pipe_friction_derivative(re,rr,rough,&df),F=re*re*f-K; + double target=sqrt(K/f); + if(fabs(target-re)*flow_per_re<=fmax(1e-13,re*flow_per_re*1e-10)) { + /* Check the returned point, not just the change of successive guesses. */ + status->relative_residual=fabs(target*target*pipe_friction_prepared(target,rr,rough)/K-1); + if(isfinite(target) && status->relative_residual<=1e-9){status->converged=1;return target;} + } + if(F>0)hi=re;else lo=re; + double next=re-F/(2*re*f+re*re*df); + if(!isfinite(next) || next<=lo || next>=hi){next=.5*(lo+hi);status->bisections++;} + re=next; + } + return NAN; +} double native_pipe_flow(const NativeMedium *m, double p1, double p2, double T, double d, double length, double rr, int kind) { + return native_pipe_flow_context(NULL,m,p1,p2,T,d,length,rr,kind); +} +double native_pipe_flow_context(NativePropertyCache *cache,const NativeMedium *m,double p1,double p2,double T, + double d,double length,double rr,int kind) { if(fabs(p1-p2)<=1e-8) return 0; double p=fmax(fmax(p1,p2),1),pd=fmin(p1,p2),sign=p1>p2?1:-1; T=fmax(T,1); - double area=PI*d*d/4,mu=native_viscosity(m,T,0),den=PI*d*mu; + NativePropertyState scratch,*up=property_pt(cache,m,p,T,&scratch); + double area=PI*d*d/4,mu=property_viscosity(up),den=PI*d*mu; /* PNL0001/2/3 share compressible flow and its near-equilibrium smoothing. PNL0003 differs in storage placement, not in the resistance law. */ - double cm,vel;medium_valve(m,p,pd,T,&cm,&vel); + double cm,vel;state_valve(cache,up,pd,&cm,&vel); if(kind==0) { double lam=pow(area*p*cm,2)/(16*PI*mu*length*T); if(4*lam/den<=1000) return sign*lam; } - double rough_limit=pipe_rough_limit(rr); - double base=area*p*cm/sqrt(T),q=sqrt(d/(length*.02))*base; - for(int i=0;i<(kind==0?64:16);i++) { - double next=sqrt(d/(length*pipe_friction_prepared(4*fabs(q)/den,rr,rough_limit)))*base; - if(fabs(next-q)<=fmax(1e-12,fabs(q)*1e-9)) return sign*next; - q=.5*(q+next); - } - return sign*q; + double base=area*p*cm/sqrt(T),K=pow(4*base/den,2)*d/length; + return sign*native_pipe_resistance(K,rr,den/4,NULL)*den/4; } double native_pipe_flow_cached(NativePipeCache *cache, const NativeMedium *m, double p1, double p2, double T, double d, double length, double rr, int kind) { + return native_pipe_flow_cached_context(NULL,cache,m,p1,p2,T,d,length,rr,kind); +} +double native_pipe_flow_cached_context(NativePropertyCache *properties,NativePipeCache *cache,const NativeMedium *m, + double p1,double p2,double T,double d,double length,double rr,int kind) { if(cache->valid && cache->p1==p1 && cache->p2==p2 && cache->T==T && cache->diameter==d && cache->length==length && cache->roughness==rr && cache->kind==kind && - cache->medium.real_helium==m->real_helium && cache->medium.R==m->R && - cache->medium.cp==m->cp && cache->medium.Tref==m->Tref && cache->medium.slope==m->slope && - cache->medium.mu==m->mu && cache->medium.muT==m->muT && cache->medium.S==m->S) + same_medium(&cache->medium,m)) return cache->flow; - double result=native_pipe_flow(m,p1,p2,T,d,length,rr,kind); + double result=properties?native_pipe_flow_context(properties,m,p1,p2,T,d,length,rr,kind): + native_pipe_flow(m,p1,p2,T,d,length,rr,kind); cache->valid=0; if(isfinite(result)) { cache->medium=*m;cache->p1=p1;cache->p2=p2;cache->T=T; @@ -295,10 +415,16 @@ double native_pipe_flow_cached(NativePipeCache *cache, const NativeMedium *m, } void native_pipe_diagnostics(const NativeMedium *m, double q, double p, double T, double d, double length, double rr, int diagnostic, double *r) { - double area=PI*d*d/4,re=4*fabs(q)/(PI*d*native_viscosity(m,T,0)),ff=pipe_friction(re,rr); - r[0]=4*fabs(q)/(PI*d*native_viscosity(m,T,diagnostic)); + native_pipe_diagnostics_context(NULL,m,q,p,T,d,length,rr,diagnostic,r); +} +void native_pipe_diagnostics_context(NativePropertyCache *cache,const NativeMedium *m,double q,double p,double T, + double d,double length,double rr,int diagnostic,double *r) { + (void)diagnostic; + NativePropertyState scratch,*state=property_pt(cache,m,p,T,&scratch); + double area=PI*d*d/4,re=4*fabs(q)/(PI*d*property_viscosity(state)),ff=pipe_friction(re,rr); + r[0]=re; r[1]=fabs(q)*sqrt(T)/fmax(sqrt(d/(length*ff))*area*p,1e-18); - r[2]=q/(fmax(native_density(m,p,T),1e-12)*area);r[3]=ff; + r[2]=q/(fmax(property_density(state),1e-12)*area);r[3]=ff; } double native_limit_force(double penetration, double velocity, double stiffness, double damping, double depth, int signed_force) { diff --git a/native/include/kernels.h b/native/include/kernels.h index c89f907..9da2fef 100644 --- a/native/include/kernels.h +++ b/native/include/kernels.h @@ -1,5 +1,6 @@ #ifndef NATIVE_KERNELS_H #define NATIVE_KERNELS_H +#include typedef struct { double p, T, rho, u, h; } NativeGas; typedef struct { int velocity_index; @@ -8,6 +9,30 @@ typedef struct { } NativeStop; /* Constants are emitted per medium instance by the model compiler. */ typedef struct { int real_helium; double R, cp, Tref, slope, mu, muT, S; } NativeMedium; +/* Caller-owned scratch for ONE model_eval. Keys use exact values and a copy of + * every medium constant. No process/thread global cache or approximate reuse. */ +enum { + NATIVE_PROPERTY_PT=1, NATIVE_PROPERTY_H=2, NATIVE_PROPERTY_RHO=4, + NATIVE_PROPERTY_MU=8, NATIVE_PROPERTY_ISENTROPIC=16 +}; +typedef struct { + NativeMedium medium; + double p, T, h, rho, mu, isentropic_factor, isentropic_exponent; + unsigned valid; +} NativePropertyState; +typedef struct { + NativePropertyState *states; + size_t count, capacity; +} NativePropertyCache; +void native_properties_init(NativePropertyCache *, NativePropertyState *, size_t capacity); +int native_gas_context(NativePropertyCache *, double m, double U, double V, NativeGas *); +int native_medium_gas_context(NativePropertyCache *, const NativeMedium *, double m, double U, double V, NativeGas *); +double native_temperature_ph_context(NativePropertyCache *, const NativeMedium *, double p, double h); +double native_density_context(NativePropertyCache *, const NativeMedium *, double p, double T); +int native_orifice_context(NativePropertyCache *, double p1, double p2, double h1, double h2, + double area, double opening, double *q, double *cm, double *v); +int native_medium_orifice_context(NativePropertyCache *, const NativeMedium *, double p1, double p2, double h1, double h2, + double area, double opening, double *q, double *cm, double *v); /* One entry per pipe branch, zero-initialized for each model_eval. Never shared * across solver trials. Exact inputs, including medium constants, form the key. */ typedef struct { @@ -15,6 +40,16 @@ typedef struct { double p1, p2, T, diameter, length, roughness, flow; int kind, valid; } NativePipeCache; +typedef struct { int converged, iterations, bisections; double relative_residual; } NativePipeSolve; +/* Solve Re^2*f(Re)=K for the shared pipe resistance law. On failure returns + * NAN, with converged=0; flow_per_re converts the local stopping test to kg/s. */ +double native_pipe_resistance(double K, double roughness, double flow_per_re, NativePipeSolve *status); +double native_pipe_flow_context(NativePropertyCache *, const NativeMedium *, double p1, double p2, double T, + double diameter, double length, double roughness, int kind); +double native_pipe_flow_cached_context(NativePropertyCache *, NativePipeCache *, const NativeMedium *, + double p1, double p2, double T, double diameter, double length, double roughness, int kind); +void native_pipe_diagnostics_context(NativePropertyCache *, const NativeMedium *, double q, double p, double T, + double diameter, double length, double roughness, int diagnostic, double *result); int native_medium_init(const NativeMedium *, double p, double T, double V, int legacy_ideal_initial, double *mU); double native_density(const NativeMedium *, double p, double T); double native_temperature_ph(const NativeMedium *, double p, double h); diff --git a/tests/test_native_catalog.py b/tests/test_native_catalog.py index e90ed2e..f9ff53c 100644 --- a/tests/test_native_catalog.py +++ b/tests/test_native_catalog.py @@ -12,11 +12,13 @@ from tests.native_reference import reference_data, reference_network def has_revised_pipe_law(case): - """The frozen Python flow laws predate the September 2026 corrections.""" - return any(c['type'] in ('amesim_pnl0001', 'amesim_pnl0003') or - (c['type'].startswith('amesim_pnl') and - c['medium']['type'] == 'AmesimHeliumPengRobinsonMedium') - for c in case['components']) + """The frozen flow values include capped pipe iterates, not checked roots. + + All PNL media now share the residual-checked scalar solver. Keep the frozen + thermodynamic oracle; the pipe suite checks the actual resistance equation + against independent bisection as well as system-level Amesim curves. + """ + return any(c['type'].startswith('amesim_pnl') for c in case['components']) class Circuit: diff --git a/tests/test_native_pipe_physics.py b/tests/test_native_pipe_physics.py index a999128..61f7988 100644 --- a/tests/test_native_pipe_physics.py +++ b/tests/test_native_pipe_physics.py @@ -68,6 +68,55 @@ int main(void) { run=subprocess.run([str(exe)],capture_output=True,text=True,timeout=15) self.assertEqual(run.returncode,0,run.stderr) + def test_scalar_pipe_solver_checks_residual_and_bracket_fallback(self): + harness = r''' +#include +#define CHECK(x) do { if(!(x)){fprintf(stderr,"line %d\n",__LINE__);return 1;} } while(0) +static double reference(double K,double rr) { + double lo=0,hi=fmax(sqrt(K/.02),1); + while(hi*hi*pipe_friction(hi,rr)K)hi=mid;else lo=mid; + } + return .5*(lo+hi); +} +int main(void) { + NativePipeSolve s;int fallbacks=0; + for(int i=0;i<=150;i++)for(int r=0;r<6;r++)for(int scale=0;scale<3;scale++) { + double re=pow(10,-6+.1*i),rr=r==0?0:pow(10,-6+r); + double K=re*re*pipe_friction(re,rr),q_per_re=pow(10,-12+6*scale); + double actual=native_pipe_resistance(K,rr,q_per_re,&s),expected=reference(K,rr); + CHECK(s.converged && isfinite(actual) && s.relative_residual<=1e-9); + CHECK(s.iterations<20 && fabs(actual-expected)<=1e-12+expected*1e-9); + CHECK(fabs(actual*actual*pipe_friction(actual,rr)/K-1)<=1e-9); + fallbacks+=s.bisections; + } + CHECK(fallbacks>0); + CHECK(native_pipe_resistance(0,0,1,&s)==0 && s.converged); + CHECK(isnan(native_pipe_resistance(-1,0,1,&s)) && !s.converged); + CHECK(isnan(native_pipe_resistance(NAN,0,1,&s)) && !s.converged); + CHECK(isnan(native_pipe_resistance(1,0,0,&s)) && !s.converged); + /* Direct laminar handling preserves the distinct PNL00R branch. */ + NativeMedium medium={0,287,1005,300,0,1.8e-5,300,110.4}; + double p=2e5,T=300,d=.01,L=1,cm,v; + medium_valve(NULL,&medium,p,p-1,T,&cm,&v); + double expected=pow(PI*d*d/4*p*cm,2)/(16*PI*native_viscosity(&medium,T,0)*L*T); + CHECK(4*expected/(PI*d*native_viscosity(&medium,T,0))<1000); + CHECK(native_pipe_flow(&medium,p,p-1,T,d,L,1e-5,0)==expected); + return 0; +} +''' + with tempfile.TemporaryDirectory(prefix='native-pipe-root-') as tmp: + directory=Path(tmp);source=directory/'check.c';exe=directory/'check.exe' + source.write_text((ROOT/'native/components/kernels.c').read_text()+harness) + build=subprocess.run([self.compiler,'-std=c11','-O3','-Wall','-Wextra','-Werror', + '-ffp-contract=off','-fno-fast-math','-static-libgcc','-I',str(ROOT/'native/include'), + str(source),'-lm','-o',str(exe)],capture_output=True,text=True,timeout=60) + self.assertEqual(build.returncode,0,build.stderr) + run=subprocess.run([str(exe)],capture_output=True,text=True,timeout=30) + self.assertEqual(run.returncode,0,run.stderr) + def test_pnl0001_uses_upstream_temperature_in_both_directions(self): for reverse in (False,True): with self.subTest(reverse=reverse): diff --git a/tests/test_native_properties.py b/tests/test_native_properties.py new file mode 100644 index 0000000..d9e9fe9 --- /dev/null +++ b/tests/test_native_properties.py @@ -0,0 +1,136 @@ +"""Property reuse must preserve state identity, trial isolation and EOS physics.""" +from pathlib import Path +import subprocess +import tempfile +import unittest + +from app.simulation.native_codegen.build import toolchain + +ROOT = Path(__file__).resolve().parents[1] + + +class NativePropertyTests(unittest.TestCase): + def test_shared_properties_and_invalidation(self): + try: + compiler = toolchain()[0] + except (OSError, RuntimeError, subprocess.SubprocessError) as exc: + self.skipTest(f"Native toolchain unavailable: {exc}") + source = (ROOT / 'native/components/kernels.c').read_text() + for signature, counter in ( + ('static double temperature_ph(double p,double h) {', 'ph_calls'), + ('static double z_factor(double p,double T) {', 'z_calls'), + ('double native_viscosity(const NativeMedium *m, double T, int diagnostic) {', 'mu_calls'), + ): + self.assertEqual(source.count(signature), 1) + source = source.replace(signature, signature + f'++{counter};') + source = 'static int ph_calls,z_calls,mu_calls;\n' + source + harness = r''' +#include +#define CHECK(x) do { if(!(x)){fprintf(stderr,"line %d\n",__LINE__);return 1;} } while(0) +static int close_to(double a,double b,double atol,double rtol) { + return isfinite(a) && isfinite(b) && fabs(a-b)<=atol+rtol*fmax(fabs(a),fabs(b)); +} +int main(void) { + NativeMedium m=helium_medium; + NativePropertyState states[32],other_states[2];NativePropertyCache cache,other; + native_properties_init(&cache,states,32);native_properties_init(&other,other_states,2); + double y[2];NativeGas gas; + CHECK(native_medium_init(&m,1.5e7,293.15,.01,0,y)); + int before=z_calls; + CHECK(native_medium_gas_context(&cache,&m,y[0],y[1],.01,&gas)); + CHECK(z_calls==before); /* h=u+p/rho does not solve the cubic again. */ + CHECK(gas.h==gas.u+gas.p/gas.rho); + before=ph_calls; + CHECK(native_temperature_ph_context(&cache,&m,gas.p,gas.h)==gas.T && ph_calls==before); + NativeMedium copy=m; + CHECK(native_temperature_ph_context(&cache,©,gas.p,gas.h)==gas.T && ph_calls==before); + double q,cm,v;before=z_calls; + CHECK(native_medium_orifice_context(&cache,&m,gas.p,.6*gas.p,gas.h,gas.h,.001,.5,&q,&cm,&v)); + CHECK(z_calls==before+1); /* Only the different downstream state needs rho. */ + CHECK(native_medium_orifice_context(&cache,&m,gas.p,.6*gas.p,gas.h,gas.h,.001,.5,&q,&cm,&v)); + CHECK(z_calls==before+1); + before=mu_calls; + q=native_pipe_flow_context(&cache,&m,gas.p,.6*gas.p,gas.T,.01,1,1e-5,1); + CHECK(isfinite(q) && mu_calls==before+1); + double diag[4];native_pipe_diagnostics_context(&cache,&m,q,gas.p,gas.T,.01,1,1e-5,1,diag); + CHECK(mu_calls==before+1); + + /* Changed pressure (even one ULP), mixed enthalpy and every medium field + must force a fresh PH evaluation. Equal h alone is not a state key. */ + double pressures[]={nextafter(gas.p,INFINITY),.2*gas.p,gas.p}; + for(int i=0;i<3;i++) { + double h=gas.h+(i==2?1000:0),expected=native_temperature_ph(&m,pressures[i],h); + before=ph_calls; + CHECK(native_temperature_ph_context(&cache,&m,pressures[i],h)==expected); + CHECK(ph_calls==before+1); + CHECK(native_temperature_ph_context(&cache,&m,pressures[i],h)==expected && ph_calls==before+1); + } + CHECK(fabs(native_temperature_ph(&m,.2*gas.p,gas.h)-gas.T)>1e-3); + double positive_h=gas.h+2e6; + (void)native_temperature_ph_context(&cache,&m,gas.p,positive_h); + for(int i=0;i<8;i++) { + copy=m; + double *fields[]={©.R,©.cp,©.Tref,©.slope,©.mu,©.muT,©.S}; + if(i==7)copy.real_helium=0;else *fields[i]+=fmax(fabs(*fields[i])*.001,.00001); + double expected=native_temperature_ph(©,gas.p,positive_h); + size_t count=cache.count; + CHECK(native_temperature_ph_context(&cache,©,gas.p,positive_h)==expected); + CHECK(cache.count>count); + } + /* Separate contexts, a new trial, and exhausted capacity remain correct. */ + before=ph_calls; + CHECK(isfinite(native_temperature_ph_context(&other,&m,gas.p,gas.h)) && ph_calls==before+1); + native_properties_init(&cache,states,1); + CHECK(native_medium_gas_context(&cache,&m,y[0],y[1],.01,&gas)); + for(int i=1;i<30;i++) { + double p=gas.p*(1+i*.01),h=gas.h+i; + CHECK(native_temperature_ph_context(&cache,&m,p,h)==native_temperature_ph(&m,p,h)); + CHECK(cache.count==1); + } + native_properties_init(&cache,states,32);before=ph_calls; + CHECK(isfinite(native_temperature_ph_context(&cache,&m,gas.p,gas.h)) && ph_calls==before+1); + size_t count=cache.count; + (void)native_temperature_ph_context(&cache,&m,NAN,gas.h); + CHECK(cache.count==count); + + /* Thermodynamic round trips use an independent PT enthalpy formula. + Ideal gases include temperature-dependent heat capacity. */ + for(int ideal=0;ideal<2;ideal++) { + NativeMedium fluid=ideal?(NativeMedium){0,287,1005,300,.2,1.8e-5,300,110.4}:m; + for(double T=100;T<=1000;T+=100)for(double p=1e4;p<=1e8;p*=10) { + native_properties_init(&cache,states,32); + CHECK(native_medium_init(&fluid,p,T,.01,0,y)); + CHECK(native_medium_gas_context(&cache,&fluid,y[0],y[1],.01,&gas)); + double ref_h=ideal?fluid.cp*T+.5*fluid.slope*pow(T-fluid.Tref,2):h_ideal(T)+h_departure(p,T); + CHECK(close_to(gas.h,ref_h,2e-5,1e-9)); + CHECK(close_to(gas.p,p,.001,1e-9)); + CHECK(close_to(gas.T,T,1e-7,1e-9)); + CHECK(close_to(native_temperature_ph(&fluid,gas.p,gas.h),gas.T,1e-7,1e-9)); + double oldq,oldcm,oldv; + CHECK(native_medium_orifice(&fluid,gas.p,.8*gas.p,gas.h,gas.h,.001,.5,&oldq,&oldcm,&oldv)); + CHECK(native_medium_orifice_context(&cache,&fluid,gas.p,.8*gas.p,gas.h,gas.h,.001,.5,&q,&cm,&v)); + CHECK(close_to(q,oldq,1e-10,2e-9) && close_to(v,oldv,1e-7,2e-9)); + } + } + /* Subcritical helium retains the pre-existing vapor-root selection. */ + native_properties_init(&cache,states,32); + CHECK(native_gas_init(1e4,3,.01,y)); + CHECK(native_medium_gas_context(&cache,&m,y[0],y[1],.01,&gas)); + CHECK(gas.T