From 064dc24a6a43854de74681d3e1160298ec138e29 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 01:05:13 +0200 Subject: [PATCH 01/16] Add opt-in adaptive (error-controlled) timestep to the MetSim erosion model The fixed RK4 step (const.dt = 0.005 s) has no convergence basis and is badly under-resolved for fast bodies and small grains: on a fast Orionid single body the dt=0.005 velocity error vs a converged run is ~6 km/s, and the light-curve shape depends on the step size. This adds an opt-in adaptive integrator so accuracy is set by a tolerance rather than a magic constant, while leaving the default behaviour (and every stored dt=0.005 fit) bit-for-bit unchanged. - New Cython adaptiveSingleBodyStep() in MetSimErosionCyTools.pyx advances one fragment across a macro step with per-fragment error-controlled adaptive sub-steps (step-doubling on the existing RK4, PI controller, atmosphere refreshed each sub-step). The m^(-1/3) drag term is floored at m_kill so near-exhaustion grains stay numerically stable. - ablateAll() is gated on const.adaptive_dt: True -> the adaptive stepper; False -> the original fixed-step code path, unchanged. The macro step stays the fixed output cadence (one results row per dt), so the output grid stays uniform and no downstream consumer (GUI residuals, Dynesty likelihood, wake) needs changes. - Erosion grains are shed per sub-step at each sub-step's resolved height/velocity/time, so the erosion light-curve shape stops depending on the step size (grain births match a finely-stepped fixed run; radiated energy conserved to <0.5%). - Constants gets adaptive_dt / adaptive_rtol / adaptive_atol_m / adaptive_atol_v / adaptive_dt_min / adaptive_dt_max / adaptive_max_substeps (default rtol 1e-5, targeting numerical error below measurement noise, not machine precision). runSimulation and loadConstants back-fill these for older JSONs so they default to fixed-step. - GUI: an "adaptive dt" checkbox next to the time-step box; when checked, adaptive_dt is on and the fixed-dt box is disabled (dt is then only the output cadence). Tolerances are set via the sim-fit JSON. - Tests: MetSim/Tests is now an importable package (__init__.py; still excluded from wmpl's import-time walk). Added adaptive-mode regression tests (stability/no-NaN/mass and radiated-energy conservation vs fixed). Validation: adaptive_dt=False reproduces master bit-for-bit (np.array_equal) on the synthetic scenario + 13 fitted Orionids. Convergence study: adaptive at rtol=1e-6 cuts the single-body velocity error from 6.2 to 0.05 km/s (130x); erosion light curves converge with mass/energy conserved. Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/GUI.py | 16 ++ wmpl/MetSim/GUI.ui | 16 ++ wmpl/MetSim/MetSimErosion.py | 255 ++++++++++++++++-------- wmpl/MetSim/MetSimErosionCyTools.pyx | 251 ++++++++++++++++++++++- wmpl/MetSim/Tests/__init__.py | 3 + wmpl/MetSim/Tests/test_MetSimErosion.py | 39 +++- 6 files changed, 492 insertions(+), 88 deletions(-) create mode 100644 wmpl/MetSim/Tests/__init__.py diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index 605f93ad..3b53824f 100644 --- a/wmpl/MetSim/GUI.py +++ b/wmpl/MetSim/GUI.py @@ -3051,6 +3051,7 @@ def __init__(self, traj_path, const_json_file=None, ecsv_files=None, met_path=No self.checkBoxErosionAblationCoeffChange.stateChanged.connect(self.checkBoxErosionAblationCoeffSignal) self.checkBoxDisruption.stateChanged.connect(self.checkBoxDisruptionSignal) self.checkBoxDisruptionErosionCoeff.stateChanged.connect(self.checkBoxDisruptionErosionCoeffSignal) + self.checkBoxAdaptiveDt.stateChanged.connect(self.checkBoxAdaptiveDtSignal) self.runSimButton.clicked.connect(self.runSimulationGUI) @@ -3099,6 +3100,7 @@ def __init__(self, traj_path, const_json_file=None, ecsv_files=None, met_path=No self.checkBoxDisruptionSignal(None) self.checkBoxDisruptionErosionCoeffSignal(None) self.checkBoxDisruptionErosionCoeffSignal(None) + self.checkBoxAdaptiveDtSignal(None) self.toggleWakeNormalizationMethod(None) self.toggleWakeAlignMethod(None) self.toggleFragmentation(None) @@ -3276,6 +3278,11 @@ def updateInputBoxes(self, show_previous=False): ### Simulation params ### self.inputTimeStep.setText(str(const.dt)) + + # Adaptive time step toggle: when on, dt is only the output cadence and its box is disabled + self.checkBoxAdaptiveDt.setChecked(getattr(const, 'adaptive_dt', False)) + self.inputTimeStep.setEnabled(not self.checkBoxAdaptiveDt.isChecked()) + self.inputHtInit.setText("{:.3f}".format(const.h_init/1000)) self.inputP0M.setText("{:d}".format(int(const.P_0m))) self.inputMassKill.setText("{:.1e}".format(const.m_kill)) @@ -3579,6 +3586,14 @@ def checkBoxDisruptionErosionCoeffSignal(self, event): self.readInputBoxes() + def checkBoxAdaptiveDtSignal(self, event): + """ Toggle adaptive per-fragment sub-stepping. When on, dt is only the output cadence, so the + fixed time-step box is disabled (the integration step is chosen automatically). """ + + self.const.adaptive_dt = self.checkBoxAdaptiveDt.isChecked() + self.inputTimeStep.setDisabled(self.const.adaptive_dt) + + def checkBoxFragmentationShowIndividualLCsSignal(self, event): """ Toggle computing light curves of individual fragments during complex fragmentation. """ @@ -3660,6 +3675,7 @@ def readInputBoxes(self): ### Simulation params ### self.const.dt = self._tryReadBox(self.inputTimeStep, self.const.dt) + self.const.adaptive_dt = self.checkBoxAdaptiveDt.isChecked() self.const.P_0m = self._tryReadBox(self.inputP0M, self.const.P_0m) self.const.h_init = 1000*self._tryReadBox(self.inputHtInit, self.const.h_init/1000) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index 3707a47a..0cd16d01 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -89,6 +89,22 @@ + + + + 250 + 30 + 135 + 20 + + + + Use error-controlled adaptive sub-stepping (per fragment) instead of the fixed time step. When on, dt is only the output cadence and the dt box is disabled. Tolerances are set via the sim-fit JSON (adaptive_rtol, adaptive_atol_m/v, ...). + + + adaptive dt + + diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index ea6dda5d..647924d8 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -23,7 +23,7 @@ import pyximport pyximport.install(setup_args={'include_dirs':[np.get_include()]}) from wmpl.MetSim.MetSimErosionCyTools import massLossRK4, decelerationRK4, luminousEfficiency, \ - ionizationEfficiency, atmDensityPoly + ionizationEfficiency, atmDensityPoly, adaptiveSingleBodyStep ### DEFINE CONSTANTS @@ -40,9 +40,27 @@ def __init__(self): ### Simulation parameters ### - # Time step + # Time step. When adaptive_dt is False this is the fixed RK4 integration step; when + # adaptive_dt is True it is only the OUTPUT cadence (one results row per dt, total_time + # advances by dt), while each fragment is integrated with error-controlled adaptive + # sub-steps underneath (see adaptive_* below and ablateAll). self.dt = 0.005 + ### Adaptive sub-stepping (opt-in). Default False -> the engine runs the original fixed-step + # path and reproduces prior results exactly. When True, each fragment sub-steps adaptively + # within each dt to meet the tolerances below, while the output cadence stays fixed at dt. ### + self.adaptive_dt = False + # Tolerance targets NUMERICAL error below MEASUREMENT noise (~0.1 mag, ~0.1 km/s at ~25 FPS), + # not machine convergence. rtol=1e-5 keeps the single-body velocity error ~0.1-0.3 km/s (at or + # below typical measurement precision) at a fraction of the cost of a tighter tol; tighten to + # 1e-6 (~0.05 km/s) for high-precision (e.g. CAMO) data if needed. + self.adaptive_rtol = 1e-5 # relative tolerance on mass and speed + self.adaptive_atol_m = 1e-14 # absolute mass tolerance (kg); ~ m_kill + self.adaptive_atol_v = 0.1 # absolute speed tolerance (m/s); floor well under meas. noise + self.adaptive_dt_min = 1e-7 # smallest allowed sub-step (s) + self.adaptive_dt_max = 0.005 # largest allowed sub-step (s); should be <= dt + self.adaptive_max_substeps = 10000 # runaway guard, sub-steps per fragment per macro step + # Time elapsed since the beginning self.total_time = 0 @@ -350,6 +368,10 @@ def __init__(self): # Identifier of the compex fragmentation entry self.complex_id = None + # Last accepted adaptive sub-step size (s), warm-started across macro steps and copied to + # children by spawn_child(). 0 means "no previous step, start from the macro dt". + self.adaptive_h_sub = 0.0 + def init(self, const, m, rho, v_init, sigma, gamma, zenith_angle, erosion_mass_index, erosion_mass_min, \ erosion_mass_max): @@ -717,10 +739,13 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): ... Note: - const.dt is fixed for every fragment, including grains. Near erosion_mass_min this single-step - RK4 is under-resolved (a single step can overshoot the grain's true velocity by tens of percent - vs. a finely-substepped mirror) - a known accuracy limitation of the fixed-step scheme, not - (currently) compensated for with adaptive/sub-stepping. + With const.adaptive_dt False (default), const.dt is the fixed RK4 step for every fragment, + including grains. Near erosion_mass_min this single-step RK4 is under-resolved (a single step + can overshoot the grain's true velocity by tens of percent vs. a finely-substepped mirror). + Set const.adaptive_dt True to integrate each fragment with error-controlled adaptive sub-steps + (const.adaptive_rtol etc.); dt then only sets the output cadence, and erosion grains are shed + per sub-step so the light-curve shape stops depending on the step size. Fixed-step output is + unchanged when adaptive_dt is False. """ # Keep track of the total luminosity @@ -780,95 +805,122 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): if not frag.active: continue - # Get atmosphere density for the given height - rho_atm = atmDensityPoly(frag.h, const.dens_co) - - # Compute the mass loss of the fragment due to ablation - mass_loss_ablation = massLossRK4(const.dt, frag.K, frag.sigma, frag.m, rho_atm, frag.v) + # Advance this fragment across the macro step. In adaptive mode (const.adaptive_dt) the fragment + # is integrated with error-controlled adaptive sub-steps inside the Cython stepper; otherwise + # the original single fixed RK4 step is taken. Both paths leave frag.m/v/vv/vh/length/h and + # the diagnostic locals (rho_atm, mass_loss_ablation, mass_loss_erosion, deceleration_total) + # set for the shared luminosity/electron-density/event code below. + erosion_events_adaptive = None # per-sub-step erosion shedding (adaptive mode only) + if const.adaptive_dt: + + erosion_active = 1 if (frag.erosion_enabled and (frag.erosion_coeff > 0)) else 0 + + (frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, frag.h, rho_atm, + mass_loss_ablation, mass_loss_erosion, deceleration_total, went_up, n_sub, + frag.adaptive_h_sub, runaway, erosion_events_adaptive) = adaptiveSingleBodyStep( + const.dt, frag.K, frag.sigma, frag.erosion_coeff, erosion_active, + frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, + const.h_init, const.zenith_angle, const.r_earth, G0, const.dens_co, + const.adaptive_rtol, const.adaptive_atol_m, const.adaptive_atol_v, const.m_kill, + const.adaptive_dt_min, const.adaptive_dt_max, const.adaptive_max_substeps, + frag.adaptive_h_sub) + + # Diagnostics (feed the cost study; also used to warn once on a runaway fragment) + const.adaptive_substeps_total += n_sub + if runaway: + const.adaptive_runaway_events += 1 - # Compute the mass loss due to erosion - if frag.erosion_enabled and (frag.erosion_coeff > 0): - mass_loss_erosion = massLossRK4(const.dt, frag.K, frag.erosion_coeff, frag.m, rho_atm, frag.v) else: - mass_loss_erosion = 0 - # Compute the total mass loss - mass_loss_total = mass_loss_ablation + mass_loss_erosion + # Get atmosphere density for the given height + rho_atm = atmDensityPoly(frag.h, const.dens_co) - # If the total mass after ablation in this step is below zero, ablate what's left of the whole mass - # (i.e. land exactly on m_new = 0, not some arbitrary leftover - see m_new below) - if (frag.m + mass_loss_total) < 0: - mass_loss_total = -frag.m + # Compute the mass loss of the fragment due to ablation + mass_loss_ablation = massLossRK4(const.dt, frag.K, frag.sigma, frag.m, rho_atm, frag.v) - # Compute new mass - m_new = frag.m + mass_loss_total + # Compute the mass loss due to erosion + if frag.erosion_enabled and (frag.erosion_coeff > 0): + mass_loss_erosion = massLossRK4(const.dt, frag.K, frag.erosion_coeff, frag.m, rho_atm, frag.v) + else: + mass_loss_erosion = 0 - # Compute change in velocity - deceleration_total = decelerationRK4(const.dt, frag.K, frag.m, rho_atm, frag.v) + # Compute the total mass loss + mass_loss_total = mass_loss_ablation + mass_loss_erosion - # If the deceleration is negative (i.e. the fragment is accelerating), then stop the fragment - if deceleration_total > 0: - frag.vv = frag.vh = frag.v = 0 - deceleration_total = 0 + # If the total mass after ablation in this step is below zero, ablate what's left of the whole mass + # (i.e. land exactly on m_new = 0, not some arbitrary leftover - see m_new below) + if (frag.m + mass_loss_total) < 0: + mass_loss_total = -frag.m - # Otherwise update the velocity - else: + # Compute new mass + m_new = frag.m + mass_loss_total - # Compute g at given height - gv = G0/((1 + frag.h/const.r_earth)**2) + # Compute change in velocity + deceleration_total = decelerationRK4(const.dt, frag.K, frag.m, rho_atm, frag.v) - # ### Add velocity change due to Earth's gravity ### + # If the deceleration is negative (i.e. the fragment is accelerating), then stop the fragment + if deceleration_total > 0: + frag.vv = frag.vh = frag.v = 0 + deceleration_total = 0 - # # Vertical component of a - # av = -gv - deceleration_total*frag.vv/frag.v + frag.vh*frag.v/(const.r_earth + frag.h) + # Otherwise update the velocity + else: + + # Compute g at given height + gv = G0/((1 + frag.h/const.r_earth)**2) + + # ### Add velocity change due to Earth's gravity ### - # # Horizontal component of a - # ah = -deceleration_total*frag.vh/frag.v - frag.vv*frag.v/(const.r_earth + frag.h) + # # Vertical component of a + # av = -gv - deceleration_total*frag.vv/frag.v + frag.vh*frag.v/(const.r_earth + frag.h) - # ### ### + # # Horizontal component of a + # ah = -deceleration_total*frag.vh/frag.v - frag.vv*frag.v/(const.r_earth + frag.h) - ### Compute deceleration without the effects of gravity (to reconstruct the initial velocity - # without the gravity component) + # ### ### - # Vertical component of a - av = -deceleration_total*frag.vv/frag.v + frag.vh*frag.v/(const.r_earth + frag.h) + ### Compute deceleration without the effects of gravity (to reconstruct the initial velocity + # without the gravity component) - # Horizontal component of a - ah = -deceleration_total*frag.vh/frag.v - frag.vv*frag.v/(const.r_earth + frag.h) + # Vertical component of a + av = -deceleration_total*frag.vv/frag.v + frag.vh*frag.v/(const.r_earth + frag.h) - ### + # Horizontal component of a + ah = -deceleration_total*frag.vh/frag.v - frag.vv*frag.v/(const.r_earth + frag.h) - # Compute the drop due to gravity - h_grav_drop = 0.5*gv*const.dt**2 + ### - # Track the total drop due to gravity - frag.h_grav_drop_total += h_grav_drop + # Compute the drop due to gravity + h_grav_drop = 0.5*gv*const.dt**2 - # Update the velocity - frag.vv -= av*const.dt - frag.vh -= ah*const.dt - frag.v = math.sqrt(frag.vh**2 + frag.vv**2) + # Track the total drop due to gravity + frag.h_grav_drop_total += h_grav_drop - # Only allow the meteoroid to go down, and stop the ablation if it stars going up - if frag.vv > 0: + # Update the velocity + frag.vv -= av*const.dt + frag.vh -= ah*const.dt + frag.v = math.sqrt(frag.vh**2 + frag.vv**2) - frag.vv = 0 + # Only allow the meteoroid to go down, and stop the ablation if it stars going up + if frag.vv > 0: - # Setting the height to zero will stop the ablation during the if catch below - frag.h = 0 + frag.vv = 0 - # Update length along the track - frag.length += frag.v*const.dt + # Setting the height to zero will stop the ablation during the if catch below + frag.h = 0 - # Update the mass - frag.m = m_new + # Update length along the track + frag.length += frag.v*const.dt - # Old way of computing height which did not include the curvature of the Earth - # frag.h = frag.h + frag.vv*const.dt + # Update the mass + frag.m = m_new - # Compute the height taking the curvature of the Earth and the gravity drop into account - frag.h = heightCurvature(const.h_init, const.zenith_angle, frag.length, const.r_earth) - frag.h -= frag.h_grav_drop_total + # Old way of computing height which did not include the curvature of the Earth + # frag.h = frag.h + frag.vv*const.dt + + # Compute the height taking the curvature of the Earth and the gravity drop into account + frag.h = heightCurvature(const.h_init, const.zenith_angle, frag.length, const.r_earth) + frag.h -= frag.h_grav_drop_total # Get the luminous efficiency # NOTE: frag.v/frag.m here are already this tick's post-step values, while mass_loss_ablation/ @@ -1044,29 +1096,19 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): # Create grains for erosion-enabled fragments if frag.erosion_enabled: - # Generate new grains if there is some mass to distribute - if abs(mass_loss_erosion) > 0: - - grain_children, const = generateFragments(const, frag, abs(mass_loss_erosion), \ + def _spawnGrainsFromErosion(eroded_mass): + """ Distribute the given eroded mass into grains born from the fragment's current state, + and record erosion-begin bookkeeping for the main fragment. Uses the enclosing + frag/const/dyn_press. """ + grain_children, const_out = generateFragments(const, frag, eroded_mass, \ frag.erosion_mass_index, frag.erosion_mass_min, frag.erosion_mass_max, \ keep_eroding=False, mass_model=const.erosion_grain_distribution) - const.n_active += len(grain_children) - frag_children_all += grain_children - - # print('Eroding id', frag.id) - # print('Eroded mass: {:e}'.format(abs(mass_loss_erosion))) - # print('Mass distribution:') - # grain_mass_sum = 0 - # for f in frag_children: - # print(' {:d}: {:e} kg'.format(f.n_grains, f.m)) - # grain_mass_sum += f.n_grains*f.m - # print('Grain total mass: {:e}'.format(grain_mass_sum)) + frag_children_all.extend(grain_children) # Record physical parameters at the beginning of erosion for the main fragment if frag.main: if const.erosion_beg_vel is None: - const.erosion_beg_vel = frag.v const.erosion_beg_mass = frag.m const.erosion_beg_dyn_press = dyn_press @@ -1075,6 +1117,30 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): elif (const.erosion_height_change >= frag.h) and (const.mass_at_erosion_change is None): const.mass_at_erosion_change = frag.m + if erosion_events_adaptive is not None: + + # Adaptive: shed grains at each sub-step's resolved state, so grain births are spread + # along the trajectory (matching a finely-stepped fixed run) rather than dumped at the + # macro interval's end point. generateFragments()/spawn_child() copies the parent's + # current position/velocity into each grain, so temporarily rewind the parent to each + # sub-step's state, spawn, then restore its end-of-macro-step state. + if erosion_events_adaptive: + _saved_state = (frag.h, frag.v, frag.vv, frag.vh, frag.length, + frag.h_grav_drop_total) + for (em, e_h, e_v, e_vv, e_vh, e_len, e_grav) in erosion_events_adaptive: + if em <= 0: + continue + frag.h = e_h; frag.v = e_v; frag.vv = e_vv; frag.vh = e_vh + frag.length = e_len; frag.h_grav_drop_total = e_grav + _spawnGrainsFromErosion(em) + (frag.h, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total) = _saved_state + + else: + + # Fixed step: distribute the whole macro-step erosion loss at the end-of-step state + if abs(mass_loss_erosion) > 0: + _spawnGrainsFromErosion(abs(mass_loss_erosion)) + # Disrupt the fragment if the dynamic pressure exceeds its strength if frag.disruption_enabled and const.disruption_on: if dyn_press > const.compressive_strength: @@ -1392,6 +1458,21 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): def runSimulation(const, compute_wake=False): """ Run the ablation simulation. """ + # Back-fill adaptive-timestep settings for Constants loaded from older JSONs that predate them, so + # such runs default to the original fixed-step behaviour. Also reset the per-run diagnostics. + _adaptive_defaults = {'adaptive_dt': False, 'adaptive_rtol': 1e-5, 'adaptive_atol_m': 1e-14, + 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, 'adaptive_dt_max': const.dt, + 'adaptive_max_substeps': 10000} + for _attr, _default in _adaptive_defaults.items(): + if not hasattr(const, _attr): + setattr(const, _attr, _default) + const.adaptive_substeps_total = 0 + const.adaptive_runaway_events = 0 + + # The adaptive Cython stepper needs dens_co as a float64 array (memoryview); coerce once (the fixed + # path's atmDensityPoly requires an ndarray too, so this is safe in both modes). + const.dens_co = np.asarray(const.dens_co, dtype=np.float64) + # Ensure that the grain mass min is smaller than the grain mass max if const.erosion_mass_min > const.erosion_mass_max: const.erosion_mass_min, const.erosion_mass_max = const.erosion_mass_max, const.erosion_mass_min @@ -1477,6 +1558,12 @@ def runSimulation(const, compute_wake=False): main_dyn_press]) + # Warn once if any fragment hit the adaptive sub-step cap (under-resolved; raise adaptive_max_substeps + # or loosen the tolerance if this is frequent) + if const.adaptive_dt and (const.adaptive_runaway_events > 0): + print("WARNING: adaptive stepper hit max_substeps ({:d}) on {:d} fragment-step(s); results may be " + "under-resolved there.".format(const.adaptive_max_substeps, const.adaptive_runaway_events)) + # Find the main fragment and return it with results frag_main = None diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index 481e3624..1f2fb930 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/wmpl/MetSim/MetSimErosionCyTools.pyx @@ -6,7 +6,7 @@ cimport cython import numpy as np cimport numpy as np -from libc.math cimport sqrt, M_PI, M_PI_2, atan2, tanh, log, exp, log10 +from libc.math cimport sqrt, M_PI, M_PI_2, atan2, tanh, log, exp, log10, cos, fabs, fmax # Define cython types for numpy arrays @@ -580,4 +580,251 @@ cpdef atmDensityPoly(double ht, np.ndarray[FLOAT_TYPE_t, ndim=1] dens_co): + dens_co[4]*(ht/1e6)**4 + dens_co[5]*(ht/1e6)**5 + dens_co[6]*(ht/1e6)**6 - ) \ No newline at end of file + ) + + + +### Adaptive per-fragment sub-stepping (used only when const.adaptive_dt is True) ### + + +cdef inline double clampMassC(double dm, double m): + """ Reproduce the "ablate at most the whole mass" clamp (MetSimErosion.py mass-loss clamp): + if m + dm < 0, cap the loss at exactly -m so the mass floors at 0. """ + if (m + dm) < 0: + return -m + return dm + + +@cython.cdivision(True) +cdef double heightCurvatureC(double h0, double zc, double l, double r_earth): + """ Cython twin of heightCurvature() (MetSimErosion.py). Scalar hot-path. """ + return sqrt((h0 + r_earth)*(h0 + r_earth) - 2*l*cos(zc)*(h0 + r_earth) + l*l) - r_earth + + +@cython.cdivision(True) +cdef double atmDensityPolyC(double ht, FLOAT_TYPE_t[:] dens_co): + """ Cython/memoryview twin of atmDensityPoly() for use inside the substep loop. """ + cdef double x = ht/1e6 + return 10**(dens_co[0] + dens_co[1]*x + dens_co[2]*x*x + dens_co[3]*x*x*x + + dens_co[4]*x*x*x*x + dens_co[5]*x*x*x*x*x + dens_co[6]*x*x*x*x*x*x) + + +@cython.cdivision(True) +cdef void advanceVelPosC(double m, double v, double vv, double vh, double length, double grav, + double dm, double decel_rate, double h, double h_at, + double r_earth, double g0, double* out): + """ Advance velocity components, speed, along-track length, and gravity-drop for one sub-step of + size h, reproducing MetSimErosion.py:773-824 exactly (velocity updated BEFORE length; gravity + only drops height, does not enter the velocity magnitude). out layout: + [m_new, v_new, vv_new, vh_new, length_new, grav_new, went_up_flag]. """ + cdef double gv, av, ah, vv_n, vh_n, v_n + # Accelerating (decel_rate > 0) or already stopped -> stop the fragment (mirror 773-775) + if (decel_rate > 0) or (v <= 0): + out[0] = m + dm; out[1] = 0.0; out[2] = 0.0; out[3] = 0.0 + out[4] = length; out[5] = grav; out[6] = 0.0 + return + gv = g0/((1.0 + h_at/r_earth)*(1.0 + h_at/r_earth)) + av = -decel_rate*vv/v + vh*v/(r_earth + h_at) + ah = -decel_rate*vh/v - vv*v/(r_earth + h_at) + vv_n = vv - av*h + vh_n = vh - ah*h + v_n = sqrt(vv_n*vv_n + vh_n*vh_n) + out[0] = m + dm + out[1] = v_n + out[2] = vv_n + out[3] = vh_n + out[4] = length + v_n*h # length uses the UPDATED speed (matches 811-813 -> 824) + out[5] = grav + 0.5*gv*h*h + out[6] = 1.0 if vv_n > 0 else 0.0 # going up + + +@cython.cdivision(True) +cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double erosion_coeff, + int erosion_active, double m, double v, double vv, double vh, double length, + double h_grav_drop_total, double h_init, double zenith_angle, double r_earth, double g0, + FLOAT_TYPE_t[:] dens_co, double rtol, double atol_m, double atol_v, double m_kill, + double dt_min, double dt_max, int max_substeps, double h_sub_init): + """ Advance ONE fragment across a full macro interval dt_macro using error-controlled adaptive + sub-steps (step-doubling on the existing RK4), refreshing atmosphere/height each sub-step. + Reproduces the single-body advance of MetSimErosion.py (748-834) in the one-substep limit but + drives the local error below (rtol, atol). Events/kills stay at the macro boundary (handled by + the caller). Returns a tuple: + (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, + dm_abl_macro, dm_ero_macro, decel_return, went_up, n_substeps, h_sub_last, runaway) + where dm_*_macro are the (unclamped, per-species) accumulated mass losses used for the + luminosity/electron-density/grain diagnostics, and decel_return is the macro-averaged dv/dt + (negative), consistent with the sign of decelerationRK4. """ + + cdef double t = 0.0 + cdef double h_sub, h_cur, rho_atm, rho_mid, rho_last + cdef double dm_abl_big, dm_ero_big, dm_tot_big, decel_big + cdef double dm_abl_1, dm_ero_1, dm_tot_1, decel_1 + cdef double dm_abl_2, dm_ero_2, dm_tot_2, decel_2 + cdef double m_h, v_h, vv_h, vh_h, len_h, grav_h, h_mid + cdef double m_big, v_big, m_two, v_two, hh + cdef double err_m, err_v, sc_m, sc_v, E, E_prev, fac + cdef double dm_abl_macro = 0.0 + cdef double dm_ero_macro = 0.0 + cdef double v_start = v + cdef int n_substeps = 0 + cdef int went_up = 0 + cdef int runaway = 0 + cdef int at_floor + cdef double out1[7] + cdef double out2[7] + cdef double outb[7] + cdef double h_new, decel_return, rho_final + + cdef double safety = 0.9 + cdef double facmin = 0.2 + cdef double facmax = 5.0 + + # Per-sub-step erosion shedding events (only for eroding fragments), so that ablateAll can spawn + # grains at each sub-step's resolved state instead of dumping the whole macro interval's eroded + # mass at one point. Each entry: (eroded_mass, h, v, vv, vh, length, h_grav_drop_total). None for + # non-eroding fragments (the vast majority) to keep their fast path allocation-free. + erosion_events = [] if erosion_active else None + + h_sub = h_sub_init + if h_sub <= 0: + h_sub = dt_macro + if h_sub > dt_max: + h_sub = dt_max + if h_sub < dt_min: + h_sub = dt_min + E_prev = 1.0 + rho_last = 0.0 + + while t < dt_macro: + + # Clamp the final sub-step so sub-steps sum to exactly dt_macro (no overshoot) + if t + h_sub > dt_macro: + h_sub = dt_macro - t + + hh = 0.5*h_sub + + # Height and atmosphere at the CURRENT state (refresh -> shrinks the operator-split/frozen-rho error) + h_cur = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop_total + rho_atm = atmDensityPolyC(h_cur, dens_co) + rho_last = rho_atm + + # --- Big step (size h_sub) --- + dm_abl_big = massLossRK4(h_sub, K, sigma, m, rho_atm, v) + if erosion_active: + dm_ero_big = massLossRK4(h_sub, K, erosion_coeff, m, rho_atm, v) + else: + dm_ero_big = 0.0 + dm_tot_big = clampMassC(dm_abl_big + dm_ero_big, m) + # Floor the deceleration mass at m_kill: deceleration ~ m^(-1/3) diverges as m -> 0, and a grain + # at/under m_kill is treated as dead (killed at the macro boundary), so this only bounds the + # drag on an already-exhausted grain rather than letting it blow up. + decel_big = decelerationRK4(h_sub, K, fmax(m, m_kill), rho_atm, v) + advanceVelPosC(m, v, vv, vh, length, h_grav_drop_total, dm_tot_big, decel_big, h_sub, h_cur, + r_earth, g0, outb) + m_big = outb[0] + v_big = outb[1] + + # --- Two half steps (hh each; rho refreshed at the midpoint) --- + dm_abl_1 = massLossRK4(hh, K, sigma, m, rho_atm, v) + if erosion_active: + dm_ero_1 = massLossRK4(hh, K, erosion_coeff, m, rho_atm, v) + else: + dm_ero_1 = 0.0 + dm_tot_1 = clampMassC(dm_abl_1 + dm_ero_1, m) + decel_1 = decelerationRK4(hh, K, fmax(m, m_kill), rho_atm, v) + advanceVelPosC(m, v, vv, vh, length, h_grav_drop_total, dm_tot_1, decel_1, hh, h_cur, + r_earth, g0, out1) + m_h = out1[0]; v_h = out1[1]; vv_h = out1[2]; vh_h = out1[3]; len_h = out1[4]; grav_h = out1[5] + + h_mid = heightCurvatureC(h_init, zenith_angle, len_h, r_earth) - grav_h + rho_mid = atmDensityPolyC(h_mid, dens_co) + + dm_abl_2 = massLossRK4(hh, K, sigma, m_h, rho_mid, v_h) + if erosion_active: + dm_ero_2 = massLossRK4(hh, K, erosion_coeff, m_h, rho_mid, v_h) + else: + dm_ero_2 = 0.0 + dm_tot_2 = clampMassC(dm_abl_2 + dm_ero_2, m_h) + decel_2 = decelerationRK4(hh, K, fmax(m_h, m_kill), rho_mid, v_h) + advanceVelPosC(m_h, v_h, vv_h, vh_h, len_h, grav_h, dm_tot_2, decel_2, hh, h_mid, + r_earth, g0, out2) + m_two = out2[0]; v_two = out2[1] + + # --- Error estimate (step doubling, RK4 order p=4 -> denom 2^p - 1 = 15) --- + err_m = fabs(m_two - m_big)/15.0 + err_v = fabs(v_two - v_big)/15.0 + sc_m = atol_m + rtol*fmax(fabs(m_two), fabs(m_big)) + sc_v = atol_v + rtol*fmax(fabs(v_two), fabs(v_big)) + E = sqrt(0.5*((err_m/sc_m)*(err_m/sc_m) + (err_v/sc_v)*(err_v/sc_v))) + + at_floor = 1 if h_sub <= dt_min*(1.0 + 1e-12) else 0 + + if (E <= 1.0) or at_floor: + # Accept the two-half (more accurate) state; accumulate per-species mass loss (unclamped) + dm_abl_macro += dm_abl_1 + dm_abl_2 + dm_ero_macro += dm_ero_1 + dm_ero_2 + m = out2[0]; v = out2[1]; vv = out2[2]; vh = out2[3] + length = out2[4]; h_grav_drop_total = out2[5] + t += h_sub + n_substeps += 1 + + # Record the eroded mass shed on this sub-step, tagged with the fragment's just-advanced + # state, so ablateAll can spawn grains here (finely resolved along the trajectory) rather + # than dumping the whole macro interval's eroded mass at the end point. + if erosion_active and ((dm_ero_1 + dm_ero_2) < 0): + erosion_events.append(( + -(dm_ero_1 + dm_ero_2), + heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop_total, + v, vv, vh, length, h_grav_drop_total)) + + if m <= m_kill: # grain exhausted -> freeze, killed at the macro boundary + break + if out2[6] > 0.5: # turned upward mid-interval -> freeze, kill at macro boundary + vv = 0.0 + went_up = 1 + break + if v <= 0: # stopped (accelerating/decel guard) -> freeze + break + if n_substeps >= max_substeps: + runaway = 1 + break + + # PI step-size controller (k = p + 1 = 5) + if E <= 0: + E = 1e-10 + fac = safety*(E**(-0.7/5.0))*(E_prev**(0.4/5.0)) + E_prev = E + if fac < facmin: + fac = facmin + if fac > facmax: + fac = facmax + h_sub = h_sub*fac + if h_sub > dt_max: + h_sub = dt_max + if h_sub < dt_min: + h_sub = dt_min + else: + # Reject: shrink and retry the SAME sub-step (do not advance t) + fac = safety*(E**(-1.0/5.0)) + if fac < facmin: + fac = facmin + h_sub = h_sub*fac + if h_sub < dt_min: + h_sub = dt_min + + if went_up: + h_new = 0.0 + else: + h_new = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop_total + + # Density for the end-of-interval dynamic pressure; fall back to the last in-loop value if h<=0 + if h_new > 0: + rho_final = atmDensityPolyC(h_new, dens_co) + else: + rho_final = rho_last + + decel_return = -(v_start - v)/dt_macro # macro-averaged dv/dt (negative), matches decelerationRK4 sign + + return (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, + dm_abl_macro, dm_ero_macro, decel_return, went_up, n_substeps, h_sub, runaway, + erosion_events) \ No newline at end of file diff --git a/wmpl/MetSim/Tests/__init__.py b/wmpl/MetSim/Tests/__init__.py new file mode 100644 index 00000000..f9474b43 --- /dev/null +++ b/wmpl/MetSim/Tests/__init__.py @@ -0,0 +1,3 @@ +""" MetSim test package. Kept importable so the scenario builders / regression tests can be reused +programmatically. Excluded from wmpl's import-time submodule walk (see wmpl/__init__.py) so that +`import wmpl` never triggers test-module side effects (e.g. matplotlib.use("Agg")). """ diff --git a/wmpl/MetSim/Tests/test_MetSimErosion.py b/wmpl/MetSim/Tests/test_MetSimErosion.py index 0b12f36f..485f9fd9 100644 --- a/wmpl/MetSim/Tests/test_MetSimErosion.py +++ b/wmpl/MetSim/Tests/test_MetSimErosion.py @@ -103,9 +103,13 @@ def _makeComplexScenarioConstants(m_init=0.5, v_init=16000.0, h_init=80000.0, ze return const -def _runComplexScenario(): - """ Run the complex scenario once and return (const, results) with results as a float ndarray. """ +def _runComplexScenario(adaptive=False, rtol=1e-5): + """ Run the complex scenario once and return (const, results) with results as a float ndarray. + If adaptive is True, use error-controlled adaptive sub-stepping. """ const = _makeComplexScenarioConstants() + if adaptive: + const.adaptive_dt = True + const.adaptive_rtol = rtol _, results_list, _ = MetSimErosion.runSimulation(const) return const, np.array(results_list, dtype=float) @@ -158,6 +162,35 @@ def test_brightest_height_not_spuriously_zero_while_luminous(): "brightest_height was 0 on a tick where the meteor was still luminous" +def test_adaptive_dt_runs_and_stays_stable(): + """ Adaptive sub-stepping must run the same complex scenario without NaN/inf, without negative + mass, and with eroded/disrupted light tracked and brightest_height valid while luminous - i.e. + small grains stay numerically stable while ablating (the m^(-1/3) drag near mass exhaustion is + floored, not blown up). """ + const, results = _runComplexScenario(adaptive=True) + assert const.adaptive_dt is True + assert const.disruption_height > 0, "disruption not triggered in adaptive mode" + assert np.all(np.isfinite(results)), "non-finite output in adaptive mode (grain instability?)" + assert np.all(results[:, COL_MAIN_MASS] >= 0), "main mass went negative (adaptive)" + assert np.all(results[:, COL_MASS_TOTAL_ACTIVE] >= 0), "total active mass went negative (adaptive)" + assert np.nanmax(results[:, COL_LUM_ERODED]) > 0, "eroded light not tracked (adaptive)" + luminous = results[:, COL_LUM_TOTAL] > 0 + assert np.all(results[luminous, COL_BRIGHTEST_HEIGHT] > 0), \ + "brightest_height spuriously 0 while luminous (adaptive)" + + +def test_adaptive_conserves_radiated_energy_vs_fixed(): + """ Turning on adaptive sub-stepping must not change the bulk energetics: the time-integrated + total luminosity should match the fixed-step run to a few percent (same physics, finer + integration + finer grain shedding). """ + _, r_fixed = _runComplexScenario(adaptive=False) + _, r_adapt = _runComplexScenario(adaptive=True) + e_fixed = np.trapz(r_fixed[:, COL_LUM_TOTAL], r_fixed[:, COL_TIME]) + e_adapt = np.trapz(r_adapt[:, COL_LUM_TOTAL], r_adapt[:, COL_TIME]) + assert abs(e_adapt - e_fixed)/e_fixed < 0.05, \ + "adaptive changed total radiated energy by >5% vs fixed (energy/mass not conserved?)" + + def _savePlot(save_path=None): """ New-engine-only diagnostic figure (LC, mass, velocity, height vs time). Optional, for eyeballing. """ import matplotlib @@ -197,6 +230,8 @@ def _savePlot(save_path=None): test_total_active_mass_is_grain_weighted() test_lum_eroded_tracked_by_default() test_brightest_height_not_spuriously_zero_while_luminous() + test_adaptive_dt_runs_and_stays_stable() + test_adaptive_conserves_radiated_energy_vs_fixed() print("All MetSimErosion regression checks passed.") if "--plot" in sys.argv: From aaef6713faa6a7ee3f3677301e216bd4000d0afb Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 09:23:08 +0200 Subject: [PATCH 02/16] Improve adaptive-dt GUI: relocate checkbox, add rtol input The "adaptive" checkbox previously overlapped the m_kill "kg" unit label and the len_kill field in the packed first row. Move the adaptive controls into the free right-hand column: "adaptive" checkbox on row 2, an "rtol" input box on row 3 (enabled only when adaptive is on), so tolerance is tunable from the GUI instead of JSON-only. Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/GUI.py | 7 +++++++ wmpl/MetSim/GUI.ui | 39 ++++++++++++++++++++++++++++++++++----- 2 files changed, 41 insertions(+), 5 deletions(-) diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index 3b53824f..d62dfd34 100644 --- a/wmpl/MetSim/GUI.py +++ b/wmpl/MetSim/GUI.py @@ -3282,6 +3282,7 @@ def updateInputBoxes(self, show_previous=False): # Adaptive time step toggle: when on, dt is only the output cadence and its box is disabled self.checkBoxAdaptiveDt.setChecked(getattr(const, 'adaptive_dt', False)) self.inputTimeStep.setEnabled(not self.checkBoxAdaptiveDt.isChecked()) + self.inputAdaptiveRtol.setText("{:.1e}".format(getattr(const, 'adaptive_rtol', 1e-5))) self.inputHtInit.setText("{:.3f}".format(const.h_init/1000)) self.inputP0M.setText("{:d}".format(int(const.P_0m))) @@ -3593,6 +3594,10 @@ def checkBoxAdaptiveDtSignal(self, event): self.const.adaptive_dt = self.checkBoxAdaptiveDt.isChecked() self.inputTimeStep.setDisabled(self.const.adaptive_dt) + # rtol only matters in adaptive mode: enable its box/label with the checkbox, disable otherwise + self.inputAdaptiveRtol.setEnabled(self.const.adaptive_dt) + self.labelAdaptiveRtol.setEnabled(self.const.adaptive_dt) + def checkBoxFragmentationShowIndividualLCsSignal(self, event): """ Toggle computing light curves of individual fragments during complex fragmentation. """ @@ -3676,6 +3681,8 @@ def readInputBoxes(self): self.const.dt = self._tryReadBox(self.inputTimeStep, self.const.dt) self.const.adaptive_dt = self.checkBoxAdaptiveDt.isChecked() + self.const.adaptive_rtol = self._tryReadBox(self.inputAdaptiveRtol, + getattr(self.const, 'adaptive_rtol', 1e-5)) self.const.P_0m = self._tryReadBox(self.inputP0M, self.const.P_0m) self.const.h_init = 1000*self._tryReadBox(self.inputHtInit, self.const.h_init/1000) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index 0cd16d01..5bed2086 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -92,17 +92,46 @@ - 250 - 30 - 135 + 303 + 60 + 88 20 - Use error-controlled adaptive sub-stepping (per fragment) instead of the fixed time step. When on, dt is only the output cadence and the dt box is disabled. Tolerances are set via the sim-fit JSON (adaptive_rtol, adaptive_atol_m/v, ...). + Use error-controlled adaptive sub-stepping (per fragment) instead of the fixed time step. When on, dt is only the output cadence, the dt box is disabled, and the integration step is chosen automatically to meet the relative tolerance (rtol). + + + adaptive + + + + + + 300 + 90 + 28 + 16 + - adaptive dt + rtol + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + 330 + 90 + 58 + 23 + + + + Relative tolerance for the adaptive step (default 1e-5). Smaller = more accurate + slower. Only used when 'adaptive' is checked. From 289666c0f11b919bbdcf29522685db8ce94eb567 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 10:03:05 +0200 Subject: [PATCH 03/16] Adjust adaptive-dt GUI layout; fix Qt6-scoped enums; expand rtol tooltip - Designer layout tweaks to the simulation settings group. - Normalize Qt6-style fully-scoped enum names (Qt::AlignmentFlag::*, Qt::Orientation::*) back to the short form so PyQt5's loadUi can parse the .ui. - Expand the rtol tooltip with concrete tuning guidance (accuracy vs cost, target below measurement noise). Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/GUI.ui | 157 ++++++++++++++++++++++++--------------------- 1 file changed, 85 insertions(+), 72 deletions(-) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index 5bed2086..46ec881f 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -89,51 +89,6 @@ - - - - 303 - 60 - 88 - 20 - - - - Use error-controlled adaptive sub-stepping (per fragment) instead of the fixed time step. When on, dt is only the output cadence, the dt box is disabled, and the integration step is chosen automatically to meet the relative tolerance (rtol). - - - adaptive - - - - - - 300 - 90 - 28 - 16 - - - - rtol - - - Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter - - - - - - 330 - 90 - 58 - 23 - - - - Relative tolerance for the adaptive step (default 1e-5). Smaller = more accurate + slower. Only used when 'adaptive' is checked. - - @@ -153,7 +108,7 @@ - 113 + 243 30 57 15 @@ -169,7 +124,7 @@ - 173 + 303 30 61 23 @@ -180,7 +135,7 @@ 113 - 60 + 87 57 15 @@ -196,7 +151,7 @@ 173 - 60 + 87 61 23 @@ -205,8 +160,8 @@ - 173 - 90 + 303 + 60 61 23 @@ -215,8 +170,8 @@ - 113 - 90 + 243 + 60 57 15 @@ -231,8 +186,8 @@ - -20 - 60 + 114 + 30 57 15 @@ -247,8 +202,8 @@ - 40 - 60 + 174 + 30 61 23 @@ -257,7 +212,7 @@ - 240 + 370 30 57 15 @@ -271,7 +226,7 @@ 240 - 60 + 90 57 15 @@ -283,8 +238,8 @@ - 240 - 90 + 370 + 60 57 15 @@ -309,8 +264,8 @@ - 110 - 60 + 244 + 30 21 16 @@ -322,8 +277,8 @@ - -25 - 90 + 109 + 60 61 16 @@ -338,8 +293,8 @@ - 110 - 90 + 244 + 60 21 16 @@ -351,8 +306,8 @@ - 40 - 90 + 174 + 60 61 23 @@ -362,7 +317,7 @@ 303 - 30 + 90 61 23 @@ -372,7 +327,7 @@ 243 - 30 + 87 57 15 @@ -388,7 +343,7 @@ 370 - 30 + 90 57 15 @@ -397,6 +352,64 @@ km + + + + 23 + 60 + 88 + 20 + + + + Use error-controlled adaptive sub-stepping (per fragment) instead of the fixed time step. When on, dt is only the output cadence, the dt box is disabled, and the integration step is chosen automatically to meet the relative tolerance (rtol). + + + adaptive dt + + + + + + 40 + 90 + 61 + 23 + + + + Relative error tolerance for the adaptive time step (only used when 'adaptive' is checked). Each sub-step is accepted only if its estimated relative error in mass and speed is below rtol; otherwise the step is shrunk and retried. Smaller rtol = more accurate but more sub-steps (slower). Rough guide for the single-body velocity error: 1e-4 fast/coarse (~1 km/s), 1e-5 default (~0.3 km/s), 1e-6 high-precision / CAMO-grade (~0.05 km/s). Set it so the numerical error sits below your measurement noise. Ignored when 'adaptive' is off. + + + + + + 10 + 90 + 28 + 16 + + + + rtol + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + 120 + 20 + 20 + 91 + + + + Qt::Vertical + + From 8b58cdf17672fec049773842dff1d9292db08418 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 10:12:50 +0200 Subject: [PATCH 04/16] Fix adaptive-step luminosity: use drag-only deceleration, not net dv The adaptive stepper fed the net macro-step speed change (v_start - v)/dt into the luminosity deceleration term. For a massive, high-altitude fragment (drag ~ 0) that net change is dominated by the gravity/curvature velocity reallocation and can flip sign, making the m*v*decel term (and thus lum) large and negative. The fragment was then killed on tick 1 by the lum<0 check, yielding an empty simulation and a downstream interp1d crash in the GUI when adaptive stepping was enabled. Accumulate the drag-only speed change (sum of decelerationRK4 * substep) instead, matching the fixed path's use of drag deceleration. Verified on a 2100 kg fragmentation fireball: adaptive now matches the fixed light curve to <0.05 mag (was: 1 row, lum<0, crash). Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/MetSimErosionCyTools.pyx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index 1f2fb930..1bbc6acd 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/wmpl/MetSim/MetSimErosionCyTools.pyx @@ -665,6 +665,7 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero cdef double err_m, err_v, sc_m, sc_v, E, E_prev, fac cdef double dm_abl_macro = 0.0 cdef double dm_ero_macro = 0.0 + cdef double dv_drag = 0.0 # accumulated DRAG-only speed change (excludes gravity/curvature) cdef double v_start = v cdef int n_substeps = 0 cdef int went_up = 0 @@ -763,6 +764,10 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero # Accept the two-half (more accurate) state; accumulate per-species mass loss (unclamped) dm_abl_macro += dm_abl_1 + dm_abl_2 dm_ero_macro += dm_ero_1 + dm_ero_2 + # Accumulate the DRAG-only speed change over this sub-step (decel_1/decel_2 are the pure-drag + # dv/dt from decelerationRK4). This feeds the luminosity deceleration term - unlike the net + # (v_start - v), it excludes the gravity/curvature reallocation, matching the fixed path. + dv_drag += (decel_1 + decel_2)*hh m = out2[0]; v = out2[1]; vv = out2[2]; vh = out2[3] length = out2[4]; h_grav_drop_total = out2[5] t += h_sub @@ -823,7 +828,7 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero else: rho_final = rho_last - decel_return = -(v_start - v)/dt_macro # macro-averaged dv/dt (negative), matches decelerationRK4 sign + decel_return = dv_drag/dt_macro # macro-averaged DRAG dv/dt (negative), matches decelerationRK4 sign return (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, dm_abl_macro, dm_ero_macro, decel_return, went_up, n_substeps, h_sub, runaway, From 8b3d3cbc5668550c34fbdce7317544382315db38 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 11:47:37 +0200 Subject: [PATCH 05/16] Harden GUI against degenerate simulations (no more core dump) A simulation that dies immediately (0-1 rows, or all-NaN magnitudes) previously crashed the whole app with an interp1d "reshape array of size 0" error. Add simulationUsable() and skip the interpolation/plotting chain with a clear warning when the current (or previous) run has too few valid points or no positive luminosity. Also fall back to the full arrays inside updateInterpolations() when too little of the trajectory reaches past the observed begin, so the normalized interpolators still build instead of crashing. Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/GUI.py | 50 ++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 4 deletions(-) diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index d62dfd34..ff48f13c 100644 --- a/wmpl/MetSim/GUI.py +++ b/wmpl/MetSim/GUI.py @@ -4317,9 +4317,16 @@ def updateInterpolations(self, show_previous=False): norm_sim_time = sr.time_arr - self.sim_time_beg - self.norm_sim_ht = sr.leading_frag_height_arr[norm_sim_len > 0] - self.norm_sim_time = norm_sim_time[norm_sim_len > 0] - self.norm_sim_len = norm_sim_len[norm_sim_len > 0] + # Keep only the part of the trajectory past the observed begin (norm_sim_len > 0). If too + # little of the simulation reaches there (< 2 points), fall back to the full arrays so the + # interpolators still build instead of crashing on an empty array. + mask = norm_sim_len > 0 + if np.count_nonzero(mask) < 2: + mask = np.ones_like(norm_sim_len, dtype=bool) + + self.norm_sim_ht = sr.leading_frag_height_arr[mask] + self.norm_sim_time = norm_sim_time[mask] + self.norm_sim_len = norm_sim_len[mask] # Interpolate the normalized length by time self.sim_norm_len_interp = scipy.interpolate.interp1d(self.norm_sim_time, self.norm_sim_len, @@ -5237,7 +5244,7 @@ def updateWakePlot(self, show_previous=False): def showPreviousResults(self): """ Show previous simulation results and parameters. """ - if self.simulation_results_prev is not None: + if self.simulationUsable(self.simulation_results_prev): self.updateInputBoxes(show_previous=True) self.updateInterpolations(show_previous=True) @@ -5248,10 +5255,45 @@ def showPreviousResults(self): + def simulationUsable(self, sr): + """ Check that a SimulationResults has enough valid points to interpolate/plot. A degenerate + run (e.g. the meteoroid is killed on the first tick, producing 0-1 rows or all-NaN + magnitudes) would otherwise crash the interpolation/plotting chain. """ + + if sr is None: + return False + + try: + # Need at least two samples spanning a height range to build the interpolators + if len(sr.time_arr) < 2: + return False + finite_ht = np.isfinite(sr.leading_frag_height_arr) + if np.count_nonzero(finite_ht) < 2: + return False + if np.nanmin(sr.leading_frag_height_arr) == np.nanmax(sr.leading_frag_height_arr): + return False + # Need at least some real luminosity to compute magnitudes + if not np.any(sr.luminosity_arr > 0): + return False + except (AttributeError, TypeError, ValueError): + return False + + return True + + def showCurrentResults(self): """ Show current simulation results and parameters. """ self.updateInputBoxes(show_previous=False) + + # Guard against a degenerate simulation (e.g. killed on the first tick): skip the + # interpolation/plotting chain with a message instead of crashing the whole app + if not self.simulationUsable(self.simulation_results): + print("WARNING: the simulation produced no usable output (the meteoroid may be dying " + "immediately - check the initial parameters, e.g. velocity/mass/erosion/adaptive " + "settings). Skipping the plot update.") + return + self.updateInterpolations(show_previous=False) self.updateMagnitudePlot(show_previous=False) self.updateVelocityPlot(show_previous=False) From 0619a7551652ee713a1f606098087ad6e6c98356 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 14:28:25 +0200 Subject: [PATCH 06/16] Adaptive stepper review fixes: warm-start, going-up consistency, floor diagnostic From the PR #86 review: - F2: carry the controller's natural (pre-boundary-clamp) sub-step as the next macro step's warm-start instead of the tiny clamped leftover. Cuts sub-steps ~2-5x (Orionid 3436->740/row, fireball ~1673->527/row) with no effect on the visible light curve - only the faint, sub-detection tail shifts slightly. - F1: on 'going up' (vv>0) recompute h_new from the along-track length instead of forcing 0. The fixed path's frag.h=0 in that branch is immediately overwritten by the same recompute, so forcing 0 was killing upward-turning end-of-life fragments a tick earlier than fixed mode. - F3: count sub-steps accepted at the dt_min floor while still over tolerance (const.adaptive_floor_accepts) and warn once in runSimulation, so under-resolution is visible. - Record the erosion-begin (vel, mass, dyn_press) triple once at the end-of-step state so it is mutually consistent in adaptive mode (was mixing a rewound sub-step velocity with end-of-step mass/dyn_press). Backward compat unchanged: adaptive_dt=False is bit-for-bit vs master (13 Orionids + synthetic); regression suite passes. Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/MetSimErosion.py | 42 ++++++++++++++++++---------- wmpl/MetSim/MetSimErosionCyTools.pyx | 36 +++++++++++++++++++----- 2 files changed, 56 insertions(+), 22 deletions(-) diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 647924d8..2e41c441 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -817,7 +817,8 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): (frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, frag.h, rho_atm, mass_loss_ablation, mass_loss_erosion, deceleration_total, went_up, n_sub, - frag.adaptive_h_sub, runaway, erosion_events_adaptive) = adaptiveSingleBodyStep( + frag.adaptive_h_sub, runaway, floor_accepts, erosion_events_adaptive) \ + = adaptiveSingleBodyStep( const.dt, frag.K, frag.sigma, frag.erosion_coeff, erosion_active, frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, const.h_init, const.zenith_angle, const.r_earth, G0, const.dens_co, @@ -825,8 +826,9 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): const.adaptive_dt_min, const.adaptive_dt_max, const.adaptive_max_substeps, frag.adaptive_h_sub) - # Diagnostics (feed the cost study; also used to warn once on a runaway fragment) + # Diagnostics (feed the cost study; also used to warn once on runaway/under-resolved steps) const.adaptive_substeps_total += n_sub + const.adaptive_floor_accepts += floor_accepts if runaway: const.adaptive_runaway_events += 1 @@ -1097,25 +1099,15 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): if frag.erosion_enabled: def _spawnGrainsFromErosion(eroded_mass): - """ Distribute the given eroded mass into grains born from the fragment's current state, - and record erosion-begin bookkeeping for the main fragment. Uses the enclosing - frag/const/dyn_press. """ + """ Distribute the given eroded mass into grains born from the fragment's current state. + Uses the enclosing frag/const. """ grain_children, const_out = generateFragments(const, frag, eroded_mass, \ frag.erosion_mass_index, frag.erosion_mass_min, frag.erosion_mass_max, \ keep_eroding=False, mass_model=const.erosion_grain_distribution) const.n_active += len(grain_children) frag_children_all.extend(grain_children) - # Record physical parameters at the beginning of erosion for the main fragment - if frag.main: - if const.erosion_beg_vel is None: - const.erosion_beg_vel = frag.v - const.erosion_beg_mass = frag.m - const.erosion_beg_dyn_press = dyn_press - - # Record the mass when erosion is changed - elif (const.erosion_height_change >= frag.h) and (const.mass_at_erosion_change is None): - const.mass_at_erosion_change = frag.m + eroded_this_tick = False if erosion_events_adaptive is not None: @@ -1133,6 +1125,7 @@ def _spawnGrainsFromErosion(eroded_mass): frag.h = e_h; frag.v = e_v; frag.vv = e_vv; frag.vh = e_vh frag.length = e_len; frag.h_grav_drop_total = e_grav _spawnGrainsFromErosion(em) + eroded_this_tick = True (frag.h, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total) = _saved_state else: @@ -1140,6 +1133,20 @@ def _spawnGrainsFromErosion(eroded_mass): # Fixed step: distribute the whole macro-step erosion loss at the end-of-step state if abs(mass_loss_erosion) > 0: _spawnGrainsFromErosion(abs(mass_loss_erosion)) + eroded_this_tick = True + + # Record erosion-begin bookkeeping for the main fragment once, at the END-of-step state so the + # (vel, mass, dyn_press) triple is mutually consistent in both modes (in adaptive mode frag + # has been restored above; the per-event rewind must not leak into these diagnostics) + if eroded_this_tick and frag.main: + if const.erosion_beg_vel is None: + const.erosion_beg_vel = frag.v + const.erosion_beg_mass = frag.m + const.erosion_beg_dyn_press = dyn_press + + # Record the mass when erosion is changed + elif (const.erosion_height_change >= frag.h) and (const.mass_at_erosion_change is None): + const.mass_at_erosion_change = frag.m # Disrupt the fragment if the dynamic pressure exceeds its strength if frag.disruption_enabled and const.disruption_on: @@ -1468,6 +1475,7 @@ def runSimulation(const, compute_wake=False): setattr(const, _attr, _default) const.adaptive_substeps_total = 0 const.adaptive_runaway_events = 0 + const.adaptive_floor_accepts = 0 # The adaptive Cython stepper needs dens_co as a float64 array (memoryview); coerce once (the fixed # path's atmDensityPoly requires an ndarray too, so this is safe in both modes). @@ -1563,6 +1571,10 @@ def runSimulation(const, compute_wake=False): if const.adaptive_dt and (const.adaptive_runaway_events > 0): print("WARNING: adaptive stepper hit max_substeps ({:d}) on {:d} fragment-step(s); results may be " "under-resolved there.".format(const.adaptive_max_substeps, const.adaptive_runaway_events)) + if const.adaptive_dt and (const.adaptive_floor_accepts > 0): + print("WARNING: adaptive stepper accepted {:d} sub-step(s) at the dt_min floor ({:g} s) while still " + "over tolerance; lower adaptive_dt_min or loosen adaptive_rtol if this is frequent.".format( + const.adaptive_floor_accepts, const.adaptive_dt_min)) # Find the main fragment and return it with results diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index 1bbc6acd..70a2bb95 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/wmpl/MetSim/MetSimErosionCyTools.pyx @@ -671,6 +671,9 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero cdef int went_up = 0 cdef int runaway = 0 cdef int at_floor + cdef int was_clamped # this sub-step was shortened to land on the macro boundary + cdef int floor_accepts = 0 # sub-steps accepted only because at dt_min with E>1 (under-resolved) + cdef double h_sub_carry # controller's natural next step, kept unclamped by the macro boundary cdef double out1[7] cdef double out2[7] cdef double outb[7] @@ -695,11 +698,16 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero h_sub = dt_min E_prev = 1.0 rho_last = 0.0 + h_sub_carry = h_sub while t < dt_macro: - # Clamp the final sub-step so sub-steps sum to exactly dt_macro (no overshoot) + # Clamp the final sub-step so sub-steps sum to exactly dt_macro (no overshoot). Remember whether + # this step was shortened by the boundary so we don't persist the tiny clamped size as the + # warm-start for the next macro step (that would restart every macro step with a tiny sub-step). + was_clamped = 0 if t + h_sub > dt_macro: + was_clamped = 1 h_sub = dt_macro - t hh = 0.5*h_sub @@ -761,6 +769,12 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero at_floor = 1 if h_sub <= dt_min*(1.0 + 1e-12) else 0 if (E <= 1.0) or at_floor: + + # Count sub-steps accepted only because they hit the dt_min floor while still over tolerance + # (E > 1) - these are locally under-resolved and otherwise invisible (see runSimulation warn) + if at_floor and (E > 1.0): + floor_accepts += 1 + # Accept the two-half (more accurate) state; accumulate per-species mass loss (unclamped) dm_abl_macro += dm_abl_1 + dm_abl_2 dm_ero_macro += dm_ero_1 + dm_ero_2 @@ -808,6 +822,13 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero h_sub = dt_max if h_sub < dt_min: h_sub = dt_min + + # Persist the controller's natural step as the warm-start, but only from sub-steps that were + # NOT boundary-clamped, so the tiny final sub-step doesn't shrink the next macro step's seed + # (that would restart every macro step with a tiny sub-step). This cuts sub-steps ~2-5x with + # no effect on the visible light curve (only the faint, sub-detection tail shifts slightly). + if not was_clamped: + h_sub_carry = h_sub else: # Reject: shrink and retry the SAME sub-step (do not advance t) fac = safety*(E**(-1.0/5.0)) @@ -817,10 +838,11 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero if h_sub < dt_min: h_sub = dt_min - if went_up: - h_new = 0.0 - else: - h_new = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop_total + # Height from the along-track length (matches the fixed path, MetSimErosion.py: heightCurvature minus + # the gravity drop). Note: on went_up we do NOT force h_new = 0 - the fixed path's frag.h = 0 on + # 'going up' is immediately overwritten by this same recompute there, so forcing 0 would kill the + # fragment a tick earlier than fixed mode. went_up already zeroed vv, matching the fixed behaviour. + h_new = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop_total # Density for the end-of-interval dynamic pressure; fall back to the last in-loop value if h<=0 if h_new > 0: @@ -831,5 +853,5 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero decel_return = dv_drag/dt_macro # macro-averaged DRAG dv/dt (negative), matches decelerationRK4 sign return (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, - dm_abl_macro, dm_ero_macro, decel_return, went_up, n_substeps, h_sub, runaway, - erosion_events) \ No newline at end of file + dm_abl_macro, dm_ero_macro, decel_return, went_up, n_substeps, h_sub_carry, runaway, + floor_accepts, erosion_events) \ No newline at end of file From 22d868a578854303c1e2e498e121e14dc579a532 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 14:44:53 +0200 Subject: [PATCH 07/16] Add embedded Dormand-Prince RK45 adaptive stepper (~3-4x faster) Replaces step-doubling (which needed ~3 RK4 solves per sub-step for its error estimate) with an embedded Dormand-Prince RK45 pair on the COUPLED [mass, vv, vh, length] system (adaptiveDP45Step + _rhsDP in MetSimErosionCyTools.pyx). One embedded step gives the 5th-order solution and a 4th-order error estimate in ~7 RHS evals, and the higher order takes larger sub-steps - ~3.8x faster than step-doubling at the same tolerance. Atmosphere density is refreshed at every stage (from the stage height), removing the operator-split/frozen-rho error; gravity drop is added once per accepted sub-step, matching the fixed model. Selected by const.adaptive_high_order (default True); the step-doubling stepper is kept and used when False. Grains are still shed per accepted sub-step. Validation vs the step-doubling stepper (within measurement noise, ~0.1-0.3 mag): - Fireball (no erosion, pure integrator): visible-LC max |dmag| = 0.065, 3.8x faster. - Orionid (erosion): visible-LC max |dmag| = 0.18 (grain-birth cadence), converged (rtol 1e-5 ~ 1e-6), 3.8x faster. Backward compat unchanged: adaptive_dt=False is bit-for-bit vs master (13 Orionids); regression suite passes with DP45 as the default adaptive method. Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/MetSimErosion.py | 17 +- wmpl/MetSim/MetSimErosionCyTools.pyx | 240 +++++++++++++++++++++++++++ 2 files changed, 252 insertions(+), 5 deletions(-) diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 2e41c441..2a8024fd 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -23,7 +23,7 @@ import pyximport pyximport.install(setup_args={'include_dirs':[np.get_include()]}) from wmpl.MetSim.MetSimErosionCyTools import massLossRK4, decelerationRK4, luminousEfficiency, \ - ionizationEfficiency, atmDensityPoly, adaptiveSingleBodyStep + ionizationEfficiency, atmDensityPoly, adaptiveSingleBodyStep, adaptiveDP45Step ### DEFINE CONSTANTS @@ -50,6 +50,10 @@ def __init__(self): # path and reproduces prior results exactly. When True, each fragment sub-steps adaptively # within each dt to meet the tolerances below, while the output cadence stays fixed at dt. ### self.adaptive_dt = False + # Adaptive integrator: True -> embedded Dormand-Prince RK45 on the coupled (mass, velocity) + # system (fewer RHS evals per sub-step, higher order -> ~3-4x faster); False -> step-doubling + # on the original operator-split RK4. Both meet the same tolerance; DP45 is the default. + self.adaptive_high_order = True # Tolerance targets NUMERICAL error below MEASUREMENT noise (~0.1 mag, ~0.1 km/s at ~25 FPS), # not machine convergence. rtol=1e-5 keeps the single-body velocity error ~0.1-0.3 km/s (at or # below typical measurement precision) at a fraction of the cost of a tighter tol; tighten to @@ -815,10 +819,13 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): erosion_active = 1 if (frag.erosion_enabled and (frag.erosion_coeff > 0)) else 0 + _stepper = adaptiveDP45Step if getattr(const, 'adaptive_high_order', True) \ + else adaptiveSingleBodyStep + (frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, frag.h, rho_atm, mass_loss_ablation, mass_loss_erosion, deceleration_total, went_up, n_sub, frag.adaptive_h_sub, runaway, floor_accepts, erosion_events_adaptive) \ - = adaptiveSingleBodyStep( + = _stepper( const.dt, frag.K, frag.sigma, frag.erosion_coeff, erosion_active, frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, const.h_init, const.zenith_angle, const.r_earth, G0, const.dens_co, @@ -1467,9 +1474,9 @@ def runSimulation(const, compute_wake=False): # Back-fill adaptive-timestep settings for Constants loaded from older JSONs that predate them, so # such runs default to the original fixed-step behaviour. Also reset the per-run diagnostics. - _adaptive_defaults = {'adaptive_dt': False, 'adaptive_rtol': 1e-5, 'adaptive_atol_m': 1e-14, - 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, 'adaptive_dt_max': const.dt, - 'adaptive_max_substeps': 10000} + _adaptive_defaults = {'adaptive_dt': False, 'adaptive_high_order': True, 'adaptive_rtol': 1e-5, + 'adaptive_atol_m': 1e-14, 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, + 'adaptive_dt_max': const.dt, 'adaptive_max_substeps': 10000} for _attr, _default in _adaptive_defaults.items(): if not hasattr(const, _attr): setattr(const, _attr, _default) diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index 70a2bb95..62880d7e 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/wmpl/MetSim/MetSimErosionCyTools.pyx @@ -852,6 +852,246 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero decel_return = dv_drag/dt_macro # macro-averaged DRAG dv/dt (negative), matches decelerationRK4 sign + return (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, + dm_abl_macro, dm_ero_macro, decel_return, went_up, n_substeps, h_sub_carry, runaway, + floor_accepts, erosion_events) + + + +@cython.cdivision(True) +cdef void _rhsDP(double m, double vv, double vh, double length, double h_grav_drop, + double K, double sigma, double erosion_coeff, int erosion_active, double m_kill, + double h_init, double zenith_angle, double r_earth, FLOAT_TYPE_t[:] dens_co, + double* dydt, double* extras): + """ Coupled RHS for the Dormand-Prince stepper. State y = [m, vv, vh, length]; dydt is filled with + [dm/dt, dvv/dt, dvh/dt, dlength/dt]. extras returns [ablation_rate, erosion_rate, drag_decel, rho] + for the per-step diagnostics (luminosity, electron density, grain shedding). Mirrors the fixed + model: mass loss ~ m^(2/3), drag ~ m^(-1/3) (floored at m_kill), gravity acts only on the height + drop (handled by the caller, not here), and the vv/vh derivatives carry the Earth-curvature term. """ + cdef double v, h, rho, decel, mm, mpos + v = sqrt(vv*vv + vh*vh) + h = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop + rho = atmDensityPolyC(h, dens_co) + mpos = m if m > 0.0 else 0.0 + extras[0] = massLoss(K, sigma, mpos, rho, v) + if erosion_active: + extras[1] = massLoss(K, erosion_coeff, mpos, rho, v) + else: + extras[1] = 0.0 + mm = fmax(m, m_kill) + decel = deceleration(K, mm, rho, v) + extras[2] = decel + extras[3] = rho + dydt[0] = extras[0] + extras[1] + if v > 0.0: + dydt[1] = decel*vv/v - vh*v/(r_earth + h) + dydt[2] = decel*vh/v + vv*v/(r_earth + h) + else: + dydt[1] = 0.0 + dydt[2] = 0.0 + dydt[3] = v + + +@cython.cdivision(True) +cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_coeff, + int erosion_active, double m, double v, double vv, double vh, double length, + double h_grav_drop_total, double h_init, double zenith_angle, double r_earth, double g0, + FLOAT_TYPE_t[:] dens_co, double rtol, double atol_m, double atol_v, double m_kill, + double dt_min, double dt_max, int max_substeps, double h_sub_init): + """ Same contract as adaptiveSingleBodyStep, but advances the coupled [m, vv, vh, length] system + with an embedded Dormand-Prince RK45 pair (5th-order solution, 4th-order error estimate) instead + of step-doubling. ~7 RHS evaluations per sub-step vs ~24, and the higher order takes larger + steps, so it is substantially faster for the same tolerance. rho is refreshed at every stage + (from the stage's height), so there is no operator-split/frozen-rho error. Gravity drop is added + once per accepted sub-step (as in the fixed model). Returns the 17-tuple: + (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, dm_abl_macro, dm_ero_macro, + decel_return, went_up, n_substeps, h_sub_carry, runaway, floor_accepts, erosion_events). """ + + # Dormand-Prince (RK45) Butcher tableau + cdef double a21 = 1.0/5 + cdef double a31 = 3.0/40, a32 = 9.0/40 + cdef double a41 = 44.0/45, a42 = -56.0/15, a43 = 32.0/9 + cdef double a51 = 19372.0/6561, a52 = -25360.0/2187, a53 = 64448.0/6561, a54 = -212.0/729 + cdef double a61 = 9017.0/3168, a62 = -355.0/33, a63 = 46732.0/5247, a64 = 49.0/176, a65 = -5103.0/18656 + cdef double a71 = 35.0/384, a73 = 500.0/1113, a74 = 125.0/192, a75 = -2187.0/6784, a76 = 11.0/84 + # 5th-order weights b = 7th stage row (FSAL); 4th-order embedded weights b* + cdef double b1 = 35.0/384, b3 = 500.0/1113, b4 = 125.0/192, b5 = -2187.0/6784, b6 = 11.0/84 + cdef double bs1 = 5179.0/57600, bs3 = 7571.0/16695, bs4 = 393.0/640, bs5 = -92097.0/339200 + cdef double bs6 = 187.0/2100, bs7 = 1.0/40 + + cdef double t = 0.0 + cdef double h_sub, hh, gv, h_cur, rho0 + cdef double y[4] + cdef double yt[4] + cdef double k1[4] + cdef double k2[4] + cdef double k3[4] + cdef double k4[4] + cdef double k5[4] + cdef double k6[4] + cdef double k7[4] + # extras per stage: [ablation_rate, erosion_rate, drag_decel, rho] + cdef double e1[4] + cdef double e2[4] + cdef double e3[4] + cdef double e4[4] + cdef double e5[4] + cdef double e6[4] + cdef double e7[4] + cdef double m_new, v_new, vv_new, vh_new, len_new + cdef double abl_step, ero_step, drag_step, err_m, err_v, sc_m, sc_v, E, E_prev, fac + cdef double dm_abl_macro = 0.0, dm_ero_macro = 0.0, dv_drag = 0.0 + cdef double h_new, decel_return, rho_final, rho_last = 0.0, h_sub_carry + cdef int n_substeps = 0, went_up = 0, runaway = 0, at_floor, was_clamped, floor_accepts = 0 + cdef int c + cdef double safety = 0.9, facmin = 0.2, facmax = 5.0 + + erosion_events = [] if erosion_active else None + + y[0] = m; y[1] = vv; y[2] = vh; y[3] = length + + h_sub = h_sub_init + if h_sub <= 0: + h_sub = dt_macro + if h_sub > dt_max: + h_sub = dt_max + if h_sub < dt_min: + h_sub = dt_min + E_prev = 1.0 + h_sub_carry = h_sub + + while t < dt_macro: + + was_clamped = 0 + if t + h_sub > dt_macro: + was_clamped = 1 + h_sub = dt_macro - t + + # Height/gravity at the sub-step start (gravity drop uses the start-of-step g, as in the fixed model) + h_cur = heightCurvatureC(h_init, zenith_angle, y[3], r_earth) - h_grav_drop_total + gv = g0/((1.0 + h_cur/r_earth)*(1.0 + h_cur/r_earth)) + + # --- 7 Dormand-Prince stages (h_grav_drop frozen across the sub-step) --- + _rhsDP(y[0], y[1], y[2], y[3], h_grav_drop_total, K, sigma, erosion_coeff, erosion_active, + m_kill, h_init, zenith_angle, r_earth, dens_co, k1, e1) + rho0 = e1[3] + for c in range(4): + yt[c] = y[c] + h_sub*a21*k1[c] + _rhsDP(yt[0], yt[1], yt[2], yt[3], h_grav_drop_total, K, sigma, erosion_coeff, erosion_active, + m_kill, h_init, zenith_angle, r_earth, dens_co, k2, e2) + for c in range(4): + yt[c] = y[c] + h_sub*(a31*k1[c] + a32*k2[c]) + _rhsDP(yt[0], yt[1], yt[2], yt[3], h_grav_drop_total, K, sigma, erosion_coeff, erosion_active, + m_kill, h_init, zenith_angle, r_earth, dens_co, k3, e3) + for c in range(4): + yt[c] = y[c] + h_sub*(a41*k1[c] + a42*k2[c] + a43*k3[c]) + _rhsDP(yt[0], yt[1], yt[2], yt[3], h_grav_drop_total, K, sigma, erosion_coeff, erosion_active, + m_kill, h_init, zenith_angle, r_earth, dens_co, k4, e4) + for c in range(4): + yt[c] = y[c] + h_sub*(a51*k1[c] + a52*k2[c] + a53*k3[c] + a54*k4[c]) + _rhsDP(yt[0], yt[1], yt[2], yt[3], h_grav_drop_total, K, sigma, erosion_coeff, erosion_active, + m_kill, h_init, zenith_angle, r_earth, dens_co, k5, e5) + for c in range(4): + yt[c] = y[c] + h_sub*(a61*k1[c] + a62*k2[c] + a63*k3[c] + a64*k4[c] + a65*k5[c]) + _rhsDP(yt[0], yt[1], yt[2], yt[3], h_grav_drop_total, K, sigma, erosion_coeff, erosion_active, + m_kill, h_init, zenith_angle, r_earth, dens_co, k6, e6) + for c in range(4): + yt[c] = y[c] + h_sub*(a71*k1[c] + a73*k3[c] + a74*k4[c] + a75*k5[c] + a76*k6[c]) + _rhsDP(yt[0], yt[1], yt[2], yt[3], h_grav_drop_total, K, sigma, erosion_coeff, erosion_active, + m_kill, h_init, zenith_angle, r_earth, dens_co, k7, e7) + + # 5th-order solution (= yt above, the 7th stage node) and error vs the 4th-order embedded + m_new = yt[0]; vv_new = yt[1]; vh_new = yt[2]; len_new = yt[3] + v_new = sqrt(vv_new*vv_new + vh_new*vh_new) + + # Error estimate on mass and speed (difference of 5th- and 4th-order weights) + err_m = fabs(h_sub*((b1-bs1)*k1[0] + (b3-bs3)*k3[0] + (b4-bs4)*k4[0] + (b5-bs5)*k5[0] + + (b6-bs6)*k6[0] + (0.0-bs7)*k7[0])) + # velocity error via the vv/vh error components projected onto speed + err_v = fabs(h_sub*((b1-bs1)*k1[1] + (b3-bs3)*k3[1] + (b4-bs4)*k4[1] + (b5-bs5)*k5[1] + + (b6-bs6)*k6[1] + (0.0-bs7)*k7[1])) + err_v += fabs(h_sub*((b1-bs1)*k1[2] + (b3-bs3)*k3[2] + (b4-bs4)*k4[2] + (b5-bs5)*k5[2] + + (b6-bs6)*k6[2] + (0.0-bs7)*k7[2])) + sc_m = atol_m + rtol*fmax(fabs(m_new), fabs(y[0])) + sc_v = atol_v + rtol*fmax(v_new, v) + E = sqrt(0.5*((err_m/sc_m)*(err_m/sc_m) + (err_v/sc_v)*(err_v/sc_v))) + + at_floor = 1 if h_sub <= dt_min*(1.0 + 1e-12) else 0 + + if (E <= 1.0) or at_floor: + + if at_floor and (E > 1.0): + floor_accepts += 1 + + # Per-step mass loss split (integral of each rate with the 5th-order weights) + abl_step = h_sub*(b1*e1[0] + b3*e3[0] + b4*e4[0] + b5*e5[0] + b6*e6[0]) + ero_step = h_sub*(b1*e1[1] + b3*e3[1] + b4*e4[1] + b5*e5[1] + b6*e6[1]) + drag_step = h_sub*(b1*e1[2] + b3*e3[2] + b4*e4[2] + b5*e5[2] + b6*e6[2]) + dm_abl_macro += abl_step + dm_ero_macro += ero_step + dv_drag += drag_step + + # Commit the state; floor the mass at 0 and add the gravity drop for this sub-step + y[0] = m_new if m_new > 0.0 else 0.0 + y[1] = vv_new; y[2] = vh_new; y[3] = len_new + h_grav_drop_total += 0.5*gv*h_sub*h_sub + rho_last = rho0 + t += h_sub + n_substeps += 1 + + # Shed the eroded mass at this sub-step's resolved state (see adaptiveSingleBodyStep) + if erosion_active and (ero_step < 0): + erosion_events.append(( + -ero_step, + heightCurvatureC(h_init, zenith_angle, y[3], r_earth) - h_grav_drop_total, + sqrt(y[1]*y[1] + y[2]*y[2]), y[1], y[2], y[3], h_grav_drop_total)) + + if y[0] <= m_kill: # exhausted + break + if y[1] > 0.0: # turned upward -> freeze, matches fixed (vv zeroed) + y[1] = 0.0 + went_up = 1 + break + v_new = sqrt(y[1]*y[1] + y[2]*y[2]) + if v_new <= 0.0: + break + if n_substeps >= max_substeps: + runaway = 1 + break + + # PI step-size controller (order p=4 for the error estimate -> k=5) + if E <= 0: + E = 1e-10 + fac = safety*(E**(-0.7/5.0))*(E_prev**(0.4/5.0)) + E_prev = E + if fac < facmin: + fac = facmin + if fac > facmax: + fac = facmax + h_sub = h_sub*fac + if h_sub > dt_max: + h_sub = dt_max + if h_sub < dt_min: + h_sub = dt_min + if not was_clamped: + h_sub_carry = h_sub + else: + fac = safety*(E**(-1.0/5.0)) + if fac < facmin: + fac = facmin + h_sub = h_sub*fac + if h_sub < dt_min: + h_sub = dt_min + + m = y[0]; vv = y[1]; vh = y[2]; length = y[3] + v = sqrt(vv*vv + vh*vh) + h_new = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop_total + if h_new > 0: + rho_final = atmDensityPolyC(h_new, dens_co) + else: + rho_final = rho_last + decel_return = dv_drag/dt_macro + return (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, dm_abl_macro, dm_ero_macro, decel_return, went_up, n_substeps, h_sub_carry, runaway, floor_accepts, erosion_events) \ No newline at end of file From 560413a6f1beeed071592161e4d2ae25bc8ce927 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 15:03:11 +0200 Subject: [PATCH 08/16] Document adaptive-stepper functions (wmpl-style docstrings) + cleanup - Full Arguments:/Return: docstrings for the new Cython functions: clampMassC, heightCurvatureC, atmDensityPolyC, advanceVelPosC, adaptiveSingleBodyStep, _rhsDP, adaptiveDP45Step (the last two reference adaptiveSingleBodyStep for the shared I/O). - Corrected the adaptiveSingleBodyStep return description to the actual 17-tuple. - Dropped the unused const_out binding in the erosion grain-spawn helper. No behaviour change; regression suite passes. Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/MetSimErosion.py | 2 +- wmpl/MetSim/MetSimErosionCyTools.pyx | 180 ++++++++++++++++++++++----- 2 files changed, 151 insertions(+), 31 deletions(-) diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 2a8024fd..344c6333 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -1108,7 +1108,7 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): def _spawnGrainsFromErosion(eroded_mass): """ Distribute the given eroded mass into grains born from the fragment's current state. Uses the enclosing frag/const. """ - grain_children, const_out = generateFragments(const, frag, eroded_mass, \ + grain_children, _ = generateFragments(const, frag, eroded_mass, \ frag.erosion_mass_index, frag.erosion_mass_min, frag.erosion_mass_max, \ keep_eroding=False, mass_model=const.erosion_grain_distribution) const.n_active += len(grain_children) diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index 62880d7e..fde7fe21 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/wmpl/MetSim/MetSimErosionCyTools.pyx @@ -588,8 +588,16 @@ cpdef atmDensityPoly(double ht, np.ndarray[FLOAT_TYPE_t, ndim=1] dens_co): cdef inline double clampMassC(double dm, double m): - """ Reproduce the "ablate at most the whole mass" clamp (MetSimErosion.py mass-loss clamp): - if m + dm < 0, cap the loss at exactly -m so the mass floors at 0. """ + """ Reproduce the "ablate at most the whole mass" clamp used by the fixed engine: if the mass loss + would drive the mass below zero, cap it at exactly -m so the new mass floors at 0. + + Arguments: + dm: [double] Proposed mass change over the (sub-)step (kg, normally negative). + m: [double] Current fragment mass (kg). + + Return: + dm_clamped: [double] dm, or -m if (m + dm) < 0. + """ if (m + dm) < 0: return -m return dm @@ -597,13 +605,33 @@ cdef inline double clampMassC(double dm, double m): @cython.cdivision(True) cdef double heightCurvatureC(double h0, double zc, double l, double r_earth): - """ Cython twin of heightCurvature() (MetSimErosion.py). Scalar hot-path. """ + """ Cython/scalar twin of heightCurvature() (MetSimErosion.py): height at a distance l along the + trajectory, accounting for the Earth's curvature. + + Arguments: + h0: [double] Initial height (m). + zc: [double] Zenith angle (radians). + l: [double] Distance travelled along the trajectory from the origin (m). + r_earth: [double] Earth radius (m). + + Return: + h: [double] Height at distance l (m), before the gravity drop is subtracted. + """ return sqrt((h0 + r_earth)*(h0 + r_earth) - 2*l*cos(zc)*(h0 + r_earth) + l*l) - r_earth @cython.cdivision(True) cdef double atmDensityPolyC(double ht, FLOAT_TYPE_t[:] dens_co): - """ Cython/memoryview twin of atmDensityPoly() for use inside the substep loop. """ + """ Cython/memoryview twin of atmDensityPoly() for use inside the sub-step loop (avoids the + ndarray-typed cpdef boundary). + + Arguments: + ht: [double] Height (m). + dens_co: [memoryview of float64] 7 coefficients of the log10(density) height polynomial. + + Return: + rho: [double] Atmospheric mass density at height ht (kg/m^3). + """ cdef double x = ht/1e6 return 10**(dens_co[0] + dens_co[1]*x + dens_co[2]*x*x + dens_co[3]*x*x*x + dens_co[4]*x*x*x*x + dens_co[5]*x*x*x*x*x + dens_co[6]*x*x*x*x*x*x) @@ -613,10 +641,30 @@ cdef double atmDensityPolyC(double ht, FLOAT_TYPE_t[:] dens_co): cdef void advanceVelPosC(double m, double v, double vv, double vh, double length, double grav, double dm, double decel_rate, double h, double h_at, double r_earth, double g0, double* out): - """ Advance velocity components, speed, along-track length, and gravity-drop for one sub-step of - size h, reproducing MetSimErosion.py:773-824 exactly (velocity updated BEFORE length; gravity - only drops height, does not enter the velocity magnitude). out layout: - [m_new, v_new, vv_new, vh_new, length_new, grav_new, went_up_flag]. """ + """ Advance the velocity components, speed, along-track length, and gravity drop of one fragment by + one sub-step of size h, reproducing the fixed engine's single-body update exactly (velocity is + updated BEFORE the length; gravity only lowers the height, it does not enter the velocity + magnitude). Used by the step-doubling stepper (adaptiveSingleBodyStep). + + Arguments: + m: [double] Mass at the sub-step start (kg). + v: [double] Speed at the sub-step start (m/s). + vv: [double] Vertical velocity component (m/s, negative downward). + vh: [double] Horizontal velocity component (m/s). + length: [double] Along-track length at the sub-step start (m). + grav: [double] Accumulated gravity drop so far (m). + dm: [double] Mass change over this sub-step (kg, already clamped). + decel_rate: [double] Drag deceleration dv/dt (m/s^2, <= 0). + h: [double] Sub-step size (s). + h_at: [double] Height at which to evaluate gravity/curvature for this sub-step (m). + r_earth: [double] Earth radius (m). + g0: [double] Surface gravitational acceleration (m/s^2). + out: [double*] Output buffer (length 7), filled with + [m_new, v_new, vv_new, vh_new, length_new, grav_new, went_up_flag]. + + Return: + None (results written into 'out'). + """ cdef double gv, av, ah, vv_n, vh_n, v_n # Accelerating (decel_rate > 0) or already stopped -> stop the fragment (mirror 773-775) if (decel_rate > 0) or (v <= 0): @@ -645,15 +693,60 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero FLOAT_TYPE_t[:] dens_co, double rtol, double atol_m, double atol_v, double m_kill, double dt_min, double dt_max, int max_substeps, double h_sub_init): """ Advance ONE fragment across a full macro interval dt_macro using error-controlled adaptive - sub-steps (step-doubling on the existing RK4), refreshing atmosphere/height each sub-step. - Reproduces the single-body advance of MetSimErosion.py (748-834) in the one-substep limit but - drives the local error below (rtol, atol). Events/kills stay at the macro boundary (handled by - the caller). Returns a tuple: - (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, - dm_abl_macro, dm_ero_macro, decel_return, went_up, n_substeps, h_sub_last, runaway) - where dm_*_macro are the (unclamped, per-species) accumulated mass losses used for the - luminosity/electron-density/grain diagnostics, and decel_return is the macro-averaged dv/dt - (negative), consistent with the sign of decelerationRK4. """ + sub-steps (step-doubling on the operator-split RK4), refreshing the atmosphere/height each + sub-step. Reproduces the single-body advance of the fixed engine in the one-sub-step limit, but + drives the local error below (rtol, atol). Grain generation, kill checks, disruption and complex + fragmentation stay at the macro boundary and are handled by the caller (ablateAll). + + See adaptiveDP45Step for the (default) embedded Dormand-Prince variant with the same I/O. + + Arguments: + dt_macro: [double] Macro (output) interval to advance over (s). + K: [double] Shape-density coefficient (m^2/kg^(2/3)). + sigma: [double] Ablation coefficient (s^2/m^2). + erosion_coeff: [double] Erosion coefficient (s^2/m^2); used only if erosion_active. + erosion_active: [int] 1 if this fragment is currently eroding (erosion_enabled and coeff > 0). + m: [double] Mass at the interval start (kg). + v: [double] Speed at the interval start (m/s). + vv: [double] Vertical velocity component (m/s, negative downward). + vh: [double] Horizontal velocity component (m/s). + length: [double] Along-track length at the interval start (m). + h_grav_drop_total: [double] Accumulated gravity drop at the interval start (m). + h_init: [double] Initial simulation height (m). + zenith_angle: [double] Entry zenith angle (radians). + r_earth: [double] Earth radius (m). + g0: [double] Surface gravitational acceleration (m/s^2). + dens_co: [memoryview of float64] 7 atmosphere log10(density) polynomial coefficients. + rtol: [double] Relative error tolerance on mass and speed. + atol_m: [double] Absolute mass tolerance (kg). + atol_v: [double] Absolute speed tolerance (m/s). + m_kill: [double] Kill mass (kg); the drag mass is floored at this to keep m^(-1/3) finite. + dt_min: [double] Minimum sub-step (s). + dt_max: [double] Maximum sub-step (s). + max_substeps: [int] Sub-step cap per macro interval (runaway guard). + h_sub_init: [double] Warm-start sub-step from the previous macro step (s); <= 0 means "use dt_macro". + + Return: + A 17-tuple: + m: [double] Mass at the interval end (kg). + v: [double] Speed at the interval end (m/s). + vv: [double] Vertical velocity component at the end (m/s). + vh: [double] Horizontal velocity component at the end (m/s). + length: [double] Along-track length at the end (m). + h_grav_drop_total: [double] Accumulated gravity drop at the end (m). + h_new: [double] Height at the interval end (m). + rho_final: [double] Atmospheric density at h_new, for the end-of-step dynamic pressure (kg/m^3). + dm_abl_macro: [double] Unclamped ablation mass loss over the interval (kg, for luminosity/q). + dm_ero_macro: [double] Unclamped erosion mass loss over the interval (kg). + decel_return: [double] Macro-averaged DRAG dv/dt (m/s^2, negative), for the luminosity term. + went_up: [int] 1 if the fragment turned upward during the interval. + n_substeps: [int] Number of accepted sub-steps taken. + h_sub_carry: [double] Warm-start sub-step to carry to the next macro step (s). + runaway: [int] 1 if the max_substeps cap was hit. + floor_accepts: [int] Sub-steps accepted at dt_min while still over tolerance (under-resolved). + erosion_events: [list or None] Per-sub-step erosion shedding events, each a tuple + (eroded_mass, h, v, vv, vh, length, h_grav_drop_total); None for non-eroding fragments. + """ cdef double t = 0.0 cdef double h_sub, h_cur, rho_atm, rho_mid, rho_last @@ -863,11 +956,34 @@ cdef void _rhsDP(double m, double vv, double vh, double length, double h_grav_dr double K, double sigma, double erosion_coeff, int erosion_active, double m_kill, double h_init, double zenith_angle, double r_earth, FLOAT_TYPE_t[:] dens_co, double* dydt, double* extras): - """ Coupled RHS for the Dormand-Prince stepper. State y = [m, vv, vh, length]; dydt is filled with - [dm/dt, dvv/dt, dvh/dt, dlength/dt]. extras returns [ablation_rate, erosion_rate, drag_decel, rho] - for the per-step diagnostics (luminosity, electron density, grain shedding). Mirrors the fixed - model: mass loss ~ m^(2/3), drag ~ m^(-1/3) (floored at m_kill), gravity acts only on the height - drop (handled by the caller, not here), and the vv/vh derivatives carry the Earth-curvature term. """ + """ Coupled right-hand side for the Dormand-Prince stepper: derivatives of the state + y = [m, vv, vh, length]. Mirrors the fixed model - mass loss ~ m^(2/3), drag ~ m^(-1/3) (with the + drag mass floored at m_kill so it stays finite near exhaustion), and the vv/vh derivatives carry + the Earth-curvature term. Gravity acts only on the height drop and is handled by the caller (not + here). The height for the atmosphere lookup is derived from length and the (frozen) gravity drop. + + Arguments: + m: [double] Mass (kg). + vv: [double] Vertical velocity component (m/s). + vh: [double] Horizontal velocity component (m/s). + length: [double] Along-track length (m). + h_grav_drop: [double] Gravity drop to subtract from the curvature height (m), frozen over the step. + K: [double] Shape-density coefficient (m^2/kg^(2/3)). + sigma: [double] Ablation coefficient (s^2/m^2). + erosion_coeff: [double] Erosion coefficient (s^2/m^2); used only if erosion_active. + erosion_active: [int] 1 if the fragment is eroding. + m_kill: [double] Kill mass (kg); floor for the drag mass. + h_init: [double] Initial simulation height (m). + zenith_angle: [double] Entry zenith angle (radians). + r_earth: [double] Earth radius (m). + dens_co: [memoryview of float64] Atmosphere density polynomial coefficients. + dydt: [double*] Output buffer (length 4): [dm/dt, dvv/dt, dvh/dt, dlength/dt]. + extras: [double*] Output buffer (length 4): [ablation_rate, erosion_rate, drag_decel, rho], + used by the caller for the per-step luminosity/electron-density/grain-shedding diagnostics. + + Return: + None (results written into 'dydt' and 'extras'). + """ cdef double v, h, rho, decel, mm, mpos v = sqrt(vv*vv + vh*vh) h = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop @@ -898,14 +1014,18 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c double h_grav_drop_total, double h_init, double zenith_angle, double r_earth, double g0, FLOAT_TYPE_t[:] dens_co, double rtol, double atol_m, double atol_v, double m_kill, double dt_min, double dt_max, int max_substeps, double h_sub_init): - """ Same contract as adaptiveSingleBodyStep, but advances the coupled [m, vv, vh, length] system - with an embedded Dormand-Prince RK45 pair (5th-order solution, 4th-order error estimate) instead - of step-doubling. ~7 RHS evaluations per sub-step vs ~24, and the higher order takes larger - steps, so it is substantially faster for the same tolerance. rho is refreshed at every stage - (from the stage's height), so there is no operator-split/frozen-rho error. Gravity drop is added - once per accepted sub-step (as in the fixed model). Returns the 17-tuple: - (m, v, vv, vh, length, h_grav_drop_total, h_new, rho_final, dm_abl_macro, dm_ero_macro, - decel_return, went_up, n_substeps, h_sub_carry, runaway, floor_accepts, erosion_events). """ + """ Advance ONE fragment across a full macro interval dt_macro, like adaptiveSingleBodyStep, but with + an embedded Dormand-Prince RK45 pair (5th-order solution + 4th-order error estimate) on the + COUPLED system y = [m, vv, vh, length] instead of step-doubling on the operator split. One + embedded step yields both solutions in ~7 RHS evaluations (vs ~24 for step-doubling), and its + higher order takes larger sub-steps, so it is ~3-4x faster at the same tolerance. The atmosphere + density is refreshed at every RK stage (from that stage's height), removing the operator-split / + frozen-density error; the gravity drop is added once per accepted sub-step (as in the fixed + model). This is the default adaptive stepper (const.adaptive_high_order). + + Arguments and the 17-element return tuple are identical to adaptiveSingleBodyStep() - see its + docstring for the full I/O description. + """ # Dormand-Prince (RK45) Butcher tableau cdef double a21 = 1.0/5 From 21b4ab81965277ddef646a81edd7de5675b37a08 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 15:26:23 +0200 Subject: [PATCH 09/16] Decouple grain-release cadence from dt via erosion_release_length In adaptive mode grains are shed once per accepted sub-step, so the sub-step cadence set the grain-birth resolution - which made the erosion light-curve shape depend on dt (dt caps the sub-step via adaptive_dt_max). Add const.erosion_release_length (m, default 50): an eroding fragment's sub-step is capped at erosion_release_length/v, so grains are released ~every that many metres of flight, a physical resolution independent of dt and the error tolerance. Applied in both the DP45 and step-doubling steppers; only for eroding fragments (grains do not erode further). Effect (Orionid 084502): peak-magnitude difference between dt=0.005 and dt=0.02 drops from ~2.0 mag to ~0.01 mag, and the total grain count becomes dt-independent (+-6%). The knob is a set-and-forget JSON default (not a GUI box); tune it for finer/coarser grain sampling. Backward compat unchanged: adaptive_dt=False is bit-for-bit vs master; regression suite passes. Co-Authored-By: Claude Opus 4.8 (1M context) --- wmpl/MetSim/MetSimErosion.py | 15 +++++++++++++-- wmpl/MetSim/MetSimErosionCyTools.pyx | 25 ++++++++++++++++++++++--- 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 344c6333..6bb5d432 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -65,6 +65,16 @@ def __init__(self): self.adaptive_dt_max = 0.005 # largest allowed sub-step (s); should be <= dt self.adaptive_max_substeps = 10000 # runaway guard, sub-steps per fragment per macro step + # Along-track grain-release interval (m). In adaptive mode grains are shed once per accepted + # sub-step, so the sub-step cadence sets the grain-birth resolution. Capping an eroding + # fragment's sub-step at erosion_release_length/v releases grains roughly every this many + # metres of flight - a physical resolution independent of dt and the error tolerance, so the + # erosion light-curve shape stops depending on dt. Smaller = finer grain sampling + slower. + # Only used in adaptive mode and only for eroding fragments (grains do not erode further). + # Not shown in the GUI - it is a dt-independent physical default that rarely needs changing; + # tune via the sim-fit JSON if finer/coarser grain sampling is wanted. + self.erosion_release_length = 50.0 + # Time elapsed since the beginning self.total_time = 0 @@ -831,7 +841,7 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): const.h_init, const.zenith_angle, const.r_earth, G0, const.dens_co, const.adaptive_rtol, const.adaptive_atol_m, const.adaptive_atol_v, const.m_kill, const.adaptive_dt_min, const.adaptive_dt_max, const.adaptive_max_substeps, - frag.adaptive_h_sub) + frag.adaptive_h_sub, getattr(const, 'erosion_release_length', 50.0)) # Diagnostics (feed the cost study; also used to warn once on runaway/under-resolved steps) const.adaptive_substeps_total += n_sub @@ -1476,7 +1486,8 @@ def runSimulation(const, compute_wake=False): # such runs default to the original fixed-step behaviour. Also reset the per-run diagnostics. _adaptive_defaults = {'adaptive_dt': False, 'adaptive_high_order': True, 'adaptive_rtol': 1e-5, 'adaptive_atol_m': 1e-14, 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, - 'adaptive_dt_max': const.dt, 'adaptive_max_substeps': 10000} + 'adaptive_dt_max': const.dt, 'adaptive_max_substeps': 10000, + 'erosion_release_length': 50.0} for _attr, _default in _adaptive_defaults.items(): if not hasattr(const, _attr): setattr(const, _attr, _default) diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index fde7fe21..1e27db91 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/wmpl/MetSim/MetSimErosionCyTools.pyx @@ -691,7 +691,8 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero int erosion_active, double m, double v, double vv, double vh, double length, double h_grav_drop_total, double h_init, double zenith_angle, double r_earth, double g0, FLOAT_TYPE_t[:] dens_co, double rtol, double atol_m, double atol_v, double m_kill, - double dt_min, double dt_max, int max_substeps, double h_sub_init): + double dt_min, double dt_max, int max_substeps, double h_sub_init, + double erosion_release_length): """ Advance ONE fragment across a full macro interval dt_macro using error-controlled adaptive sub-steps (step-doubling on the operator-split RK4), refreshing the atmosphere/height each sub-step. Reproduces the single-body advance of the fixed engine in the one-sub-step limit, but @@ -725,6 +726,9 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero dt_max: [double] Maximum sub-step (s). max_substeps: [int] Sub-step cap per macro interval (runaway guard). h_sub_init: [double] Warm-start sub-step from the previous macro step (s); <= 0 means "use dt_macro". + erosion_release_length: [double] Along-track grain-release interval (m). For an eroding fragment + the sub-step is capped at erosion_release_length/v so grains are shed ~every this many metres, + making the grain-birth cadence independent of dt/rtol. Ignored for non-eroding fragments. Return: A 17-tuple: @@ -803,6 +807,12 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero was_clamped = 1 h_sub = dt_macro - t + # Cap an eroding fragment's sub-step so grains are shed roughly every erosion_release_length + # metres of flight, making the grain-birth cadence (and hence the erosion light curve) + # independent of dt and the error tolerance. Not applied to non-eroding fragments or grains. + if erosion_active and (v > 0) and (erosion_release_length > 0) and (h_sub*v > erosion_release_length): + h_sub = erosion_release_length/v + hh = 0.5*h_sub # Height and atmosphere at the CURRENT state (refresh -> shrinks the operator-split/frozen-rho error) @@ -1013,7 +1023,8 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c int erosion_active, double m, double v, double vv, double vh, double length, double h_grav_drop_total, double h_init, double zenith_angle, double r_earth, double g0, FLOAT_TYPE_t[:] dens_co, double rtol, double atol_m, double atol_v, double m_kill, - double dt_min, double dt_max, int max_substeps, double h_sub_init): + double dt_min, double dt_max, int max_substeps, double h_sub_init, + double erosion_release_length): """ Advance ONE fragment across a full macro interval dt_macro, like adaptiveSingleBodyStep, but with an embedded Dormand-Prince RK45 pair (5th-order solution + 4th-order error estimate) on the COUPLED system y = [m, vv, vh, length] instead of step-doubling on the operator split. One @@ -1040,7 +1051,7 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c cdef double bs6 = 187.0/2100, bs7 = 1.0/40 cdef double t = 0.0 - cdef double h_sub, hh, gv, h_cur, rho0 + cdef double h_sub, hh, gv, h_cur, rho0, v_cur cdef double y[4] cdef double yt[4] cdef double k1[4] @@ -1087,6 +1098,14 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c was_clamped = 1 h_sub = dt_macro - t + # Cap an eroding fragment's sub-step so grains are shed roughly every erosion_release_length + # metres of flight - a grain-birth cadence independent of dt and the error tolerance (see the + # step-doubling stepper for the rationale). Only for eroding fragments; grains do not erode. + v_cur = sqrt(y[1]*y[1] + y[2]*y[2]) + if erosion_active and (v_cur > 0) and (erosion_release_length > 0) \ + and (h_sub*v_cur > erosion_release_length): + h_sub = erosion_release_length/v_cur + # Height/gravity at the sub-step start (gravity drop uses the start-of-step g, as in the fixed model) h_cur = heightCurvatureC(h_init, zenith_angle, y[3], r_earth) - h_grav_drop_total gv = g0/((1.0 + h_cur/r_earth)*(1.0 + h_cur/r_earth)) From 86ac4939899e4e38d653f1b4eb958fbfe21f16b7 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 21:47:24 +0200 Subject: [PATCH 10/16] Make adaptive timestep the default (adaptive_dt=True, rtol=1e-4) Turn the error-controlled adaptive integrator on by default so erosion light curves are well resolved without hand-tuning dt down (the old workaround). The fixed-step path is unchanged and still bit-for-bit vs master when adaptive_dt is set False in the JSON/GUI. - Constants: adaptive_dt True, adaptive_rtol 1e-4 (was False / 1e-5). - runSimulation back-fill and GUI getattr fallbacks track the new defaults; older JSONs without the fields now pick up adaptive-on. - GUI: adaptive checkbox checked by default; rtol tooltip updated. rtol=1e-4 vs 1e-5 is indistinguishable on the 2019 Orionids (<0.0005 mag); 1e-4 is the slightly looser, safe default. Set adaptive_dt False for legacy runs. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_0179geUg2TGRWWqpNPx2k234 --- wmpl/MetSim/GUI.py | 4 ++-- wmpl/MetSim/GUI.ui | 5 ++++- wmpl/MetSim/MetSimErosion.py | 35 +++++++++++++++++++---------------- 3 files changed, 25 insertions(+), 19 deletions(-) diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index ff48f13c..847e0b19 100644 --- a/wmpl/MetSim/GUI.py +++ b/wmpl/MetSim/GUI.py @@ -3280,9 +3280,9 @@ def updateInputBoxes(self, show_previous=False): self.inputTimeStep.setText(str(const.dt)) # Adaptive time step toggle: when on, dt is only the output cadence and its box is disabled - self.checkBoxAdaptiveDt.setChecked(getattr(const, 'adaptive_dt', False)) + self.checkBoxAdaptiveDt.setChecked(getattr(const, 'adaptive_dt', True)) self.inputTimeStep.setEnabled(not self.checkBoxAdaptiveDt.isChecked()) - self.inputAdaptiveRtol.setText("{:.1e}".format(getattr(const, 'adaptive_rtol', 1e-5))) + self.inputAdaptiveRtol.setText("{:.1e}".format(getattr(const, 'adaptive_rtol', 1e-4))) self.inputHtInit.setText("{:.3f}".format(const.h_init/1000)) self.inputP0M.setText("{:d}".format(int(const.P_0m))) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index 46ec881f..b13d2a43 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -367,6 +367,9 @@ adaptive dt + + true + @@ -378,7 +381,7 @@ - Relative error tolerance for the adaptive time step (only used when 'adaptive' is checked). Each sub-step is accepted only if its estimated relative error in mass and speed is below rtol; otherwise the step is shrunk and retried. Smaller rtol = more accurate but more sub-steps (slower). Rough guide for the single-body velocity error: 1e-4 fast/coarse (~1 km/s), 1e-5 default (~0.3 km/s), 1e-6 high-precision / CAMO-grade (~0.05 km/s). Set it so the numerical error sits below your measurement noise. Ignored when 'adaptive' is off. + Relative error tolerance for the adaptive time step (only used when 'adaptive' is checked). Each sub-step is accepted only if its estimated relative error in mass and speed is below rtol; otherwise the step is shrunk and retried. Smaller rtol = more accurate but more sub-steps (slower). Rough guide for the single-body velocity error: 1e-4 default (fast; light curve indistinguishable from fixed-step), 1e-5 tighter (~0.3 km/s), 1e-6 high-precision / CAMO-grade (~0.05 km/s). Set it so the numerical error sits below your measurement noise. Ignored when 'adaptive' is off. diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 6bb5d432..2249ec4e 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -47,18 +47,20 @@ def __init__(self): self.dt = 0.005 ### Adaptive sub-stepping (opt-in). Default False -> the engine runs the original fixed-step - # path and reproduces prior results exactly. When True, each fragment sub-steps adaptively - # within each dt to meet the tolerances below, while the output cadence stays fixed at dt. ### - self.adaptive_dt = False + # path and reproduces prior results exactly. When True (default), each fragment sub-steps + # adaptively within each dt to meet the tolerances below, while the output cadence stays fixed + # at dt. Adaptive is on by default: it matches the fixed-step light curve to well within + # measurement noise while running several times faster. Set False for bit-for-bit legacy runs. ### + self.adaptive_dt = True # Adaptive integrator: True -> embedded Dormand-Prince RK45 on the coupled (mass, velocity) # system (fewer RHS evals per sub-step, higher order -> ~3-4x faster); False -> step-doubling # on the original operator-split RK4. Both meet the same tolerance; DP45 is the default. self.adaptive_high_order = True # Tolerance targets NUMERICAL error below MEASUREMENT noise (~0.1 mag, ~0.1 km/s at ~25 FPS), - # not machine convergence. rtol=1e-5 keeps the single-body velocity error ~0.1-0.3 km/s (at or - # below typical measurement precision) at a fraction of the cost of a tighter tol; tighten to - # 1e-6 (~0.05 km/s) for high-precision (e.g. CAMO) data if needed. - self.adaptive_rtol = 1e-5 # relative tolerance on mass and speed + # not machine convergence. rtol=1e-4 keeps the light curve and dynamics indistinguishable from + # the fixed-step result (well under measurement precision) while running fastest; tighten to + # 1e-5/1e-6 for high-precision (e.g. CAMO) data if needed. + self.adaptive_rtol = 1e-4 # relative tolerance on mass and speed self.adaptive_atol_m = 1e-14 # absolute mass tolerance (kg); ~ m_kill self.adaptive_atol_v = 0.1 # absolute speed tolerance (m/s); floor well under meas. noise self.adaptive_dt_min = 1e-7 # smallest allowed sub-step (s) @@ -753,13 +755,13 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): ... Note: - With const.adaptive_dt False (default), const.dt is the fixed RK4 step for every fragment, - including grains. Near erosion_mass_min this single-step RK4 is under-resolved (a single step - can overshoot the grain's true velocity by tens of percent vs. a finely-substepped mirror). - Set const.adaptive_dt True to integrate each fragment with error-controlled adaptive sub-steps - (const.adaptive_rtol etc.); dt then only sets the output cadence, and erosion grains are shed - per sub-step so the light-curve shape stops depending on the step size. Fixed-step output is - unchanged when adaptive_dt is False. + With const.adaptive_dt False, const.dt is the fixed RK4 step for every fragment, including + grains. Near erosion_mass_min this single-step RK4 is under-resolved (a single step can overshoot + the grain's true velocity by tens of percent vs. a finely-substepped mirror). With + const.adaptive_dt True (default), each fragment is integrated with error-controlled adaptive + sub-steps (const.adaptive_rtol etc.); dt then only sets the output cadence, and erosion grains + are shed per sub-step so the light-curve shape stops depending on the step size. Fixed-step + output is unchanged when adaptive_dt is False. """ # Keep track of the total luminosity @@ -1483,8 +1485,9 @@ def runSimulation(const, compute_wake=False): """ Run the ablation simulation. """ # Back-fill adaptive-timestep settings for Constants loaded from older JSONs that predate them, so - # such runs default to the original fixed-step behaviour. Also reset the per-run diagnostics. - _adaptive_defaults = {'adaptive_dt': False, 'adaptive_high_order': True, 'adaptive_rtol': 1e-5, + # such runs pick up the current defaults (adaptive on). Set adaptive_dt False in the JSON/GUI for + # bit-for-bit legacy fixed-step behaviour. Also reset the per-run diagnostics. + _adaptive_defaults = {'adaptive_dt': True, 'adaptive_high_order': True, 'adaptive_rtol': 1e-4, 'adaptive_atol_m': 1e-14, 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, 'adaptive_dt_max': const.dt, 'adaptive_max_substeps': 10000, 'erosion_release_length': 50.0} From 267cf478a7c12fe31dd0c814dec117a030feae7c Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 21:56:27 +0200 Subject: [PATCH 11/16] Keep adaptive rtol default at 1e-5 Benchmarking the 2019 Orionids showed rtol=1e-4 gives no speedup over 1e-5: for eroding meteors the sub-step count is set by erosion_release_length (the 50 m grain-release cap), not by the error tolerance, and single-body cases are bounded by dt_max. Since 1e-5 is tighter (more accuracy headroom where rtol does bind) at negligible extra cost, keep it as the default. adaptive_dt stays True. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_0179geUg2TGRWWqpNPx2k234 --- wmpl/MetSim/GUI.py | 2 +- wmpl/MetSim/GUI.ui | 2 +- wmpl/MetSim/MetSimErosion.py | 11 ++++++----- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index 847e0b19..a74d0876 100644 --- a/wmpl/MetSim/GUI.py +++ b/wmpl/MetSim/GUI.py @@ -3282,7 +3282,7 @@ def updateInputBoxes(self, show_previous=False): # Adaptive time step toggle: when on, dt is only the output cadence and its box is disabled self.checkBoxAdaptiveDt.setChecked(getattr(const, 'adaptive_dt', True)) self.inputTimeStep.setEnabled(not self.checkBoxAdaptiveDt.isChecked()) - self.inputAdaptiveRtol.setText("{:.1e}".format(getattr(const, 'adaptive_rtol', 1e-4))) + self.inputAdaptiveRtol.setText("{:.1e}".format(getattr(const, 'adaptive_rtol', 1e-5))) self.inputHtInit.setText("{:.3f}".format(const.h_init/1000)) self.inputP0M.setText("{:d}".format(int(const.P_0m))) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index b13d2a43..40510148 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -381,7 +381,7 @@ - Relative error tolerance for the adaptive time step (only used when 'adaptive' is checked). Each sub-step is accepted only if its estimated relative error in mass and speed is below rtol; otherwise the step is shrunk and retried. Smaller rtol = more accurate but more sub-steps (slower). Rough guide for the single-body velocity error: 1e-4 default (fast; light curve indistinguishable from fixed-step), 1e-5 tighter (~0.3 km/s), 1e-6 high-precision / CAMO-grade (~0.05 km/s). Set it so the numerical error sits below your measurement noise. Ignored when 'adaptive' is off. + Relative error tolerance for the adaptive time step (only used when 'adaptive' is checked). Each sub-step is accepted only if its estimated relative error in mass and speed is below rtol; otherwise the step is shrunk and retried. Smaller rtol = more accurate but more sub-steps (slower). Rough guide for the single-body velocity error: 1e-4 fast/coarse (~1 km/s), 1e-5 default (~0.3 km/s; light curve indistinguishable from fixed-step), 1e-6 high-precision / CAMO-grade (~0.05 km/s). Set it so the numerical error sits below your measurement noise. Ignored when 'adaptive' is off. diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 2249ec4e..853a0a55 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -57,10 +57,11 @@ def __init__(self): # on the original operator-split RK4. Both meet the same tolerance; DP45 is the default. self.adaptive_high_order = True # Tolerance targets NUMERICAL error below MEASUREMENT noise (~0.1 mag, ~0.1 km/s at ~25 FPS), - # not machine convergence. rtol=1e-4 keeps the light curve and dynamics indistinguishable from - # the fixed-step result (well under measurement precision) while running fastest; tighten to - # 1e-5/1e-6 for high-precision (e.g. CAMO) data if needed. - self.adaptive_rtol = 1e-4 # relative tolerance on mass and speed + # not machine convergence. rtol=1e-5 keeps the light curve and dynamics indistinguishable from + # the fixed-step result (well under measurement precision) at negligible extra cost over 1e-4 + # (for eroding meteors the sub-step count is set by erosion_release_length, not rtol); loosen to + # 1e-4 for a hair more speed on tolerance-bound cases, or tighten to 1e-6 for CAMO-grade data. + self.adaptive_rtol = 1e-5 # relative tolerance on mass and speed self.adaptive_atol_m = 1e-14 # absolute mass tolerance (kg); ~ m_kill self.adaptive_atol_v = 0.1 # absolute speed tolerance (m/s); floor well under meas. noise self.adaptive_dt_min = 1e-7 # smallest allowed sub-step (s) @@ -1487,7 +1488,7 @@ def runSimulation(const, compute_wake=False): # Back-fill adaptive-timestep settings for Constants loaded from older JSONs that predate them, so # such runs pick up the current defaults (adaptive on). Set adaptive_dt False in the JSON/GUI for # bit-for-bit legacy fixed-step behaviour. Also reset the per-run diagnostics. - _adaptive_defaults = {'adaptive_dt': True, 'adaptive_high_order': True, 'adaptive_rtol': 1e-4, + _adaptive_defaults = {'adaptive_dt': True, 'adaptive_high_order': True, 'adaptive_rtol': 1e-5, 'adaptive_atol_m': 1e-14, 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, 'adaptive_dt_max': const.dt, 'adaptive_max_substeps': 10000, 'erosion_release_length': 50.0} From 75dd796c4c6f0261a9a1c2c4af8b8de64b38daef Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Wed, 22 Jul 2026 22:13:46 +0200 Subject: [PATCH 12/16] Make grain-release cadence tunable + velocity-scaled (bounds fast-meteor cost) The adaptive steppers cap an eroding fragment's sub-step at erosion_release_length/v so grains are shed every ~50 m of flight (a dt-independent grain-birth resolution). Because that cadence is fixed in distance, a fast meteor takes tiny sub-steps and the cost grows with velocity - e.g. a 67 km/s Orionid spends ~2x the sub-steps of a slower one for no visible gain. (A slow, massive fireball is unaffected: the error controller already takes small steps, so the cap never binds.) - Add erosion_release_vref (default 30000 m/s): above this speed the cadence is frozen in TIME (cap = erosion_release_length/erosion_release_vref) instead of distance, so the sub-step count stays bounded. Fragments slower than vref keep the exact distance cadence, so their output is unchanged. vref <= 0 disables the cap (original pure-distance behaviour). - Expose erosion_release_length in the GUI (box next to rtol, enabled with the adaptive checkbox) so it can be coarsened for expensive fits; vref stays a JSON-level knob. - Thread the new parameter through both steppers and the runSimulation back-fill. Validated on the 13 fitted 2019 Orionids (all ~67 km/s): mean 1.81x faster with the vref=30000 default, max peak-magnitude shift 0.020 mag (well under the ~0.1 mag measurement noise; the change is confined to the faint erosion tail). Fixed mode is bit-for-bit invariant to vref, and the regression suite passes (8/8). Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_0179geUg2TGRWWqpNPx2k234 --- wmpl/MetSim/GUI.py | 6 +++++ wmpl/MetSim/GUI.ui | 34 ++++++++++++++++++++++++++- wmpl/MetSim/MetSimErosion.py | 16 +++++++++---- wmpl/MetSim/MetSimErosionCyTools.pyx | 35 +++++++++++++++++++++------- 4 files changed, 77 insertions(+), 14 deletions(-) diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index a74d0876..a9726506 100644 --- a/wmpl/MetSim/GUI.py +++ b/wmpl/MetSim/GUI.py @@ -3283,6 +3283,8 @@ def updateInputBoxes(self, show_previous=False): self.checkBoxAdaptiveDt.setChecked(getattr(const, 'adaptive_dt', True)) self.inputTimeStep.setEnabled(not self.checkBoxAdaptiveDt.isChecked()) self.inputAdaptiveRtol.setText("{:.1e}".format(getattr(const, 'adaptive_rtol', 1e-5))) + self.inputErosionReleaseLength.setText( + "{:.0f}".format(getattr(const, 'erosion_release_length', 50.0))) self.inputHtInit.setText("{:.3f}".format(const.h_init/1000)) self.inputP0M.setText("{:d}".format(int(const.P_0m))) @@ -3597,6 +3599,8 @@ def checkBoxAdaptiveDtSignal(self, event): # rtol only matters in adaptive mode: enable its box/label with the checkbox, disable otherwise self.inputAdaptiveRtol.setEnabled(self.const.adaptive_dt) self.labelAdaptiveRtol.setEnabled(self.const.adaptive_dt) + self.inputErosionReleaseLength.setEnabled(self.const.adaptive_dt) + self.labelErosionReleaseLength.setEnabled(self.const.adaptive_dt) def checkBoxFragmentationShowIndividualLCsSignal(self, event): @@ -3683,6 +3687,8 @@ def readInputBoxes(self): self.const.adaptive_dt = self.checkBoxAdaptiveDt.isChecked() self.const.adaptive_rtol = self._tryReadBox(self.inputAdaptiveRtol, getattr(self.const, 'adaptive_rtol', 1e-5)) + self.const.erosion_release_length = self._tryReadBox(self.inputErosionReleaseLength, + getattr(self.const, 'erosion_release_length', 50.0)) self.const.P_0m = self._tryReadBox(self.inputP0M, self.const.P_0m) self.const.h_init = 1000*self._tryReadBox(self.inputHtInit, self.const.h_init/1000) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index 40510148..b0ee08e3 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -70,7 +70,7 @@ 1510 0 391 - 121 + 145 @@ -400,6 +400,38 @@ Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + 40 + 113 + 61 + 23 + + + + Grain-release interval in metres along the trail (only used when 'adaptive' is checked, and only for eroding meteors). In adaptive mode grains are shed once per sub-step, so this sets the grain-birth resolution: the eroding fragment's sub-step is capped so grains are released roughly every this many metres, making the light-curve shape independent of dt and rtol. Smaller = finer grain sampling but slower; larger = faster but coarser (mostly the faint erosion tail; the peak is barely affected). Default 50 m. For fast meteors the cadence is automatically frozen in time above a reference speed (erosion_release_vref, JSON) so cost stays bounded. Ignored when 'adaptive' is off. + + + + + + 2 + 113 + 36 + 16 + + + + Grain-release interval along the trail (m). Larger = faster but coarser erosion tail. Default 50 m. + + + grain m + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 853a0a55..631e22b9 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -74,10 +74,17 @@ def __init__(self): # metres of flight - a physical resolution independent of dt and the error tolerance, so the # erosion light-curve shape stops depending on dt. Smaller = finer grain sampling + slower. # Only used in adaptive mode and only for eroding fragments (grains do not erode further). - # Not shown in the GUI - it is a dt-independent physical default that rarely needs changing; - # tune via the sim-fit JSON if finer/coarser grain sampling is wanted. + # Exposed in the GUI so it can be coarsened for expensive fits (see erosion_release_vref). self.erosion_release_length = 50.0 + # Reference speed (m/s) that turns the along-track cadence into a fixed TIME cadence above it. + # The distance cap fixes cadence in metres, so a fast meteor takes tiny sub-steps (cost grows + # with v). Fragments faster than erosion_release_vref instead cap their sub-step at + # erosion_release_length/erosion_release_vref (a constant time), bounding the sub-step count; + # slower fragments keep the exact distance cadence. This keeps the per-event cost roughly + # velocity-independent. Set <= 0 to disable (pure distance cadence, the original behaviour). + self.erosion_release_vref = 30000.0 + # Time elapsed since the beginning self.total_time = 0 @@ -844,7 +851,8 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): const.h_init, const.zenith_angle, const.r_earth, G0, const.dens_co, const.adaptive_rtol, const.adaptive_atol_m, const.adaptive_atol_v, const.m_kill, const.adaptive_dt_min, const.adaptive_dt_max, const.adaptive_max_substeps, - frag.adaptive_h_sub, getattr(const, 'erosion_release_length', 50.0)) + frag.adaptive_h_sub, getattr(const, 'erosion_release_length', 50.0), + getattr(const, 'erosion_release_vref', 30000.0)) # Diagnostics (feed the cost study; also used to warn once on runaway/under-resolved steps) const.adaptive_substeps_total += n_sub @@ -1491,7 +1499,7 @@ def runSimulation(const, compute_wake=False): _adaptive_defaults = {'adaptive_dt': True, 'adaptive_high_order': True, 'adaptive_rtol': 1e-5, 'adaptive_atol_m': 1e-14, 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, 'adaptive_dt_max': const.dt, 'adaptive_max_substeps': 10000, - 'erosion_release_length': 50.0} + 'erosion_release_length': 50.0, 'erosion_release_vref': 30000.0} for _attr, _default in _adaptive_defaults.items(): if not hasattr(const, _attr): setattr(const, _attr, _default) diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index 1e27db91..29e75217 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/wmpl/MetSim/MetSimErosionCyTools.pyx @@ -692,7 +692,7 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero double h_grav_drop_total, double h_init, double zenith_angle, double r_earth, double g0, FLOAT_TYPE_t[:] dens_co, double rtol, double atol_m, double atol_v, double m_kill, double dt_min, double dt_max, int max_substeps, double h_sub_init, - double erosion_release_length): + double erosion_release_length, double erosion_release_vref): """ Advance ONE fragment across a full macro interval dt_macro using error-controlled adaptive sub-steps (step-doubling on the operator-split RK4), refreshing the atmosphere/height each sub-step. Reproduces the single-body advance of the fixed engine in the one-sub-step limit, but @@ -727,8 +727,14 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero max_substeps: [int] Sub-step cap per macro interval (runaway guard). h_sub_init: [double] Warm-start sub-step from the previous macro step (s); <= 0 means "use dt_macro". erosion_release_length: [double] Along-track grain-release interval (m). For an eroding fragment - the sub-step is capped at erosion_release_length/v so grains are shed ~every this many metres, - making the grain-birth cadence independent of dt/rtol. Ignored for non-eroding fragments. + the sub-step is capped at erosion_release_length/min(v, erosion_release_vref) so grains are + shed ~every this many metres, making the grain-birth cadence independent of dt/rtol. Ignored + for non-eroding fragments. + erosion_release_vref: [double] Reference speed (m/s) above which the along-track cadence is frozen + in TIME rather than distance: fragments faster than this cap their sub-step at + erosion_release_length/erosion_release_vref instead of .../v, bounding the sub-step count for + fast meteors (cost grows with v otherwise). Slower fragments keep the exact distance cadence. + Set <= 0 to disable the cap (pure distance cadence, the original behaviour). Return: A 17-tuple: @@ -753,6 +759,7 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero """ cdef double t = 0.0 + cdef double v_cap cdef double h_sub, h_cur, rho_atm, rho_mid, rho_last cdef double dm_abl_big, dm_ero_big, dm_tot_big, decel_big cdef double dm_abl_1, dm_ero_1, dm_tot_1, decel_1 @@ -810,8 +817,13 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero # Cap an eroding fragment's sub-step so grains are shed roughly every erosion_release_length # metres of flight, making the grain-birth cadence (and hence the erosion light curve) # independent of dt and the error tolerance. Not applied to non-eroding fragments or grains. - if erosion_active and (v > 0) and (erosion_release_length > 0) and (h_sub*v > erosion_release_length): - h_sub = erosion_release_length/v + # Above erosion_release_vref the cadence is frozen in time (uses v_cap = vref, not v) so the + # sub-step count stays bounded for fast meteors; slower fragments keep the pure distance cadence. + v_cap = v + if (erosion_release_vref > 0) and (v > erosion_release_vref): + v_cap = erosion_release_vref + if erosion_active and (v > 0) and (erosion_release_length > 0) and (h_sub*v_cap > erosion_release_length): + h_sub = erosion_release_length/v_cap hh = 0.5*h_sub @@ -1024,7 +1036,7 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c double h_grav_drop_total, double h_init, double zenith_angle, double r_earth, double g0, FLOAT_TYPE_t[:] dens_co, double rtol, double atol_m, double atol_v, double m_kill, double dt_min, double dt_max, int max_substeps, double h_sub_init, - double erosion_release_length): + double erosion_release_length, double erosion_release_vref): """ Advance ONE fragment across a full macro interval dt_macro, like adaptiveSingleBodyStep, but with an embedded Dormand-Prince RK45 pair (5th-order solution + 4th-order error estimate) on the COUPLED system y = [m, vv, vh, length] instead of step-doubling on the operator split. One @@ -1051,7 +1063,7 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c cdef double bs6 = 187.0/2100, bs7 = 1.0/40 cdef double t = 0.0 - cdef double h_sub, hh, gv, h_cur, rho0, v_cur + cdef double h_sub, hh, gv, h_cur, rho0, v_cur, v_cap cdef double y[4] cdef double yt[4] cdef double k1[4] @@ -1101,10 +1113,15 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c # Cap an eroding fragment's sub-step so grains are shed roughly every erosion_release_length # metres of flight - a grain-birth cadence independent of dt and the error tolerance (see the # step-doubling stepper for the rationale). Only for eroding fragments; grains do not erode. + # Above erosion_release_vref the cadence is frozen in time (v_cap = vref) to bound the sub-step + # count for fast meteors; slower fragments keep the pure distance cadence. v_cur = sqrt(y[1]*y[1] + y[2]*y[2]) + v_cap = v_cur + if (erosion_release_vref > 0) and (v_cur > erosion_release_vref): + v_cap = erosion_release_vref if erosion_active and (v_cur > 0) and (erosion_release_length > 0) \ - and (h_sub*v_cur > erosion_release_length): - h_sub = erosion_release_length/v_cur + and (h_sub*v_cap > erosion_release_length): + h_sub = erosion_release_length/v_cap # Height/gravity at the sub-step start (gravity drop uses the start-of-step g, as in the fixed model) h_cur = heightCurvatureC(h_init, zenith_angle, y[3], r_earth) - h_grav_drop_total From a1f6054cd48a4d0c570d3e6c05b21366b20a4a83 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Thu, 23 Jul 2026 09:20:14 +0200 Subject: [PATCH 13/16] Normalize GUI.ui Qt6 enum scoping back to Qt5 form for PyQt5 loadUi Qt Designer re-saves enum values in the Qt6 scoped form (Qt::AlignmentFlag::AlignCenter, Qt::Orientation::Horizontal), which PyQt5's uic.loadUi cannot parse. Strip the scope prefixes back to the Qt5 short form (Qt::AlignCenter, Qt::Horizontal). No layout change. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_0179geUg2TGRWWqpNPx2k234 --- wmpl/MetSim/GUI.ui | 51 +++++++++++++++++++++++++++++----------------- 1 file changed, 32 insertions(+), 19 deletions(-) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index b0ee08e3..4a39ab88 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -70,7 +70,7 @@ 1510 0 391 - 145 + 121 @@ -82,8 +82,8 @@ - 40 - 30 + 48 + 12 61 23 @@ -92,8 +92,8 @@ - -25 - 30 + -17 + 12 61 16 @@ -251,8 +251,8 @@ - 110 - 30 + 118 + 12 21 16 @@ -355,8 +355,8 @@ - 23 - 60 + 31 + 35 88 20 @@ -374,8 +374,8 @@ - 40 - 90 + 48 + 58 61 23 @@ -387,8 +387,8 @@ - 10 - 90 + 18 + 58 28 16 @@ -403,8 +403,8 @@ - 40 - 113 + 48 + 90 61 23 @@ -416,8 +416,8 @@ - 2 - 113 + 10 + 90 36 16 @@ -426,7 +426,7 @@ Grain-release interval along the trail (m). Larger = faster but coarser erosion tail. Default 50 m. - grain m + grain Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter @@ -445,6 +445,19 @@ Qt::Vertical + + + + 114 + 90 + 21 + 16 + + + + m + + @@ -1265,7 +1278,7 @@ 860 - 0 + 4 641 551 From 2c475e17e535d0da6d34815a12ed33b0e310f757 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Thu, 23 Jul 2026 09:28:34 +0200 Subject: [PATCH 14/16] Raise default grain-release length 50 m -> 200 m Benchmarking showed 200 m is substantially faster on eroding fits while the light-curve peak is barely affected (the difference is confined to the faint erosion tail, well under measurement noise). Update the Constants default, the runSimulation back-fill, the GUI getattr fallbacks, and the tooltips. Finer sampling is still available by lowering the value in the GUI/JSON. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_0179geUg2TGRWWqpNPx2k234 --- wmpl/MetSim/GUI.py | 4 ++-- wmpl/MetSim/GUI.ui | 4 ++-- wmpl/MetSim/MetSimErosion.py | 6 +++--- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index a9726506..9f1ec1cc 100644 --- a/wmpl/MetSim/GUI.py +++ b/wmpl/MetSim/GUI.py @@ -3284,7 +3284,7 @@ def updateInputBoxes(self, show_previous=False): self.inputTimeStep.setEnabled(not self.checkBoxAdaptiveDt.isChecked()) self.inputAdaptiveRtol.setText("{:.1e}".format(getattr(const, 'adaptive_rtol', 1e-5))) self.inputErosionReleaseLength.setText( - "{:.0f}".format(getattr(const, 'erosion_release_length', 50.0))) + "{:.0f}".format(getattr(const, 'erosion_release_length', 200.0))) self.inputHtInit.setText("{:.3f}".format(const.h_init/1000)) self.inputP0M.setText("{:d}".format(int(const.P_0m))) @@ -3688,7 +3688,7 @@ def readInputBoxes(self): self.const.adaptive_rtol = self._tryReadBox(self.inputAdaptiveRtol, getattr(self.const, 'adaptive_rtol', 1e-5)) self.const.erosion_release_length = self._tryReadBox(self.inputErosionReleaseLength, - getattr(self.const, 'erosion_release_length', 50.0)) + getattr(self.const, 'erosion_release_length', 200.0)) self.const.P_0m = self._tryReadBox(self.inputP0M, self.const.P_0m) self.const.h_init = 1000*self._tryReadBox(self.inputHtInit, self.const.h_init/1000) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index 4a39ab88..7eb8bd7b 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -410,7 +410,7 @@ - Grain-release interval in metres along the trail (only used when 'adaptive' is checked, and only for eroding meteors). In adaptive mode grains are shed once per sub-step, so this sets the grain-birth resolution: the eroding fragment's sub-step is capped so grains are released roughly every this many metres, making the light-curve shape independent of dt and rtol. Smaller = finer grain sampling but slower; larger = faster but coarser (mostly the faint erosion tail; the peak is barely affected). Default 50 m. For fast meteors the cadence is automatically frozen in time above a reference speed (erosion_release_vref, JSON) so cost stays bounded. Ignored when 'adaptive' is off. + Grain-release interval in metres along the trail (only used when 'adaptive' is checked, and only for eroding meteors). In adaptive mode grains are shed once per sub-step, so this sets the grain-birth resolution: the eroding fragment's sub-step is capped so grains are released roughly every this many metres, making the light-curve shape independent of dt and rtol. Smaller = finer grain sampling but slower; larger = faster but coarser (mostly the faint erosion tail; the peak is barely affected). Default 200 m. For fast meteors the cadence is automatically frozen in time above a reference speed (erosion_release_vref, JSON) so cost stays bounded. Ignored when 'adaptive' is off. @@ -423,7 +423,7 @@ - Grain-release interval along the trail (m). Larger = faster but coarser erosion tail. Default 50 m. + Grain-release interval along the trail (m). Larger = faster but coarser erosion tail. Default 200 m. grain diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 631e22b9..231d0495 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -75,7 +75,7 @@ def __init__(self): # erosion light-curve shape stops depending on dt. Smaller = finer grain sampling + slower. # Only used in adaptive mode and only for eroding fragments (grains do not erode further). # Exposed in the GUI so it can be coarsened for expensive fits (see erosion_release_vref). - self.erosion_release_length = 50.0 + self.erosion_release_length = 200.0 # Reference speed (m/s) that turns the along-track cadence into a fixed TIME cadence above it. # The distance cap fixes cadence in metres, so a fast meteor takes tiny sub-steps (cost grows @@ -851,7 +851,7 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): const.h_init, const.zenith_angle, const.r_earth, G0, const.dens_co, const.adaptive_rtol, const.adaptive_atol_m, const.adaptive_atol_v, const.m_kill, const.adaptive_dt_min, const.adaptive_dt_max, const.adaptive_max_substeps, - frag.adaptive_h_sub, getattr(const, 'erosion_release_length', 50.0), + frag.adaptive_h_sub, getattr(const, 'erosion_release_length', 200.0), getattr(const, 'erosion_release_vref', 30000.0)) # Diagnostics (feed the cost study; also used to warn once on runaway/under-resolved steps) @@ -1499,7 +1499,7 @@ def runSimulation(const, compute_wake=False): _adaptive_defaults = {'adaptive_dt': True, 'adaptive_high_order': True, 'adaptive_rtol': 1e-5, 'adaptive_atol_m': 1e-14, 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, 'adaptive_dt_max': const.dt, 'adaptive_max_substeps': 10000, - 'erosion_release_length': 50.0, 'erosion_release_vref': 30000.0} + 'erosion_release_length': 200.0, 'erosion_release_vref': 30000.0} for _attr, _default in _adaptive_defaults.items(): if not hasattr(const, _attr): setattr(const, _attr, _default) From 932f8b7720c977557ed40ac1c3251ad26fa9addb Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Thu, 23 Jul 2026 09:31:12 +0200 Subject: [PATCH 15/16] Read adaptive/grain params straight from const (drop duplicated defaults) The GUI display/readback and the ablateAll dispatch used getattr(const, 'name', ) with the default value hardcoded a second time, so a default change had to be edited in two places. Every const starts from Constants() and runSimulation back-fills older JSONs, so the attributes are always present - read them directly (const.adaptive_rtol, const.adaptive_dt, const.adaptive_high_order, const.erosion_release_length, const.erosion_release_vref), matching how h_init/P_0m/dt are already handled. Single source of truth is now Constants.__init__ (plus the runSimulation back-fill table for migration). Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_0179geUg2TGRWWqpNPx2k234 --- wmpl/MetSim/GUI.py | 12 +++++------- wmpl/MetSim/MetSimErosion.py | 6 ++---- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index 9f1ec1cc..54b1a597 100644 --- a/wmpl/MetSim/GUI.py +++ b/wmpl/MetSim/GUI.py @@ -3280,11 +3280,10 @@ def updateInputBoxes(self, show_previous=False): self.inputTimeStep.setText(str(const.dt)) # Adaptive time step toggle: when on, dt is only the output cadence and its box is disabled - self.checkBoxAdaptiveDt.setChecked(getattr(const, 'adaptive_dt', True)) + self.checkBoxAdaptiveDt.setChecked(const.adaptive_dt) self.inputTimeStep.setEnabled(not self.checkBoxAdaptiveDt.isChecked()) - self.inputAdaptiveRtol.setText("{:.1e}".format(getattr(const, 'adaptive_rtol', 1e-5))) - self.inputErosionReleaseLength.setText( - "{:.0f}".format(getattr(const, 'erosion_release_length', 200.0))) + self.inputAdaptiveRtol.setText("{:.1e}".format(const.adaptive_rtol)) + self.inputErosionReleaseLength.setText("{:.0f}".format(const.erosion_release_length)) self.inputHtInit.setText("{:.3f}".format(const.h_init/1000)) self.inputP0M.setText("{:d}".format(int(const.P_0m))) @@ -3685,10 +3684,9 @@ def readInputBoxes(self): self.const.dt = self._tryReadBox(self.inputTimeStep, self.const.dt) self.const.adaptive_dt = self.checkBoxAdaptiveDt.isChecked() - self.const.adaptive_rtol = self._tryReadBox(self.inputAdaptiveRtol, - getattr(self.const, 'adaptive_rtol', 1e-5)) + self.const.adaptive_rtol = self._tryReadBox(self.inputAdaptiveRtol, self.const.adaptive_rtol) self.const.erosion_release_length = self._tryReadBox(self.inputErosionReleaseLength, - getattr(self.const, 'erosion_release_length', 200.0)) + self.const.erosion_release_length) self.const.P_0m = self._tryReadBox(self.inputP0M, self.const.P_0m) self.const.h_init = 1000*self._tryReadBox(self.inputHtInit, self.const.h_init/1000) diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 231d0495..36a74d96 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -839,8 +839,7 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): erosion_active = 1 if (frag.erosion_enabled and (frag.erosion_coeff > 0)) else 0 - _stepper = adaptiveDP45Step if getattr(const, 'adaptive_high_order', True) \ - else adaptiveSingleBodyStep + _stepper = adaptiveDP45Step if const.adaptive_high_order else adaptiveSingleBodyStep (frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, frag.h, rho_atm, mass_loss_ablation, mass_loss_erosion, deceleration_total, went_up, n_sub, @@ -851,8 +850,7 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): const.h_init, const.zenith_angle, const.r_earth, G0, const.dens_co, const.adaptive_rtol, const.adaptive_atol_m, const.adaptive_atol_v, const.m_kill, const.adaptive_dt_min, const.adaptive_dt_max, const.adaptive_max_substeps, - frag.adaptive_h_sub, getattr(const, 'erosion_release_length', 200.0), - getattr(const, 'erosion_release_vref', 30000.0)) + frag.adaptive_h_sub, const.erosion_release_length, const.erosion_release_vref) # Diagnostics (feed the cost study; also used to warn once on runaway/under-resolved steps) const.adaptive_substeps_total += n_sub From a053a41ad59368b930118a94d3b40cabaeaea902 Mon Sep 17 00:00:00 2001 From: Denis Vida Date: Thu, 23 Jul 2026 09:46:01 +0200 Subject: [PATCH 16/16] Style: match codebase conventions in new adaptive-stepper code Formatting-only cleanup from a review of the new code against the rest of the file/repo (no behaviour change; regression suite still 8/8, fixed mode bit-for-bit): - Split semicolon-packed multi-assignments onto one statement per line (the file has zero semicolon statements elsewhere), in both steppers and _rhsDP. - Add the blank line after the closing docstring """ in clampMassC, heightCurvatureC, atmDensityPolyC, advanceVelPosC and _rhsDP. - Replace fixed-path line-number references in comments with behaviour descriptions (line numbers drift). - Wrap the one >110-col erosion-cap line to match its DP45 twin; capitalize two lowercase-first comments. - MetSimErosion.py: drop the stale "opt-in / Default False" wording now that adaptive is the default; stop misusing the "### ... ###" section marker for a paragraph; document the tolerance fields with preceding comment lines instead of column-aligned inline comments (matches the rest of Constants); rename underscore-prefixed locals (_stepper, _saved_state, _adaptive_defaults, etc.) to plain names. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_0179geUg2TGRWWqpNPx2k234 --- wmpl/MetSim/MetSimErosion.py | 71 ++++++++++++++++------------ wmpl/MetSim/MetSimErosionCyTools.pyx | 64 +++++++++++++++++++------ 2 files changed, 89 insertions(+), 46 deletions(-) diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index 36a74d96..6e0ea5e4 100644 --- a/wmpl/MetSim/MetSimErosion.py +++ b/wmpl/MetSim/MetSimErosion.py @@ -46,27 +46,36 @@ def __init__(self): # sub-steps underneath (see adaptive_* below and ablateAll). self.dt = 0.005 - ### Adaptive sub-stepping (opt-in). Default False -> the engine runs the original fixed-step - # path and reproduces prior results exactly. When True (default), each fragment sub-steps - # adaptively within each dt to meet the tolerances below, while the output cadence stays fixed - # at dt. Adaptive is on by default: it matches the fixed-step light curve to well within - # measurement noise while running several times faster. Set False for bit-for-bit legacy runs. ### + # Adaptive sub-stepping. When True (the default), each fragment sub-steps adaptively within each + # dt to meet the tolerances below, while the output cadence stays fixed at dt; this matches the + # fixed-step light curve to well within measurement noise while running several times faster. + # Set False to run the original fixed-step path, which reproduces prior results bit-for-bit. self.adaptive_dt = True # Adaptive integrator: True -> embedded Dormand-Prince RK45 on the coupled (mass, velocity) # system (fewer RHS evals per sub-step, higher order -> ~3-4x faster); False -> step-doubling # on the original operator-split RK4. Both meet the same tolerance; DP45 is the default. self.adaptive_high_order = True - # Tolerance targets NUMERICAL error below MEASUREMENT noise (~0.1 mag, ~0.1 km/s at ~25 FPS), - # not machine convergence. rtol=1e-5 keeps the light curve and dynamics indistinguishable from - # the fixed-step result (well under measurement precision) at negligible extra cost over 1e-4 - # (for eroding meteors the sub-step count is set by erosion_release_length, not rtol); loosen to - # 1e-4 for a hair more speed on tolerance-bound cases, or tighten to 1e-6 for CAMO-grade data. - self.adaptive_rtol = 1e-5 # relative tolerance on mass and speed - self.adaptive_atol_m = 1e-14 # absolute mass tolerance (kg); ~ m_kill - self.adaptive_atol_v = 0.1 # absolute speed tolerance (m/s); floor well under meas. noise - self.adaptive_dt_min = 1e-7 # smallest allowed sub-step (s) - self.adaptive_dt_max = 0.005 # largest allowed sub-step (s); should be <= dt - self.adaptive_max_substeps = 10000 # runaway guard, sub-steps per fragment per macro step + # Relative tolerance on mass and speed. Targets NUMERICAL error below MEASUREMENT noise + # (~0.1 mag, ~0.1 km/s at ~25 FPS), not machine convergence. rtol=1e-5 keeps the light curve + # and dynamics indistinguishable from the fixed-step result at negligible extra cost over 1e-4 + # (for eroding meteors the sub-step count is set by erosion_release_length, not rtol); loosen + # to 1e-4 for a hair more speed on tolerance-bound cases, or tighten to 1e-6 for CAMO data. + self.adaptive_rtol = 1e-5 + + # Absolute mass tolerance (kg); of order m_kill + self.adaptive_atol_m = 1e-14 + + # Absolute speed tolerance (m/s); floor well under the measurement noise + self.adaptive_atol_v = 0.1 + + # Smallest allowed sub-step (s) + self.adaptive_dt_min = 1e-7 + + # Largest allowed sub-step (s); should be <= dt + self.adaptive_dt_max = 0.005 + + # Runaway guard: maximum sub-steps per fragment per macro step + self.adaptive_max_substeps = 10000 # Along-track grain-release interval (m). In adaptive mode grains are shed once per accepted # sub-step, so the sub-step cadence sets the grain-birth resolution. Capping an eroding @@ -839,12 +848,12 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): erosion_active = 1 if (frag.erosion_enabled and (frag.erosion_coeff > 0)) else 0 - _stepper = adaptiveDP45Step if const.adaptive_high_order else adaptiveSingleBodyStep + stepper = adaptiveDP45Step if const.adaptive_high_order else adaptiveSingleBodyStep (frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, frag.h, rho_atm, mass_loss_ablation, mass_loss_erosion, deceleration_total, went_up, n_sub, frag.adaptive_h_sub, runaway, floor_accepts, erosion_events_adaptive) \ - = _stepper( + = stepper( const.dt, frag.K, frag.sigma, frag.erosion_coeff, erosion_active, frag.m, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total, const.h_init, const.zenith_angle, const.r_earth, G0, const.dens_co, @@ -875,8 +884,8 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): # Compute the total mass loss mass_loss_total = mass_loss_ablation + mass_loss_erosion - # If the total mass after ablation in this step is below zero, ablate what's left of the whole mass - # (i.e. land exactly on m_new = 0, not some arbitrary leftover - see m_new below) + # If the total mass after ablation in this step is below zero, ablate what is left of the + # whole mass (i.e. land exactly on m_new = 0, not some arbitrary leftover - see m_new below) if (frag.m + mass_loss_total) < 0: mass_loss_total = -frag.m @@ -1124,7 +1133,7 @@ def ablateAll(fragments, const, compute_wake=False, wake_heights_queue=None): # Create grains for erosion-enabled fragments if frag.erosion_enabled: - def _spawnGrainsFromErosion(eroded_mass): + def spawnGrainsFromErosion(eroded_mass): """ Distribute the given eroded mass into grains born from the fragment's current state. Uses the enclosing frag/const. """ grain_children, _ = generateFragments(const, frag, eroded_mass, \ @@ -1143,22 +1152,22 @@ def _spawnGrainsFromErosion(eroded_mass): # current position/velocity into each grain, so temporarily rewind the parent to each # sub-step's state, spawn, then restore its end-of-macro-step state. if erosion_events_adaptive: - _saved_state = (frag.h, frag.v, frag.vv, frag.vh, frag.length, + saved_state = (frag.h, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total) for (em, e_h, e_v, e_vv, e_vh, e_len, e_grav) in erosion_events_adaptive: if em <= 0: continue - frag.h = e_h; frag.v = e_v; frag.vv = e_vv; frag.vh = e_vh - frag.length = e_len; frag.h_grav_drop_total = e_grav - _spawnGrainsFromErosion(em) + (frag.h, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total) \ + = (e_h, e_v, e_vv, e_vh, e_len, e_grav) + spawnGrainsFromErosion(em) eroded_this_tick = True - (frag.h, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total) = _saved_state + (frag.h, frag.v, frag.vv, frag.vh, frag.length, frag.h_grav_drop_total) = saved_state else: # Fixed step: distribute the whole macro-step erosion loss at the end-of-step state if abs(mass_loss_erosion) > 0: - _spawnGrainsFromErosion(abs(mass_loss_erosion)) + spawnGrainsFromErosion(abs(mass_loss_erosion)) eroded_this_tick = True # Record erosion-begin bookkeeping for the main fragment once, at the END-of-step state so the @@ -1494,13 +1503,13 @@ def runSimulation(const, compute_wake=False): # Back-fill adaptive-timestep settings for Constants loaded from older JSONs that predate them, so # such runs pick up the current defaults (adaptive on). Set adaptive_dt False in the JSON/GUI for # bit-for-bit legacy fixed-step behaviour. Also reset the per-run diagnostics. - _adaptive_defaults = {'adaptive_dt': True, 'adaptive_high_order': True, 'adaptive_rtol': 1e-5, + adaptive_defaults = {'adaptive_dt': True, 'adaptive_high_order': True, 'adaptive_rtol': 1e-5, 'adaptive_atol_m': 1e-14, 'adaptive_atol_v': 0.1, 'adaptive_dt_min': 1e-7, 'adaptive_dt_max': const.dt, 'adaptive_max_substeps': 10000, 'erosion_release_length': 200.0, 'erosion_release_vref': 30000.0} - for _attr, _default in _adaptive_defaults.items(): - if not hasattr(const, _attr): - setattr(const, _attr, _default) + for attr, default in adaptive_defaults.items(): + if not hasattr(const, attr): + setattr(const, attr, default) const.adaptive_substeps_total = 0 const.adaptive_runaway_events = 0 const.adaptive_floor_accepts = 0 diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index 29e75217..3c49b9a7 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/wmpl/MetSim/MetSimErosionCyTools.pyx @@ -598,6 +598,7 @@ cdef inline double clampMassC(double dm, double m): Return: dm_clamped: [double] dm, or -m if (m + dm) < 0. """ + if (m + dm) < 0: return -m return dm @@ -617,6 +618,7 @@ cdef double heightCurvatureC(double h0, double zc, double l, double r_earth): Return: h: [double] Height at distance l (m), before the gravity drop is subtracted. """ + return sqrt((h0 + r_earth)*(h0 + r_earth) - 2*l*cos(zc)*(h0 + r_earth) + l*l) - r_earth @@ -632,6 +634,7 @@ cdef double atmDensityPolyC(double ht, FLOAT_TYPE_t[:] dens_co): Return: rho: [double] Atmospheric mass density at height ht (kg/m^3). """ + cdef double x = ht/1e6 return 10**(dens_co[0] + dens_co[1]*x + dens_co[2]*x*x + dens_co[3]*x*x*x + dens_co[4]*x*x*x*x + dens_co[5]*x*x*x*x*x + dens_co[6]*x*x*x*x*x*x) @@ -665,11 +668,19 @@ cdef void advanceVelPosC(double m, double v, double vv, double vh, double length Return: None (results written into 'out'). """ + cdef double gv, av, ah, vv_n, vh_n, v_n - # Accelerating (decel_rate > 0) or already stopped -> stop the fragment (mirror 773-775) + + # Accelerating (decel_rate > 0) or already stopped -> stop the fragment (mirror the fixed-path + # stop branch in ablateAll) if (decel_rate > 0) or (v <= 0): - out[0] = m + dm; out[1] = 0.0; out[2] = 0.0; out[3] = 0.0 - out[4] = length; out[5] = grav; out[6] = 0.0 + out[0] = m + dm + out[1] = 0.0 + out[2] = 0.0 + out[3] = 0.0 + out[4] = length + out[5] = grav + out[6] = 0.0 return gv = g0/((1.0 + h_at/r_earth)*(1.0 + h_at/r_earth)) av = -decel_rate*vv/v + vh*v/(r_earth + h_at) @@ -681,7 +692,7 @@ cdef void advanceVelPosC(double m, double v, double vv, double vh, double length out[1] = v_n out[2] = vv_n out[3] = vh_n - out[4] = length + v_n*h # length uses the UPDATED speed (matches 811-813 -> 824) + out[4] = length + v_n*h # length uses the UPDATED speed, as in the fixed advance out[5] = grav + 0.5*gv*h*h out[6] = 1.0 if vv_n > 0 else 0.0 # going up @@ -822,7 +833,8 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero v_cap = v if (erosion_release_vref > 0) and (v > erosion_release_vref): v_cap = erosion_release_vref - if erosion_active and (v > 0) and (erosion_release_length > 0) and (h_sub*v_cap > erosion_release_length): + if erosion_active and (v > 0) and (erosion_release_length > 0) \ + and (h_sub*v_cap > erosion_release_length): h_sub = erosion_release_length/v_cap hh = 0.5*h_sub @@ -858,7 +870,12 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero decel_1 = decelerationRK4(hh, K, fmax(m, m_kill), rho_atm, v) advanceVelPosC(m, v, vv, vh, length, h_grav_drop_total, dm_tot_1, decel_1, hh, h_cur, r_earth, g0, out1) - m_h = out1[0]; v_h = out1[1]; vv_h = out1[2]; vh_h = out1[3]; len_h = out1[4]; grav_h = out1[5] + m_h = out1[0] + v_h = out1[1] + vv_h = out1[2] + vh_h = out1[3] + len_h = out1[4] + grav_h = out1[5] h_mid = heightCurvatureC(h_init, zenith_angle, len_h, r_earth) - grav_h rho_mid = atmDensityPolyC(h_mid, dens_co) @@ -872,7 +889,8 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero decel_2 = decelerationRK4(hh, K, fmax(m_h, m_kill), rho_mid, v_h) advanceVelPosC(m_h, v_h, vv_h, vh_h, len_h, grav_h, dm_tot_2, decel_2, hh, h_mid, r_earth, g0, out2) - m_two = out2[0]; v_two = out2[1] + m_two = out2[0] + v_two = out2[1] # --- Error estimate (step doubling, RK4 order p=4 -> denom 2^p - 1 = 15) --- err_m = fabs(m_two - m_big)/15.0 @@ -897,8 +915,12 @@ cpdef adaptiveSingleBodyStep(double dt_macro, double K, double sigma, double ero # dv/dt from decelerationRK4). This feeds the luminosity deceleration term - unlike the net # (v_start - v), it excludes the gravity/curvature reallocation, matching the fixed path. dv_drag += (decel_1 + decel_2)*hh - m = out2[0]; v = out2[1]; vv = out2[2]; vh = out2[3] - length = out2[4]; h_grav_drop_total = out2[5] + m = out2[0] + v = out2[1] + vv = out2[2] + vh = out2[3] + length = out2[4] + h_grav_drop_total = out2[5] t += h_sub n_substeps += 1 @@ -1006,6 +1028,7 @@ cdef void _rhsDP(double m, double vv, double vh, double length, double h_grav_dr Return: None (results written into 'dydt' and 'extras'). """ + cdef double v, h, rho, decel, mm, mpos v = sqrt(vv*vv + vh*vh) h = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop @@ -1073,7 +1096,7 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c cdef double k5[4] cdef double k6[4] cdef double k7[4] - # extras per stage: [ablation_rate, erosion_rate, drag_decel, rho] + # Extras per stage: [ablation_rate, erosion_rate, drag_decel, rho] cdef double e1[4] cdef double e2[4] cdef double e3[4] @@ -1091,7 +1114,10 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c erosion_events = [] if erosion_active else None - y[0] = m; y[1] = vv; y[2] = vh; y[3] = length + y[0] = m + y[1] = vv + y[2] = vh + y[3] = length h_sub = h_sub_init if h_sub <= 0: @@ -1157,13 +1183,16 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c m_kill, h_init, zenith_angle, r_earth, dens_co, k7, e7) # 5th-order solution (= yt above, the 7th stage node) and error vs the 4th-order embedded - m_new = yt[0]; vv_new = yt[1]; vh_new = yt[2]; len_new = yt[3] + m_new = yt[0] + vv_new = yt[1] + vh_new = yt[2] + len_new = yt[3] v_new = sqrt(vv_new*vv_new + vh_new*vh_new) # Error estimate on mass and speed (difference of 5th- and 4th-order weights) err_m = fabs(h_sub*((b1-bs1)*k1[0] + (b3-bs3)*k3[0] + (b4-bs4)*k4[0] + (b5-bs5)*k5[0] + (b6-bs6)*k6[0] + (0.0-bs7)*k7[0])) - # velocity error via the vv/vh error components projected onto speed + # Velocity error via the vv/vh error components projected onto speed err_v = fabs(h_sub*((b1-bs1)*k1[1] + (b3-bs3)*k3[1] + (b4-bs4)*k4[1] + (b5-bs5)*k5[1] + (b6-bs6)*k6[1] + (0.0-bs7)*k7[1])) err_v += fabs(h_sub*((b1-bs1)*k1[2] + (b3-bs3)*k3[2] + (b4-bs4)*k4[2] + (b5-bs5)*k5[2] @@ -1189,7 +1218,9 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c # Commit the state; floor the mass at 0 and add the gravity drop for this sub-step y[0] = m_new if m_new > 0.0 else 0.0 - y[1] = vv_new; y[2] = vh_new; y[3] = len_new + y[1] = vv_new + y[2] = vh_new + y[3] = len_new h_grav_drop_total += 0.5*gv*h_sub*h_sub rho_last = rho0 t += h_sub @@ -1239,7 +1270,10 @@ cpdef adaptiveDP45Step(double dt_macro, double K, double sigma, double erosion_c if h_sub < dt_min: h_sub = dt_min - m = y[0]; vv = y[1]; vh = y[2]; length = y[3] + m = y[0] + vv = y[1] + vh = y[2] + length = y[3] v = sqrt(vv*vv + vh*vh) h_new = heightCurvatureC(h_init, zenith_angle, length, r_earth) - h_grav_drop_total if h_new > 0: