diff options
| -rw-r--r-- | src/fluid_simulator.adb | 155 |
1 files changed, 87 insertions, 68 deletions
diff --git a/src/fluid_simulator.adb b/src/fluid_simulator.adb index 01d3772..9d9186b 100644 --- a/src/fluid_simulator.adb +++ b/src/fluid_simulator.adb @@ -11,6 +11,8 @@ procedure Fluid_Simulator is package Latin renames Ada.Characters.Latin_1; package IO renames Ada.Text_IO; + + procedure Clear_Screen is begin -- ANSI control sequence Erase in Display @@ -72,6 +74,28 @@ procedure Fluid_Simulator is + procedure Read_Input + (Store : out Particle_Vectors.Vector) + is + Input : Character; + X, Y : Quantity := 1.0; + begin + while not IO.End_Of_File (IO.Standard_Input) loop + IO.Get_Immediate (Input); + if Input = Latin.LF then + X := 1.0; + Y := Y + 1.0; + else + if Input > Latin.Space and Input < Latin.DEL then + Store.Append (Create (X, Y, (Input = '#'))); + end if; + X := X + 1.0; + end if; + end loop; + end Read_Input; + + + function Marching_Squares (Input : in Particle_Vectors.Vector) return String @@ -101,95 +125,90 @@ procedure Fluid_Simulator is return Output; end Marching_Squares; -begin - declare - Input : Character; - X, Y : Quantity := 1.0; + + procedure Calculate_Density + (Store : in out Particle_Vectors.Vector) + is + Rij, W : Quantity; begin - while not IO.End_Of_File (IO.Standard_Input) loop - IO.Get_Immediate (Input); - if Input = Latin.LF then - X := 1.0; - Y := Y + 1.0; - else - if Input > Latin.Space and Input < Latin.DEL then - Particles.Append (Create (X, Y, (Input = '#'))); + for P of Store loop + P.Density := (if P.Solid then 9.0 else 0.0); + for Q of Store loop + Rij := Fixed.Modulus (P.Place - Q.Place); + W := (Rij / Particle_Radius - 1.0) ** 2; + if Rij < Particle_Radius then + P.Density := P.Density + Particle_Mass * W; end if; - X := X + 1.0; - end if; - end loop; - end; - - loop - -- Calculate density - for P of Particles loop - if P.Solid then - P.Density := 9.0; - else - P.Density := 0.0; - end if; - for Q of Particles loop - declare - Rij : Quantity := Fixed.Modulus (P.Place - Q.Place); - W : Quantity := (Rij / Particle_Radius - 1.0) ** 2; - begin - if Rij < Particle_Radius then - P.Density := P.Density + Particle_Mass * W; - end if; - end; end loop; end loop; + end Calculate_Density; + + - -- Calculate interaction forces - for P of Particles loop + procedure Calculate_Interaction + (Store : in out Particle_Vectors.Vector) + is + Displacement, Pressure, Viscosity : Fixed.Complex; + Rij : Quantity; + begin + for P of Store loop P.Acceleration := Gravity_Factor; - for Q of Particles loop - declare - Displacement : Fixed.Complex := P.Place - Q.Place; - Rij : Quantity := Fixed.Modulus (Displacement); - begin - if Rij < Particle_Radius then - declare - Pressure : Fixed.Complex := - (P.Density + Q.Density - 2.0 * P0) * Pressure_Factor * Displacement; - Viscosity : Fixed.Complex := - (P.Velocity - Q.Velocity) * Viscosity_Factor; - begin - P.Acceleration := P.Acceleration + - Fixed.Compose_From_Cartesian (1.0 - Rij / Particle_Radius) / - P.Density * (Pressure - Viscosity); - end; - end if; - end; + for Q of Store loop + Displacement := P.Place - Q.Place; + Rij := Fixed.Modulus (Displacement); + if Rij < Particle_Radius then + Pressure := (P.Density + Q.Density - 2.0 * P0) * Pressure_Factor * Displacement; + Viscosity := (P.Velocity - Q.Velocity) * Viscosity_Factor; + P.Acceleration := P.Acceleration + + Fixed.Compose_From_Cartesian (1.0 - Rij / Particle_Radius) / + P.Density * (Pressure - Viscosity); + end if; end loop; end loop; + end Calculate_Interaction; + - -- Render using marching squares - Clear_Screen; - Reset_Cursor; - IO.Put (Marching_Squares (Particles)); - -- Update acceleration, velocity, and position - for P of Particles loop + procedure Update_Position + (Store : in out Particle_Vectors.Vector) is + begin + for P of Store loop if not P.Solid then P.Velocity := P.Velocity + P.Acceleration / 10.0; P.Place := P.Place + P.Velocity; end if; end loop; + end Update_Position; + + - -- Cull particles going out of bounds before any overflows happen - for C in reverse Particles.First_Index .. Particles.Last_Index loop - if Fixed.Re (Particles (C).Place) < -50.0 or - Fixed.Re (Particles (C).Place) > 130.0 or - Fixed.Im (Particles (C).Place) < -50.0 or - Fixed.Im (Particles (C).Place) > 75.0 + procedure Cull_Outside_Bounds + (Store : in out Particle_Vectors.Vector; + Threshold : in Quantity) is + begin + for C in reverse Store.First_Index .. Store.Last_Index loop + if Fixed.Re (Store (C).Place) < 0.0 - Threshold or + Fixed.Re (Store (C).Place) > 80.0 + Threshold or + Fixed.Im (Store (C).Place) < 0.0 - Threshold or + Fixed.Im (Store (C).Place) > 25.0 + Threshold then - Particles.Delete (C); + Store.Delete (C); end if; end loop; + end Cull_Outside_Bounds; + +begin - -- Sleep for a bit + Read_Input (Particles); + loop + Clear_Screen; + Reset_Cursor; + IO.Put (Marching_Squares (Particles)); + Calculate_Density (Particles); + Calculate_Interaction (Particles); + Update_Position (Particles); + Cull_Outside_Bounds (Particles, 50.0); delay 0.012321; end loop; |