"""Conservation regressions for reference nodes and prescribed/moving volumes.""" import json import math from itertools import product from pathlib import Path import subprocess import tempfile import unittest from app.simulation.components.amesim.media.mediums import AmesimIdealAirMedium from app.simulation.config import SolveIVPConfig from app.simulation.native_codegen.build import build_native, toolchain from app.simulation.native_codegen.compiler import compile_native_program from app.simulation.native_codegen.extended import compile_extended_program from app.simulation.native_codegen.runner import execute_native from tests.test_native_catalog import Circuit def probe(program, build, state=None, time=0): if state is None: state = json.loads(subprocess.check_output( [str(build.executable), '--init'], text=True, timeout=15)) return (*probes(program, build, [state], time)[0], state) def probes(program, build, states, time=0): run = subprocess.run([str(build.executable), '--probe'], input=''.join(' '.join(format(x, '.17g') for x in (time, *state))+'\n' for state in states), capture_output=True, text=True, check=True, timeout=30) rows = [json.loads(line) for line in run.stdout.splitlines()] if len(rows) != len(states) or not all(row['success'] for row in rows): raise AssertionError(run) return [(dict(zip((v.key for v in program.variables), row['outputs'])), dict(zip(program.state_keys, row['rhs']))) for row in rows] def flow_circuit(node_kind, consumer_kind='amesim_pnch023', *, nested=False, K=1e-6): b = Circuit(AmesimIdealAirMedium()) params = (dict(p1_0=2e5, p2_0=2e5, T1_0=400, T2_0=400, mode=2) if consumer_kind == 'amesim_pnl0003' else dict(p0=2e5, T0=400)) if consumer_kind == 'amesim_pnl0001': params['mode'] = 2 consumer = b.add(consumer_kind, 'consumer', **params) node = b.add(node_kind, 'node') reference_port = ('port_2' if consumer_kind == 'amesim_pnl0001' else next(p.name for p in consumer.active_port_definitions if p.domain == 'pneumatic')) b.connect(node, 'port_2', consumer, reference_port) if nested: child = b.add('amesim_pn3node2', 'child') b.connect(child, 'port_2', node, 'port_1') branch_ports = [(child, 'port_1'), (child, 'port_3'), (node, 'port_3')] else: branch_ports = [(node, p.name) for p in node.active_port_definitions if p.name != 'port_2'] for i, (junction, port) in enumerate(branch_ports): gas = b.chamber('branch'+str(i), p0=3e5, T0=(700, 300, 500)[i]) resistance = b.add('orifice', 'flow'+str(i), K=K) b.connect(gas, 'port_1', resistance, 'port_a') b.connect(resistance, 'port_b', junction, port) return b.seal(), len(branch_ports) def piston_circuit(route, *, prescribed=0, count=1): b = Circuit() gas = b.add('amesim_pnch012', 'gas', cvol0=.01, p0=2e5, T0=400, vol1=.001, dvol1=prescribed) targets = [(gas, 'port_'+str(i+1)) for i in range(count)] if route != 'direct': node = b.add('amesim_p4node2' if route == 'nested' else route, 'node') b.connect(node, 'port_2', gas, 'port_1') targets = [(node, port) for port in ('port_1', 'port_3', 'port_4')[:count]] if route == 'nested': child = b.add('amesim_pn3node2', 'child') b.connect(child, 'port_2', node, 'port_1') targets[0] = child, 'port_1' for index, (target, port) in enumerate(targets): piston = b.add('amesim_pnrp17', 'piston'+str(index), dp=.1, dr=0, x0=.2) b.connect(piston, 'port_1', target, port) for side, a, z, speed in [('moving', 'port_2', 'port_5', .1), ('body', 'port_3', 'port_4', 0)]: name = side+str(index) mass = b.add('amesim_mecmas21', name, mass=1, stoptype=4, useFriction=1, v0=speed, x0=0) free = b.add('amesim_f000', name+'free') end = b.add('amesim_f000', name+'end') b.connect(piston, a, mass, 'port_1') b.connect(mass, 'port_2', free, 'port_1') b.connect(piston, z, end, 'port_1') return b.seal() class NativeNodeConservationTests(unittest.TestCase): @classmethod def setUpClass(cls): toolchain() def assert_balance(self, rates): self.assertLessEqual(abs(sum(rates)), 1e-20+1e-10*sum(map(abs, rates))) def check_flows(self, net, count): program = compile_native_program(net) build = build_native(program) self.addCleanup(build.close) _, _, initial = probe(program, build) cases = [(1, -.5), (.5, -1), (1, -1), (1, -1+1e-8)] if count == 2 else [ (1, -.5, -.5), (.5, .5, -1), (1, -.5, -.5+1e-8)] cases += list(product((-1, 0, 1), repeat=count)) states = [] for flows in cases: state = list(initial) for i, flow in enumerate(flows): pressure = 2e5+math.copysign(1e5*flow*flow, flow) for field in ('m', 'U'): j = program.state_keys.index('branch'+str(i)+'.'+field) state[j] *= pressure/3e5 states.append(state) for flows, (outputs, rhs) in zip(cases, probes(program, build, states)): with self.subTest(flows=flows): for fields in (('m', 'm1', 'm2'), ('U', 'U1', 'U2')): self.assert_balance([v for k, v in rhs.items() if k.rsplit('.', 1)[1] in fields]) ref_h = outputs.get('consumer.h', outputs.get('consumer.h1')) expected = 0 for i in range(count): q = outputs['flow'+str(i)+'.port_a.m_flow'] expected += q*(outputs['branch'+str(i)+'.h'] if q > 0 else ref_h) actual = sum(v for k, v in rhs.items() if k.startswith('consumer.U')) self.assertAlmostEqual(actual, expected, delta=1e-20+abs(expected)*1e-10) if count == 2 and flows == (1, -1): self.assertLess(abs(outputs['node.port_2.m_flow']), 1e-12) self.assertGreater(abs(actual), 1) def test_reference_energy_for_both_nodes_and_all_storage_consumers(self): for kind in ('amesim_pn3node2', 'amesim_p4node2'): for consumer in ('amesim_pnch023', 'amesim_pnch012', 'amesim_pnl0001', 'amesim_pnl0003', 'cylinder', 'tank'): with self.subTest(node=kind, consumer=consumer): self.check_flows(*flow_circuit(kind, consumer)) def test_nested_reference_nodes_and_reversed_component_order(self): net, count = flow_circuit('amesim_pn3node2', nested=True) net.components = dict(reversed(list(net.components.items()))) net.connections.reverse() self.check_flows(net, count) def test_sub_picogram_flows_still_conserve_energy(self): net, count = flow_circuit('amesim_p4node2', K=1e-16) self.check_flows(net, count) def test_piston_volume_and_work_pass_through_nodes_once(self): for route in ('direct', 'amesim_pn3node2', 'amesim_p4node2', 'nested'): net = piston_circuit(route, prescribed=2e-4, count=2) for compiler in (compile_native_program, compile_extended_program): with self.subTest(route=route, compiler=compiler.__name__): program = compiler(net) build = build_native(program) self.addCleanup(build.close) outputs, rhs, state = probe(program, build) expected_v = .011+sum(outputs['piston'+str(i)+'.volume'] for i in range(2)) expected_dv = 2e-4+sum(outputs['piston'+str(i)+'.volume_flow'] for i in range(2)) self.assertAlmostEqual(outputs['gas.vol'], expected_v, places=13) self.assertAlmostEqual(outputs['gas.dvol'], expected_dv, places=13) self.assertAlmostEqual(outputs['gas.p'], 2e5, delta=1e-7) self.assertAlmostEqual(rhs['gas.U'], -outputs['gas.p']*expected_dv, places=8) # Move only one piston: its geometry already contains the # displacement and must not be integrated a second time. state[program.state_keys.index('moving0.x')] += .03 later, _, _ = probe(program, build, state, time=7) self.assertAlmostEqual(later['gas.vol']-outputs['gas.vol'], math.pi*.1**2/4*.03, places=13) def test_prescribed_rates_integrate_from_simulation_start_in_both_compilers(self): for compiler in (compile_native_program, compile_extended_program): for rate in (.002, -.002): b = Circuit() b.add('amesim_pnch012', 'gas', cvol0=.01, p0=2e5, T0=400, vol1=.002, vol3=.001, dvol1=rate*.75, dvol4=rate*.25) program = compiler(b.seal()) build = build_native(program) self.addCleanup(build.close) with tempfile.TemporaryDirectory() as directory: for method in ('RK45', 'BDF'): with self.subTest(compiler=compiler.__name__, rate=rate, method=method): cfg = SolveIVPConfig(t_start=5, t_stop=5.1, max_step=.01, rtol=1e-9, method=method) result = execute_native(build, cfg, .02, run_dir=Path(directory)/method) self.assertTrue(result['success'], result) series = result['series'] for t, volume, dvol in zip(series['time'], series['gas.vol'], series['gas.dvol']): self.assertAlmostEqual(volume, .013+rate*(t-5), places=12) self.assertAlmostEqual(dvol, rate, places=13) self.assertEqual(series['gas.m'][0], series['gas.m'][-1]) self.assertLess((series['gas.U'][-1]-series['gas.U'][0])*rate, 0) def test_volume_lower_limit_stops_work_and_allows_expansion_at_boundary(self): for compiler in (compile_native_program, compile_extended_program): for rate in (-1, 1): b = Circuit() # Exactly representable boundary: cvol0/100 == cvol0+vol1. b.add('amesim_pnch012', 'gas', cvol0=100, vol1=-99, dvol1=rate) program = compiler(b.seal()) build = build_native(program) self.addCleanup(build.close) outputs, rhs, state = probe(program, build) self.assertEqual(outputs['gas.vol'], 1) self.assertEqual(outputs['gas.dvol'], max(rate, 0)) self.assertAlmostEqual(rhs['gas.U'], -outputs['gas.p']*max(rate, 0)) state[program.state_keys.index('gas._volume_displacement')] = -2 below, below_rhs, _ = probe(program, build, state) self.assertEqual(below['gas.vol'], 1) self.assertEqual(below['gas.dvol'], 0) self.assertEqual(below_rhs['gas.U'], 0) def test_prescribed_ideal_gas_volume_obeys_adiabatic_law(self): medium = AmesimIdealAirMedium() gamma = medium.cp_ref/(medium.cp_ref-medium.R_gas) b = Circuit(medium) b.add('amesim_pnch012', 'gas', cvol0=.01, p0=2e5, T0=400, dvol1=.002) build = build_native(compile_native_program(b.seal())) self.addCleanup(build.close) with tempfile.TemporaryDirectory() as directory: for method in ('RK45', 'BDF'): cfg = SolveIVPConfig(t_stop=.2, max_step=.01, rtol=1e-9, method=method) result = execute_native(build, cfg, .02, run_dir=Path(directory)/method) self.assertTrue(result['success'], result) series = result['series'] for volume, pressure, temperature in zip(series['gas.vol'], series['gas.p'], series['gas.T']): self.assertAlmostEqual(pressure, 2e5*(.01/volume)**gamma, delta=.005) self.assertAlmostEqual(temperature, 400*(.01/volume)**(gamma-1), delta=1e-5) if __name__ == '__main__': unittest.main()