Use non-physical mode for liquid drag

common/src/comp/fluid_dynamics.rs

@@ -211,6 +211,83 @@ impl Body {
         }
     }
 
+    /// Physically incorrect (but relatively dt-independent) way to calculate drag coefficients for liquids.
+    // TODO: Remove this in favour of `aerodynamic_forces` (see: `integrate_forces`)
+    pub fn drag_coefficient_liquid(
+        &self,
+        rel_flow: &Vel,
+        fluid_density: f32,
+        wings: Option<&Wings>,
+        scale: f32,
+    ) -> Vec3<f32> {
+        let v_sq = rel_flow.0.magnitude_squared();
+        let rel_flow_dir = Dir::new(rel_flow.0 / v_sq.sqrt());
+        0.5 * fluid_density
+            * match wings {
+                Some(&Wings {
+                    aspect_ratio,
+                    planform_area,
+                    ori,
+                }) => {
+                    if aspect_ratio > 25.0 {
+                        tracing::warn!(
+                            "Calculating lift for wings with an aspect ratio of {}. The \
+                              formulas are only valid for aspect ratios below 25.",
+                            aspect_ratio
+                        )
+                    };
+                    let ar = aspect_ratio.min(24.0);
+                    // We have an elliptical wing; proceed to calculate its lift and drag
+
+                    // aoa will be positive when we're pitched up and negative otherwise
+                    let aoa = angle_of_attack(&ori, &rel_flow_dir);
+                    // c_l will be positive when aoa is positive (we have positive lift,
+                    // producing an upward force) and negative otherwise
+                    let c_l = lift_coefficient(ar, aoa);
+
+                    // lift dir will be orthogonal to the local relative flow vector.
+                    // Local relative flow is the resulting vector of (relative) freestream
+                    // flow + downwash (created by the vortices
+                    // of the wing tips)
+                    let lift_dir: Dir = {
+                        // induced angle of attack
+                        let aoa_i = c_l / (PI * ar);
+                        // effective angle of attack; the aoa as seen by aerofoil after
+                        // downwash
+                        let aoa_eff = aoa - aoa_i;
+                        // Angle between chord line and local relative wind is aoa_eff
+                        // radians. Direction of lift is
+                        // perpendicular to local relative wind.
+                        // At positive lift, local relative wind will be below our cord line
+                        // at an angle of aoa_eff. Thus if
+                        // we pitch down by aoa_eff radians then
+                        // our chord line will be colinear with local relative wind vector
+                        // and our up will be the direction
+                        // of lift.
+                        ori.pitched_down(aoa_eff).up()
+                    };
+
+                    // induced drag coefficient (drag due to lift)
+                    let cdi = {
+                        // Oswald's efficiency factor (empirically derived--very magical)
+                        // (this definition should not be used for aspect ratios > 25)
+                        let e = 1.78 * (1.0 - 0.045 * ar.powf(0.68)) - 0.64;
+                        c_l.powi(2) / (PI * e * ar)
+                    };
+
+                    // drag coefficient
+                    let c_d = zero_lift_drag_coefficient() + cdi;
+                    debug_assert!(c_d.is_sign_positive());
+                    debug_assert!(c_l.is_sign_positive() || aoa.is_sign_negative());
+
+                    planform_area * scale.powf(2.0) * (c_l * *lift_dir + c_d * *rel_flow_dir)
+                        + self.parasite_drag(scale).powf(0.75) * *rel_flow_dir
+                },
+
+                _ => self.parasite_drag(scale).powf(0.75) * *rel_flow_dir,
+            }
+    }
+
     /// Parasite drag is the sum of pressure drag and skin friction.
     /// Skin friction is the drag arising from the shear forces between a fluid
     /// and a surface, while pressure drag is due to flow separation. Both are

common/src/states/utils.rs

@@ -220,7 +220,7 @@ impl Body {
                 quadruped_low::Species::Mossdrake => 1.7,
                 _ => 2.0,
             },
-            Body::Ship(ship) if ship.has_water_thrust() => 0.2,
+            Body::Ship(ship) if ship.has_water_thrust() => 0.35,
             Body::Ship(_) => 0.12,
             Body::Arthropod(_) => 3.5,
         }
