diff --git a/app/simulation/components/amesim/boundary/sources.py b/app/simulation/components/amesim/boundary/sources.py
index 4eb81aa..00c569a 100644
--- a/app/simulation/components/amesim/boundary/sources.py
+++ b/app/simulation/components/amesim/boundary/sources.py
@@ -46,6 +46,9 @@ class AmesimPnpl01(AlgebraicComponent):
     ) -> AmesimPnpl01:
         return cls(name=name)
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        return (self.port_1.m_flow,)
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         return (
             EquationResidual(
diff --git a/app/simulation/components/amesim/flow/orifices.py b/app/simulation/components/amesim/flow/orifices.py
index 2d6e0d1..5538774 100644
--- a/app/simulation/components/amesim/flow/orifices.py
+++ b/app/simulation/components/amesim/flow/orifices.py
@@ -410,6 +410,12 @@ class AmesimPnor001(AlgebraicComponent):
             "gasvel": flow_direction * gas_velocity,
         }
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        return (
+            self.port_1.m_flow + self.port_2.m_flow,
+            self.port_1.m_flow - self.mass_flow(self.port_1.p, self.port_2.p),
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         return (
             EquationResidual(
@@ -832,6 +838,12 @@ class AmesimPnvo001FixedOpening(AlgebraicComponent):
             "gasvel": flow_direction * gas_velocity,
         }
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        return (
+            self.port_2.m_flow + self.port_3.m_flow,
+            self.port_2.m_flow - self.mass_flow(self.port_2.p, self.port_3.p),
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         return (
             EquationResidual(
diff --git a/app/simulation/components/amesim/flow/pipes.py b/app/simulation/components/amesim/flow/pipes.py
index ce8c377..cff5720 100644
--- a/app/simulation/components/amesim/flow/pipes.py
+++ b/app/simulation/components/amesim/flow/pipes.py
@@ -317,6 +317,12 @@ class AmesimPnl00r(AlgebraicComponent):
             "ff": self.friction_factor(reynolds),
         }
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        return (
+            self.port_1.m_flow + self.port_2.m_flow,
+            self.port_1.m_flow - self.mass_flow(self.port_1.p, self.port_2.p),
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         return (
             EquationResidual(
@@ -838,6 +844,13 @@ class AmesimPnl0001(ThermodynamicVolumeComponent):
             "ff": self.friction_factor(reynolds),
         }
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume)
+        return (
+            self.port_2.p - props.p,
+            self.port_1.m_flow - self.mass_flow(self.port_1.p, props.p, props.T),
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume)
         return (
@@ -1080,6 +1093,25 @@ class AmesimPnl0002(AmesimPnl0001):
             "ff": friction,
         }
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume)
+        return (
+            self.port_1.m_flow
+            - self.port_mass_flow(
+                self.port_1.p,
+                props.p,
+                props.T,
+                port_name="port_1",
+            ),
+            self.port_2.m_flow
+            - self.port_mass_flow(
+                self.port_2.p,
+                props.p,
+                props.T,
+                port_name="port_2",
+            ),
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         props = self.medium.properties_from_mU(self.state.m, self.state.U, self.volume)
         return (
@@ -1411,6 +1443,14 @@ class AmesimPnl0003(DynamicComponent):
             "ff": self.friction_factor(reynolds),
         }
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        port_1 = self._properties(self.state_1)
+        port_2 = self._properties(self.state_2)
+        return (
+            self.port_1.p - port_1.p,
+            self.port_2.p - port_2.p,
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         port_1 = self._properties(self.state_1)
         port_2 = self._properties(self.state_2)
diff --git a/app/simulation/components/amesim/junctions/nodes.py b/app/simulation/components/amesim/junctions/nodes.py
index b9414dd..125f81a 100644
--- a/app/simulation/components/amesim/junctions/nodes.py
+++ b/app/simulation/components/amesim/junctions/nodes.py
@@ -27,6 +27,19 @@ class _AmesimPneumaticNode(AlgebraicComponent):
         for definition in self.PORTS:
             setattr(self, definition.name, self.register_declared_port(definition.name))
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        reference = self.get_port(self.REFERENCE_PORT)
+        return tuple(
+            self.get_port(definition.name).p - reference.p
+            for definition in self.PORTS
+            if definition.name != self.REFERENCE_PORT
+        ) + (
+            sum(
+                self.get_port(definition.name).m_flow
+                for definition in self.PORTS
+            ),
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         reference = self.get_port(self.REFERENCE_PORT)
         residuals: list[EquationResidual] = []
diff --git a/app/simulation/components/amesim/mechanical/pistons.py b/app/simulation/components/amesim/mechanical/pistons.py
index 0bd19d1..0d2456b 100644
--- a/app/simulation/components/amesim/mechanical/pistons.py
+++ b/app/simulation/components/amesim/mechanical/pistons.py
@@ -154,6 +154,22 @@ class AmesimPnrp17(AlgebraicComponent):
     def pressure_force(self) -> float:
         return (self.port_1.p - AMESIM_REFERENCE_PRESSURE_PA) * self.effective_area
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        values = [self.port_1.m_flow]
+        effort_pairs = (("port_2", "port_5"), ("port_3", "port_4"))
+        for first_name, second_name in effort_pairs:
+            first = self.get_port(first_name)
+            second = self.get_port(second_name)
+            values.extend((first.x - second.x, first.v - second.v))
+        force = self.pressure_force
+        values.extend(
+            (
+                self.port_2.f + self.port_5.f + force,
+                self.port_3.f + self.port_4.f - force,
+            )
+        )
+        return tuple(values)
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         effort_pairs = (("port_2", "port_5"), ("port_3", "port_4"))
         residuals: list[EquationResidual] = [
diff --git a/app/simulation/components/amesim/mechanical/translational.py b/app/simulation/components/amesim/mechanical/translational.py
index aa51dab..711c47a 100644
--- a/app/simulation/components/amesim/mechanical/translational.py
+++ b/app/simulation/components/amesim/mechanical/translational.py
@@ -63,6 +63,9 @@ class AmesimF000(AlgebraicComponent):
     ) -> "AmesimF000":
         return cls(name=name)
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        return (self.port_1.f,)
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         return (
             EquationResidual(
@@ -140,6 +143,11 @@ class AmesimForc(AlgebraicComponent):
     def output_force(self) -> float:
         return self.direction * float(self.res.signal)
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        return (
+            self.port_2.f + self.output_force,
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         return (
             EquationResidual(
@@ -591,6 +599,14 @@ class AmesimMecmas21(DynamicComponent):
             port.x = self.x
             port.v = self.v
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        return (
+            self.port_1.x - self.x,
+            self.port_1.v - self.v,
+            self.port_2.x - self.x,
+            self.port_2.v - self.v,
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         return (
             self._state_residual("port_1", "x", self.port_1.x - self.x),
@@ -992,6 +1008,13 @@ class AmesimLstp00a(AlgebraicComponent):
         self._causal_port_1_v = float(self.port_1.v)
         self._causal_port_2_v = float(self.port_2.v)
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        force = self.contact_force
+        return (
+            self.port_1.f - force,
+            self.port_2.f + force,
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         force = self.contact_force
         return (
@@ -1110,6 +1133,10 @@ class AmesimLmechn1(AlgebraicComponent):
     def active_ports(self) -> tuple[str, ...]:
         return tuple(f"port_{index}" for index in range(1, self.v1 + 2))
 
+    @property
+    def active_port_definitions(self) -> tuple[PortDefinition, ...]:
+        return self.PORTS[: self.v1 + 1]
+
     @property
     def reference_port_name(self) -> str:
         return f"port_{self.v1 + 1}"
@@ -1129,6 +1156,15 @@ class AmesimLmechn1(AlgebraicComponent):
     def force_balance(self) -> float:
         return self.total_force + self.get_port(self.reference_port_name).f
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        reference = self.get_port(self.reference_port_name)
+        values: list[float] = []
+        for port_name in self.active_ports[:-1]:
+            port = self.get_port(port_name)
+            values.extend((port.x - reference.x, port.v - reference.v))
+        values.append(self.force_balance)
+        return tuple(values)
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         reference_name = self.reference_port_name
         reference = self.get_port(reference_name)
@@ -1168,39 +1204,6 @@ class AmesimLmechn1(AlgebraicComponent):
                 value=self.force_balance,
             )
         )
-        for definition in self.PORTS[self.v1 + 1 :]:
-            port = self.get_port(definition.name)
-            residuals.extend(
-                (
-                    EquationResidual(
-                        id=f"{self.name}:{definition.name}_inactive_x",
-                        owner="component",
-                        owner_id=self.name,
-                        relation="constitutive",
-                        variables=(f"{self.name}.{definition.name}.x",),
-                        role="effort",
-                        value=port.x,
-                    ),
-                    EquationResidual(
-                        id=f"{self.name}:{definition.name}_inactive_v",
-                        owner="component",
-                        owner_id=self.name,
-                        relation="constitutive",
-                        variables=(f"{self.name}.{definition.name}.v",),
-                        role="effort",
-                        value=port.v,
-                    ),
-                    EquationResidual(
-                        id=f"{self.name}:{definition.name}_inactive_force",
-                        owner="component",
-                        owner_id=self.name,
-                        relation="constitutive",
-                        variables=(f"{self.name}.{definition.name}.f",),
-                        role="flow",
-                        value=port.f,
-                    ),
-                )
-            )
         return tuple(residuals)
 
     def component_result_values(self) -> Mapping[str, float]:
diff --git a/app/simulation/components/amesim/storage/chambers.py b/app/simulation/components/amesim/storage/chambers.py
index 0be73d0..15af67f 100644
--- a/app/simulation/components/amesim/storage/chambers.py
+++ b/app/simulation/components/amesim/storage/chambers.py
@@ -217,6 +217,17 @@ class AmesimPnch023(ThermodynamicVolumeComponent):
         )
         return derivative.as_vector()
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        pressure = self.medium.properties_from_mU(
+            self.state.m,
+            self.state.U,
+            self.cvol,
+        ).p
+        return (
+            self.port_1.p - pressure,
+            self.port_2.p - pressure,
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         pressure = self.medium.properties_from_mU(
             self.state.m,
@@ -507,6 +518,17 @@ class AmesimPnch012(ThermodynamicVolumeComponent):
         energy_derivative -= props.p * self.total_volume_rate()
         return VolumeState(m=mass_derivative, U=energy_derivative).as_vector()
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        pressure = self.medium.properties_from_mU(
+            self.state.m,
+            self.state.U,
+            self.total_volume(),
+        ).p
+        return tuple(
+            self.get_port(port_name).p - pressure
+            for port_name in ("port_1", "port_2", "port_3", "port_4")
+        )
+
     def pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
         pressure = self.medium.properties_from_mU(
             self.state.m,
diff --git a/app/simulation/core/base.py b/app/simulation/core/base.py
index dcdad2f..974d626 100644
--- a/app/simulation/core/base.py
+++ b/app/simulation/core/base.py
@@ -44,13 +44,19 @@ class Component(ABC):
             if port.definition is not None
         )
 
+    @property
+    def active_port_definitions(self) -> tuple[PortDefinition, ...]:
+        """Instance ports that participate in execution and result reporting."""
+
+        return self.port_definitions
+
     @property
     def required_connection_ports(self) -> tuple[str, ...]:
         """Physical ports that must have an external connection before simulation."""
 
         return tuple(
             definition.name
-            for definition in self.port_definitions
+            for definition in self.active_port_definitions
             if definition.kind == "physical"
         )
 
@@ -134,7 +140,7 @@ class Component(ABC):
                 )
             values[name] = float(component_values[name])
 
-        for port_definition in self.port_definitions:
+        for port_definition in self.active_port_definitions:
             port = self.get_port(port_definition.name)
             for variable in port_definition.variables:
                 if not variable.result_visible:
@@ -161,7 +167,7 @@ class Component(ABC):
             for definition in self.RESULT_VARIABLES
             if definition.visible
         ]
-        for port_definition in self.port_definitions:
+        for port_definition in self.active_port_definitions:
             for variable in port_definition.variables:
                 if not variable.result_visible:
                     continue
@@ -209,6 +215,20 @@ class Component(ABC):
 
         return ()
 
+    def pressure_flow_equation_values(self) -> tuple[float, ...]:
+        """Return live residual values in the declared equation order.
+
+        Components with frequently evaluated equations can override this
+        method to avoid rebuilding immutable equation metadata during closure.
+        The default keeps third-party components compatible with the public
+        residual API.
+        """
+
+        return tuple(
+            float(equation.value)
+            for equation in self.pressure_flow_equation_residuals()
+        )
+
     def update_stream_outflows(self, connected_h: Mapping[str, float]) -> None:
         """Update connector outflow properties from current flow directions."""
 
diff --git a/app/simulation/solvers/algebraic.py b/app/simulation/solvers/algebraic.py
index bd65286..db3b9d5 100644
--- a/app/simulation/solvers/algebraic.py
+++ b/app/simulation/solvers/algebraic.py
@@ -54,6 +54,23 @@ class ExplicitFlowAssignment:
 @dataclass(frozen=True)
 class ExplicitFlowStage:
     assignments: tuple[ExplicitFlowAssignment, ...]
+    direct_evaluations: tuple[tuple[int, Callable[[], float]], ...]
+    component_evaluations: tuple["ExplicitFlowComponentEvaluation", ...]
+
+
+@dataclass(frozen=True)
+class ExplicitFlowComponentEvaluation:
+    component: object
+    evaluate: Callable[[], tuple[float, ...]]
+    assignment_indices: tuple[int, ...]
+    equation_indices: tuple[int, ...]
+    equation_ids: tuple[str, ...]
+
+
+@dataclass(frozen=True)
+class EquationScalePlan:
+    variable_names: tuple[str, ...]
+    force_unknowns: tuple[AlgebraicUnknown, ...]
 
 
 @dataclass(frozen=True)
@@ -68,6 +85,13 @@ class ConnectionEquationEvaluation:
     evaluate: Callable[[], float]
 
 
+@dataclass(frozen=True)
+class ComponentEquationEvaluation:
+    component: object
+    evaluate: Callable[[], tuple[float, ...]]
+    templates: tuple[EquationResidual, ...]
+
+
 @dataclass(frozen=True)
 class PnorPnl0001SeriesBinding:
     orifice: AmesimPnor001
@@ -147,6 +171,18 @@ class PressureFlowSolver:
             for variable in ("p", "m_flow", "x", "v", "f")
         }
         self._component_equation_owners = tuple(network.components.values())
+        self._component_equation_plan = tuple(
+            ComponentEquationEvaluation(
+                component=component,
+                evaluate=component.pressure_flow_equation_values,
+                templates=component.pressure_flow_equation_residuals(),
+            )
+            for component in self._component_equation_owners
+        )
+        self._component_equation_plans_by_id = {
+            item.component.name: item
+            for item in self._component_equation_plan
+        }
         self._estimated_flow_components = tuple(
             component
             for component in self._component_equation_owners
@@ -164,6 +200,11 @@ class PressureFlowSolver:
             )
             for equation in network.connection_equation_residuals()
         )
+        self._equation_templates = tuple(
+            equation
+            for item in self._component_equation_plan
+            for equation in item.templates
+        ) + tuple(item.template for item in self._connection_equation_plan)
         self._effort_groups = {
             variable: self._build_effort_equality_groups(variable)
             for variable in ("p", "x", "v")
@@ -176,12 +217,19 @@ class PressureFlowSolver:
         )
         self._unilateral_contact_plan = self._build_unilateral_contact_plan()
         self._explicit_flow_plan = self._build_explicit_flow_plan()
+        self._explicit_flow_plans_by_variables = {
+            selected: self._filter_explicit_flow_plan(selected)
+            for selected in (frozenset(("f", "m_flow")), frozenset(("f",)))
+        }
+        self._equation_scale_plans = self._build_equation_scale_plans(
+            self._equation_templates
+        )
         self.last_diagnostics: AlgebraicSolveDiagnostics | None = None
 
     def _build_unknowns(self) -> tuple[AlgebraicUnknown, ...]:
         unknowns: list[AlgebraicUnknown] = []
         for component in self.network.components.values():
-            for definition in component.port_definitions:
+            for definition in component.active_port_definitions:
                 if definition.kind != "physical":
                     continue
                 state = component.get_port(definition.name)
@@ -277,8 +325,8 @@ class PressureFlowSolver:
                 union(first, second)
 
         component_equations = {
-            component.name: component.pressure_flow_equation_residuals()
-            for component in self.network.components.values()
+            item.component.name: item.templates
+            for item in self._component_equation_plan
         }
         for equations in component_equations.values():
             for equation in equations:
@@ -720,12 +768,10 @@ class PressureFlowSolver:
                 f"Unsupported connection equation relation: {equation.relation}."
             )
 
-        component = self.network.components[equation.owner_id]
+        evaluation = self._component_equation_plans_by_id[equation.owner_id]
+        component = evaluation.component
         equation_id = equation.id
-        equation_ids = tuple(
-            current.id
-            for current in component.pressure_flow_equation_residuals()
-        )
+        equation_ids = tuple(current.id for current in evaluation.templates)
         try:
             equation_index = equation_ids.index(equation_id)
         except ValueError as exc:
@@ -734,38 +780,49 @@ class PressureFlowSolver:
             ) from exc
 
         def read_component_equation() -> float:
-            current_equations = component.pressure_flow_equation_residuals()
-            if (
-                equation_index >= len(current_equations)
-                or current_equations[equation_index].id != equation_id
-            ):
+            current_values = evaluation.evaluate()
+            if equation_index >= len(current_values):
                 raise RuntimeError(
-                    f"Compiled algebraic equation disappeared at runtime: {equation_id}."
+                    "Compiled algebraic equation disappeared at runtime: "
+                    f"{equation_id}."
                 )
-            return float(current_equations[equation_index].value)
+            return float(current_values[equation_index])
 
         return read_component_equation
-    def _pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
-        """Evaluate live values through a precompiled connector topology."""
-        component_residuals = tuple(
-            residual
-            for component in self._component_equation_owners
-            for residual in component.pressure_flow_equation_residuals()
-        )
-        connection_residuals = tuple(
-            EquationResidual(
-                id=item.template.id,
-                owner=item.template.owner,
-                owner_id=item.template.owner_id,
-                relation=item.template.relation,
-                variables=item.template.variables,
-                value=item.evaluate(),
-                role=item.template.role,
-            )
-            for item in self._connection_equation_plan
-        )
-        return component_residuals + connection_residuals
 
+    def _pressure_flow_equation_values(self) -> tuple[float, ...]:
+        """Evaluate live equation values through the compiled topology."""
+
+        values: list[float] = []
+        for item in self._component_equation_plan:
+            current_values = item.evaluate()
+            if len(current_values) != len(item.templates):
+                raise RuntimeError(
+                    "Compiled algebraic equation count changed at runtime for "
+                    f"{item.component.name}."
+                )
+            values.extend(current_values)
+        values.extend(item.evaluate() for item in self._connection_equation_plan)
+        return tuple(values)
+
+    def _pressure_flow_equation_residuals(self) -> tuple[EquationResidual, ...]:
+        """Materialize public residual objects only when explicitly requested."""
+
+        return tuple(
+            EquationResidual(
+                id=template.id,
+                owner=template.owner,
+                owner_id=template.owner_id,
+                relation=template.relation,
+                variables=template.variables,
+                value=value,
+                role=template.role,
+            )
+            for template, value in zip(
+                self._equation_templates,
+                self._pressure_flow_equation_values(),
+            )
+        )
 
     def _explicit_flow_assignment(
         self,
@@ -779,8 +836,9 @@ class PressureFlowSolver:
                 evaluate=self._equation_value_reader(equation),
             )
 
-        component = self.network.components[equation.owner_id]
-        equations = component.pressure_flow_equation_residuals()
+        evaluation = self._component_equation_plans_by_id[equation.owner_id]
+        component = evaluation.component
+        equations = evaluation.templates
         equation_ids = tuple(current.id for current in equations)
         try:
             equation_index = equation_ids.index(equation.id)
@@ -803,8 +861,8 @@ class PressureFlowSolver:
         seeded_ids: set[str] = set()
 
         initial_assignments: list[ExplicitFlowAssignment] = []
-        for component in self.network.components.values():
-            for equation in component.pressure_flow_equation_residuals():
+        for evaluation in self._component_equation_plan:
+            for equation in evaluation.templates:
                 if equation.relation != "constitutive" or equation.role != "flow":
                     continue
                 flow_unknowns = self._flow_unknowns_for_equation(equation)
@@ -818,9 +876,9 @@ class PressureFlowSolver:
                 )
                 seeded_ids.add(unknown.id)
         if initial_assignments:
-            stages.append(ExplicitFlowStage(tuple(initial_assignments)))
+            stages.append(self._compile_explicit_flow_stage(tuple(initial_assignments)))
 
-        equations = self._pressure_flow_equation_residuals()
+        equations = self._equation_templates
         while True:
             stage_assignments: list[ExplicitFlowAssignment] = []
             stage_unknown_ids: set[str] = set()
@@ -850,7 +908,9 @@ class PressureFlowSolver:
                 )
                 stage_unknown_ids.add(unknown.id)
             if stage_assignments:
-                stages.append(ExplicitFlowStage(tuple(stage_assignments)))
+                stages.append(
+                    self._compile_explicit_flow_stage(tuple(stage_assignments))
+                )
                 seeded_ids.update(stage_unknown_ids)
                 continue
 
@@ -880,43 +940,105 @@ class PressureFlowSolver:
                 break
             if fallback_assignment is None:
                 break
-            stages.append(ExplicitFlowStage((fallback_assignment,)))
+            stages.append(self._compile_explicit_flow_stage((fallback_assignment,)))
 
         return tuple(stages)
 
+    @staticmethod
+    def _compile_explicit_flow_stage(
+        assignments: tuple[ExplicitFlowAssignment, ...],
+    ) -> ExplicitFlowStage:
+        direct_evaluations: list[tuple[int, Callable[[], float]]] = []
+        assignments_by_component: dict[
+            object,
+            list[tuple[int, ExplicitFlowAssignment]],
+        ] = {}
+        for assignment_index, assignment in enumerate(assignments):
+            if assignment.component is None:
+                if assignment.evaluate is None:
+                    raise RuntimeError(
+                        "Explicit flow assignment has no compiled evaluator: "
+                        f"{assignment.equation_id}."
+                    )
+                direct_evaluations.append((assignment_index, assignment.evaluate))
+                continue
+            assignments_by_component.setdefault(assignment.component, []).append(
+                (assignment_index, assignment)
+            )
+
+        component_evaluations: list[ExplicitFlowComponentEvaluation] = []
+        for component, component_assignments in assignments_by_component.items():
+            if any(
+                assignment.equation_index is None
+                for _assignment_index, assignment in component_assignments
+            ):
+                raise RuntimeError(
+                    "Explicit component flow assignment has no equation index."
+                )
+            component_evaluations.append(
+                ExplicitFlowComponentEvaluation(
+                    component=component,
+                    evaluate=component.pressure_flow_equation_values,
+                    assignment_indices=tuple(
+                        assignment_index
+                        for assignment_index, _assignment in component_assignments
+                    ),
+                    equation_indices=tuple(
+                        int(assignment.equation_index)
+                        for _assignment_index, assignment in component_assignments
+                    ),
+                    equation_ids=tuple(
+                        assignment.equation_id
+                        for _assignment_index, assignment in component_assignments
+                    ),
+                )
+            )
+        return ExplicitFlowStage(
+            assignments=assignments,
+            direct_evaluations=tuple(direct_evaluations),
+            component_evaluations=tuple(component_evaluations),
+        )
+
+    def _filter_explicit_flow_plan(
+        self,
+        selected: frozenset[str],
+    ) -> tuple[ExplicitFlowStage, ...]:
+        return tuple(
+            self._compile_explicit_flow_stage(
+                tuple(
+                    assignment
+                    for assignment in stage.assignments
+                    if assignment.unknown.variable in selected
+                )
+            )
+            for stage in self._explicit_flow_plan
+        )
 
     @staticmethod
     def _evaluate_explicit_flow_stage(
-        assignments: tuple[ExplicitFlowAssignment, ...],
-    ) -> dict[str, float]:
-        values: dict[str, float] = {}
-        assignments_by_component: dict[object, list[ExplicitFlowAssignment]] = {}
-        for assignment in assignments:
-            if assignment.component is None:
-                assert assignment.evaluate is not None
-                values[assignment.equation_id] = assignment.evaluate()
-                continue
-            assignments_by_component.setdefault(assignment.component, []).append(
-                assignment
-            )
+        stage: ExplicitFlowStage,
+    ) -> tuple[float, ...]:
+        assignments = stage.assignments
+        values: list[float | None] = [None] * len(assignments)
+        for assignment_index, evaluate in stage.direct_evaluations:
+            values[assignment_index] = evaluate()
 
-        for component, component_assignments in assignments_by_component.items():
-            equations = component.pressure_flow_equation_residuals()
-            for assignment in component_assignments:
-                assert assignment.equation_index is not None
-                equation_index = assignment.equation_index
-                if (
-                    equation_index >= len(equations)
-                    or equations[equation_index].id != assignment.equation_id
-                ):
+        for evaluation in stage.component_evaluations:
+            equation_values = evaluation.evaluate()
+            for assignment_index, equation_index, equation_id in zip(
+                evaluation.assignment_indices,
+                evaluation.equation_indices,
+                evaluation.equation_ids,
+            ):
+                if equation_index >= len(equation_values):
                     raise RuntimeError(
                         "Compiled algebraic equation disappeared at runtime: "
-                        f"{assignment.equation_id}."
+                        f"{equation_id}."
                     )
-                values[assignment.equation_id] = float(
-                    equations[equation_index].value
-                )
-        return values
+                values[assignment_index] = float(equation_values[equation_index])
+        if any(value is None for value in values):
+            raise RuntimeError("Explicit flow evaluation plan returned no value.")
+        return tuple(float(value) for value in values)
 
     def _solve_explicit_flow_unknowns(
         self,
@@ -930,19 +1052,18 @@ class PressureFlowSolver:
                 unknown.write(0.0)
 
         seeded_ids: set[str] = set()
-        for stage in self._explicit_flow_plan:
-            assignments = tuple(
-                assignment
-                for assignment in stage.assignments
-                if assignment.unknown.variable in selected
-            )
-            values = self._evaluate_explicit_flow_stage(assignments)
+        plan = self._explicit_flow_plans_by_variables.get(selected)
+        if plan is None:
+            plan = self._filter_explicit_flow_plan(selected)
+        for stage in plan:
+            assignments = stage.assignments
+            values = self._evaluate_explicit_flow_stage(stage)
             targets = tuple(
                 (
                     assignment,
-                    assignment.unknown.read() - values[assignment.equation_id],
+                    assignment.unknown.read() - value,
                 )
-                for assignment in assignments
+                for assignment, value in zip(assignments, values)
             )
             for assignment, target_value in targets:
                 if not isfinite(target_value):
@@ -1101,12 +1222,38 @@ class PressureFlowSolver:
             orifice.get_port(orifice_port).p = pressure
             pipe.get_port(pipe_port).p = pressure
 
+    def _equation_scale_plan(self, equation) -> EquationScalePlan:
+        return EquationScalePlan(
+            variable_names=tuple(
+                variable.rsplit(".", 1)[-1]
+                for variable in equation.variables
+            ),
+            force_unknowns=tuple(
+                unknown
+                for variable in equation.variables
+                if (unknown := self._unknowns_by_id.get(variable)) is not None
+                and unknown.variable == "f"
+            ),
+        )
+
+    def _build_equation_scale_plans(
+        self,
+        equations: tuple[EquationResidual, ...],
+    ) -> tuple[EquationScalePlan, ...]:
+        return tuple(
+            self._equation_scale_plan(equation)
+            for equation in equations
+        )
+
     def _scales(self) -> dict[str, float]:
+        pressure_values = [
+            unknown.read() for unknown in self._unknowns_by_variable["p"]
+        ]
         pressure_scale = max(
             [
-                abs(unknown.read())
-                for unknown in self._unknowns_by_variable["p"]
-                if unknown.read() > 0.0
+                abs(value)
+                for value in pressure_values
+                if value > 0.0
             ]
             + [1e5]
         )
@@ -1166,36 +1313,24 @@ class PressureFlowSolver:
         scales = self._scales()
         pressure_scale = scales["p"]
         flow_scale = scales["m_flow"]
-        unknown_scales = {
-            unknown.id: (
-                max(abs(unknown.read()), 1.0)
-                if unknown.variable == "f"
-                else scales.get(unknown.variable, max(abs(unknown.read()), 1.0))
-            )
-            for unknown in self.unknowns
-        }
-        def variable_scale(unknown: AlgebraicUnknown) -> float:
-            return unknown_scales[unknown.id]
+        seeded_values = self._pressure_flow_equation_values()
 
-        seeded_equations = self._pressure_flow_equation_residuals()
-
-        def initial_equation_scale(equation) -> float:
-            variable_names = [
-                variable.rsplit(".", 1)[-1]
-                for variable in equation.variables
-            ]
+        def initial_equation_scale(
+            equation: EquationResidual,
+            value: float,
+            scale_plan: EquationScalePlan,
+        ) -> float:
+            variable_names = scale_plan.variable_names
             if equation.role == "flow":
                 force_scales = [
-                    unknown_scales[variable]
-                    for variable in equation.variables
-                    if variable in self._unknowns_by_id
-                    and self._unknowns_by_id[variable].variable == "f"
+                    max(abs(unknown.read()), 1.0)
+                    for unknown in scale_plan.force_unknowns
                 ]
                 if force_scales:
                     # Freeze force scaling per equation.  A 1e17 N source must
                     # not hide an unrelated 40 N piston/contact imbalance in a
                     # different mechanical branch.
-                    return max(force_scales + [abs(float(equation.value)), 1.0])
+                    return max(force_scales + [abs(value), 1.0])
                 return flow_scale
             if equation.role == "effort":
                 if "x" in variable_names:
@@ -1205,20 +1340,18 @@ class PressureFlowSolver:
                 return pressure_scale
             return max([scales.get(name, 1.0) for name in variable_names] + [1.0])
 
-        equation_scales = {
-            equation.id: initial_equation_scale(equation)
-            for equation in seeded_equations
-        }
-
-        def equation_scale(equation) -> float:
-            cached = equation_scales.get(equation.id)
-            if cached is not None:
-                return cached
-            return initial_equation_scale(equation)
+        equation_scales = tuple(
+            initial_equation_scale(equation, value, scale_plan)
+            for equation, value, scale_plan in zip(
+                self._equation_templates,
+                seeded_values,
+                self._equation_scale_plans,
+            )
+        )
 
         seeded_scaled = [
-            abs(equation.value / equation_scale(equation))
-            for equation in seeded_equations
+            abs(value / scale)
+            for value, scale in zip(seeded_values, equation_scales)
         ]
         seeded_max_scaled_residual = max(seeded_scaled, default=0.0)
         seeded_unknown_values = [
@@ -1242,13 +1375,25 @@ class PressureFlowSolver:
                 flow_scale=flow_scale,
                 max_scaled_residual=seeded_max_scaled_residual,
                 max_raw_residual=max(
-                    (abs(item.value) for item in seeded_equations),
+                    (abs(value) for value in seeded_values),
                     default=0.0,
                 ),
             )
             self.last_diagnostics = diagnostics
             return diagnostics
 
+        unknown_scales = {
+            unknown.id: (
+                max(abs(unknown.read()), 1.0)
+                if unknown.variable == "f"
+                else scales[unknown.variable]
+            )
+            for unknown in self.unknowns
+        }
+
+        def variable_scale(unknown: AlgebraicUnknown) -> float:
+            return unknown_scales[unknown.id]
+
         positive_pressures = [
             unknown.read()
             for unknown in self._unknowns_by_variable["p"]
@@ -1293,11 +1438,11 @@ class PressureFlowSolver:
         def scaled_residuals(values):
             assign(values)
             self._refresh_unilateral_contacts(contact_bindings)
-            equations = self._pressure_flow_equation_residuals()
+            equation_values = self._pressure_flow_equation_values()
             return np.asarray(
                 [
-                    equation.value / equation_scale(equation)
-                    for equation in equations
+                    value / scale
+                    for value, scale in zip(equation_values, equation_scales)
                 ],
                 dtype=float,
             )
@@ -1314,12 +1459,10 @@ class PressureFlowSolver:
         )
         assign(result.x)
         self._refresh_unilateral_contacts(contact_bindings)
-        equations = self._pressure_flow_equation_residuals()
+        equation_values = self._pressure_flow_equation_values()
         scaled = [
-            abs(
-                equation.value / equation_scale(equation)
-            )
-            for equation in equations
+            abs(value / scale)
+            for value, scale in zip(equation_values, equation_scales)
         ]
         max_scaled_residual = max(scaled, default=0.0)
         residuals_converged = (
@@ -1335,7 +1478,10 @@ class PressureFlowSolver:
             pressure_scale=pressure_scale,
             flow_scale=flow_scale,
             max_scaled_residual=max_scaled_residual,
-            max_raw_residual=max((abs(item.value) for item in equations), default=0.0),
+            max_raw_residual=max(
+                (abs(value) for value in equation_values),
+                default=0.0,
+            ),
         )
         self.last_diagnostics = diagnostics
         if not success:
diff --git a/app/simulation/solvers/mechanical.py b/app/simulation/solvers/mechanical.py
index 3827185..4d683c0 100644
--- a/app/simulation/solvers/mechanical.py
+++ b/app/simulation/solvers/mechanical.py
@@ -279,7 +279,7 @@ class MechanicalStateReducer:
         mechanical_ports = {
             (component.name, definition.name)
             for component in self.network.components.values()
-            for definition in component.port_definitions
+            for definition in component.active_port_definitions
             if definition.kind == "physical" and definition.domain == "mechanical"
         }
         parents = {
@@ -339,7 +339,7 @@ class MechanicalStateReducer:
         for component in masses:
             ports = [
                 (component.name, definition.name)
-                for definition in component.port_definitions
+                for definition in component.active_port_definitions
                 if definition.kind == "physical"
                 and definition.domain == "mechanical"
             ]
@@ -354,7 +354,7 @@ class MechanicalStateReducer:
         for component in masses:
             first_port = next(
                 (component.name, definition.name)
-                for definition in component.port_definitions
+                for definition in component.active_port_definitions
                 if definition.kind == "physical"
                 and definition.domain == "mechanical"
             )
diff --git a/app/simulation/solvers/pneumatic_storage.py b/app/simulation/solvers/pneumatic_storage.py
index 33621a9..15968da 100644
--- a/app/simulation/solvers/pneumatic_storage.py
+++ b/app/simulation/solvers/pneumatic_storage.py
@@ -99,7 +99,7 @@ def _pressure_storage_endpoint_groups(
     pneumatic_endpoints = {
         Endpoint(component.name, definition.name)
         for component in network.components.values()
-        for definition in component.port_definitions
+        for definition in component.active_port_definitions
         if definition.kind == "physical" and definition.domain == "pneumatic"
     }
     parent = {endpoint: endpoint for endpoint in pneumatic_endpoints}
diff --git a/app/simulation/solvers/pneumatic_volume.py b/app/simulation/solvers/pneumatic_volume.py
index b2d0d22..8c74551 100644
--- a/app/simulation/solvers/pneumatic_volume.py
+++ b/app/simulation/solvers/pneumatic_volume.py
@@ -38,7 +38,7 @@ class PneumaticVolumeResolver:
 
     def solve(self) -> PneumaticVolumeDiagnostics:
         for component in self.network.components.values():
-            for definition in component.port_definitions:
+            for definition in component.active_port_definitions:
                 if definition.kind == "physical" and definition.domain == "pneumatic":
                     port = component.get_port(definition.name)
                     port.volume = 0.0
diff --git a/app/simulation/solvers/solver.py b/app/simulation/solvers/solver.py
index 5e23d4d..3332de7 100644
--- a/app/simulation/solvers/solver.py
+++ b/app/simulation/solvers/solver.py
@@ -44,6 +44,36 @@ class SolveIVPConfig:
     first_step: float | None = None
 
 
+@dataclass(frozen=True)
+class SolverSegmentDiagnostics:
+    """Work performed by implicit solver instances inside one event segment."""
+
+    start_time: float
+    requested_stop_time: float
+    simulated_until: float
+    nfev: int = 0
+    njev: int = 0
+    nlu: int = 0
+    accepted_step_count: int = 0
+    solver_start_count: int = 0
+    state_transition_count: int = 0
+    recoverable_retry_count: int = 0
+
+    def as_dict(self) -> dict[str, float | int]:
+        return {
+            "startTime": self.start_time,
+            "requestedStopTime": self.requested_stop_time,
+            "simulatedUntil": self.simulated_until,
+            "nfev": self.nfev,
+            "njev": self.njev,
+            "nlu": self.nlu,
+            "acceptedStepCount": self.accepted_step_count,
+            "solverStartCount": self.solver_start_count,
+            "stateTransitionCount": self.state_transition_count,
+            "recoverableRetryCount": self.recoverable_retry_count,
+        }
+
+
 @dataclass(frozen=True)
 class ODESolution:
     t: list[float]
@@ -52,6 +82,7 @@ class ODESolution:
     message: str
     status: IntegrationStatus = "completed"
     error: Exception | None = None
+    solver_segments: tuple[SolverSegmentDiagnostics, ...] = ()
 
 
 def _vector_add(a: list[float], b: list[float], scale: float = 1.0) -> list[float]:
@@ -610,6 +641,7 @@ def _integrate_scipy_stepwise(
     same_time_transition_count = 0
     integration_progressed = False
     last_reported_step: float | None = None
+    solver_segments: list[SolverSegmentDiagnostics] = []
 
     def cancellation_message() -> str:
         return (
@@ -639,9 +671,17 @@ def _integrate_scipy_stepwise(
             math.nextafter(segment_end, -math.inf) if is_breakpoint else segment_end
         )
         has_integration_interval = integration_end > last_accepted_time
+        segment_start_time = last_accepted_time
         segment_max_step = float(config.max_step)
         recoverable_retry_count = 0
         last_recoverable_error: RecoverableTrialStateError | None = None
+        segment_nfev = 0
+        segment_njev = 0
+        segment_nlu = 0
+        segment_accepted_steps = 0
+        segment_solver_starts = 0
+        segment_state_transitions = 0
+        segment_recoverable_retries = 0
 
         while has_integration_interval and last_accepted_time < integration_end:
             if cancel_check():
@@ -681,6 +721,7 @@ def _integrate_scipy_stepwise(
                 break
             except RecoverableTrialStateError as exc:
                 recoverable_retry_count += 1
+                segment_recoverable_retries += 1
                 last_recoverable_error = exc
                 next_step = 0.5 * segment_max_step
                 minimum_step = 64.0 * math.ulp(max(abs(last_accepted_time), 1.0))
@@ -697,6 +738,7 @@ def _integrate_scipy_stepwise(
                 error = exc
                 break
 
+            segment_solver_starts += 1
 
             restart_at_transition = False
             restart_after_recoverable = False
@@ -720,6 +762,7 @@ def _integrate_scipy_stepwise(
                     break
                 except RecoverableTrialStateError as exc:
                     recoverable_retry_count += 1
+                    segment_recoverable_retries += 1
                     last_recoverable_error = exc
                     attempted_step = segment_max_step
                     next_step = 0.5 * attempted_step
@@ -744,6 +787,7 @@ def _integrate_scipy_stepwise(
                 if solver.status == "failed":
                     if last_recoverable_error is not None:
                         recoverable_retry_count += 1
+                        segment_recoverable_retries += 1
                         next_step = 0.5 * segment_max_step
                         minimum_step = 64.0 * math.ulp(
                             max(abs(last_accepted_time), 1.0)
@@ -759,6 +803,7 @@ def _integrate_scipy_stepwise(
                     message = str(step_message or "Integration step failed.")
                     break
 
+                segment_accepted_steps += 1
                 step_end_time = float(solver.t)
                 step_end_state = [float(value) for value in solver.y]
                 dense_output = (
@@ -807,6 +852,7 @@ def _integrate_scipy_stepwise(
                         break
 
                 if transition is not None:
+                    segment_state_transitions += 1
                     try:
                         same_time_transition_count = (
                             _next_same_time_transition_count(
@@ -896,6 +942,9 @@ def _integrate_scipy_stepwise(
                     )
                 report_step(reported_time)
 
+            segment_nfev += int(getattr(solver, "nfev", 0))
+            segment_njev += int(getattr(solver, "njev", 0))
+            segment_nlu += int(getattr(solver, "nlu", 0))
             if status != "completed":
                 break
             if restart_after_recoverable:
@@ -903,6 +952,22 @@ def _integrate_scipy_stepwise(
             if not restart_at_transition:
                 break
 
+        solver_segments.append(
+            SolverSegmentDiagnostics(
+                start_time=float(segment_start_time),
+                requested_stop_time=float(segment_end),
+                simulated_until=float(
+                    segment_end if status == "completed" else last_accepted_time
+                ),
+                nfev=segment_nfev,
+                njev=segment_njev,
+                nlu=segment_nlu,
+                accepted_step_count=segment_accepted_steps,
+                solver_start_count=segment_solver_starts,
+                state_transition_count=segment_state_transitions,
+                recoverable_retry_count=segment_recoverable_retries,
+            )
+        )
         if status != "completed":
             break
 
@@ -949,6 +1014,7 @@ def _integrate_scipy_stepwise(
         message=message,
         status=status,
         error=error,
+        solver_segments=tuple(solver_segments),
     )
 
 
diff --git a/app/simulation/systems/generic.py b/app/simulation/systems/generic.py
index a447689..05cd2ae 100644
--- a/app/simulation/systems/generic.py
+++ b/app/simulation/systems/generic.py
@@ -160,7 +160,7 @@ def simulation_preparation_issues(
         for component in network.components.values()
         if any(
             definition.kind == "physical"
-            for definition in component.port_definitions
+            for definition in component.active_port_definitions
         )
     }
     adjacency = {name: set() for name in physical_component_names}
@@ -449,6 +449,22 @@ class GenericFluidSystem:
             self._jacobian_sparsity = self._build_jacobian_sparsity()
         return self._jacobian_sparsity
 
+    def jacobian_sparsity_diagnostics(self) -> dict[str, float | int]:
+        from scipy.optimize._numdiff import group_columns
+
+        sparsity = self.jacobian_sparsity()
+        group_count = int(group_columns(sparsity).max(initial=-1)) + 1
+        state_count = int(sparsity.shape[0])
+        return {
+            "nonzeroCount": int(sparsity.nnz),
+            "density": (
+                float(sparsity.nnz) / float(state_count * state_count)
+                if state_count
+                else 0.0
+            ),
+            "colorGroupCount": group_count,
+        }
+
     def _close_current_state(self, time: float) -> dict[str, dict[str, float]]:
         signal = self.signal_resolver.solve(time)
         self.signal_propagation_count += signal.propagated
@@ -471,7 +487,7 @@ class GenericFluidSystem:
         physical_ports = tuple(
             port
             for component in self.network.components.values()
-            for definition in component.port_definitions
+            for definition in component.active_port_definitions
             if definition.kind == "physical"
             for port in (component.get_port(definition.name),)
         )
@@ -655,6 +671,54 @@ class GenericFluidSystem:
         )
         report_progress(postprocess_progress, "postprocessing", force=True)
         times = [float(value) for value in solution.t]
+        if isinstance(solution, ODESolution):
+            solver_segment_diagnostics = [
+                segment.as_dict() for segment in solution.solver_segments
+            ]
+        else:
+            solver_segment_diagnostics = [
+                {
+                    "startTime": float(config.t_start),
+                    "requestedStopTime": float(config.t_stop),
+                    "simulatedUntil": times[-1] if times else float(config.t_start),
+                    "nfev": int(getattr(solution, "nfev", 0)),
+                    "njev": int(getattr(solution, "njev", 0)),
+                    "nlu": int(getattr(solution, "nlu", 0)),
+                    "acceptedStepCount": 0,
+                    "solverStartCount": 1,
+                    "stateTransitionCount": 0,
+                    "recoverableRetryCount": 0,
+                }
+            ]
+        solver_total_keys = (
+            "nfev",
+            "njev",
+            "nlu",
+            "acceptedStepCount",
+            "solverStartCount",
+            "stateTransitionCount",
+            "recoverableRetryCount",
+        )
+        solver_totals = {
+            key: sum(int(segment[key]) for segment in solver_segment_diagnostics)
+            for key in solver_total_keys
+        }
+        jacobian_diagnostics = (
+            self.jacobian_sparsity_diagnostics()
+            if integration_config.method in {"BDF", "Radau"}
+            else None
+        )
+        if jacobian_diagnostics is not None:
+            color_group_count = int(jacobian_diagnostics["colorGroupCount"])
+            for segment in solver_segment_diagnostics:
+                segment["finiteDifferenceRhsEstimate"] = (
+                    int(segment["njev"]) * color_group_count
+                )
+            solver_totals["finiteDifferenceRhsEstimate"] = sum(
+                int(segment["finiteDifferenceRhsEstimate"])
+                for segment in solver_segment_diagnostics
+            )
+
         series: dict[str, list[float]] = {"time": []}
         postprocessing_error: Exception | None = None
         self.mechanical_state_reducer.reset_constraint_modes()
@@ -693,6 +757,13 @@ class GenericFluidSystem:
             if key != "time" and values
         }
         diagnostics = {
+            "integration": {
+                "method": integration_config.method,
+                "jacobianSparsity": jacobian_diagnostics,
+                "segmentCount": len(solver_segment_diagnostics),
+                "segments": solver_segment_diagnostics,
+                "totals": solver_totals,
+            },
             "pressureFlow": {
                 "solveCount": self.algebraic_solve_count,
                 "maxScaledResidual": self.max_algebraic_residual,
diff --git a/app/simulation/systems/network.py b/app/simulation/systems/network.py
index 25cf2c0..bd4d825 100644
--- a/app/simulation/systems/network.py
+++ b/app/simulation/systems/network.py
@@ -229,7 +229,7 @@ class SimulationNetwork:
         return tuple(
             f"{component.name}.{definition.name}.{variable.name}"
             for component in self.components.values()
-            for definition in component.port_definitions
+            for definition in component.active_port_definitions
             if definition.kind == "physical"
             for variable in definition.variables
             if variable.role in {"effort", "flow"}
@@ -304,7 +304,7 @@ class SimulationNetwork:
                     "parameters": component.parameter_interface_dicts(),
                     "ports": [
                         definition.as_interface_dict()
-                        for definition in component.port_definitions
+                        for definition in component.active_port_definitions
                     ],
                     "resultVariables": [
                         variable.as_dict()
@@ -320,7 +320,7 @@ class SimulationNetwork:
             "unconnectedPorts": [
                 {"component": component.name, "port": definition.name}
                 for component in self.components.values()
-                for definition in component.port_definitions
+                for definition in component.active_port_definitions
                 if (component.name, definition.name) not in connected_endpoints
             ],
         }
diff --git a/tests/test_amesim_mechanical_public_components.py b/tests/test_amesim_mechanical_public_components.py
index e0d833c..bf815d6 100644
--- a/tests/test_amesim_mechanical_public_components.py
+++ b/tests/test_amesim_mechanical_public_components.py
@@ -41,6 +41,10 @@ class AmesimMechanicalPublicComponentTests(unittest.TestCase):
         residuals = converter.pressure_flow_equation_residuals()
 
         self.assertAlmostEqual(converter.output_force, 20.0)
+        self.assertEqual(
+            converter.pressure_flow_equation_values(),
+            tuple(residual.value for residual in residuals),
+        )
         self.assertAlmostEqual(residuals[0].value, 0.0)
         self.assertEqual(converter.component_result_values(), {"force": 20.0})
         self.assertEqual(converter.parameter_values, {"direction": 1.0})
@@ -54,10 +58,13 @@ class AmesimMechanicalPublicComponentTests(unittest.TestCase):
         converter.res.signal = 20.0
         converter.port_2.f = 20.0
 
-        self.assertAlmostEqual(
-            converter.pressure_flow_equation_residuals()[0].value,
-            0.0,
+        residuals = converter.pressure_flow_equation_residuals()
+
+        self.assertEqual(
+            converter.pressure_flow_equation_values(),
+            tuple(residual.value for residual in residuals),
         )
+        self.assertAlmostEqual(residuals[0].value, 0.0)
         self.assertEqual(converter.parameter_values, {"direction": -1.0})
         self.assertEqual(converter.component_result_values(), {"force": -20.0})
 
@@ -231,9 +238,17 @@ class AmesimMechanicalPublicComponentTests(unittest.TestCase):
         residuals = node.pressure_flow_equation_residuals()
 
         self.assertEqual(node.active_ports, ("port_1", "port_2", "port_3"))
+        self.assertEqual(
+            tuple(definition.name for definition in node.active_port_definitions),
+            node.active_ports,
+        )
         self.assertAlmostEqual(node.total_force, 10.0)
         self.assertAlmostEqual(node.force_balance, 0.0)
-        self.assertEqual(len(residuals), 5 + 18 * 3)
+        self.assertEqual(len(residuals), 5)
+        self.assertEqual(
+            node.pressure_flow_equation_values(),
+            tuple(residual.value for residual in residuals),
+        )
         self.assertTrue(all(abs(residual.value) <= 1.0e-12 for residual in residuals))
         self.assertEqual(node.component_result_values(), {"tforce": 10.0})
 
diff --git a/tests/test_amesim_mechanical_xml.py b/tests/test_amesim_mechanical_xml.py
index 4b98527..6bc3555 100644
--- a/tests/test_amesim_mechanical_xml.py
+++ b/tests/test_amesim_mechanical_xml.py
@@ -423,6 +423,8 @@ class AmesimMechanicalXmlTests(unittest.TestCase):
         result = run_system_xml_simulation(xml)
 
         self.assertTrue(result["success"], result["message"])
+        self.assertNotIn("node_1.port_4.x", result["series"])
+        self.assertNotIn("node_1.port_21.f", result["series"])
         self.assertAlmostEqual(result["series"]["mass_1.a"][0], 5.0)
 
     def test_elastic_contact_project_compiles_and_simulates(self) -> None:
diff --git a/tests/test_core_solver.py b/tests/test_core_solver.py
index 585a4db..bfdfa39 100644
--- a/tests/test_core_solver.py
+++ b/tests/test_core_solver.py
@@ -241,6 +241,77 @@ class IntegrateOdeTests(unittest.TestCase):
                     )
                 )
 
+    def test_segmented_solver_reports_implicit_work_by_event_segment(self) -> None:
+        import numpy as np
+        import scipy.integrate
+
+        class CountingBDF:
+            def __init__(self, _fun, t0, y0, t_bound, **_kwargs):
+                self.t = float(t0)
+                self.y = np.asarray(y0, dtype=float)
+                self.t_bound = float(t_bound)
+                self.status = "running"
+                self.nfev = 2
+                self.njev = 1
+                self.nlu = 0
+
+            def step(self):
+                self.t = self.t_bound
+                self.nfev += 3
+                self.nlu += 2
+                self.status = "finished"
+                return None
+
+            def dense_output(self):
+                state = self.y.copy()
+                return lambda _time: state.copy()
+
+        with patch.object(scipy.integrate, "BDF", CountingBDF):
+            result = integrate_ode(
+                rhs=lambda _time, _state: [0.0],
+                initial_state=[1.0],
+                config=SolveIVPConfig(
+                    t_start=0.0,
+                    t_stop=1.0,
+                    method="BDF",
+                    max_step=1.0,
+                ),
+                t_eval=[0.0, 0.4, 1.0],
+                breakpoints=[0.4],
+            )
+
+        self.assertTrue(result.success, result.message)
+        self.assertEqual(len(result.solver_segments), 2)
+        self.assertEqual(
+            [segment.as_dict() for segment in result.solver_segments],
+            [
+                {
+                    "startTime": 0.0,
+                    "requestedStopTime": 0.4,
+                    "simulatedUntil": 0.4,
+                    "nfev": 5,
+                    "njev": 1,
+                    "nlu": 2,
+                    "acceptedStepCount": 1,
+                    "solverStartCount": 1,
+                    "stateTransitionCount": 0,
+                    "recoverableRetryCount": 0,
+                },
+                {
+                    "startTime": 0.4,
+                    "requestedStopTime": 1.0,
+                    "simulatedUntil": 1.0,
+                    "nfev": 5,
+                    "njev": 1,
+                    "nlu": 2,
+                    "acceptedStepCount": 1,
+                    "solverStartCount": 1,
+                    "stateTransitionCount": 0,
+                    "recoverableRetryCount": 0,
+                },
+            ],
+        )
+
     def test_segmented_solver_can_cancel_after_crossing_a_breakpoint(self) -> None:
         callback_times: list[float] = []
         cancellation_requested = False
diff --git a/tests/test_generic_system_xml_simulation.py b/tests/test_generic_system_xml_simulation.py
index d8d746a..85a3519 100644
--- a/tests/test_generic_system_xml_simulation.py
+++ b/tests/test_generic_system_xml_simulation.py
@@ -322,6 +322,23 @@ class GenericSystemXmlSimulationTests(unittest.TestCase):
         self.assertFalse(result.success)
         self.assertEqual(result.status, "cancelled")
         self.assertGreaterEqual(result.diagnostics["sampleCount"], 2)
+        integration = result.diagnostics["integration"]
+        self.assertEqual(integration["method"], "BDF")
+        self.assertEqual(integration["segmentCount"], len(integration["segments"]))
+        self.assertEqual(
+            integration["totals"]["nfev"],
+            sum(segment["nfev"] for segment in integration["segments"]),
+        )
+        sparsity = integration["jacobianSparsity"]
+        self.assertGreater(sparsity["nonzeroCount"], 0)
+        self.assertGreater(sparsity["colorGroupCount"], 0)
+        self.assertEqual(
+            integration["totals"]["finiteDifferenceRhsEstimate"],
+            sum(
+                segment["finiteDifferenceRhsEstimate"]
+                for segment in integration["segments"]
+            ),
+        )
         self.assertGreater(result.simulated_until, 0.0)
         self.assertLess(result.simulated_until, 0.05)
         self.assertEqual(
diff --git a/tests/test_pressure_flow_solver_initialization.py b/tests/test_pressure_flow_solver_initialization.py
index adacb56..5b3fa6a 100644
--- a/tests/test_pressure_flow_solver_initialization.py
+++ b/tests/test_pressure_flow_solver_initialization.py
@@ -219,6 +219,25 @@ class PressureFlowSolverInitializationTests(unittest.TestCase):
         self.assertAlmostEqual(high.port_b.m_flow, -valve.port_2.m_flow, places=10)
         self.assertAlmostEqual(low.port_a.m_flow, -valve.port_3.m_flow, places=10)
 
+    def test_compiled_plan_does_not_rebuild_equation_metadata_during_solve(self) -> None:
+        boundary = AmesimPnpl01("compiled_closed")
+        boundary.port_1.p = 100_000.0
+        boundary.port_1.m_flow = 1.0
+        network = SimulationNetwork("compiled-closed-boundary")
+        network.add_component(boundary)
+        solver = PressureFlowSolver(network)
+
+        with patch.object(
+            boundary,
+            "pressure_flow_equation_residuals",
+            side_effect=AssertionError("equation metadata rebuilt at runtime"),
+        ):
+            diagnostics = solver.solve()
+
+        self.assertTrue(diagnostics.success)
+        self.assertEqual(diagnostics.evaluations, 0)
+        self.assertEqual(boundary.port_1.m_flow, 0.0)
+
     @staticmethod
     def _closed_boundary_solver() -> PressureFlowSolver:
         boundary = AmesimPnpl01("closed")