562 lines
32 KiB
Python
562 lines
32 KiB
Python
"""Static C lowering for the complete built-in component catalog.
|
|
|
|
Python constructs the graph and eliminates constant linear constraints once.
|
|
All thermodynamics, flow/stream closure and derivatives execute in the EXE.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import json
|
|
from importlib import import_module
|
|
|
|
from .compiler import NativeCapabilityError, NativeProgram, _Groups, _number as num
|
|
from .contracts import SUPPORTED_VERSIONS
|
|
|
|
|
|
GAS_TYPES = {'amesim_pnch023', 'amesim_pnch012', 'amesim_pnl0001',
|
|
'amesim_pnl0002', 'amesim_pnl0003', 'cylinder', 'tank'}
|
|
NODES = {'amesim_pn3node2', 'amesim_p4node2', 'tee'}
|
|
RESISTORS = {'amesim_pnor001', 'amesim_pnvo001_fixed', 'amesim_pnvo001',
|
|
'amesim_pnl00r', 'pipe', 'orifice'}
|
|
|
|
|
|
def catalog_contracts():
|
|
from app.simulation.components.amesim.library import LIBRARY as a
|
|
from app.simulation.components.experimental.library import LIBRARY as e
|
|
result = {}
|
|
for entry in (*a.models, *e.models):
|
|
module, name = entry.split(':')
|
|
cls = getattr(import_module(module), name)
|
|
result[cls.MODEL_TYPE] = cls
|
|
return result
|
|
|
|
|
|
def linear_schedule(equations, unknowns, free_flows=()):
|
|
"""Eliminate constant coefficients; retain RHS expressions as C temporaries.
|
|
|
|
Exact pivot elimination avoids a numeric pseudoinverse and its tiny spurious
|
|
dependencies. Redundant rows are left to the nodal pressure closure.
|
|
"""
|
|
rows = [[dict(a), {i: 1.0}] for i, (a, _) in enumerate(equations)]
|
|
pivots = []
|
|
for key in unknowns:
|
|
found = next((i for i in range(len(pivots), len(rows)) if abs(rows[i][0].get(key, 0)) > 1e-14), None)
|
|
if found is None and key in free_flows:
|
|
# Ideal coupled pipe compliances have a redundant internal flow.
|
|
# Their total m/U derivative is distributed by physical volume.
|
|
index = len(equations)
|
|
equations.append(({key: 1.0}, '0.0'))
|
|
rows.append([{key: 1.0}, {index: 1.0}])
|
|
found = len(rows)-1
|
|
if found is None:
|
|
raise NativeCapabilityError(f'Underdetermined native connection constraint: {key}')
|
|
j = len(pivots)
|
|
rows[j], rows[found] = rows[found], rows[j]
|
|
a, b = rows[j]
|
|
pivot = a[key]
|
|
rows[j] = [{k: v/pivot for k, v in a.items()}, {k: v/pivot for k, v in b.items()}]
|
|
for i, (a, b) in enumerate(rows):
|
|
if i == j:
|
|
continue
|
|
scale = a.get(key, 0)
|
|
if not scale:
|
|
continue
|
|
for target, source in zip((a, b), rows[j]):
|
|
for k, v in source.items():
|
|
value = target.get(k, 0) - scale*v
|
|
if abs(value) < 1e-14:
|
|
target.pop(k, None)
|
|
else:
|
|
target[k] = value
|
|
pivots.append(key)
|
|
lines = [f'double b{i} = {expr};' for i, (_, expr) in enumerate(equations)
|
|
if any(i in rows[j][1] for j in range(len(pivots)))]
|
|
for j, key in enumerate(pivots):
|
|
expr = ' + '.join(f'({num(v)})*b{i}' for i, v in rows[j][1].items()) or '0.0'
|
|
lines.append(f'{key} = {expr};')
|
|
return lines
|
|
|
|
|
|
def compile_extended_program(network):
|
|
components = list(network.components.values())
|
|
contracts = catalog_contracts()
|
|
for c in components:
|
|
if type(c) is not contracts.get(c.model_type):
|
|
raise NativeCapabilityError(f'{c.name}: no native contract for {c.model_type} / {type(c).__name__}')
|
|
if c.MODEL_VERSION != SUPPORTED_VERSIONS.get(c.model_type):
|
|
raise NativeCapabilityError(f'{c.name}: native kernel version does not match the component contract')
|
|
variables = tuple(v for c in components for v in c.result_variable_metadata())
|
|
if not variables:
|
|
raise NativeCapabilityError('Simulation requires a runtime component; medium definitions alone have no outputs')
|
|
slots = {v.key: i for i, v in enumerate(variables)}
|
|
ports = {(c.name, p.name): p for c in components for p in c.active_port_definitions}
|
|
groups = _Groups(ports)
|
|
adjacent = {}
|
|
for edge in network.connections:
|
|
a, b = (p.key for p in edge.endpoints)
|
|
if edge.kind == 'physical':
|
|
adjacent[a], adjacent[b] = b, a
|
|
groups.union(a, b)
|
|
else:
|
|
source, target = (a, b) if ports[a].nominal_role == 'output' else (b, a)
|
|
adjacent[target] = source
|
|
for c in components:
|
|
for port in c.required_connection_ports:
|
|
if (c.name, port) not in adjacent:
|
|
raise NativeCapabilityError(f'{c.name}.{port}: unconnected required port')
|
|
names = [p.name for p in c.active_port_definitions]
|
|
if c.model_type in NODES | {'amesim_pnch023', 'amesim_pnch012', 'amesim_mecmas21', 'amesim_lmechn1'} or (c.model_type=='pipe' and c.lambda_darcy==0):
|
|
for p in names[1:]:
|
|
groups.union((c.name, names[0]), (c.name, p))
|
|
if c.model_type == 'amesim_pnrp17':
|
|
for a, b in [('port_2', 'port_5'), ('port_3', 'port_4')]:
|
|
groups.union((c.name, a), (c.name, b))
|
|
|
|
state_keys, initial, gas_initializers = [], [], []
|
|
states = {}
|
|
def add_states(c, fields, values):
|
|
for field, value in zip(fields, values):
|
|
states[c.name, field] = len(state_keys)
|
|
state_keys.append(f'{c.name}.{field}')
|
|
initial.append(float(value))
|
|
mass_groups = {}
|
|
for c in components:
|
|
if c.model_type in GAS_TYPES:
|
|
fields = ('m1', 'U1', 'm2', 'U2') if c.model_type == 'amesim_pnl0003' else ('m', 'U')
|
|
add_states(c, fields, [0.0]*len(fields))
|
|
elif c.model_type == 'amesim_mecmas21':
|
|
root = groups.find((c.name, 'port_1'))
|
|
mass_groups.setdefault(root, []).append(c)
|
|
if len(mass_groups[root]) == 1:
|
|
add_states(c, ('v', 'x'), (c.v0, c.x0))
|
|
else:
|
|
ref = mass_groups[root][0]
|
|
if abs(c.v0-ref.v0)>1e-10*max(abs(c.v0),1) or abs(c.x0-ref.x0)>1e-10*max(abs(c.x0),1):
|
|
raise NativeCapabilityError('Rigidly connected masses require consistent initial x/v')
|
|
for field in ('v', 'x'):
|
|
states[c.name, field] = states[ref.name, field]
|
|
if len(state_keys)>1024 or len(variables)>16384:
|
|
raise NativeCapabilityError('Native model exceeds the 1024-state / 16384-output resource limit')
|
|
# Algebraic models use an internal constant state; the public state map stays empty.
|
|
nstates = max(len(state_keys), 1)
|
|
if not initial:
|
|
initial = [0.0]
|
|
lines, declarations, breaks = [], [], []
|
|
assigned = set()
|
|
def w(c, field):
|
|
name = c if isinstance(c, str) else c.name
|
|
return f'w[{slots[name+"."+field]}]'
|
|
def put(c, field, expr, dest=None):
|
|
(lines if dest is None else dest).append(f'{w(c,field)} = {expr};')
|
|
assigned.add((c if isinstance(c,str) else c.name)+'.'+field)
|
|
def y(c, field):
|
|
return f'y[{states[c.name,field]}]'
|
|
def ep(c, port):
|
|
return c.name, port
|
|
def pnames(c):
|
|
return [p.name for p in c.active_port_definitions if p.domain == 'pneumatic']
|
|
def mnames(c):
|
|
return [p.name for p in c.active_port_definitions if p.domain == 'mechanical']
|
|
for c in components:
|
|
if c.model_type == 'amesim_mecmas21':
|
|
for f in ('v', 'x'):
|
|
put(c, f, y(c, f))
|
|
for name in mnames(c):
|
|
root = groups.find(ep(c, name))
|
|
if root not in mass_groups:
|
|
raise NativeCapabilityError(f'{c.name}.{name}: mechanical group has no inertia anchor')
|
|
mass = mass_groups[root][0]
|
|
for f in ('v', 'x'):
|
|
put(c, name+'.'+f, y(mass, f))
|
|
if c.model_type == 'amesim_step0':
|
|
put(c, 'y', f't < {num(c.time)} ? {num(c.initial)} : {num(c.final)}')
|
|
breaks.append(f'if(t < {num(c.time)}) result=fmin(result,{num(c.time)});')
|
|
elif c.model_type == 'amesim_ud00':
|
|
ident = 'signal_'+str(len(declarations))
|
|
declarations.append(f'static const double {ident}[24] = {{'+','.join(num(v) for vs in (c.starts,c.ends,c.durations) for v in vs)+'};')
|
|
args=f'{num(c.tstart)}, {c.nstages}, {int(c.iscyclic)}, {ident}'
|
|
put(c, 'y', f'native_signal(t, {args})')
|
|
breaks.append(f'result=fmin(result,native_signal_break(t,end,{args}));')
|
|
else:
|
|
continue
|
|
put(c, 'out.signal', w(c, 'y'))
|
|
for c in components:
|
|
for port in c.active_port_definitions:
|
|
if port.kind == 'signal' and port.nominal_role == 'input':
|
|
target = adjacent.get(ep(c, port.name))
|
|
if target:
|
|
expr = w(target[0], target[1]+'.signal')
|
|
elif c.model_type == 'amesim_pnvo001':
|
|
expr = num(c.opening0)
|
|
else:
|
|
raise NativeCapabilityError(f'{c.name}.{port.name}: signal input missing')
|
|
put(c, port.name+'.signal', expr)
|
|
|
|
media = {}
|
|
def medium(c):
|
|
m = c.medium
|
|
key = (getattr(m, 'SUBSTANCE_ID', None), getattr(m, 'PROPERTY_METHOD_ID', None),
|
|
m.R_gas, m.cp_ref, m.T_ref, m.cp_slope, m.viscosity_ref, m.viscosity_T_ref, m.sutherland_constant)
|
|
if key not in media:
|
|
if key[1] not in (None, 'ideal_gas', 'peng_robinson') or (key[1]=='peng_robinson' and key[0]!='helium'):
|
|
raise NativeCapabilityError(f'{c.name}: unsupported medium contract {key[:2]}')
|
|
ident = f'medium_{len(media)}'
|
|
declarations.append(f'static const NativeMedium {ident} = {{'+str(int(key[1]=='peng_robinson'))+','+','.join(num(v) for v in key[2:])+'};')
|
|
media[key]=ident
|
|
return '&'+media[key]
|
|
|
|
pistons = [c for c in components if c.model_type == 'amesim_pnrp17']
|
|
for c in pistons:
|
|
put(c, 'length', f'{num(c.x0)}+{w(c,"port_5.x")}-{w(c,"port_4.x")}')
|
|
put(c, 'volume', f'{num(c.effective_area)}*{w(c,"length")}')
|
|
put(c, 'volume_flow', f'{num(c.effective_area)}*({w(c,"port_5.v")}-{w(c,"port_4.v")})')
|
|
|
|
pneu = [endpoint for endpoint, port in ports.items() if port.domain == 'pneumatic']
|
|
pi = {endpoint: i for i, endpoint in enumerate(pneu)}
|
|
pgroups = list(dict.fromkeys(groups.find(endpoint) for endpoint in pneu))
|
|
pgi = {root: i for i, root in enumerate(pgroups)}
|
|
def p(c, name): return f'p[{pgi[groups.find(ep(c,name))]}]'
|
|
def q(c, name): return f'q[{pi[ep(c,name)]}]'
|
|
def h(c, name): return f'h[{pi[ep(c,name)]}]'
|
|
def hin(c, name, temperature=False):
|
|
other = adjacent[ep(c,name)]
|
|
obj = network.components[other[0]]
|
|
if temperature and obj.model_type in NODES-{'tee'}:
|
|
other = adjacent[(obj.name, 'port_2')]
|
|
return f'h[{pi[other]}]'
|
|
gases, anchor, port_gas, volume_rate = {}, {}, {}, {}
|
|
anchor_partitions = {}
|
|
gas_count = 0
|
|
for c in components:
|
|
if c.model_type not in GAS_TYPES:
|
|
continue
|
|
kind = c.model_type
|
|
if kind == 'amesim_pnch023':
|
|
volume, rate = num(c.cvol), '0.0'
|
|
elif kind == 'amesim_pnch012':
|
|
attached = [d for d in pistons if adjacent[ep(d,'port_1')][0] == c.name]
|
|
volume = num(c.cvol0+sum(c.external_volumes.values()))+''.join('+'+w(d,'volume') for d in attached)
|
|
put(c, 'vol', f'fmax({num(c.cvol0/100)}, {volume})')
|
|
rate = num(sum(c.external_volume_rates.values()))+''.join('+'+w(d,'volume_flow') for d in attached)
|
|
put(c, 'dvol', f'{w(c,"vol")} <= {num(c.cvol0/100)} ? 0.0 : ({rate})')
|
|
volume, rate = w(c,'vol'), w(c,'dvol')
|
|
elif kind in ('cylinder', 'tank'):
|
|
volume, rate = num(c.V), '0.0'
|
|
else:
|
|
volume, rate = num(c.volume), '0.0'
|
|
volume_rate[c.name] = rate
|
|
halves = (1,2) if kind == 'amesim_pnl0003' else (0,)
|
|
for half in halves:
|
|
suffix = str(half) if half else ''
|
|
gas = f'g[{gas_count}]';gas_count += 1
|
|
gases[c.name,half] = gas
|
|
V = num(c.compliance_volume) if half else volume
|
|
if half:
|
|
p0,T0=c.parameter_values[f'p{half}_0'],c.parameter_values[f'T{half}_0']
|
|
initial_volume=c.compliance_volume
|
|
else:
|
|
p0,T0=c.parameter_values['p0'],c.parameter_values['T0']
|
|
initial_volume=(max(c.cvol0+sum(c.external_volumes.values()),c.cvol0/100)
|
|
if kind=='amesim_pnch012' else
|
|
c.cvol if kind=='amesim_pnch023' else c.V if kind in ('cylinder','tank') else c.volume)
|
|
gas_initializers.append(f'if(!native_medium_init({medium(c)},{num(p0)},{num(T0)},{num(initial_volume)},{int(kind in ("cylinder","tank"))},&y[{states[c.name,"m"+suffix]}])) return 0;')
|
|
lines.append(f'if(!native_medium_gas({medium(c)}, {y(c,"m"+suffix)}, {y(c,"U"+suffix)}, {V}, &{gas})) return 0;')
|
|
for field in ('m', 'U'):
|
|
put(c, field+suffix, y(c, field+suffix))
|
|
for field in ('p','T','rho','u','h'):
|
|
put(c, field+suffix, gas+'.'+field)
|
|
anchored = pnames(c)
|
|
if kind == 'amesim_pnl0001': anchored=['port_2']
|
|
if kind == 'amesim_pnl0002': anchored=[]
|
|
if half: anchored=['port_'+suffix]
|
|
for name in anchored:
|
|
root=groups.find(ep(c,name))
|
|
anchor[root]=gas
|
|
partition=(c.name,half)
|
|
anchor_partitions.setdefault(root, {})[partition]=(c, float(c.compliance_volume if half else c.volume) if kind in ('amesim_pnl0001','amesim_pnl0003') else None)
|
|
for name in (['port_'+suffix] if half else pnames(c)):
|
|
port_gas[ep(c,name)]=gas
|
|
coupled=[]
|
|
project=[]
|
|
for root,partitions in anchor_partitions.items():
|
|
if len(partitions)<2: continue
|
|
if any(volume is None for _,volume in partitions.values()) or len({id(c.medium) for c,_ in partitions.values()})!=1:
|
|
raise NativeCapabilityError('Direct gas-storage coupling requires compatible fixed pipe compliances; insert a resistance between independent chambers')
|
|
offsets=[];volumes=[]
|
|
for (name,half),(c,volume) in partitions.items():
|
|
offsets.append(states[name,'m'+(str(half) if half else '')]);volumes.append(volume)
|
|
for prefix in ('p','T'):
|
|
values=[c.parameter_values[f'{prefix}{half}_0'] if half else getattr(c,prefix+'0')
|
|
for (_,half),(c,_) in partitions.items()]
|
|
if max(values)-min(values)>1e-9*max(1,*map(abs,values)):
|
|
raise NativeCapabilityError('Ideally coupled pipe compliances require consistent initial pressure and temperature')
|
|
total=sum(volumes)
|
|
project.append('{ double mass='+ '+'.join(f'y[{i}]' for i in offsets)+',energy='+ '+'.join(f'y[{i+1}]' for i in offsets)+';')
|
|
for i,volume in zip(offsets,volumes):
|
|
project += [f'projected[{i}]=mass*{num(volume/total)};projected[{i+1}]=energy*{num(volume/total)};']
|
|
project.append('}')
|
|
coupled.append((root,offsets,volumes))
|
|
for root, gas in anchor.items():
|
|
lines.append(f'p[{pgi[root]}]={gas}.p;')
|
|
default_h = next(iter(gases.values()))+'.h' if gases else '0.0'
|
|
for endpoint, idx in pi.items():
|
|
lines.append(f'h[{idx}]={port_gas.get(endpoint,default_h.removesuffix(".h"))}.h;' if endpoint in port_gas else f'h[{idx}]={default_h};')
|
|
for c in components:
|
|
if c.model_type in {'amesim_pnpl01','amesim_pnrp17'}:
|
|
# Closed-end ports retain their declared zero outflow enthalpy.
|
|
lines.append(f'{h(c,"port_1")}=0.0;')
|
|
|
|
flow_lines, flow_known, flow_eq = [], set(), []
|
|
def flow(c, name, expr):
|
|
target=q(c,name);flow_lines.append(f'{target}={expr};');flow_known.add(target)
|
|
def flow_equation(terms):
|
|
flow_eq.append(({q(c,name):coef for c,name,coef in terms},'0.0'))
|
|
for c in components:
|
|
kind=c.model_type
|
|
names=pnames(c)
|
|
if kind in RESISTORS or kind in NODES:
|
|
flow_equation([(c,name,1) for name in names])
|
|
if kind in {'amesim_pnpl01','amesim_pnrp17'}:
|
|
flow(c,'port_1','0.0')
|
|
if kind in RESISTORS:
|
|
a,b=names
|
|
pa,pb=p(c,a),p(c,b)
|
|
if kind=='pipe' and c.lambda_darcy==0:
|
|
# A zero-loss pipe is an ideal connection, whose flow is
|
|
# determined by the neighbouring constitutive equations.
|
|
continue
|
|
if kind=='orifice':
|
|
flow(c,a,f'{num(c.K_eff)}*copysign(sqrt(fabs({pa}-{pb})),{pa}-{pb})')
|
|
elif kind=='pipe':
|
|
resistance=c.lambda_darcy*c.L/c.D
|
|
flow(c,a,f'copysign(sqrt(fabs({pa}-{pb})*2*fmax(native_density({medium(c)},fmax(.5*({pa}+{pb}),1),{num(c.T0)}),1e-12)*{num(c.area*c.area/resistance)}),{pa}-{pb})')
|
|
elif kind=='amesim_pnl00r':
|
|
T=f'native_temperature_ph({medium(c)},fmax(fmax({pa},{pb}),1),{pa}>={pb}?{hin(c,a,True)}:{hin(c,b,True)})'
|
|
flow(c,a,f'native_pipe_flow({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'
|
|
if kind!='amesim_pnor001':
|
|
value = f'fmax(0,fmin(1,{w(c,"res.signal")}))' if kind=='amesim_pnvo001' else num(c.opening)
|
|
put(c,'xv',value)
|
|
area = c.effective_cq*(c.effective_area if kind=='amesim_pnor001' else c.maximum_area)
|
|
flow_lines.append(f'if(!native_medium_orifice({medium(c)},{pa},{pb},{hin(c,a)},{hin(c,b)},{num(area)},{opening},&{q(c,a)},&{w(c,"cm")},&{w(c,"gasvel")})) return 0;')
|
|
flow_known.add(q(c,a));assigned.update((c.name+'.cm',c.name+'.gasvel'))
|
|
flow(c,b,f'-{q(c,a)}')
|
|
elif kind in ('amesim_pnl0001','amesim_pnl0002'):
|
|
gas=gases[c.name,0]
|
|
for name in (['port_1'] if kind.endswith('1') else names):
|
|
T=gas+'.T'
|
|
if kind.endswith('2'):
|
|
T=f'({p(c,name)}>{gas}.p?native_temperature_ph({medium(c)},fmax({p(c,name)},1),{hin(c,name,True)}):{gas}.T)'
|
|
flow(c,name,f'native_pipe_flow({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':
|
|
flow_eq.append(({f'q[{pi[e.key]}]':1 for e in edge.endpoints},'0.0'))
|
|
unknownq=[f'q[{i}]' for i in range(len(pneu)) if f'q[{i}]' not in flow_known]
|
|
reduced=[]
|
|
for terms,rhs in flow_eq:
|
|
known=''.join(f'-({num(v)})*{k}' for k,v in terms.items() if k in flow_known)
|
|
reduced.append(({k:v for k,v in terms.items() if k not in flow_known},rhs+known))
|
|
free_flows={f'q[{pi[e]}]' for root,_,_ in coupled for e in pneu if groups.find(e)==root}
|
|
flow_lines += linear_schedule(reduced,unknownq,free_flows)
|
|
|
|
unknownp=[root for root in pgroups if root not in anchor]
|
|
residuals=[]
|
|
for root in unknownp:
|
|
terms=[f'q[{pi[e]}]' for e in pneu if groups.find(e)==root and f'q[{pi[e]}]' in flow_known]
|
|
if not terms:
|
|
raise NativeCapabilityError('Unanchored pneumatic pressure group has no constitutive flow relation')
|
|
residuals.append(' + '.join(terms))
|
|
# A monotone nodal mass-balance solve. This is an algebraic connection
|
|
# closure; the integration solver's Jacobian policy is unchanged.
|
|
pressure_lines=[]
|
|
if unknownp:
|
|
if not anchor:
|
|
raise NativeCapabilityError('Pneumatic pressure network has no storage pressure anchor')
|
|
pressure_lines += ['double plo=INFINITY,phi=0;', *[f'plo=fmin(plo,{g}.p);phi=fmax(phi,{g}.p);' for g in anchor.values()]]
|
|
pressure_lines += [f'p[{pgi[root]}]=.5*(plo+phi);' for root in unknownp]
|
|
pressure_lines += ['int pressure_ok=0;', 'for(int sweep=0;sweep<256;sweep++) {']
|
|
for root,expr in zip(unknownp,residuals):
|
|
pressure_lines += ['{ double lo=plo,hi=phi;', 'for(int bisect=0;bisect<48;bisect++) {',f'p[{pgi[root]}]=.5*(lo+hi);',
|
|
'if(!model_flows(p,h,g,w,q)) return 0;', f'if(({expr})>0) hi=p[{pgi[root]}];else lo=p[{pgi[root]}];','}}']
|
|
pressure_lines += ['if(!model_flows(p,h,g,w,q)) return 0;', 'double residual=0;', *[f'residual=fmax(residual,fabs({expr}));' for expr in residuals],
|
|
'if(residual<1e-11) { pressure_ok=1;break; }','}', 'if(!pressure_ok) return 0;']
|
|
else:
|
|
pressure_lines=['if(!model_flows(p,h,g,w,q)) return 0;']
|
|
|
|
stream_lines=[]
|
|
for c in components:
|
|
names=pnames(c);kind=c.model_type
|
|
if kind in {'amesim_pnpl01','amesim_pnrp17'}:
|
|
stream_lines.append(f'{h(c,"port_1")}={hin(c,"port_1")};')
|
|
elif kind in RESISTORS:
|
|
a,b=names
|
|
stream_lines += [f'{h(c,a)}={hin(c,b)};',f'{h(c,b)}={hin(c,a)};']
|
|
elif kind in NODES:
|
|
stream_lines += ['{ double total=0,energy=0,average=0;', *[f'if({q(c,name)}>1e-12) {{total+={q(c,name)};energy+={q(c,name)}*{hin(c,name)};}} average+={hin(c,name)};' for name in names]]
|
|
if kind=='tee':
|
|
stream_lines += [f'double mixed=total>1e-12?energy/total:average/{len(names)};', *[f'{h(c,name)}=mixed;' for name in names]]
|
|
else:
|
|
stream_lines += [f'double ref={hin(c,"port_2")},mixed=total>1e-12?energy/total:ref;', *[f'{h(c,name)}=ref;' for name in names if name!='port_2'],f'{h(c,"port_2")}=mixed;',f'if({q(c,"port_2")}<0) {{ double e=0,scale=0;',
|
|
*[f'e+={q(c,name)}*({q(c,name)}>1e-12?{hin(c,name)}:ref);scale+=fabs({q(c,name)});' for name in names if name!='port_2'],
|
|
f'double flow={q(c,"port_2")},transition=fmax(.05*scale,1e-12);',
|
|
'double inv=-flow>=transition?1/flow:flow*(2*transition*transition-flow*flow)/pow(transition,4);',
|
|
f'{h(c,"port_2")}=mixed-(e+flow*mixed)*inv;','}']
|
|
# Nodes expose their reference temperature as an output too.
|
|
if c.name+'.T' in slots:
|
|
raise NativeCapabilityError('Unexpected node temperature output contract')
|
|
stream_lines += ['}']
|
|
if pneu:
|
|
lines += ['int closure_ok=0;',f'for(int closure=0;closure<{max(64,4*len(pneu))};closure++) {{',f'double previous[{len(pneu)}];',f'for(int i=0;i<{len(pneu)};i++) previous[i]=h[i];',*pressure_lines,*stream_lines,
|
|
'double change=0;',f'for(int i=0;i<{len(pneu)};i++) change=fmax(change,fabs(h[i]-previous[i])/fmax(1,fabs(h[i])));',
|
|
'if(change<1e-12) {closure_ok=1;break;}','}', 'if(!closure_ok) return 0;', 'if(!model_flows(p,h,g,w,q)) return 0;']
|
|
for c in components:
|
|
for name in pnames(c):
|
|
for field,expr in [('p',p(c,name)),('m_flow',q(c,name)),('h_outflow',h(c,name))]:
|
|
put(c,name+'.'+field,expr)
|
|
|
|
# Mechanical force balance includes shared accelerations for rigid groups.
|
|
feq=[]
|
|
for edge in network.connections:
|
|
if edge.domain=='mechanical':
|
|
feq.append(({w(e.component,e.port+'.f'):1 for e in edge.endpoints},'0.0'))
|
|
for c in components:
|
|
kind=c.model_type
|
|
if kind=='amesim_forc':
|
|
put(c,'force',f'{num(c.direction)}*{w(c,"res.signal")}')
|
|
feq.append(({w(c,'port_2.f'):1},'-'+w(c,'force')))
|
|
elif kind=='amesim_f000':
|
|
feq.append(({w(c,'port_1.f'):1},'0.0'))
|
|
elif kind=='amesim_lstp00a':
|
|
put(c,'gap',f'{num(c.gap0)}+{w(c,"port_2.x")}-{w(c,"port_1.x")}')
|
|
put(c,'penetration',f'fmax(-{w(c,"gap")},0)')
|
|
put(c,'force',f'native_contact({w(c,"penetration")},{w(c,"port_1.v")}-{w(c,"port_2.v")},{num(c.kcont)},{num(c.rcont)},{num(c.Pdis)},{int(c.discContactOption)})')
|
|
feq += [({w(c,'port_1.f'):1},w(c,'force')),({w(c,'port_2.f'):1},'-'+w(c,'force'))]
|
|
elif kind=='amesim_pnrp17':
|
|
put(c,'pressure_force',f'({w(c,"port_1.p")}-101300)*{num(c.effective_area)}')
|
|
feq += [({w(c,'port_2.f'):1,w(c,'port_5.f'):1},'-'+w(c,'pressure_force')),({w(c,'port_3.f'):1,w(c,'port_4.f'):1},w(c,'pressure_force'))]
|
|
elif kind=='amesim_lmechn1':
|
|
feq.append(({w(c,name+'.f'):1 for name in mnames(c)},'0.0'))
|
|
elif kind=='amesim_mecmas21':
|
|
v,x=w(c,'v'),w(c,'x')
|
|
put(c,'Fvisc',f'-{num(c.rvisc)}*{v}' if c.use_friction else '0.0')
|
|
put(c,'Ffric',f'{v}>0?-{num(c.fcoul)}:({v}<0?{num(c.fcoul)}:0)' if c.use_friction else '0.0')
|
|
for field,penetration,velocity,suffix in [('Fmin',num(c.xmin)+'-'+x,'-'+v,'min'),('Fmax',x+'-'+num(c.xmax),v,'max')]:
|
|
expr=f'native_limit_force({penetration},{velocity},{num(getattr(c,"Kb"+suffix))},{num(getattr(c,"Db"+suffix))},{num(getattr(c,"Pd"+suffix))},{int(c.discContactOption)})' if int(c.stoptype)==2 else '0.0'
|
|
put(c,field,expr)
|
|
ref=mass_groups[groups.find(ep(c,'port_1'))][0]
|
|
extra=f'{w(c,"Fvisc")}+{w(c,"Ffric")}+{w(c,"Fmin")}-{w(c,"Fmax")}'
|
|
if c.use_friction: extra+=f'-{num(c.wind)}*{v}*fabs({v})'
|
|
feq.append(({w(c,'port_1.f'):1,w(c,'port_2.f'):1,w(ref,'a'):-c.mass},f'-({extra})'))
|
|
unknownf=[w(c,name+'.f') for c in components for name in mnames(c)]+[w(group[0],'a') for group in mass_groups.values()]
|
|
lines += linear_schedule(feq,unknownf)
|
|
for c in components:
|
|
for name in mnames(c): assigned.add(c.name+'.'+name+'.f')
|
|
if c.model_type=='amesim_lmechn1':
|
|
put(c,'tforce',' + '.join(w(c,name+'.f') for name in mnames(c)[:-1]))
|
|
stops=[]
|
|
for group in mass_groups.values():
|
|
ref=group[0];vi=states[ref.name,'v'];xi=states[ref.name,'x']
|
|
limits=[c for c in group if int(c.stoptype) in (1,3)]
|
|
lines += [f'dy[{vi}]={w(ref,"a")};dy[{xi}]={y(ref,"v")};']
|
|
if limits:
|
|
lower=max(c.xmin for c in limits);upper=min(c.xmax for c in limits)
|
|
if lower>upper or ref.x0<lower-1e-12 or ref.x0>upper+1e-12:
|
|
raise NativeCapabilityError('Inconsistent discrete endstop bounds or initial position')
|
|
restitution, thresholds = [], []
|
|
for parameter, bound in [('xmin',lower),('xmax',upper)]:
|
|
active=[c for c in limits if abs(getattr(c,parameter)-bound)<=1e-12*max(abs(bound),1)]
|
|
restitution.append(0 if any(int(c.stoptype)==1 for c in active) else min(c.restcoeff for c in active))
|
|
thresholds.append(max((c.restdvel for c in active if int(c.stoptype)==3),default=0))
|
|
stops.append((vi,lower,upper,*restitution,*thresholds))
|
|
lines.append(f'native_stop_motion({y(ref,"x")},{y(ref,"v")},{num(lower)},{num(upper)},&dy[{vi}],&dy[{xi}]);')
|
|
for c in group: put(c,'a',f'dy[{vi}]')
|
|
|
|
for c in components:
|
|
kind=c.model_type
|
|
if kind not in GAS_TYPES and kind!='amesim_pnl00r': continue
|
|
if kind in GAS_TYPES:
|
|
halves=(1,2) if kind=='amesim_pnl0003' else (0,)
|
|
center='0.0'
|
|
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)')
|
|
for half in halves:
|
|
gas=gases[c.name,half];suffix=str(half) if half else ''
|
|
names=['port_'+suffix] if half else pnames(c)
|
|
mass=' + '.join(q(c,name) for name in names)
|
|
energy=' + '.join(f'{q(c,name)}*({q(c,name)}>0?{hin(c,name)}:{gas}.h)' for name in names)
|
|
if half:
|
|
sign='-' if half==1 else '+'
|
|
mass+=sign+center
|
|
energy+=f'{sign}{center}*({center}>0?{gases[c.name,1]}.h:{gases[c.name,2]}.h)'
|
|
heat='0.0'
|
|
if kind.startswith('amesim_pnch'):
|
|
heat=f'{num(c.kth*c.sth)}*({num(c.extemp)}-{gas}.T)-{gas}.p*({volume_rate[c.name]})'
|
|
elif kind.startswith('amesim_pnl') and int(c.mode)!=1:
|
|
temp=f'.5*({gases[c.name,1]}.T+{gases[c.name,2]}.T)' if half else gas+'.T'
|
|
heat=f'{num(c.kth*c.exchange_area/(2 if half else 1))}*({num(c.extemp)}-({temp}))'
|
|
lines += [f'dy[{states[c.name,"m"+suffix]}]={mass};',f'dy[{states[c.name,"U"+suffix]}]={energy}+({heat});']
|
|
if kind.startswith('amesim_pnl'):
|
|
diag=[]
|
|
if kind=='amesim_pnl0002':
|
|
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)'
|
|
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')
|
|
diag=[(flow,f'fmax(fmax({a}.p,{b}.p),1)',f'({flow}>=0?{a}.T:{b}.T)',c.le,1)]
|
|
else:
|
|
pa,pb=p(c,'port_1'),p(c,'port_2');pp=f'fmax(fmax({pa},{pb}),1)'
|
|
temp=gases[c.name,0]+'.T' if kind=='amesim_pnl0001' else f'fmax(native_temperature_ph({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);')
|
|
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')):
|
|
expr=f'acc[{i}]/{len(diag)}'
|
|
put(c,field,f'fmin({expr},64000000)' if field=='ff' else expr)
|
|
lines.append('}')
|
|
for _,offsets,volumes in coupled:
|
|
lines.append('{ double mass='+ '+'.join(f'dy[{i}]' for i in offsets)+',energy='+ '+'.join(f'dy[{i+1}]' for i in offsets)+';')
|
|
for i,volume in zip(offsets,volumes):
|
|
lines.append(f'dy[{i}]=mass*{num(volume/sum(volumes))};dy[{i+1}]=energy*{num(volume/sum(volumes))};')
|
|
lines.append('}')
|
|
missing=set(slots)-assigned
|
|
if missing: raise NativeCapabilityError(f'Native output mapping incomplete: {sorted(missing)}')
|
|
if not state_keys: lines.append('dy[0]=0;')
|
|
np,ng,nq=max(1,len(pgroups)),max(1,gas_count),max(1,len(pneu))
|
|
source='\n'.join(['#include "model.h"','#include <math.h>',*declarations,
|
|
f'const NativeStop model_stops[{max(1,len(stops))}] = {{'+(','.join('{'+str(s[0])+','+','.join(num(v) for v in s[1:])+'}' for s in stops) or '{0,0,0,0,0,0,0}')+'};',
|
|
'const double model_atol[NSTATES] = {'+','.join('1e-12' if k.endswith(('.v','.x')) else '1e-8' for k in state_keys or ['dummy'])+'};',
|
|
'const char *const model_output_keys[NOUTPUTS] = {'+(','.join(json.dumps(v.key,ensure_ascii=True) for v in variables) or '""')+'};',
|
|
'static int model_flows(const double *p,const double *h,const NativeGas *g,double *w,double *q) {',
|
|
'(void)p;(void)h;(void)g;(void)w;(void)q;',*flow_lines,'return 1;}',
|
|
'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) {',
|
|
*(['double projected[NSTATES];for(int i=0;i<NSTATES;i++) projected[i]=y[i];',*project,'y=projected;'] if project else []),
|
|
f'double p[{np}]={{0}},h[{nq}]={{0}},q[{nq}]={{0}};NativeGas g[{ng}];',
|
|
'(void)t;(void)y;(void)w;(void)p;(void)h;(void)q;(void)g;(void)model_flows;',*lines,
|
|
'for(int i=0;i<NSTATES;i++) if(!isfinite(dy[i])) return 0;',
|
|
'for(int i=0;i<NOUTPUTS;i++) if(!isfinite(w[i])) return 0;','return 1;}',
|
|
'double model_next_break(double t,double end) { double result=end;(void)t;',*breaks,'return result;}',''])
|
|
header=f'''#ifndef GENERATED_NATIVE_MODEL_H
|
|
#define GENERATED_NATIVE_MODEL_H
|
|
#include "kernels.h"
|
|
#define NSTATES {nstates}
|
|
#define NOUTPUTS {max(1,len(variables))}
|
|
#define NSTOPS {len(stops)}
|
|
extern const NativeStop model_stops[{max(1,len(stops))}];
|
|
extern const double model_atol[NSTATES];
|
|
extern const char *const model_output_keys[NOUTPUTS];
|
|
int model_init(double *y);
|
|
int model_eval(double t,const double *y,double *dy,double *w);
|
|
double model_next_break(double t,double end);
|
|
#endif
|
|
'''
|
|
return NativeProgram(source,header,tuple(state_keys),variables,tuple(sorted({c.model_type for c in components})))
|