diff options
Diffstat (limited to 'src/fluid_simulator.adb')
| -rw-r--r-- | src/fluid_simulator.adb | 197 |
1 files changed, 197 insertions, 0 deletions
diff --git a/src/fluid_simulator.adb b/src/fluid_simulator.adb new file mode 100644 index 0000000..01d3772 --- /dev/null +++ b/src/fluid_simulator.adb @@ -0,0 +1,197 @@ + +with + + Ada.Numerics.Generic_Complex_Types, + Ada.Containers.Vectors, + Ada.Characters.Latin_1, + Ada.Text_IO; + +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 + -- Variant to clear entire screen + IO.Put (Latin.ESC & "[2J"); + end Clear_Screen; + + procedure Reset_Cursor is + begin + -- ANSI control sequence Cursor Position + -- Parameters to move cursor to top left corner + IO.Put (Latin.ESC & "[1;1H"); + end Reset_Cursor; + + + + type Quantity is digits 18; + + package Fixed is new Ada.Numerics.Generic_Complex_Types (Real => Quantity); + use type Fixed.Complex; + + type Particle is record + Place : Fixed.Complex; + Solid : Boolean; + Density : Quantity := 0.0; + Acceleration : Fixed.Complex := Fixed.Compose_From_Cartesian (0.0, 0.0); + Velocity : Fixed.Complex := Fixed.Compose_From_Cartesian (0.0, 0.0); + end record; + + function Create + (X, Y : Quantity; + Solid : in Boolean) + return Particle is + begin + return (Place => Fixed.Compose_From_Cartesian (X, Y), Solid => Solid, others => <>); + end Create; + + + + -- Constant properties of particles + Particle_Radius : constant Quantity := 2.0; + Particle_Mass : constant Quantity := 1.0; + + -- Constant used in calculating fluid interaction forces + P0 : constant Quantity := 1.5; + + -- Other constant force factors + Gravity_Factor : constant Fixed.Complex := Fixed.Compose_From_Cartesian (0.0, 1.0); + Pressure_Factor : constant Quantity := 4.0; + Viscosity_Factor : constant Quantity := 8.0; + + + + package Particle_Vectors is new Ada.Containers.Vectors + (Index_Type => Positive, + Element_Type => Particle); + + Particles : Particle_Vectors.Vector := Particle_Vectors.Empty_Vector; + + + + function Marching_Squares + (Input : in Particle_Vectors.Vector) + return String + is + Liquid_Chars : String (1 .. 16) := " .,_`/[/']\\-/\#"; + Grid : array (0 .. 80, 0 .. 25) of Boolean := (others => (others => False)); + Output : String (1 .. 2025); + X, Y : Integer; + begin + for P of Input loop + X := Integer (Fixed.Re (P.Place)); + Y := Integer (Fixed.Im (P.Place)); + if X >= 0 and X <= 80 and Y >= 0 and Y <= 25 then + Grid (X, Y) := True; + end if; + end loop; + for J in Integer range 1 .. 25 loop + for I in Integer range 1 .. 80 loop + Output ((J - 1) * 81 + I) := Liquid_Chars + (Boolean'Pos (Grid (I - 1, J)) + 1 + + Boolean'Pos (Grid (I, J)) * 2 + + Boolean'Pos (Grid (I, J - 1)) * 4 + + Boolean'Pos (Grid (I - 1, J - 1)) * 8); + end loop; + Output (J * 81) := Latin.LF; + end loop; + return Output; + end Marching_Squares; + +begin + + declare + 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 + Particles.Append (Create (X, Y, (Input = '#'))); + 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; + + -- Calculate interaction forces + for P of Particles 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; + end loop; + end loop; + + -- Render using marching squares + Clear_Screen; + Reset_Cursor; + IO.Put (Marching_Squares (Particles)); + + -- Update acceleration, velocity, and position + for P of Particles loop + if not P.Solid then + P.Velocity := P.Velocity + P.Acceleration / 10.0; + P.Place := P.Place + P.Velocity; + end if; + end loop; + + -- 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 + then + Particles.Delete (C); + end if; + end loop; + + -- Sleep for a bit + delay 0.012321; + end loop; + +end Fluid_Simulator; + |