summaryrefslogtreecommitdiff
path: root/src/fluid_simulator.adb
blob: e9e86b9fc39ba7d76bd173fbee1c9c8554124eb0 (plain)
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
153
154

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
            if not P.Solid then
                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 if;
        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;