Merge pull request #2464 from JoePfeiffer/euler-time-step

Refine time step in AbstractEulerStepper
This commit is contained in:
Joe Pfeiffer 2024-03-16 08:19:26 -06:00 committed by GitHub
commit 8cb9ba3e8f
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
5 changed files with 158 additions and 101 deletions

View File

@ -4,6 +4,7 @@ import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import net.sf.openrocket.l10n.Translator;
import net.sf.openrocket.logging.SimulationAbort;
import net.sf.openrocket.models.atmosphere.AtmosphericConditions;
import net.sf.openrocket.rocketcomponent.InstanceMap;
import net.sf.openrocket.rocketcomponent.RecoveryDevice;
@ -38,51 +39,55 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
public void step(SimulationStatus status, double maxTimeStep) throws SimulationException {
// Get the atmospheric conditions
AtmosphericConditions atmosphere = modelAtmosphericConditions(status);
final AtmosphericConditions atmosphere = modelAtmosphericConditions(status);
//// Local wind speed and direction
Coordinate windSpeed = modelWindVelocity(status);
final Coordinate windSpeed = modelWindVelocity(status);
Coordinate airSpeed = status.getRocketVelocity().add(windSpeed);
// Compute drag force
double mach = airSpeed.length() / atmosphere.getMachSpeed();
double dynP = (0.5 * atmosphere.getDensity() * airSpeed.length2());
double dragForce = getCD() * dynP * status.getConfiguration().getReferenceArea();
final double mach = airSpeed.length() / atmosphere.getMachSpeed();
final double CdA = getCD() * status.getConfiguration().getReferenceArea();
final double dragForce = 0.5 * CdA * atmosphere.getDensity() * airSpeed.length2();
double rocketMass = calculateStructureMass(status).getMass();
double motorMass = calculateMotorMass(status).getMass();
final double rocketMass = calculateStructureMass(status).getMass();
final double motorMass = calculateMotorMass(status).getMass();
double mass = rocketMass + motorMass;
final double mass = rocketMass + motorMass;
if (mass < MathUtil.EPSILON) {
throw new SimulationException(trans.get("SimulationStepper.error.totalMassZero"));
status.abortSimulation(SimulationAbort.Cause.ACTIVE_MASS_ZERO);
}
// Compute drag acceleration
Coordinate linearAcceleration;
if (airSpeed.length() > 0.001) {
linearAcceleration = airSpeed.normalize().multiply(-dragForce / mass);
} else {
linearAcceleration = Coordinate.NUL;
}
Coordinate linearAcceleration = airSpeed.normalize().multiply(-dragForce / mass);
// Add effect of gravity
double gravity = modelGravity(status);
final double gravity = modelGravity(status);
linearAcceleration = linearAcceleration.sub(0, 0, gravity);
// Add coriolis acceleration
Coordinate coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation().getCoriolisAcceleration(
final Coordinate coriolisAcceleration = status.getSimulationConditions().getGeodeticComputation().getCoriolisAcceleration(
status.getRocketWorldPosition(), status.getRocketVelocity());
linearAcceleration = linearAcceleration.add(coriolisAcceleration);
// Select tentative time step
double timeStep = RECOVERY_TIME_STEP;
// adjust based on change in acceleration (ie jerk)
final double jerk = Math.abs(linearAcceleration.sub(status.getRocketAcceleration()).multiply(1.0/status.getPreviousTimeStep()).length());
if (jerk > MathUtil.EPSILON) {
timeStep = Math.min(timeStep, 1.0/jerk);
// adjust based on acceleration
final double absAccel = linearAcceleration.length();
if (absAccel > MathUtil.EPSILON) {
timeStep = Math.min(timeStep, 1.0/absAccel);
}
// Honor max step size passed in. If the time to next time step is greater than our minimum
// we'll set our next step to just before it in order to better capture discontinuities in things like chute opening
if (maxTimeStep < timeStep) {
if (maxTimeStep > MIN_TIME_STEP) {
timeStep = maxTimeStep - MIN_TIME_STEP;
} else {
timeStep = maxTimeStep;
}
}
// but don't let it get *too* small
@ -90,33 +95,72 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
log.trace("timeStep is " + timeStep);
// Perform Euler integration
Coordinate newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)).
add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2));
EulerValues newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep);
// If I've hit the ground, recalculate time step and position
if (newPosition.z < 0) {
// Check to see if z or either of its first two derivatives have changed sign and recalculate
// time step to point of change if so
// z -- ground hit
// z' -- apogee
// z'' -- possible oscillation building up in descent rate
// Note that it's virtually impossible for apogee to occur on the same
// step as either ground hit or descent rate inflection, and if we get a ground hit
// any descent rate inflection won't matter
final double a = linearAcceleration.z;
final double v = status.getRocketVelocity().z;
final double z = status.getRocketPosition().z;
double t = timeStep;
if (newVals.pos.z < 0) {
// If I've hit the ground, the new timestep is the solution of
// 1/2 at^2 + vt + z = 0
t = (-v - Math.sqrt(v*v - 2*a*z))/a;
log.trace("ground hit changes timeStep to " + t);
} else if (v * newVals.vel.z < 0) {
// If I've got apogee, the new timestep is the solution of
// v + at = 0
t = Math.abs(v / a);
log.trace("apogee changes timeStep to " + t);
} else {
// Use jerk to estimate accleration at end of time step. Don't really need to redo all the atmospheric
// calculations to get it "right"; this will be close enough for our purposes.
// use chain rule to compute jerk
// dA/dT = dA/dV * dV/dT
final double dFdV = CdA * atmosphere.getDensity() * airSpeed.length();
final Coordinate dAdV = airSpeed.normalize().multiply(dFdV / mass);
final Coordinate jerk = linearAcceleration.multiply(dAdV);
final Coordinate newAcceleration = linearAcceleration.add(jerk.multiply(timeStep));
final double a = linearAcceleration.z;
final double v = status.getRocketVelocity().z;
final double z0 = status.getRocketPosition().z;
// Only do this one if acceleration is appreciably different from 0
if (newAcceleration.z * linearAcceleration.z < -MathUtil.EPSILON) {
// If acceleration oscillation is building up, the new timestep is the solution of
// a + j*t = 0
t = Math.abs(a / jerk.z);
log.trace("oscillation avoidance changes timeStep to " + t);
}
}
// once again, make sure new timestep isn't *too* small
t = Math.max(t, MIN_TIME_STEP);
// The new timestep is the solution of
// 1/2 at^2 + vt + z0 = 0
timeStep = (-v - Math.sqrt(v*v - 2*a*z0))/a;
log.trace("ground hit changes timeStep to " + timeStep);
// recalculate Euler integration for position and velocity if necessary.
if (Math.abs(t - timeStep) > MathUtil.EPSILON) {
timeStep = t;
newPosition = status.getRocketPosition().add(status.getRocketVelocity().multiply(timeStep)).
add(linearAcceleration.multiply(MathUtil.pow2(timeStep) / 2));
if (maxTimeStep - timeStep < MIN_TIME_STEP) {
timeStep = maxTimeStep;
}
// avoid rounding error in new altitude
newPosition = newPosition.setZ(0);
newVals = eulerIntegrate(status.getRocketPosition(), status.getRocketVelocity(), linearAcceleration, timeStep);
// If we just landed chop off rounding error
if (Math.abs(newVals.pos.z) < MathUtil.EPSILON) {
newVals.pos = newVals.pos.setZ(0);
}
}
status.setSimulationTime(status.getSimulationTime() + timeStep);
status.setPreviousTimeStep(timeStep);
status.setRocketPosition(newPosition);
status.setRocketVelocity(status.getRocketVelocity().add(linearAcceleration.multiply(timeStep)));
status.setRocketPosition(newVals.pos);
status.setRocketVelocity(newVals.vel);
status.setRocketAcceleration(linearAcceleration);
// Update the world coordinate
@ -125,64 +169,84 @@ public abstract class AbstractEulerStepper extends AbstractSimulationStepper {
status.setRocketWorldPosition(w);
// Store data
FlightDataBranch data = status.getFlightData();
data.addPoint();
final FlightDataBranch data = status.getFlightData();
// Values looked up or calculated at start of time step
data.setValue(FlightDataType.TYPE_REFERENCE_LENGTH, status.getConfiguration().getReferenceLength());
data.setValue(FlightDataType.TYPE_REFERENCE_AREA, status.getConfiguration().getReferenceArea());
data.setValue(FlightDataType.TYPE_WIND_VELOCITY, windSpeed.length());
data.setValue(FlightDataType.TYPE_AIR_TEMPERATURE, atmosphere.getTemperature());
data.setValue(FlightDataType.TYPE_AIR_PRESSURE, atmosphere.getPressure());
data.setValue(FlightDataType.TYPE_SPEED_OF_SOUND, atmosphere.getMachSpeed());
data.setValue(FlightDataType.TYPE_MACH_NUMBER, mach);
if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.FLAT) {
data.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, coriolisAcceleration.length());
}
data.setValue(FlightDataType.TYPE_GRAVITY, gravity);
data.setValue(FlightDataType.TYPE_DRAG_COEFF, getCD());
data.setValue(FlightDataType.TYPE_PRESSURE_DRAG_COEFF, getCD());
data.setValue(FlightDataType.TYPE_FRICTION_DRAG_COEFF, 0);
data.setValue(FlightDataType.TYPE_BASE_DRAG_COEFF, 0);
data.setValue(FlightDataType.TYPE_AXIAL_DRAG_COEFF, getCD());
data.setValue(FlightDataType.TYPE_THRUST_FORCE, 0);
data.setValue(FlightDataType.TYPE_DRAG_FORCE, dragForce);
data.setValue(FlightDataType.TYPE_MASS, mass);
data.setValue(FlightDataType.TYPE_MOTOR_MASS, motorMass);
data.setValue(FlightDataType.TYPE_THRUST_WEIGHT_RATIO, 0);
data.setValue(FlightDataType.TYPE_ACCELERATION_XY,
MathUtil.hypot(linearAcceleration.x, linearAcceleration.y));
data.setValue(FlightDataType.TYPE_ACCELERATION_Z, linearAcceleration.z);
data.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, linearAcceleration.length());
data.setValue(FlightDataType.TYPE_TIME_STEP, timeStep);
// Values calculated on this step
data.addPoint();
data.setValue(FlightDataType.TYPE_TIME, status.getSimulationTime());
data.setValue(FlightDataType.TYPE_ALTITUDE, status.getRocketPosition().z);
data.setValue(FlightDataType.TYPE_POSITION_X, status.getRocketPosition().x);
data.setValue(FlightDataType.TYPE_POSITION_Y, status.getRocketPosition().y);
airSpeed = status.getRocketVelocity().add(windSpeed);
data.setValue(FlightDataType.TYPE_POSITION_XY,
MathUtil.hypot(status.getRocketPosition().x, status.getRocketPosition().y));
data.setValue(FlightDataType.TYPE_POSITION_DIRECTION,
Math.atan2(status.getRocketPosition().y, status.getRocketPosition().x));
data.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad());
data.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad());
data.setValue(FlightDataType.TYPE_VELOCITY_XY,
MathUtil.hypot(status.getRocketVelocity().x, status.getRocketVelocity().y));
data.setValue(FlightDataType.TYPE_ACCELERATION_XY,
MathUtil.hypot(linearAcceleration.x, linearAcceleration.y));
data.setValue(FlightDataType.TYPE_VELOCITY_Z, status.getRocketVelocity().z);
data.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, airSpeed.length());
data.setValue(FlightDataType.TYPE_ACCELERATION_TOTAL, linearAcceleration.length());
double Re = airSpeed.length() *
airSpeed = status.getRocketVelocity().add(windSpeed);
final double Re = airSpeed.length() *
status.getConfiguration().getLengthAerodynamic() /
atmosphere.getKinematicViscosity();
data.setValue(FlightDataType.TYPE_REYNOLDS_NUMBER, Re);
data.setValue(FlightDataType.TYPE_LATITUDE, status.getRocketWorldPosition().getLatitudeRad());
data.setValue(FlightDataType.TYPE_LONGITUDE, status.getRocketWorldPosition().getLongitudeRad());
data.setValue(FlightDataType.TYPE_GRAVITY, gravity);
if (status.getSimulationConditions().getGeodeticComputation() != GeodeticComputationStrategy.FLAT) {
data.setValue(FlightDataType.TYPE_CORIOLIS_ACCELERATION, coriolisAcceleration.length());
}
data.setValue(FlightDataType.TYPE_VELOCITY_Z, status.getRocketVelocity().z);
data.setValue(FlightDataType.TYPE_ACCELERATION_Z, linearAcceleration.z);
data.setValue(FlightDataType.TYPE_VELOCITY_TOTAL, airSpeed.length());
data.setValue(FlightDataType.TYPE_MACH_NUMBER, mach);
data.setValue(FlightDataType.TYPE_MASS, mass);
data.setValue(FlightDataType.TYPE_MOTOR_MASS, motorMass);
data.setValue(FlightDataType.TYPE_THRUST_FORCE, 0);
data.setValue(FlightDataType.TYPE_DRAG_FORCE, dragForce);
data.setValue(FlightDataType.TYPE_WIND_VELOCITY, windSpeed.length());
data.setValue(FlightDataType.TYPE_AIR_TEMPERATURE, atmosphere.getTemperature());
data.setValue(FlightDataType.TYPE_AIR_PRESSURE, atmosphere.getPressure());
data.setValue(FlightDataType.TYPE_SPEED_OF_SOUND, atmosphere.getMachSpeed());
data.setValue(FlightDataType.TYPE_TIME_STEP, timeStep);
data.setValue(FlightDataType.TYPE_COMPUTATION_TIME,
(System.nanoTime() - status.getSimulationStartWallTime()) / 1000000000.0);
log.trace("time " + data.getLast(FlightDataType.TYPE_TIME) + ", altitude " + data.getLast(FlightDataType.TYPE_ALTITUDE) + ", velocity " + data.getLast(FlightDataType.TYPE_VELOCITY_Z));
}
private static class EulerValues {
/** linear velocity */
public Coordinate vel;
/** position */
public Coordinate pos;
}
private EulerValues eulerIntegrate (Coordinate pos, Coordinate v, Coordinate a, double timeStep) {
EulerValues result = new EulerValues();
result.vel = v.add(a.multiply(timeStep));
result.pos = pos.add(v.multiply(timeStep)).add(a.multiply(MathUtil.pow2(timeStep) / 2.0));
return result;
}
}

View File

@ -706,7 +706,6 @@ public class BasicEventSimulationEngine implements SimulationEngine {
double d = 0;
boolean b = false;
d += currentStatus.getSimulationTime();
d += currentStatus.getPreviousTimeStep();
b |= currentStatus.getRocketPosition().isNaN();
b |= currentStatus.getRocketVelocity().isNaN();
b |= currentStatus.getRocketOrientationQuaternion().isNaN();
@ -716,7 +715,6 @@ public class BasicEventSimulationEngine implements SimulationEngine {
if (Double.isNaN(d) || b) {
log.error("Simulation resulted in NaN value:" +
" simulationTime=" + currentStatus.getSimulationTime() +
" previousTimeStep=" + currentStatus.getPreviousTimeStep() +
" rocketPosition=" + currentStatus.getRocketPosition() +
" rocketVelocity=" + currentStatus.getRocketVelocity() +
" rocketOrientationQuaternion=" + currentStatus.getRocketOrientationQuaternion() +

View File

@ -44,8 +44,6 @@ public class SimulationStatus implements Monitorable {
private double time;
private double previousTimeStep;
private Coordinate position;
private WorldCoordinate worldPosition;
private Coordinate velocity;
@ -105,7 +103,6 @@ public class SimulationStatus implements Monitorable {
this.configuration = configuration;
this.time = 0;
this.previousTimeStep = this.simulationConditions.getTimeStep();
this.position = this.simulationConditions.getLaunchPosition();
this.velocity = this.simulationConditions.getLaunchVelocity();
this.worldPosition = this.simulationConditions.getLaunchSite();
@ -178,7 +175,6 @@ public class SimulationStatus implements Monitorable {
// FlightData is not cloned.
this.flightData = orig.flightData;
this.time = orig.time;
this.previousTimeStep = orig.previousTimeStep;
this.position = orig.position;
this.acceleration = orig.acceleration;
this.worldPosition = orig.worldPosition;
@ -267,17 +263,6 @@ public class SimulationStatus implements Monitorable {
return flightData;
}
public double getPreviousTimeStep() {
return previousTimeStep;
}
public void setPreviousTimeStep(double previousTimeStep) {
this.previousTimeStep = previousTimeStep;
this.modID++;
}
public void setRocketPosition(Coordinate position) {
this.position = position;
this.modID++;

View File

@ -190,6 +190,16 @@ public final class Coordinate implements Cloneable, Serializable {
return new Coordinate(this.x * m, this.y * m, this.z * m, this.weight * m);
}
/**
* Multiply the <code>Coordinate</code> with another <Coordinate> component-by-component
*
* @param other the other Coordinate
* @return The product.
*/
public Coordinate multiply(Coordinate other) {
return new Coordinate(this.x * other.x, this.y * other.y, this.z * other.z, this.weight * other.weight);
}
/**
* Dot product of two Coordinates, taken as vectors. Equal to
* x1*x2+y1*y2+z1*z2

View File

@ -55,9 +55,9 @@ public class FlightEventsTest extends BaseTestCase {
new FlightEvent(FlightEvent.Type.BURNOUT, 2.0, motorMountTube),
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.076, null),
new FlightEvent(FlightEvent.Type.SIMULATION_END, 43.076, null)
new FlightEvent(FlightEvent.Type.APOGEE, 2.48, rocket),
new FlightEvent(FlightEvent.Type.GROUND_HIT, 42.97, null),
new FlightEvent(FlightEvent.Type.SIMULATION_END, 42.97, null)
};
checkEvents(expectedEvents, sim, 0);
@ -138,9 +138,9 @@ public class FlightEventsTest extends BaseTestCase {
new FlightEvent(FlightEvent.Type.EJECTION_CHARGE, 2.01, centerBooster),
new FlightEvent(FlightEvent.Type.STAGE_SEPARATION, 2.01, centerBooster),
new FlightEvent(FlightEvent.Type.TUMBLE, 2.85, null),
new FlightEvent(FlightEvent.Type.APOGEE, 1200, rocket),
new FlightEvent(FlightEvent.Type.GROUND_HIT, 1200, null),
new FlightEvent(FlightEvent.Type.SIMULATION_END, 1200, null)
new FlightEvent(FlightEvent.Type.APOGEE, 3.78, rocket),
new FlightEvent(FlightEvent.Type.GROUND_HIT, 9.0, null),
new FlightEvent(FlightEvent.Type.SIMULATION_END, 9.0, null)
};
break;