@@ -236,16 +236,16 @@ impl Body {
             match self {
                 Body::Object(_) => return None,
                 Body::ItemDrop(_) => return None,
-                Body::Ship(ship) if ship.has_water_thrust() => 0.2 * self.mass().0,
+                Body::Ship(ship) if ship.has_water_thrust() => 350.0 * self.mass().0,
                 Body::Ship(_) => return None,
                 Body::BipedLarge(_) => 120.0 * self.mass().0,
-                Body::Golem(_) => 0.5 * self.mass().0,
+                Body::Golem(_) => 100.0 * self.mass().0,
                 Body::BipedSmall(_) => 1000.0 * self.mass().0,
-                Body::BirdMedium(_) => 1500.0 * self.mass().0,
+                Body::BirdMedium(_) => 400.0 * self.mass().0,
-                Body::BirdLarge(_) => 35.0 * self.mass().0,
+                Body::BirdLarge(_) => 400.0 * self.mass().0,
-                Body::FishMedium(_) => 75.0 * self.mass().0,
+                Body::FishMedium(_) => 900.0 * self.mass().0,
-                Body::FishSmall(_) => 120.0 * self.mass().0,
+                Body::FishSmall(_) => 1000.0 * self.mass().0,
-                Body::Dragon(_) => 0.5 * self.mass().0,
+                Body::Dragon(_) => 50.0 * self.mass().0,
                 // Humanoids are a bit different: we try to give them thrusts that result in similar
                 // speeds for gameplay reasons
                 Body::Humanoid(_) => 4_000_000.0 / self.mass().0,
@@ -256,16 +256,16 @@ impl Body {
                     | theropod::Species::Woodraptor
                     | theropod::Species::Dodarock
                     | theropod::Species::Axebeak
-                    | theropod::Species::Yale => 700.0 * self.mass().0,
+                    | theropod::Species::Yale => 500.0 * self.mass().0,
-                    _ => 3.0 * self.mass().0,
+                    _ => 150.0 * self.mass().0,
                 },
-                Body::QuadrupedLow(_) => 14.0 * self.mass().0,
+                Body::QuadrupedLow(_) => 1200.0 * self.mass().0,
                 Body::QuadrupedMedium(body) => match body.species {
-                    quadruped_medium::Species::Mammoth => 75.0 * self.mass().0,
+                    quadruped_medium::Species::Mammoth => 150.0 * self.mass().0,
-                    _ => 500.0 * self.mass().0,
+                    _ => 1000.0 * self.mass().0,
                 },
                 Body::QuadrupedSmall(_) => 1500.0 * self.mass().0,
-                Body::Arthropod(_) => 300.0 * self.mass().0,
+                Body::Arthropod(_) => 500.0 * self.mass().0,
             } * front_profile,
         )
     }
@@ -671,14 +671,15 @@ fn swim_move(
         // Assume that feet/flippers get less efficient as we become less submerged
         let move_z = move_z.min((submersion * 1.5 - 0.5).clamp(0.0, 1.0).powi(2));
 
-        update.vel.0 += Vec3::broadcast(data.dt.0)
+        update.vel.0 += Vec3::new(dir.x, dir.y, move_z)
-            * Vec3::new(dir.x, dir.y, move_z)
                 // TODO: Should probably be normalised, but creates odd discrepancies when treading water
                 // .try_normalized()
                 // .unwrap_or_default()
             * water_accel
             // Gives a good balance between submerged and surface speed
-            * (submersion * 3.0).clamp(0.0, 1.0);
+            * submersion.clamp(0.0, 1.0).sqrt()
+            // Good approximate compensation for dt-dependent effects
+            * data.dt.0 * 0.04;
 
         true
     } else {

common/systems/src/phys.rs

@@ -62,6 +62,21 @@ fn integrate_forces(
     // Aerodynamic/hydrodynamic forces
     if !rel_flow.0.is_approx_zero() {
         debug_assert!(!rel_flow.0.map(|a| a.is_nan()).reduce_or());
+        // HACK: We should really use the latter logic (i.e: `aerodynamic_forces`) for liquids, but it results in
+        // pretty serious dt-dependent problems that are extremely difficult to resolve. This is a compromise: for
+        // liquids only, we calculate drag using an incorrect (but still visually plausible) model that is much more
+        // resistant to differences in dt. Because it's physically incorrect anyway, there are magic coefficients that
+        // exist simply to get us closer to what water 'should' feel like.
+        if fluid.is_liquid() {
+            let fric = body.drag_coefficient_liquid(
+                &rel_flow,
+                fluid_density.0,
+                wings,
+                scale.map_or(1.0, |s| s.0),
+            ).magnitude() * 0.02;
+
+            vel.0 *= (1.0 / (1.0 + fric)).powf(dt.0 * 10.0);
+        } else {
             let impulse = dt.0
                 * body.aerodynamic_forces(
                     &rel_flow,
@@ -84,7 +99,8 @@ fn integrate_forces(
                 } else {
                     vel.0 = new_v;
                 }
-        };
+            }
+        }
         debug_assert!(!vel.0.map(|a| a.is_nan()).reduce_or());
     };
 
@@ -1732,6 +1748,8 @@ fn box_voxel_collision<T: BaseVol<Vox = Block> + ReadVol>(
                 (old_depth + old_pos.z - pos.0.z).max(depth)
             });
 
+            physics_state.ground_vel = ground_vel;
+
             Fluid::Liquid {
                 kind,
                 depth: new_depth,