From 27f3bd569ff8f2b96dcdf1c36381ff9769d8458c Mon Sep 17 00:00:00 2001 From: Kyle Gorkowski Date: Wed, 4 Mar 2026 19:41:45 -0700 Subject: [PATCH 1/5] feat(condensation): cover binned latent heat cases Normalize first-order mass transport inputs for binned pressure arrays and add discrete/PDF latent heat tests for mass conservation, isothermal parity, energy tracking, and shapes. Closes #1142 ADW-ID: d6b71842 --- .../condensation/condensation_strategies.py | 26 + .../tests/condensation_strategies_test.py | 467 +++++++++++++++++- 2 files changed, 492 insertions(+), 1 deletion(-) diff --git a/particula/dynamics/condensation/condensation_strategies.py b/particula/dynamics/condensation/condensation_strategies.py index 243a9ccb4..3f550ae77 100644 --- a/particula/dynamics/condensation/condensation_strategies.py +++ b/particula/dynamics/condensation/condensation_strategies.py @@ -876,6 +876,19 @@ def mass_transfer_rate( pressure_delta, posinf=0.0, neginf=0.0, nan=0.0 ) + first_order_mass_transport = np.asarray( + first_order_mass_transport, dtype=np.float64 + ) + if pressure_delta.ndim == 2: + if first_order_mass_transport.ndim == 0: + first_order_mass_transport = np.full( + (pressure_delta.shape[0], 1), + float(first_order_mass_transport), + dtype=np.float64, + ) + elif first_order_mass_transport.ndim == 1: + first_order_mass_transport = first_order_mass_transport[:, None] + return get_mass_transfer_rate( pressure_delta=pressure_delta, first_order_mass_transport=first_order_mass_transport, @@ -1959,6 +1972,19 @@ def mass_transfer_rate( pressure_delta, posinf=0.0, neginf=0.0, nan=0.0 ) + first_order_mass_transport = np.asarray( + first_order_mass_transport, dtype=np.float64 + ) + if pressure_delta.ndim == 2: + if first_order_mass_transport.ndim == 0: + first_order_mass_transport = np.full( + (pressure_delta.shape[0], 1), + float(first_order_mass_transport), + dtype=np.float64, + ) + elif first_order_mass_transport.ndim == 1: + first_order_mass_transport = first_order_mass_transport[:, None] + if self._latent_heat_strategy is None: return get_mass_transfer_rate( pressure_delta=pressure_delta, diff --git a/particula/dynamics/condensation/tests/condensation_strategies_test.py b/particula/dynamics/condensation/tests/condensation_strategies_test.py index dff15bd2f..e1d84ba39 100644 --- a/particula/dynamics/condensation/tests/condensation_strategies_test.py +++ b/particula/dynamics/condensation/tests/condensation_strategies_test.py @@ -22,7 +22,10 @@ _unwrap_gas, _unwrap_particle, ) -from particula.dynamics.condensation.mass_transfer import get_mass_transfer +from particula.dynamics.condensation.mass_transfer import ( + get_latent_heat_energy_released, + get_mass_transfer, +) from particula.gas.gas_data import GasData, from_species from particula.gas.latent_heat_strategies import ConstantLatentHeat from particula.gas.species import GasSpecies @@ -1831,6 +1834,94 @@ def _build_single_species_fixture( return gas_species, particle, activity, surface + def _build_binned_representation( + self, + distribution_type: str, + n_bins: int, + densities: np.ndarray | None = None, + ) -> tuple[ParticleRepresentation, ParticleData, GasSpecies, GasData]: + """Build a binned particle/gas fixture for discrete or PDF bins.""" + if densities is None: + densities = np.array([1.0e3, 1.77e3, 1.2e3]) + + molar_mass = np.array([18.015e-3, 132.14e-3, 150.0e-3]) + vp_water = par.gas.VaporPressureFactory().get_strategy("water_buck") + vp_constant = par.gas.VaporPressureFactory().get_strategy( + "constant", {"vapor_pressure": 1e-24, "vapor_pressure_units": "Pa"} + ) + water_sat = vp_water.saturation_concentration( + molar_mass[0], self.temperature + ) + water_conc = water_sat * 1.02 + gas_species = ( + par.gas.GasSpeciesBuilder() + .set_name(np.array(["H2O", "NH4HSO4", "ORG"])) + .set_molar_mass(molar_mass, "kg/mol") + .set_vapor_pressure_strategy([vp_water, vp_constant, vp_constant]) + .set_concentration(np.array([water_conc, 1e-30, 1e-30]), "kg/m^3") + .set_partitioning(True) + .build() + ) + + activity = ( + par.particles.ActivityKappaParameterBuilder() + .set_density(densities, "kg/m^3") + .set_kappa(np.array([0.0, 0.61, 0.2])) + .set_molar_mass(molar_mass, "kg/mol") + .set_water_index(0) + .build() + ) + surface = ( + par.particles.SurfaceStrategyVolumeBuilder() + .set_density(densities, "kg/m^3") + .set_surface_tension(np.array([0.072, 0.092, 0.08]), "N/m") + .build() + ) + + radius_bins = np.logspace(-9, -6, n_bins) + mode = np.array([100e-9]) + gsd = np.array([1.4]) + number_concentration = np.array([1e4]) * 1e6 + base_particle = ( + par.particles.PresetParticleRadiusBuilder() + .set_distribution_type(distribution_type) + .set_radius_bins(radius_bins, "m") + .set_mode(mode, "m") + .set_geometric_standard_deviation(gsd) + .set_number_concentration(number_concentration, "1/m^3") + .set_charge(np.zeros_like(radius_bins)) + .build() + ) + concentration = base_particle.get_concentration() + volume_bins = 4 / 3 * np.pi * radius_bins**3 + fractions = np.array([0.4, 0.3, 0.3]) + masses = volume_bins[:, np.newaxis] * densities * fractions + particle = ParticleRepresentation( + strategy=par.particles.SpeciatedMassMovingBin(), + activity=activity, + surface=surface, + distribution=masses, + density=densities, + concentration=concentration, + charge=np.zeros_like(concentration), + volume=1.0, + ) + return ( + particle, + from_representation(particle), + gas_species, + from_species(gas_species), + ) + + def _total_mass_binned( + self, particle_data: ParticleData, gas_data: GasData + ) -> float: + """Return total mass using concentration-weighted binned masses.""" + masses = particle_data.masses[0] + concentration = particle_data.concentration[0] + particle_mass = np.sum(masses * concentration[:, np.newaxis], axis=0) + return float(particle_mass.sum() + gas_data.concentration[0].sum()) + def test_instantiation_with_strategy(self): """Latent heat strategy is used when provided.""" strategy = ConstantLatentHeat(latent_heat_ref=2.26e6) @@ -2352,6 +2443,380 @@ def test_step_mass_conservation_data_only(self): relative_error = abs(final_mass - initial_mass) / initial_mass self.assertLess(relative_error, 1e-12) + def test_step_discrete_10_bins_mass_conservation(self): + """Discrete binned step conserves mass for 10 bins.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pmf", 10) + ) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6), + ) + + initial_mass = self._total_mass_binned(particle_data, gas_data) + particle_new, gas_new = cond.step( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + final_mass = self._total_mass_binned(particle_new, gas_new) + np.testing.assert_allclose(final_mass, initial_mass, rtol=1e-12) + + def test_step_discrete_50_bins_mass_conservation(self): + """Discrete binned step conserves mass for 50 bins.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pmf", 50) + ) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6), + ) + + initial_mass = self._total_mass_binned(particle_data, gas_data) + particle_new, gas_new = cond.step( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + final_mass = self._total_mass_binned(particle_new, gas_new) + np.testing.assert_allclose(final_mass, initial_mass, rtol=1e-12) + + def test_step_discrete_100_bins_mass_conservation(self): + """Discrete binned step conserves mass for 100 bins.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pmf", 100) + ) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6), + ) + + initial_mass = self._total_mass_binned(particle_data, gas_data) + particle_new, gas_new = cond.step( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + final_mass = self._total_mass_binned(particle_new, gas_new) + np.testing.assert_allclose(final_mass, initial_mass, rtol=1e-12) + + def test_step_continuous_distribution_mass_conservation(self): + """Continuous binned step conserves total mass.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pdf", 50) + ) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6), + ) + + initial_mass = self._total_mass_binned(particle_data, gas_data) + particle_new, gas_new = cond.step( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + final_mass = self._total_mass_binned(particle_new, gas_new) + np.testing.assert_allclose(final_mass, initial_mass, rtol=1e-12) + + def test_step_discrete_isothermal_parity(self): + """Discrete binned step matches isothermal when latent heat is zero.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pmf", 50) + ) + latent_cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + latent_heat=0.0, + ) + iso_cond = CondensationIsothermal( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + ) + + latent_particle, latent_gas = latent_cond.step( + particle=copy.deepcopy(particle_data), + gas_species=copy.deepcopy(gas_data), + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + iso_particle, iso_gas = iso_cond.step( + particle=copy.deepcopy(particle_data), + gas_species=copy.deepcopy(gas_data), + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + + np.testing.assert_allclose( + latent_particle.masses[0], iso_particle.masses[0], rtol=1e-15 + ) + np.testing.assert_allclose( + latent_gas.concentration[0], iso_gas.concentration[0], rtol=1e-15 + ) + + def test_step_continuous_isothermal_parity(self): + """Continuous binned step matches isothermal when latent heat is zero.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pdf", 50) + ) + latent_cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + latent_heat=0.0, + ) + iso_cond = CondensationIsothermal( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + ) + + latent_particle, latent_gas = latent_cond.step( + particle=copy.deepcopy(particle_data), + gas_species=copy.deepcopy(gas_data), + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + iso_particle, iso_gas = iso_cond.step( + particle=copy.deepcopy(particle_data), + gas_species=copy.deepcopy(gas_data), + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + + np.testing.assert_allclose( + latent_particle.masses[0], iso_particle.masses[0], rtol=1e-15 + ) + np.testing.assert_allclose( + latent_gas.concentration[0], iso_gas.concentration[0], rtol=1e-15 + ) + + def test_step_discrete_energy_tracking(self): + """Latent heat energy tracking matches mass transfer for binned data.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pmf", 50) + ) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + latent_heat_strategy=ConstantLatentHeat(latent_heat_ref=2.26e6), + ) + + mass_rate = cond.mass_transfer_rate( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + ) + volume = particle_data.volume[0] + norm_conc = particle_data.concentration[0] / volume + mass_transfer = get_mass_transfer( + mass_rate=mass_rate, + time_step=1.0, + gas_mass=gas_data.concentration[0], + particle_mass=particle_data.masses[0], + particle_concentration=norm_conc, + ) + latent_heat = cond._latent_heat_strategy.latent_heat(self.temperature) + expected_energy = np.sum( + get_latent_heat_energy_released(mass_transfer, latent_heat) + ) + + cond.step( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + + np.testing.assert_allclose( + cond.last_latent_heat_energy, expected_energy, rtol=1e-14 + ) + + def test_step_discrete_2d_mass_shapes(self): + """Binned discrete step preserves 2D mass shapes.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pmf", 10) + ) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + ) + + particle_new, _ = cond.step( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + self.assertEqual( + particle_new.masses[0].shape, + ( + particle_data.masses[0].shape[0], + particle_data.masses[0].shape[1], + ), + ) + + def test_rate_discrete_shape(self): + """rate() returns binned shape for discrete and continuous inputs.""" + for distribution_type in ("pmf", "pdf"): + with self.subTest(distribution_type=distribution_type): + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation(distribution_type, 10) + ) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=( + gas_species.pure_vapor_pressure_strategy + ), + ) + + rate = cond.rate( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + ) + self.assertEqual(rate.shape, particle_data.masses[0].shape) + + def test_mass_transfer_rate_discrete_shape(self): + """mass_transfer_rate returns binned shape for discrete/continuous.""" + for distribution_type in ("pmf", "pdf"): + with self.subTest(distribution_type=distribution_type): + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation(distribution_type, 10) + ) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=( + gas_species.pure_vapor_pressure_strategy + ), + ) + + mass_rate = cond.mass_transfer_rate( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + ) + self.assertEqual(mass_rate.shape, particle_data.masses[0].shape) + + def test_mass_transfer_rate_rejects_nonfinite_latent_heat(self): + """mass_transfer_rate rejects non-finite latent heat values.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pmf", 10) + ) + + class NonFiniteLatentHeat: + """Latent heat strategy returning non-finite values.""" + + def latent_heat(self, temperature: float) -> float: + return float("nan") + + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + latent_heat_strategy=NonFiniteLatentHeat(), + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + ) + + with self.assertRaisesRegex(ValueError, "latent_heat"): + cond.mass_transfer_rate( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + ) + + def test_step_rejects_nonpositive_concentration(self): + """Binned step rejects nonpositive concentrations.""" + particle, particle_data, gas_species, gas_data = ( + self._build_binned_representation("pmf", 10) + ) + particle_data.concentration[0][0] = -1.0 + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + activity_strategy=particle.activity, + surface_strategy=particle.surface, + vapor_pressure_strategy=gas_species.pure_vapor_pressure_strategy, + ) + + with self.assertRaisesRegex(ValueError, "concentration must be finite"): + cond.step( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + time_step=1.0, + ) + def test_step_rejects_nonpositive_volume(self): """Step rejects nonpositive volumes for norm_conc.""" particle_data, gas_data = self._make_data_inputs() From a7d53f9d548863c086f1cc1b28387d244d34bba4 Mon Sep 17 00:00:00 2001 From: Kyle Gorkowski Date: Wed, 4 Mar 2026 19:44:17 -0700 Subject: [PATCH 2/5] fix(validate): address validation gaps for #1142 Successfully fixed: - Clarified continuous-bin isothermal parity test docstring to state latent-heat-zero parity expectation Still failing: - None Closes #1142 ADW-ID: d6b71842 --- .../dynamics/condensation/tests/condensation_strategies_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/particula/dynamics/condensation/tests/condensation_strategies_test.py b/particula/dynamics/condensation/tests/condensation_strategies_test.py index e1d84ba39..5dddc400f 100644 --- a/particula/dynamics/condensation/tests/condensation_strategies_test.py +++ b/particula/dynamics/condensation/tests/condensation_strategies_test.py @@ -2593,7 +2593,7 @@ def test_step_discrete_isothermal_parity(self): ) def test_step_continuous_isothermal_parity(self): - """Continuous binned step matches isothermal when latent heat is zero.""" + """Continuous bins match isothermal when latent heat is zero.""" particle, particle_data, gas_species, gas_data = ( self._build_binned_representation("pdf", 50) ) From e9afac770e1ff6796d5892184f11953eca879121 Mon Sep 17 00:00:00 2001 From: Kyle Gorkowski Date: Wed, 4 Mar 2026 20:11:41 -0700 Subject: [PATCH 3/5] docs(dev-plans): link E5-F3-P4 to issue 1142 Record the P4 issue assignment and update timestamps in the epic and feature plans to reflect the in-progress status. Closes #1142 ADW-ID: d6b71842 --- adw-docs/dev-plans/README.md | 2 +- adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md | 4 +++- .../features/E5-F3-condensation-latent-heat-strategy.md | 5 +++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/adw-docs/dev-plans/README.md b/adw-docs/dev-plans/README.md index 969118530..9673b7856 100644 --- a/adw-docs/dev-plans/README.md +++ b/adw-docs/dev-plans/README.md @@ -77,7 +77,7 @@ and rollout. - Scope: Thermal resistance factor and non-isothermal mass transfer rate pure functions with energy tracking. - [E5-F3: CondensationLatentHeat Strategy Class][e5-f3] — Status: In Progress - (P1, #1139) + (P1, #1139; P4, #1142) - Scope: New condensation strategy with latent heat correction and energy diagnostics. - [E5-F4: Builder, Factory, and Exports][e5-f4] — Status: Planning diff --git a/adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md b/adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md index 71dc38072..2b408be08 100644 --- a/adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md +++ b/adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md @@ -5,7 +5,7 @@ **Owners**: @Gorkowski **Start Date**: 2026-03-02 **Target Date**: TBD -**Last Updated**: 2026-03-04 +**Last Updated**: 2026-03-05 **Size**: Medium (7 features, ~22 phases) ## Vision @@ -306,6 +306,7 @@ L -> 0 as T -> T_c. Used in engineering thermodynamics and EOS-based models. single-species and multi-species particle-resolved cases - [ ] **E5-F3-P4**: Add discrete and continuous distribution support with tests + - Issue: #1142 | Status: In Progress - Extend `step()` to handle discrete (binned) and continuous (PDF) distribution types - These use the same `get_mass_transfer()` routing (single vs multiple @@ -623,3 +624,4 @@ class CondensationLatentHeat(CondensationStrategy): | 2026-03-02 | Initial epic creation | ADW | | 2026-03-02 | Split E5-F1-P3 into P3 (builders) + P4 (factory+exports); split E5-F3-P3 into P3 (particle-resolved step) + P4 (discrete+continuous) + P5 (data-only parity); added missing details: function signatures, file references, thermal conductivity source, vapor_pressure_surface parameter, test tolerances, literature targets | ADW | | 2026-03-04 | Noted E5-F3-P1 issue #1139 and logging expectations | ADW | +| 2026-03-05 | Linked E5-F3-P4 to issue #1142 | ADW | diff --git a/adw-docs/dev-plans/features/E5-F3-condensation-latent-heat-strategy.md b/adw-docs/dev-plans/features/E5-F3-condensation-latent-heat-strategy.md index 97bf96004..b0f200f5c 100644 --- a/adw-docs/dev-plans/features/E5-F3-condensation-latent-heat-strategy.md +++ b/adw-docs/dev-plans/features/E5-F3-condensation-latent-heat-strategy.md @@ -6,7 +6,7 @@ **Owners**: @Gorkowski **Start Date**: 2026-03-04 **Target Date**: TBD -**Last Updated**: 2026-03-04 +**Last Updated**: 2026-03-05 **Size**: Large (5 phases) ## Summary @@ -162,7 +162,7 @@ additions: * L) to machine precision, sign conventions, single and multi-species - [ ] **E5-F3-P4**: Add discrete and continuous distribution support with tests - - Issue: TBD | Size: S (~60 LOC) | Status: Not Started + - Issue: #1142 | Size: S (~60 LOC) | Status: In Progress - Extend `step()` to handle discrete (binned) and continuous (PDF) distribution types - Uses same `get_mass_transfer()` routing as isothermal -- the only @@ -241,3 +241,4 @@ additions: |------|--------|--------| | 2026-03-02 | Initial feature document created from E5 epic | ADW | | 2026-03-04 | Marked P1 in progress for issue #1139 | ADW | +| 2026-03-05 | Linked P4 to issue #1142 | ADW | From 08f1ca52d59d77cf568ee74efe15336bb1291f84 Mon Sep 17 00:00:00 2001 From: Kyle Gorkowski Date: Thu, 5 Mar 2026 06:05:07 -0700 Subject: [PATCH 4/5] fix(validate): address validation gaps for #1149 Successfully fixed: - Added time_step validation and tests for isothermal and latent heat steps. - Guarded non-finite Kelvin terms in vapor pressure surface handling. - Clamped small negative concentrations in norm_conc with tolerance. - Deduplicated transport normalization and added coverage for 2D shapes. - Improved staggered batch updates to reduce per-particle overhead. - Adjusted parity test tolerances and renamed negative concentration test. - Updated condensation planning docs to reflect current issue status. Still failing (if any): - No remaining failures. Notes: - Wrapped long lines in condensation_strategies.py to satisfy ruff E501. Closes #1149 ADW-ID: 86d29772 --- adw-docs/dev-plans/README.md | 2 +- .../epics/E5-non-isothermal-condensation.md | 4 +- ...E5-F3-condensation-latent-heat-strategy.md | 5 +- .../condensation/condensation_strategies.py | 305 +++++++++++++++--- .../tests/condensation_strategies_test.py | 93 +++++- 5 files changed, 346 insertions(+), 63 deletions(-) diff --git a/adw-docs/dev-plans/README.md b/adw-docs/dev-plans/README.md index 9673b7856..969118530 100644 --- a/adw-docs/dev-plans/README.md +++ b/adw-docs/dev-plans/README.md @@ -77,7 +77,7 @@ and rollout. - Scope: Thermal resistance factor and non-isothermal mass transfer rate pure functions with energy tracking. - [E5-F3: CondensationLatentHeat Strategy Class][e5-f3] — Status: In Progress - (P1, #1139; P4, #1142) + (P1, #1139) - Scope: New condensation strategy with latent heat correction and energy diagnostics. - [E5-F4: Builder, Factory, and Exports][e5-f4] — Status: Planning diff --git a/adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md b/adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md index 2b408be08..71dc38072 100644 --- a/adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md +++ b/adw-docs/dev-plans/epics/E5-non-isothermal-condensation.md @@ -5,7 +5,7 @@ **Owners**: @Gorkowski **Start Date**: 2026-03-02 **Target Date**: TBD -**Last Updated**: 2026-03-05 +**Last Updated**: 2026-03-04 **Size**: Medium (7 features, ~22 phases) ## Vision @@ -306,7 +306,6 @@ L -> 0 as T -> T_c. Used in engineering thermodynamics and EOS-based models. single-species and multi-species particle-resolved cases - [ ] **E5-F3-P4**: Add discrete and continuous distribution support with tests - - Issue: #1142 | Status: In Progress - Extend `step()` to handle discrete (binned) and continuous (PDF) distribution types - These use the same `get_mass_transfer()` routing (single vs multiple @@ -624,4 +623,3 @@ class CondensationLatentHeat(CondensationStrategy): | 2026-03-02 | Initial epic creation | ADW | | 2026-03-02 | Split E5-F1-P3 into P3 (builders) + P4 (factory+exports); split E5-F3-P3 into P3 (particle-resolved step) + P4 (discrete+continuous) + P5 (data-only parity); added missing details: function signatures, file references, thermal conductivity source, vapor_pressure_surface parameter, test tolerances, literature targets | ADW | | 2026-03-04 | Noted E5-F3-P1 issue #1139 and logging expectations | ADW | -| 2026-03-05 | Linked E5-F3-P4 to issue #1142 | ADW | diff --git a/adw-docs/dev-plans/features/E5-F3-condensation-latent-heat-strategy.md b/adw-docs/dev-plans/features/E5-F3-condensation-latent-heat-strategy.md index b0f200f5c..97bf96004 100644 --- a/adw-docs/dev-plans/features/E5-F3-condensation-latent-heat-strategy.md +++ b/adw-docs/dev-plans/features/E5-F3-condensation-latent-heat-strategy.md @@ -6,7 +6,7 @@ **Owners**: @Gorkowski **Start Date**: 2026-03-04 **Target Date**: TBD -**Last Updated**: 2026-03-05 +**Last Updated**: 2026-03-04 **Size**: Large (5 phases) ## Summary @@ -162,7 +162,7 @@ additions: * L) to machine precision, sign conventions, single and multi-species - [ ] **E5-F3-P4**: Add discrete and continuous distribution support with tests - - Issue: #1142 | Size: S (~60 LOC) | Status: In Progress + - Issue: TBD | Size: S (~60 LOC) | Status: Not Started - Extend `step()` to handle discrete (binned) and continuous (PDF) distribution types - Uses same `get_mass_transfer()` routing as isothermal -- the only @@ -241,4 +241,3 @@ additions: |------|--------|--------| | 2026-03-02 | Initial feature document created from E5 epic | ADW | | 2026-03-04 | Marked P1 in progress for issue #1139 | ADW | -| 2026-03-05 | Linked P4 to issue #1142 | ADW | diff --git a/particula/dynamics/condensation/condensation_strategies.py b/particula/dynamics/condensation/condensation_strategies.py index 3f550ae77..62eb08b22 100644 --- a/particula/dynamics/condensation/condensation_strategies.py +++ b/particula/dynamics/condensation/condensation_strategies.py @@ -124,6 +124,12 @@ def _require_single_box(n_boxes: int, label: str) -> None: ) +def _validate_time_step(time_step: float) -> None: + """Validate that the timestep is finite and nonnegative.""" + if not np.isfinite(time_step) or time_step < 0.0: + raise ValueError("time_step must be finite and nonnegative.") + + def _pure_vapor_pressure_from_strategy( vapor_pressure_strategy: VaporPressureStrategy | Sequence[VaporPressureStrategy], @@ -174,6 +180,24 @@ def _partial_pressure_from_strategy( ) +def _normalize_first_order_mass_transport( + first_order_mass_transport: float | NDArray[np.float64], + pressure_delta: NDArray[np.float64], +) -> NDArray[np.float64]: + """Normalize transport shape to match pressure-delta dimensionality.""" + mass_transport = np.asarray(first_order_mass_transport, dtype=np.float64) + if pressure_delta.ndim == 2: + if mass_transport.ndim == 0: + return np.full( + (pressure_delta.shape[0], 1), + float(mass_transport), + dtype=np.float64, + ) + if mass_transport.ndim == 1: + return mass_transport[:, None] + return mass_transport + + # mass transfer abstract class class CondensationStrategy(ABC): """Abstract base class for condensation strategies. @@ -579,9 +603,10 @@ def _get_particle_partial_pressure_and_kelvin( mass_concentration=mass_concentration_in_particle, temperature=temperature, ) - return partial_pressure_particle, np.asarray( - kelvin_term, dtype=np.float64 - ) + kelvin_term_array = np.asarray(kelvin_term, dtype=np.float64) + if not np.all(np.isfinite(kelvin_term_array)): + raise ValueError("kelvin_term must be finite for vapor pressure.") + return partial_pressure_particle, kelvin_term_array def _apply_skip_partitioning( self, array: NDArray[np.float64] @@ -707,6 +732,7 @@ def step( temperature: float, pressure: float, time_step: float, + dynamic_viscosity: Optional[float] = None, ) -> Tuple[ParticleRepresentation | ParticleData, GasSpecies | GasData]: """Perform one timestep of isothermal condensation on the particle. @@ -876,18 +902,10 @@ def mass_transfer_rate( pressure_delta, posinf=0.0, neginf=0.0, nan=0.0 ) - first_order_mass_transport = np.asarray( - first_order_mass_transport, dtype=np.float64 - ) - if pressure_delta.ndim == 2: - if first_order_mass_transport.ndim == 0: - first_order_mass_transport = np.full( - (pressure_delta.shape[0], 1), - float(first_order_mass_transport), - dtype=np.float64, - ) - elif first_order_mass_transport.ndim == 1: - first_order_mass_transport = first_order_mass_transport[:, None] + first_order_mass_transport = _normalize_first_order_mass_transport( + first_order_mass_transport, + pressure_delta, + ) return get_mass_transfer_rate( pressure_delta=pressure_delta, @@ -971,6 +989,7 @@ def step( temperature: float, pressure: float, time_step: float, + dynamic_viscosity: Optional[float] = None, ) -> Tuple[ParticleRepresentation, GasSpecies]: ... @overload @@ -981,6 +1000,7 @@ def step( temperature: float, pressure: float, time_step: float, + dynamic_viscosity: Optional[float] = None, ) -> Tuple[ParticleData, GasData]: ... def step( @@ -990,6 +1010,7 @@ def step( temperature: float, pressure: float, time_step: float, + dynamic_viscosity: Optional[float] = None, ) -> Tuple[ParticleRepresentation | ParticleData, GasSpecies | GasData]: """Advance the simulation one timestep using isothermal condensation. @@ -1014,6 +1035,8 @@ def step( ) ``` """ + _validate_time_step(time_step) + # Calculate the mass transfer rate particle_data, particle_is_legacy = _unwrap_particle(particle) gas_data, gas_is_legacy = _unwrap_gas(gas_species) @@ -1026,6 +1049,7 @@ def step( gas_species=gas_species, temperature=temperature, pressure=pressure, + dynamic_viscosity=dynamic_viscosity, ) # Type guard: ensure mass_rate is an array @@ -1452,6 +1476,7 @@ def _calculate_single_particle_transfer( dt_local: float, radii: Optional[NDArray[np.float64]] = None, first_order_mass_transport: Optional[NDArray[np.float64]] = None, + dynamic_viscosity: Optional[float] = None, ) -> NDArray[np.float64]: """Calculate single-particle mass change without mutating inputs. @@ -1515,6 +1540,7 @@ def _calculate_single_particle_transfer( particle_radius=radius, temperature=temperature, pressure=pressure, + dynamic_viscosity=dynamic_viscosity, ) elif np.ndim(first_order_mass_transport) == 0: mass_transport = first_order_mass_transport @@ -1581,6 +1607,108 @@ def _calculate_single_particle_transfer( # Squeeze back to 1D (n_species,) for single particle return mass_to_change.squeeze() + def _calculate_batch_particle_transfer( + self, + particle_data: ParticleData, + gas_data: GasData, + activity_strategy: ActivityStrategy, + surface_strategy: SurfaceStrategy, + vapor_pressure_strategy: VaporPressureStrategy + | Sequence[VaporPressureStrategy], + particle_indices: NDArray[np.int64], + gas_concentration: NDArray[np.float64], + temperature: float, + pressure: float, + dt_local: float, + radii: NDArray[np.float64], + first_order_mass_transport: Optional[NDArray[np.float64]] = None, + dynamic_viscosity: Optional[float] = None, + ) -> NDArray[np.float64]: + """Calculate batch mass changes without mutating inputs.""" + particle_mass = particle_data.masses[0][particle_indices] + particle_concentration = np.asarray( + particle_data.concentration[0][particle_indices] + ) + gas_mass = np.asarray(gas_concentration, dtype=np.float64) + radius_batch = np.maximum( + radii[particle_indices], MIN_PARTICLE_RADIUS_M + ) + + if first_order_mass_transport is None: + mass_transport = self.first_order_mass_transport( + particle_radius=radius_batch, + temperature=temperature, + pressure=pressure, + dynamic_viscosity=dynamic_viscosity, + ) + elif np.ndim(first_order_mass_transport) == 0: + mass_transport = first_order_mass_transport + else: + mass_transport = first_order_mass_transport[particle_indices] + + pure_vapor_pressure = _pure_vapor_pressure_from_strategy( + vapor_pressure_strategy, temperature + ) + partial_pressure_particle = np.asarray( + activity_strategy.partial_pressure( + pure_vapor_pressure=pure_vapor_pressure, + mass_concentration=particle_mass, + ) + ) + if ( + partial_pressure_particle.ndim == 2 + and partial_pressure_particle.shape[1] == 1 + ): + partial_pressure_particle = partial_pressure_particle[:, 0] + + molar_mass = np.asarray(gas_data.molar_mass) + partial_pressure_gas = _partial_pressure_from_strategy( + vapor_pressure_strategy, + concentration=gas_mass, + molar_mass=molar_mass, + temperature=temperature, + ) + kelvin_term = surface_strategy.kelvin_term( + radius=radius_batch, + molar_mass=self.molar_mass, + mass_concentration=particle_mass, + temperature=temperature, + ) + kelvin_term_array = np.asarray(kelvin_term, dtype=np.float64) + if not np.all(np.isfinite(kelvin_term_array)): + raise ValueError("kelvin_term must be finite for vapor pressure.") + + pressure_delta = get_partial_pressure_delta( + partial_pressure_gas=partial_pressure_gas, + partial_pressure_particle=partial_pressure_particle, + kelvin_term=kelvin_term_array, + ) + pressure_delta = np.nan_to_num( + pressure_delta, posinf=0.0, neginf=0.0, nan=0.0 + ) + + mass_transport = _normalize_first_order_mass_transport( + mass_transport, + np.asarray(pressure_delta), + ) + mass_rate = get_mass_transfer_rate( + pressure_delta=pressure_delta, + first_order_mass_transport=mass_transport, + temperature=temperature, + molar_mass=self.molar_mass, + ) + mass_rate_array = np.asarray(mass_rate, dtype=np.float64) + if particle_mass.ndim == 2 and mass_rate_array.ndim == 1: + mass_rate_array = mass_rate_array[:, None] + + return get_mass_transfer( + mass_rate=mass_rate_array, + time_step=dt_local, + gas_mass=np.atleast_1d(gas_mass), + particle_mass=particle_mass, + particle_concentration=particle_concentration, + ) + # pylint: disable=too-many-positional-arguments, too-many-arguments @overload # type: ignore[override] def step( @@ -1590,6 +1718,7 @@ def step( temperature: float, pressure: float, time_step: float, + dynamic_viscosity: Optional[float] = None, ) -> Tuple[ParticleRepresentation, GasSpecies]: ... @overload @@ -1600,6 +1729,7 @@ def step( temperature: float, pressure: float, time_step: float, + dynamic_viscosity: Optional[float] = None, ) -> Tuple[ParticleData, GasData]: ... def step( # noqa: C901 @@ -1609,6 +1739,7 @@ def step( # noqa: C901 temperature: float, pressure: float, time_step: float, + dynamic_viscosity: Optional[float] = None, ) -> Tuple[ParticleRepresentation | ParticleData, GasSpecies | GasData]: """Perform two-pass staggered condensation update. @@ -1637,6 +1768,8 @@ def step( # noqa: C901 skip-partitioning, and gas updates run only when ``update_gases`` is True. """ + _validate_time_step(time_step) + particle_data, particle_is_legacy = _unwrap_particle(particle) gas_data, gas_is_legacy = _unwrap_gas(gas_species) _require_matching_types(particle_is_legacy, gas_is_legacy) @@ -1661,6 +1794,13 @@ def step( # noqa: C901 particle_radius=radii, temperature=temperature, pressure=pressure, + dynamic_viscosity=dynamic_viscosity, + ) + ) + + activity_strategy, surface_strategy, vapor_pressure_strategy = ( + self._resolve_strategies( + particle, gas_species, particle_is_legacy, gas_is_legacy ) ) @@ -1678,23 +1818,54 @@ def step( # noqa: C901 for batch in batches: batch_dm_total.fill(0.0) # Gauss-Seidel: update gas after each batch in this pass. - for particle_index in batch: - dt_local = theta_dt_first[particle_index] - if dt_local <= 0.0: - continue - mass_change = self._calculate_single_particle_transfer( - particle=particle, - particle_index=int(particle_index), - gas_species=gas_species, + batch_indices = np.asarray(batch, dtype=int) + dt_local = theta_dt_first[batch_indices] + active_mask = dt_local > 0.0 + if not np.any(active_mask): + continue + batch_indices = batch_indices[active_mask] + dt_local = dt_local[active_mask] + use_vectorized = np.allclose(dt_local, dt_local[0]) and getattr( + self._calculate_single_particle_transfer, "__func__", None + ) is ( + CondensationIsothermalStaggered._calculate_single_particle_transfer + ) + if use_vectorized: + mass_change = self._calculate_batch_particle_transfer( + particle_data=particle_data, + gas_data=gas_data, + activity_strategy=activity_strategy, + surface_strategy=surface_strategy, + vapor_pressure_strategy=vapor_pressure_strategy, + particle_indices=batch_indices, gas_concentration=working_gas_concentration, temperature=temperature, pressure=pressure, - dt_local=float(dt_local), + dt_local=float(dt_local[0]), radii=radii, first_order_mass_transport=first_order_mass_transport, + dynamic_viscosity=dynamic_viscosity, ) - mass_changes[particle_index] += mass_change - batch_dm_total += mass_change + mass_changes[batch_indices] += mass_change + batch_dm_total += mass_change.sum(axis=0) + else: + for particle_index, dt_value in zip( + batch_indices, dt_local, strict=False + ): + mass_change = self._calculate_single_particle_transfer( + particle=particle, + particle_index=int(particle_index), + gas_species=gas_species, + gas_concentration=working_gas_concentration, + temperature=temperature, + pressure=pressure, + dt_local=float(dt_value), + radii=radii, + first_order_mass_transport=first_order_mass_transport, + dynamic_viscosity=dynamic_viscosity, + ) + mass_changes[particle_index] += mass_change + batch_dm_total += mass_change working_gas_concentration = np.maximum( working_gas_concentration - batch_dm_total, 0.0 ) @@ -1703,23 +1874,54 @@ def step( # noqa: C901 batch_dm_total.fill(0.0) # Second pass updates gas after each batch for parity with # num_batches=1 behavior. - for particle_index in batch: - dt_local = theta_dt_second[particle_index] - if dt_local <= 0.0: - continue - mass_change = self._calculate_single_particle_transfer( - particle=particle, - particle_index=int(particle_index), - gas_species=gas_species, + batch_indices = np.asarray(batch, dtype=int) + dt_local = theta_dt_second[batch_indices] + active_mask = dt_local > 0.0 + if not np.any(active_mask): + continue + batch_indices = batch_indices[active_mask] + dt_local = dt_local[active_mask] + use_vectorized = np.allclose(dt_local, dt_local[0]) and getattr( + self._calculate_single_particle_transfer, "__func__", None + ) is ( + CondensationIsothermalStaggered._calculate_single_particle_transfer + ) + if use_vectorized: + mass_change = self._calculate_batch_particle_transfer( + particle_data=particle_data, + gas_data=gas_data, + activity_strategy=activity_strategy, + surface_strategy=surface_strategy, + vapor_pressure_strategy=vapor_pressure_strategy, + particle_indices=batch_indices, gas_concentration=working_gas_concentration, temperature=temperature, pressure=pressure, - dt_local=float(dt_local), + dt_local=float(dt_local[0]), radii=radii, first_order_mass_transport=first_order_mass_transport, + dynamic_viscosity=dynamic_viscosity, ) - mass_changes[particle_index] += mass_change - batch_dm_total += mass_change + mass_changes[batch_indices] += mass_change + batch_dm_total += mass_change.sum(axis=0) + else: + for particle_index, dt_value in zip( + batch_indices, dt_local, strict=False + ): + mass_change = self._calculate_single_particle_transfer( + particle=particle, + particle_index=int(particle_index), + gas_species=gas_species, + gas_concentration=working_gas_concentration, + temperature=temperature, + pressure=pressure, + dt_local=float(dt_value), + radii=radii, + first_order_mass_transport=first_order_mass_transport, + dynamic_viscosity=dynamic_viscosity, + ) + mass_changes[particle_index] += mass_change + batch_dm_total += mass_change working_gas_concentration = np.maximum( working_gas_concentration - batch_dm_total, 0.0 ) @@ -1972,18 +2174,10 @@ def mass_transfer_rate( pressure_delta, posinf=0.0, neginf=0.0, nan=0.0 ) - first_order_mass_transport = np.asarray( - first_order_mass_transport, dtype=np.float64 - ) - if pressure_delta.ndim == 2: - if first_order_mass_transport.ndim == 0: - first_order_mass_transport = np.full( - (pressure_delta.shape[0], 1), - float(first_order_mass_transport), - dtype=np.float64, - ) - elif first_order_mass_transport.ndim == 1: - first_order_mass_transport = first_order_mass_transport[:, None] + first_order_mass_transport = _normalize_first_order_mass_transport( + first_order_mass_transport, + pressure_delta, + ) if self._latent_heat_strategy is None: return get_mass_transfer_rate( @@ -2103,12 +2297,17 @@ def _calculate_norm_conc( """ if not np.isfinite(volume) or volume <= 0.0: raise ValueError("volume must be finite and positive.") - if not np.all(np.isfinite(concentration)) or np.any( - concentration < 0.0 - ): + if not np.all(np.isfinite(concentration)): raise ValueError( "concentration must be finite and nonnegative for norm_conc." ) + tolerance = 1e-12 + if np.any(concentration < -tolerance): + raise ValueError( + "concentration must be finite and nonnegative for norm_conc." + ) + if np.any(concentration < 0.0): + concentration = np.maximum(concentration, 0.0) return concentration / volume def _normalize_mass_transfer_shape( @@ -2232,6 +2431,8 @@ def step( Returns: Tuple containing updated particle and gas species objects. """ + _validate_time_step(time_step) + particle_data, particle_is_legacy = _unwrap_particle(particle) gas_data, gas_is_legacy = _unwrap_gas(gas_species) _require_matching_types(particle_is_legacy, gas_is_legacy) diff --git a/particula/dynamics/condensation/tests/condensation_strategies_test.py b/particula/dynamics/condensation/tests/condensation_strategies_test.py index 5dddc400f..1424ca27b 100644 --- a/particula/dynamics/condensation/tests/condensation_strategies_test.py +++ b/particula/dynamics/condensation/tests/condensation_strategies_test.py @@ -15,12 +15,14 @@ CondensationIsothermal, CondensationIsothermalStaggered, CondensationLatentHeat, + _normalize_first_order_mass_transport, _partial_pressure_from_strategy, _pure_vapor_pressure_from_strategy, _require_matching_types, _require_single_box, _unwrap_gas, _unwrap_particle, + _validate_time_step, ) from particula.dynamics.condensation.mass_transfer import ( get_latent_heat_energy_released, @@ -134,6 +136,12 @@ def _make_data_inputs(self) -> tuple[ParticleData, GasData]: from_species(self.gas_species), ) + def test_validate_time_step_allows_zero_and_rejects_invalid(self): + """_validate_time_step accepts zero and rejects invalid values.""" + _validate_time_step(0.0) + with self.assertRaisesRegex(ValueError, "time_step must be finite"): + _validate_time_step(-0.5) + def test_mean_free_path(self): """Test the mean free path call.""" result = self.strategy.mean_free_path( @@ -268,6 +276,28 @@ def test_first_order_mass_transport(self): ) self.assertIsNotNone(result) + def test_normalize_first_order_mass_transport_2d(self): + """Normalize transport shape when pressure delta is 2D.""" + pressure_delta = np.ones((3, 2)) + mass_transport = np.array([1.0, 2.0, 3.0]) + + normalized = _normalize_first_order_mass_transport( + mass_transport, pressure_delta + ) + + self.assertEqual(normalized.shape, (3, 1)) + + def test_step_rejects_invalid_time_step(self): + """step() rejects non-finite or negative timesteps.""" + with self.assertRaisesRegex(ValueError, "time_step must be finite"): + self.strategy.step( + particle=self.particle, + gas_species=self.gas_species, + temperature=self.temperature, + pressure=self.pressure, + time_step=-1.0, + ) + def test_fill_zero_radius(self): """_fill_zero_radius changes zeros to max radius.""" radii = np.array([0.0, 1e-9, 2e-9]) @@ -2631,10 +2661,16 @@ def test_step_continuous_isothermal_parity(self): ) np.testing.assert_allclose( - latent_particle.masses[0], iso_particle.masses[0], rtol=1e-15 + latent_particle.masses[0], + iso_particle.masses[0], + rtol=1e-15, + atol=1e-17, ) np.testing.assert_allclose( - latent_gas.concentration[0], iso_gas.concentration[0], rtol=1e-15 + latent_gas.concentration[0], + iso_gas.concentration[0], + rtol=1e-15, + atol=1e-17, ) def test_step_discrete_energy_tracking(self): @@ -2793,8 +2829,8 @@ def latent_heat(self, temperature: float) -> float: pressure=self.pressure, ) - def test_step_rejects_nonpositive_concentration(self): - """Binned step rejects nonpositive concentrations.""" + def test_step_rejects_negative_concentration(self): + """Binned step rejects negative concentrations.""" particle, particle_data, gas_species, gas_data = ( self._build_binned_representation("pmf", 10) ) @@ -3278,6 +3314,55 @@ def test_step_clamps_gas_concentration_nonnegative(self): np.testing.assert_allclose(gas_new.concentration[0], 0.0, rtol=1e-15) + def test_calculate_norm_conc_clamps_small_negative(self): + """_calculate_norm_conc clamps small negative concentrations.""" + volume = 1.0 + concentration = np.array([1.0, -1.0e-13]) + + norm_conc = self.data_strategy._calculate_norm_conc( + volume, concentration + ) + + self.assertTrue(np.all(norm_conc >= 0.0)) + self.assertEqual(norm_conc[1], 0.0) + + def test_step_rejects_invalid_time_step(self): + """step() rejects non-finite or negative timesteps.""" + particle_data, gas_data = self._make_data_inputs() + with self.assertRaisesRegex(ValueError, "time_step must be finite"): + self.data_strategy.step( + particle=particle_data, + gas_species=gas_data, + temperature=self.temperature, + pressure=self.pressure, + time_step=-0.5, + ) + + def test_get_vapor_pressure_surface_rejects_nonfinite_kelvin(self): + """Non-finite Kelvin terms raise a clear error.""" + particle = copy.deepcopy(self.particle) + gas_species = copy.deepcopy(self.gas_species) + cond = CondensationLatentHeat( + molar_mass=self.molar_mass, + diffusion_coefficient=self.diffusion_coefficient, + accommodation_coefficient=self.accommodation_coefficient, + ) + + with patch.object( + particle.surface, + "kelvin_term", + return_value=np.array([np.nan, np.nan]), + ): + with self.assertRaisesRegex( + ValueError, "kelvin_term must be finite" + ): + cond._get_vapor_pressure_surface( + particle=particle, + gas_species=gas_species, + temperature=self.temperature, + radius=particle.data.radii[0], + ) + def test_step_dynamic_viscosity_pass_through(self): """Dynamic viscosity override is forwarded to mass_transfer_rate.""" particle = copy.deepcopy(self.particle) From 42c24a6091803264720419a713c353c95c99688a Mon Sep 17 00:00:00 2001 From: Kyle Gorkowski Date: Thu, 5 Mar 2026 06:16:43 -0700 Subject: [PATCH 5/5] fix(condensation): support per-particle time steps Handle vectorized time steps in staggered condensation by routing array inputs through shared mass-transfer limiting utilities. This prevents test failures when dt varies per particle and keeps limits consistent across scalar and vectorized paths. --- .../condensation/condensation_strategies.py | 57 +++++++++++++++++-- .../condensation/mass_transfer_utils.py | 4 +- 2 files changed, 54 insertions(+), 7 deletions(-) diff --git a/particula/dynamics/condensation/condensation_strategies.py b/particula/dynamics/condensation/condensation_strategies.py index 62eb08b22..89bdd3b73 100644 --- a/particula/dynamics/condensation/condensation_strategies.py +++ b/particula/dynamics/condensation/condensation_strategies.py @@ -38,6 +38,12 @@ get_mass_transfer_rate, get_mass_transfer_rate_latent_heat, ) +from particula.dynamics.condensation.mass_transfer_utils import ( + apply_condensation_limit, + apply_evaporation_limit, + apply_per_bin_limit, + calc_mass_to_change, +) from particula.gas import get_molecule_mean_free_path from particula.gas.gas_data import GasData from particula.gas.latent_heat_strategies import ( @@ -130,6 +136,37 @@ def _validate_time_step(time_step: float) -> None: raise ValueError("time_step must be finite and nonnegative.") +def _get_mass_transfer_variable_time_step( + mass_rate: NDArray[np.float64], + time_step: NDArray[np.float64], + gas_mass: NDArray[np.float64], + particle_mass: NDArray[np.float64], + particle_concentration: NDArray[np.float64], +) -> NDArray[np.float64]: + """Compute mass transfer for vectorized per-particle time steps.""" + mass_to_change = calc_mass_to_change( + mass_rate=mass_rate, + time_step=time_step, + particle_concentration=particle_concentration, + ) + mass_to_change, evap_sum, neg_mask = apply_condensation_limit( + mass_to_change=mass_to_change, + gas_mass=gas_mass, + ) + mass_to_change = apply_evaporation_limit( + mass_to_change=mass_to_change, + particle_mass=particle_mass, + particle_concentration=particle_concentration, + evap_sum=evap_sum, + neg_mask=neg_mask, + ) + return apply_per_bin_limit( + mass_to_change=mass_to_change, + particle_mass=particle_mass, + particle_concentration=particle_concentration, + ) + + def _pure_vapor_pressure_from_strategy( vapor_pressure_strategy: VaporPressureStrategy | Sequence[VaporPressureStrategy], @@ -1619,7 +1656,7 @@ def _calculate_batch_particle_transfer( gas_concentration: NDArray[np.float64], temperature: float, pressure: float, - dt_local: float, + dt_local: float | NDArray[np.float64], radii: NDArray[np.float64], first_order_mass_transport: Optional[NDArray[np.float64]] = None, dynamic_viscosity: Optional[float] = None, @@ -1701,9 +1738,19 @@ def _calculate_batch_particle_transfer( if particle_mass.ndim == 2 and mass_rate_array.ndim == 1: mass_rate_array = mass_rate_array[:, None] - return get_mass_transfer( + dt_local_array = np.asarray(dt_local, dtype=np.float64) + if dt_local_array.shape == (): + return get_mass_transfer( + mass_rate=mass_rate_array, + time_step=float(dt_local_array), + gas_mass=np.atleast_1d(gas_mass), + particle_mass=particle_mass, + particle_concentration=particle_concentration, + ) + + return _get_mass_transfer_variable_time_step( mass_rate=mass_rate_array, - time_step=dt_local, + time_step=dt_local_array, gas_mass=np.atleast_1d(gas_mass), particle_mass=particle_mass, particle_concentration=particle_concentration, @@ -1825,7 +1872,7 @@ def step( # noqa: C901 continue batch_indices = batch_indices[active_mask] dt_local = dt_local[active_mask] - use_vectorized = np.allclose(dt_local, dt_local[0]) and getattr( + use_vectorized = getattr( self._calculate_single_particle_transfer, "__func__", None ) is ( CondensationIsothermalStaggered._calculate_single_particle_transfer @@ -1881,7 +1928,7 @@ def step( # noqa: C901 continue batch_indices = batch_indices[active_mask] dt_local = dt_local[active_mask] - use_vectorized = np.allclose(dt_local, dt_local[0]) and getattr( + use_vectorized = getattr( self._calculate_single_particle_transfer, "__func__", None ) is ( CondensationIsothermalStaggered._calculate_single_particle_transfer diff --git a/particula/dynamics/condensation/mass_transfer_utils.py b/particula/dynamics/condensation/mass_transfer_utils.py index 56e34f2d0..6f7baa1b7 100644 --- a/particula/dynamics/condensation/mass_transfer_utils.py +++ b/particula/dynamics/condensation/mass_transfer_utils.py @@ -16,7 +16,7 @@ def calc_mass_to_change( mass_rate: NDArray[np.float64], - time_step: float, + time_step: float | NDArray[np.float64], particle_concentration: NDArray[np.float64], ) -> NDArray[np.float64]: """Calculate the requested mass change for every particle/species pair. @@ -30,7 +30,7 @@ def calc_mass_to_change( Args: mass_rate: Mass transfer rate (ṁ) for each particle or ``(particle, species)`` pair in kg s⁻¹. - time_step: Time step Δt in seconds. + time_step: Time step Δt in seconds (scalar or per-particle array). particle_concentration: Number concentration C in m⁻³. Returns: