1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
with
Datatypes,
ANSI_Terminal,
Ada.Characters.Latin_1,
Ada.Text_IO;
use
Datatypes;
use type
Datatypes.Plane.Complex;
procedure Fluid_Simulator is
package ANSI renames ANSI_Terminal;
package Latin renames Ada.Characters.Latin_1;
package IO renames Ada.Text_IO;
Particles : Particle_Vector := Particle_Vectors.Empty_Vector;
-- 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 Plane.Complex := Plane.Compose_From_Cartesian (0.0, 1.0);
Pressure_Factor : constant Quantity := 4.0;
Viscosity_Factor : constant Quantity := 8.0;
procedure Read_Input
(Store : out Particle_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 + 2.0;
else
if Input > Latin.Space and Input < Latin.DEL then
Store.Append (Create (X, Y, (Input = '#')));
Store.Append (Create (X, Y + 1.0, (Input = '#')));
end if;
X := X + 1.0;
end if;
end loop;
end Read_Input;
procedure Calculate_Density
(Store : in out Particle_Vector)
is
Rij, W : Quantity;
begin
for P of Store loop
P.Density := (if P.Solid then 9.0 else 0.0);
for Q of Store loop
Rij := Plane.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;
end loop;
end loop;
end Calculate_Density;
procedure Calculate_Interaction
(Store : in out Particle_Vector)
is
Displacement, Pressure, Viscosity : Plane.Complex;
Rij : Quantity;
begin
for P of Store loop
P.Acceleration := Gravity_Factor;
for Q of Store loop
Displacement := P.Place - Q.Place;
Rij := Plane.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 +
Plane.Compose_From_Cartesian (1.0 - Rij / Particle_Radius) /
P.Density * (Pressure - Viscosity);
end if;
end loop;
end loop;
end Calculate_Interaction;
procedure Update_Position
(Store : in out Particle_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;
procedure Cull_Outside_Bounds
(Store : in out Particle_Vector;
Threshold : in Quantity) is
begin
for C in reverse Store.First_Index .. Store.Last_Index loop
if Plane.Re (Store (C).Place) < 1.0 - Threshold or
Plane.Re (Store (C).Place) > 80.0 + Threshold or
Plane.Im (Store (C).Place) < 1.0 - Threshold or
Plane.Im (Store (C).Place) > 50.0 + Threshold
then
Store.Delete (C);
end if;
end loop;
end Cull_Outside_Bounds;
begin
Read_Input (Particles);
loop
Calculate_Density (Particles);
ANSI.Clear_Screen;
ANSI.Reset_Cursor;
IO.Put (ANSI.Marching_Squares (Particles));
Calculate_Interaction (Particles);
Update_Position (Particles);
Cull_Outside_Bounds (Particles, 50.0);
delay 0.012321;
end loop;
end Fluid_Simulator;
|
