From e51b76e1769cc0828806ec5bdfa93e56a86f4951 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Thu, 4 Jan 2024 09:52:01 -0700 Subject: [PATCH 1/2] Prevent RK4 simulation timestep from extending past next scheduled event. Calculate thrust at sample points instead of averaged across timestep --- .../simulation/AbstractSimulationStepper.java | 14 ++-- .../BasicEventSimulationEngine.java | 9 +++ .../simulation/RK4SimulationStepper.java | 77 +++++++------------ 3 files changed, 41 insertions(+), 59 deletions(-) diff --git a/core/src/net/sf/openrocket/simulation/AbstractSimulationStepper.java b/core/src/net/sf/openrocket/simulation/AbstractSimulationStepper.java index 637d5c6e3..58beae1e6 100644 --- a/core/src/net/sf/openrocket/simulation/AbstractSimulationStepper.java +++ b/core/src/net/sf/openrocket/simulation/AbstractSimulationStepper.java @@ -157,20 +157,18 @@ public abstract class AbstractSimulationStepper implements SimulationStepper { /** - * Calculate the average thrust produced by the motors in the current configuration, allowing - * listeners to override. The average is taken between status.time and - * status.time + timestep. + * Calculate the thrust produced by the motors in the current configuration, allowing + * listeners to override. The thrust is calculated at time status.time. *

* TODO: HIGH: This method does not take into account any moments generated by off-center motors. * * @param status the current simulation status. - * @param timestep the time step of the current iteration. * @param acceleration the current (approximate) acceleration * @param atmosphericConditions the current atmospheric conditions * @param stepMotors whether to step the motors forward or work on a clone object - * @return the average thrust during the time step. + * @return the thrust at the specified time */ - protected double calculateAverageThrust(SimulationStatus status, double timestep, + protected double calculateThrust(SimulationStatus status, double acceleration, AtmosphericConditions atmosphericConditions, boolean stepMotors) throws SimulationException { double thrust; @@ -182,11 +180,9 @@ public abstract class AbstractSimulationStepper implements SimulationStepper { } thrust = 0; - final double currentTime = status.getSimulationTime() + timestep; Collection activeMotorList = status.getActiveMotors(); for (MotorClusterState currentMotorState : activeMotorList ) { - thrust += currentMotorState.getAverageThrust( status.getSimulationTime(), currentTime ); - //thrust += currentMotorState.getThrust( currentTime ); + thrust += currentMotorState.getThrust( status.getSimulationTime() ); } // Post-listeners diff --git a/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java b/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java index 8f6c0ee10..b7f46b7bf 100644 --- a/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java +++ b/core/src/net/sf/openrocket/simulation/BasicEventSimulationEngine.java @@ -12,6 +12,7 @@ import net.sf.openrocket.logging.WarningSet; import net.sf.openrocket.l10n.Translator; import net.sf.openrocket.motor.MotorConfiguration; import net.sf.openrocket.motor.MotorConfigurationId; +import net.sf.openrocket.motor.ThrustCurveMotor; import net.sf.openrocket.rocketcomponent.AxialStage; import net.sf.openrocket.rocketcomponent.DeploymentConfiguration; import net.sf.openrocket.rocketcomponent.FlightConfiguration; @@ -412,6 +413,14 @@ public class BasicEventSimulationEngine implements SimulationEngine { continue; } + // Queue an altitude event for every point in the thrust curve to set the RK4 simulation time + // steps + ThrustCurveMotor motor = (ThrustCurveMotor) motorState.getMotor(); + double[] timePoints = motor.getTimePoints(); + for (double point : timePoints) { + addEvent(new FlightEvent(FlightEvent.Type.ALTITUDE, point, event.getSource(), motorState)); + } + // and queue up the burnout for this motor, as well. double duration = motorState.getBurnTime(); double burnout = currentStatus.getSimulationTime() + duration; diff --git a/core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java b/core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java index 0f32bb142..60e281517 100644 --- a/core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java +++ b/core/src/net/sf/openrocket/simulation/RK4SimulationStepper.java @@ -93,26 +93,12 @@ public class RK4SimulationStepper extends AbstractSimulationStepper { RK4SimulationStatus status2; RK4Parameters k1, k2, k3, k4; - - - /* - * Start with previous time step which is used to compute the initial thrust estimate. - * Don't make it longer than maxTimeStep, but at least MIN_TIME_STEP. - */ - store.timestep = MathUtil.max(MathUtil.min(store.timestep, maxTimeStep), MIN_TIME_STEP); - checkNaN(store.timestep); /* * Get the current atmospheric conditions */ calculateFlightConditions(status, store); store.atmosphericConditions = store.flightConditions.getAtmosphericConditions(); - - /* - * Compute the initial thrust estimate. This is used for the first time step computation. - */ - store.thrustForce = calculateAverageThrust(status, store.timestep, store.longitudinalAcceleration, - store.atmosphericConditions, false); /* * Perform RK4 integration. Decide the time step length after the first step. @@ -162,48 +148,38 @@ public class RK4SimulationStepper extends AbstractSimulationStepper { limitingValue = i; } } - + + log.trace("Selected time step " + store.timestep + " (limiting factor " + limitingValue + ")"); + + // If we have a scheduled event coming up before the end of our timestep, truncate step + // else if the time from the end of our timestep to the next scheduled event time is less than + // minTimeStep, stretch it double minTimeStep = status.getSimulationConditions().getTimeStep() / 20; + FlightEvent nextEvent = status.getEventQueue().peek(); + if (nextEvent != null) { + double nextEventTime = nextEvent.getTime(); + if (status.getSimulationTime() + store.timestep > nextEventTime) { + store.timestep = nextEventTime - status.getSimulationTime(); + log.trace("scheduled event at " + nextEventTime + " truncates timestep to " + store.timestep); + } else if ((status.getSimulationTime() + store.timestep < nextEventTime) && + (status.getSimulationTime() + store.timestep + minTimeStep > nextEventTime)) { + store.timestep = nextEventTime - status.getSimulationTime(); + log.trace("Scheduled event at " + nextEventTime + " stretches timestep to " + store.timestep); + } + } + + // If we've wound up with a too-small timestep, increase it avoid numerical instability even at the + // cost of not being *quite* on an event if (store.timestep < minTimeStep) { log.trace("Too small time step " + store.timestep + " (limiting factor " + limitingValue + "), using " + minTimeStep + " instead."); store.timestep = minTimeStep; - } else { - log.trace("Selected time step " + store.timestep + " (limiting factor " + limitingValue + ")"); } + checkNaN(store.timestep); - - /* - * Compute the correct thrust for this time step. If the original thrust estimate differs more - * than 10% from the true value then recompute the RK4 step 1. The 10% error in step 1 is - * diminished by it affecting only 1/6th of the total, so it's an acceptable error. - */ - double thrustEstimate = store.thrustForce; - store.thrustForce = calculateAverageThrust(status, store.timestep, store.longitudinalAcceleration, - store.atmosphericConditions, true); - log.trace("Thrust = " + store.thrustForce); - double thrustDiff = Math.abs(store.thrustForce - thrustEstimate); - // Log if difference over 1%, recompute if over 10% - if (thrustDiff > 0.01 * thrustEstimate) { - if (thrustDiff > 0.1 * thrustEstimate + 0.001) { - log.debug("Thrust estimate differs from correct value by " + - (Math.rint(1000 * (thrustDiff + 0.000001) / thrustEstimate) / 10.0) + "%," + - " estimate=" + thrustEstimate + - " correct=" + store.thrustForce + - " timestep=" + store.timestep + - ", recomputing k1 parameters"); - k1 = computeParameters(status, store); - } else { - log.trace("Thrust estimate differs from correct value by " + - (Math.rint(1000 * (thrustDiff + 0.000001) / thrustEstimate) / 10.0) + "%," + - " estimate=" + thrustEstimate + - " correct=" + store.thrustForce + - " timestep=" + store.timestep + - ", error acceptable"); - } - } - + + //// Second position, k2 = f(t + h/2, y + k1*h/2) status2 = status.clone(); @@ -345,8 +321,9 @@ public class RK4SimulationStepper extends AbstractSimulationStepper { store.dragForce = store.forces.getCDaxial() * dynP * refArea; double fN = store.forces.getCN() * dynP * refArea; double fSide = store.forces.getCside() * dynP * refArea; - - double forceZ = store.thrustForce - store.dragForce; + + store.thrustForce = calculateThrust(status, store.longitudinalAcceleration, store.flightConditions.getAtmosphericConditions(), false); + double forceZ = store.thrustForce - store.dragForce; store.linearAcceleration = new Coordinate(-fN / store.rocketMass.getMass(), -fSide / store.rocketMass.getMass(), From 9e63380f0660157697e27669fc2567993dc770e5 Mon Sep 17 00:00:00 2001 From: JoePfeiffer Date: Thu, 4 Jan 2024 14:02:17 -0700 Subject: [PATCH 2/2] update unit tests --- .../simulation/FlightEventsTest.java | 31 ++++++++++--------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/core/test/net/sf/openrocket/simulation/FlightEventsTest.java b/core/test/net/sf/openrocket/simulation/FlightEventsTest.java index 0e21a9834..5e097d5d8 100644 --- a/core/test/net/sf/openrocket/simulation/FlightEventsTest.java +++ b/core/test/net/sf/openrocket/simulation/FlightEventsTest.java @@ -56,8 +56,8 @@ public class FlightEventsTest extends BaseTestCase { new FlightEvent(FlightEvent.Type.EJECTION_CHARGE, 2.0, stage), new FlightEvent(FlightEvent.Type.RECOVERY_DEVICE_DEPLOYMENT, 2.001, parachute), new FlightEvent(FlightEvent.Type.APOGEE, 2.48338, rocket), - new FlightEvent(FlightEvent.Type.GROUND_HIT, 43.1, null), - new FlightEvent(FlightEvent.Type.SIMULATION_END, 43.1, null) + new FlightEvent(FlightEvent.Type.GROUND_HIT, 43.076, null), + new FlightEvent(FlightEvent.Type.SIMULATION_END, 43.076, null) }; checkEvents(expectedEvents, sim, 0); @@ -99,18 +99,18 @@ public class FlightEventsTest extends BaseTestCase { new FlightEvent(FlightEvent.Type.LAUNCH, 0.0, rocket), new FlightEvent(FlightEvent.Type.IGNITION, 0.0, boosterMotorTubes), new FlightEvent(FlightEvent.Type.IGNITION, 0.0, coreBody), - new FlightEvent(FlightEvent.Type.LIFTOFF, 0.1225, null), - new FlightEvent(FlightEvent.Type.LAUNCHROD, 0.125, null), - new FlightEvent(FlightEvent.Type.APOGEE, 1.867, rocket), + new FlightEvent(FlightEvent.Type.LIFTOFF, 0.1275, null), + new FlightEvent(FlightEvent.Type.LAUNCHROD, 0.130, null), new FlightEvent(FlightEvent.Type.BURNOUT, 2.0, boosterMotorTubes), new FlightEvent(FlightEvent.Type.EJECTION_CHARGE, 2.0, boosterStage), new FlightEvent(FlightEvent.Type.STAGE_SEPARATION, 2.0, boosterStage), new FlightEvent(FlightEvent.Type.BURNOUT, 2.0, coreBody), new FlightEvent(FlightEvent.Type.EJECTION_CHARGE, 2.0, coreStage), new FlightEvent(FlightEvent.Type.STAGE_SEPARATION, 2.0, coreStage), - new FlightEvent(FlightEvent.Type.TUMBLE, 2.38, null), - new FlightEvent(FlightEvent.Type.GROUND_HIT, 13.2, null), - new FlightEvent(FlightEvent.Type.SIMULATION_END, 13.2, null) + new FlightEvent(FlightEvent.Type.APOGEE, 2.12, rocket), + new FlightEvent(FlightEvent.Type.TUMBLE, 2.202, null), + new FlightEvent(FlightEvent.Type.GROUND_HIT, 14.12, null), + new FlightEvent(FlightEvent.Type.SIMULATION_END, 14.12, null) }; break; // Core stage @@ -120,8 +120,10 @@ public class FlightEventsTest extends BaseTestCase { new FlightEvent(FlightEvent.Type.BURNOUT, 2.0, coreBody), new FlightEvent(FlightEvent.Type.EJECTION_CHARGE, 2.0, coreStage), new FlightEvent(FlightEvent.Type.STAGE_SEPARATION, 2.0, coreStage), - new FlightEvent(FlightEvent.Type.GROUND_HIT, 5.9, null), - new FlightEvent(FlightEvent.Type.SIMULATION_END, 5.9, null) + new FlightEvent(FlightEvent.Type.APOGEE, 2.1, rocket), + new FlightEvent(FlightEvent.Type.TUMBLE, 2.15, null), + new FlightEvent(FlightEvent.Type.GROUND_HIT, 7.26, null), + new FlightEvent(FlightEvent.Type.SIMULATION_END, 7.26, null) }; break; // Booster stage @@ -131,9 +133,10 @@ public class FlightEventsTest extends BaseTestCase { new FlightEvent(FlightEvent.Type.BURNOUT, 2.0, boosterMotorTubes), new FlightEvent(FlightEvent.Type.EJECTION_CHARGE, 2.0, boosterStage), new FlightEvent(FlightEvent.Type.STAGE_SEPARATION, 2.0, boosterStage), - new FlightEvent(FlightEvent.Type.TUMBLE, 3.428, null), - new FlightEvent(FlightEvent.Type.GROUND_HIT, 10.3, null), - new FlightEvent(FlightEvent.Type.SIMULATION_END, 10.3, null) + new FlightEvent(FlightEvent.Type.TUMBLE, 2.05, null), + new FlightEvent(FlightEvent.Type.APOGEE, 2.06, rocket), + new FlightEvent(FlightEvent.Type.GROUND_HIT, 12.0, null), + new FlightEvent(FlightEvent.Type.SIMULATION_END, 12.0, null) }; break; default: @@ -147,7 +150,7 @@ public class FlightEventsTest extends BaseTestCase { private void checkEvents(FlightEvent[] expectedEvents, Simulation sim, int branchNo) { FlightEvent[] actualEvents = sim.getSimulatedData().getBranch(branchNo).getEvents().toArray(new FlightEvent[0]); - + // Test event count assertEquals("Branch " + branchNo + " invalid number of events ", expectedEvents.length, actualEvents.length);