diff --git a/wmpl/MetSim/GUI.py b/wmpl/MetSim/GUI.py index 605f93ad..54b1a597 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,13 @@ 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(const.adaptive_dt) + self.inputTimeStep.setEnabled(not self.checkBoxAdaptiveDt.isChecked()) + 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))) self.inputMassKill.setText("{:.1e}".format(const.m_kill)) @@ -3579,6 +3588,20 @@ 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) + + # 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): """ Toggle computing light curves of individual fragments during complex fragmentation. """ @@ -3660,6 +3683,10 @@ def readInputBoxes(self): ### Simulation params ### 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, self.const.adaptive_rtol) + self.const.erosion_release_length = self._tryReadBox(self.inputErosionReleaseLength, + 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) @@ -4294,9 +4321,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, @@ -5214,7 +5248,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) @@ -5225,10 +5259,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) diff --git a/wmpl/MetSim/GUI.ui b/wmpl/MetSim/GUI.ui index 3707a47a..7eb8bd7b 100644 --- a/wmpl/MetSim/GUI.ui +++ b/wmpl/MetSim/GUI.ui @@ -82,8 +82,8 @@ - 40 - 30 + 48 + 12 61 23 @@ -92,8 +92,8 @@ - -25 - 30 + -17 + 12 61 16 @@ -108,7 +108,7 @@ - 113 + 243 30 57 15 @@ -124,7 +124,7 @@ - 173 + 303 30 61 23 @@ -135,7 +135,7 @@ 113 - 60 + 87 57 15 @@ -151,7 +151,7 @@ 173 - 60 + 87 61 23 @@ -160,8 +160,8 @@ - 173 - 90 + 303 + 60 61 23 @@ -170,8 +170,8 @@ - 113 - 90 + 243 + 60 57 15 @@ -186,8 +186,8 @@ - -20 - 60 + 114 + 30 57 15 @@ -202,8 +202,8 @@ - 40 - 60 + 174 + 30 61 23 @@ -212,7 +212,7 @@ - 240 + 370 30 57 15 @@ -226,7 +226,7 @@ 240 - 60 + 90 57 15 @@ -238,8 +238,8 @@ - 240 - 90 + 370 + 60 57 15 @@ -251,8 +251,8 @@ - 110 - 30 + 118 + 12 21 16 @@ -264,8 +264,8 @@ - 110 - 60 + 244 + 30 21 16 @@ -277,8 +277,8 @@ - -25 - 90 + 109 + 60 61 16 @@ -293,8 +293,8 @@ - 110 - 90 + 244 + 60 21 16 @@ -306,8 +306,8 @@ - 40 - 90 + 174 + 60 61 23 @@ -317,7 +317,7 @@ 303 - 30 + 90 61 23 @@ -327,7 +327,7 @@ 243 - 30 + 87 57 15 @@ -343,7 +343,7 @@ 370 - 30 + 90 57 15 @@ -352,6 +352,112 @@ km + + + + 31 + 35 + 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 + + + true + + + + + + 48 + 58 + 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; 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. + + + + + + 18 + 58 + 28 + 16 + + + + rtol + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + 48 + 90 + 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 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. + + + + + + 10 + 90 + 36 + 16 + + + + Grain-release interval along the trail (m). Larger = faster but coarser erosion tail. Default 200 m. + + + grain + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + 120 + 20 + 20 + 91 + + + + Qt::Vertical + + + + + + 114 + 90 + 21 + 16 + + + + m + + @@ -1172,7 +1278,7 @@ 860 - 0 + 4 641 551 diff --git a/wmpl/MetSim/MetSimErosion.py b/wmpl/MetSim/MetSimErosion.py index ea6dda5d..6e0ea5e4 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, adaptiveDP45Step ### DEFINE CONSTANTS @@ -40,9 +40,60 @@ 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. 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 + # 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 + # 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). + # Exposed in the GUI so it can be coarsened for expensive fits (see erosion_release_vref). + 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 + # 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 @@ -350,6 +401,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 +772,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, 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 @@ -780,95 +838,126 @@ 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) + # 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 + + 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( + 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, 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 + const.adaptive_floor_accepts += floor_accepts + if runaway: + const.adaptive_runaway_events += 1 - # 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 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 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 - # Otherwise update the velocity - else: + # Compute new mass + m_new = frag.m + mass_loss_total + + # Compute change in velocity + deceleration_total = decelerationRK4(const.dt, frag.K, frag.m, rho_atm, frag.v) + + # 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 - # Compute g at given height - gv = G0/((1 + frag.h/const.r_earth)**2) + # 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 ### + # ### Add velocity change due to Earth's gravity ### - # # Vertical component of a - # av = -gv - deceleration_total*frag.vv/frag.v + frag.vh*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) + # # 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) + ### 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) + # 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) + # 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 + # 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 + # Track the total drop due to gravity + frag.h_grav_drop_total += h_grav_drop - # Update the velocity - frag.vv -= av*const.dt - frag.vh -= ah*const.dt - frag.v = math.sqrt(frag.vh**2 + frag.vv**2) + # Update the velocity + frag.vv -= av*const.dt + frag.vh -= ah*const.dt + frag.v = math.sqrt(frag.vh**2 + frag.vv**2) - # Only allow the meteoroid to go down, and stop the ablation if it stars going up - if frag.vv > 0: + # Only allow the meteoroid to go down, and stop the ablation if it stars going up + if frag.vv > 0: - frag.vv = 0 + frag.vv = 0 - # Setting the height to zero will stop the ablation during the if catch below - frag.h = 0 + # Setting the height to zero will stop the ablation during the if catch below + frag.h = 0 - # Update length along the track - frag.length += frag.v*const.dt + # Update length along the track + frag.length += frag.v*const.dt - # Update the mass - frag.m = m_new + # Update the mass + frag.m = m_new - # Old way of computing height which did not include the curvature of the Earth - # frag.h = frag.h + frag.vv*const.dt + # 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 + # 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,36 +1133,55 @@ 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. + Uses the enclosing frag/const. """ + 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) - frag_children_all += grain_children + frag_children_all.extend(grain_children) + + eroded_this_tick = False + + 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, 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 - # 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)) + else: - # Record physical parameters at the beginning of erosion for the main fragment - if frag.main: - if const.erosion_beg_vel is None: + # 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 - const.erosion_beg_vel = frag.v - const.erosion_beg_mass = frag.m - const.erosion_beg_dyn_press = dyn_press + # 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 + # 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: @@ -1392,6 +1500,24 @@ 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 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_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) + 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). + 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 +1603,16 @@ 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)) + 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 frag_main = None diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index 481e3624..3c49b9a7 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,708 @@ 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 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 + + +@cython.cdivision(True) +cdef double heightCurvatureC(double h0, double zc, double l, double r_earth): + """ 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 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) + + +@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 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 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 + 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, 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 + + +@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, + 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 + 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". + erosion_release_length: [double] Along-track grain-release interval (m). For an eroding fragment + 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: + 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 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 + 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 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 + 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] + 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 + h_sub_carry = h_sub + + while t < dt_macro: + + # 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 + + # 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. + # 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 + + # 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: + + # 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 + # 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 + 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 + + # 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)) + if fac < facmin: + fac = facmin + h_sub = h_sub*fac + if h_sub < dt_min: + h_sub = dt_min + + # 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: + rho_final = atmDensityPolyC(h_new, dens_co) + else: + rho_final = rho_last + + 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 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 + 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, + 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 + 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 + 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, v_cur, v_cap + 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 + + # 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_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 + 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 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: