summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/fluid_simulator.adb155
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;