Skip to content

Commit

Permalink
First attempts at fixing the simulation. Pretty bugged right now
Browse files Browse the repository at this point in the history
  • Loading branch information
@Zang3th
Zang3th committed May 1, 2024
1 parent 3226673 commit e89be3c
Show file tree
Hide file tree
Showing 10 changed files with 185 additions and 122 deletions.
10 changes: 5 additions & 5 deletions Apps/Liquefied/LiquefiedApp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,11 +158,11 @@ namespace Liq
if(_physicsTimer->CheckElapsedAndReset())
{
//Add a horizontal turbine (initial velocity)
_fluidSimulator->AddHorizonalTurbine(1, 48, (float)Engine::LiquiefiedParams::turbinePower);
_fluidSimulator->AddHorizonalTurbine(1, 49, (float)Engine::LiquiefiedParams::turbinePower);
_fluidSimulator->AddHorizonalTurbine(1, 50, (float)Engine::LiquiefiedParams::turbinePower);
_fluidSimulator->AddHorizonalTurbine(1, 51, (float)Engine::LiquiefiedParams::turbinePower);
_fluidSimulator->AddHorizonalTurbine(1, 52, (float)Engine::LiquiefiedParams::turbinePower);
_fluidSimulator->AddHorizonalTurbine(1, 48, 0.1f);
_fluidSimulator->AddHorizonalTurbine(1, 49, 0.1f);
_fluidSimulator->AddHorizonalTurbine(1, 50, 0.1f);
_fluidSimulator->AddHorizonalTurbine(1, 51, 0.1f);
_fluidSimulator->AddHorizonalTurbine(1, 52, 0.1f);

//Run simulation timestep and visualize result
_fluidSimulator->TimeStep();
Expand Down
8 changes: 4 additions & 4 deletions Apps/Liquefied/LiquefiedInterface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,11 +38,11 @@ namespace Liq
CenterText("Numerical value monitoring");
ImGui::Separator();
ImGui::NewLine();
for(auto const& entry : Engine::Monitor::values)
for(auto const& entry : Engine::Monitoring::values)
{
ImGui::Text("\t%s (min): %5.3f", entry.first, entry.second.min);
ImGui::Text("\t%s (max): %5.3f", entry.first, entry.second.max);
ImGui::Text("\t%s (val): %5.3f", entry.first, entry.second.val);
ImGui::Text("\t%s (min): %5.3f at (%d, %d)", entry.first, entry.second.min, entry.second.x, entry.second.y);
ImGui::Text("\t%s (max): %5.3f at (%d, %d)", entry.first, entry.second.max, entry.second.x, entry.second.y);
ImGui::Text("\t%s (val): %5.3f at (%d, %d)", entry.first, entry.second.val, entry.second.x, entry.second.y);
ImGui::NewLine();
}
ImGui::Separator();
Expand Down
6 changes: 3 additions & 3 deletions Engine/Application/GlobalParams.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,9 +164,9 @@ namespace Engine
inline static constexpr uint32 SIMULATION_WIDTH = 150;
inline static constexpr uint32 SIMULATION_HEIGHT = 100;
inline static constexpr uint32 LIQUID_NUM_CELLS = SIMULATION_WIDTH * SIMULATION_HEIGHT;
inline static constexpr uint32 GAUSS_SEIDEL_ITERATIONS = 20;
inline static constexpr float GAUSS_SEIDEL_OVERRELAXATION = 1.9f;
inline static uint32 turbinePower = 50;
inline static constexpr uint32 GAUSS_SEIDEL_ITERATIONS = 25;
inline static constexpr float GAUSS_SEIDEL_OVERRELAXATION = 1.95f;
inline static uint32 turbinePower = 1;
inline static bool visualizeSmoke = true;
inline static bool scientificColorScheme = false;
inline static bool pauseSimulation = true;
Expand Down
15 changes: 9 additions & 6 deletions Engine/Core/Monitor.cpp → Engine/Core/Monitoring.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
#include "Monitor.hpp"
#include "Monitoring.hpp"

namespace Engine
{
// ----- Public -----

void Monitor::MinMaxAvg(const char* name, const float value)
void Monitoring::MinMaxAvgAt(const char* name, const float value, const uint32 x, const uint32 y)
{
//Add new value if it's not already getting monitored
if(values.find(name) == values.end())
{
values[name] = {value, value, value};
values[name] = {value, value, value, x, y};
return;
}

Expand All @@ -23,15 +23,18 @@ namespace Engine
else if(valueStruct.max < value)
valueStruct.max = value;

//Add value
//Add value and position
valueStruct.val = value;
valueStruct.x = x;
valueStruct.y = y;

//Assign new values
values[name] = valueStruct;
}

void Monitor::Reset()
void Monitoring::Reset()
{
values.clear();
for(auto& val : values)
val.second = {0.0f, 0.0f, 0.0f, UINT32_MAX, UINT32_MAX};
}
}
6 changes: 3 additions & 3 deletions Engine/Core/Monitor.hpp → Engine/Core/Monitoring.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,12 @@

namespace Engine
{
class Monitor
class Monitoring
{
public:
Monitor() = delete;
Monitoring() = delete;

static void MinMaxAvg(const char* name, float value);
static void MinMaxAvgAt(const char* name, float value, uint32 x, uint32 y);
static void Reset();

inline static std::map<const char*, MinMaxAvg_t> values = std::map<const char*, MinMaxAvg_t>();
Expand Down
2 changes: 2 additions & 0 deletions Engine/Core/Types.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,5 +45,7 @@ namespace Engine
float min;
float max;
float val;
uint32 x;
uint32 y;
} MinMaxAvg_t;
}
2 changes: 1 addition & 1 deletion Engine/Engine.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
#include "App.hpp"
#include "Utility.hpp"
#include "Timer.hpp"
#include "Monitor.hpp"
#include "Monitoring.hpp"

// ----- Resources and renderables -----
#include "ResourceManager.hpp"
Expand Down
132 changes: 92 additions & 40 deletions Engine/Physics/FluidSimulation/FluidSimulator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,18 @@ namespace Engine
{
// ----- Private -----

float FluidSimulator::ForwardEuler(const float dt, const float dx, const float vel_field, const float q, const float q_right, const float q_left)
float FluidSimulator::ForwardEuler(const float dt, const float u, const float q, const float q_right, const float q_left) const
{
return q - (dt * vel_field * ((q_right - q_left) / (2 * dx)));
/* The update step looks something like this: *
* q[i] = q[i] - dt * u[i] * (q[i+1] - q[i-1]) / (2 * dx); *
* - q := The quantity we want to advect. *
* - q_right := q[i+1] - The right/upper neighbor. *
* - q_left := q[i-1] - The left/lower neighbor . *
* - dt := The timestep size. *
* - u := The velocity (horizonal or vertical). *
* - dx := The cell size (1 in our simple case). */

return q - (dt * u * ((q_right - q_left) / (2 * _grid.dx)));
}

float FluidSimulator::BackwardEuler()
Expand All @@ -24,18 +33,38 @@ namespace Engine
for(uint32 y = 0; y < _grid.height; y++)
{
//Check for border cell
if(_grid.s_At(x, y) != 0.0f)
if(_grid.s_At(x, y) != 0.0f && _grid.s_At(x, y-1) != 0.0f)
{
_grid.v_At(x, y) += GRAVITY * dt;
}
}
}
}

/* There shouldn't be a need for this function to exist, if Project() would work correctly. *
* Furthermore, this function doesn't fix the problem. */
void FluidSimulator::CorrectForces()
{
for(uint32 x = 0; x < _grid.width; x++)
{
for(uint32 y = 0; y < _grid.height; y++)
{
if((x == 0) || (y == 0) || (x == _grid.width-1) || (y == _grid.height-1))
{
_grid.u_At(x, y) = 0.0f;
_grid.v_At(x, y) = 0.0f;
}
}
}
}

/* Projection calculates and applies the right amount of pressure to make the *
* velocity field divergence-free and also enforces solid wall boundary conditions. */
void FluidSimulator::Project(const float dt)
{
//ToDo: Project() is currently bugged. It does not correctly enforce solid wall boundary conditions.
//ToDo: Fix it. Divergence bugged at (148, 98).

//Simple solution using Gauss-Seidel
for(uint32 it = 0; it < ITER; it++)
{
Expand All @@ -44,34 +73,44 @@ namespace Engine
for(uint32 y = 1; y < _grid.height-1; y++)
{
//Skip border cells
if(_grid.s_At(x, y) == 0)
if(_grid.s_At(x, y) == 0.0f)
{
continue;
}

//Get the amount of fluid that is entering or leaving the cell
float divergence = _grid.u_At(x+1, y) - _grid.u_At(x, y) +
_grid.v_At(x, y+1) - _grid.v_At(x, y);

Monitor::MinMaxAvg("Divergence", divergence);

//Apply overrelaxation to speed up convergence
divergence *= OVERRELAX;

//Get the amount of border cells in the area
const float rightNeighbor = _grid.s_At(x+1, y);
const float leftNeighbor = _grid.s_At(x-1, y);
const float upperNeigbor = _grid.s_At(x, y+1);
const float lowerNeighbor = _grid.s_At(x, y-1);

// Sum them up to later divide the divergence by the correct amount
//Sum them up to divide the divergence by the correct amount
const float s_sum = rightNeighbor + leftNeighbor + upperNeigbor + lowerNeighbor;

//Skip if there are only solid cells
if(s_sum == 0.0f)
{
continue;
}

//Get the amount of fluid that is entering or leaving the cell
float divergence = _grid.u_At(x+1, y) - _grid.u_At(x, y) + _grid.v_At(x, y+1) - _grid.v_At(x, y);

//Calculate pressure
float pressure = divergence / s_sum;
Monitoring::MinMaxAvgAt("Pressure", pressure, x, y);

//Apply overrelaxation to speed up convergence
pressure *= OVERRELAX;

//Push all velocities out by the same amout to force incompressibility
_grid.u_At(x, y) += divergence * (leftNeighbor / s_sum);
_grid.u_At(x+1, y) -= divergence * (rightNeighbor / s_sum);
_grid.v_At(x, y) += divergence * (lowerNeighbor / s_sum);
_grid.v_At(x, y+1) -= divergence * (upperNeigbor / s_sum);
_grid.u_At(x, y) += pressure * leftNeighbor;
_grid.u_At(x+1, y) -= pressure * rightNeighbor;
_grid.v_At(x, y) += pressure * lowerNeighbor;
_grid.v_At(x, y+1) -= pressure * upperNeigbor;

//Get the amount of fluid that is entering or leaving the cell
divergence = _grid.u_At(x+1, y) - _grid.u_At(x, y) + _grid.v_At(x, y+1) - _grid.v_At(x, y);
}
}
}
Expand All @@ -80,33 +119,41 @@ namespace Engine
/* Advection moves the quantity, or molecules, or particles, along the velocity field. */
void FluidSimulator::AdvectVelocity(const float dt)
{
/* Iterate over grid and update all velocities according to our numerical integrator. *
* For now we use Forward-Euler (unstable!) for the time derivative and a central spatial derivative. *
* *
* Our Update-Algorithm will look something like this: *
* q[i] = q[i] - delta_t * u[i] * (q[i+1] - q[i-1]) / (2 * delta_x); *
* - q := The quantity we want to advect. *
* - delta_t := The timestep size. *
* - u := The horizontal velocity (we need to advect the vertical velocity v as well). *
* - delta_x := The cell size (1 in our simple case). */

/* Iterate over grid and update all velocities according to the chosen numerical integrator. */
for(uint32 x = 1; x < _grid.width-1; x++)
{
for(uint32 y = 1; y < _grid.height-1; y++)
{
//ToDo: Add more extensive debugging capabilities, like a window to scroll through all the values at runtime.

//Monitor the divergence to make sure that the fluid is incompressible
const float divergence = _grid.u_At(x+1, y) - _grid.u_At(x, y) + _grid.v_At(x, y+1) - _grid.v_At(x, y);
Monitoring::MinMaxAvgAt("Div", divergence, x, y);
Logger::Print("Divergence: " + std::to_string(divergence) + " at (" + std::to_string(x) + ", " + std::to_string(y) + ")");

//Skip border cells
if(_grid.s_At(x, y) == 0)
{
continue;
}

//u-component (horizontal advection)
_grid.u_temp_At(x, y) = ForwardEuler(dt, _grid.deltaX, _grid.u_Avg_At(x, y), _grid.u_At(x, y), _grid.u_At(x+1, y), _grid.u_At(x-1, y));
Monitor::MinMaxAvg("u-component", _grid.u_temp_At(x, y));
if(LiquiefiedParams::integratorChoice == Integrator::ForwardEuler)
{
//ToDo: Fix Project() first, after that look at the unstability of Forward-Euler.

//v-component (vertical advection)
_grid.v_temp_At(x, y) = ForwardEuler(dt, _grid.deltaY, _grid.v_Avg_At(x, y), _grid.v_At(x, y), _grid.v_At(x, y+1), _grid.v_At(x, y-1));
Monitor::MinMaxAvg("v-component", _grid.v_temp_At(x, y));
/*const float newU = ForwardEuler(dt, _grid.u_Avg_At(x, y), _grid.u_At(x, y), _grid.u_At(x+1, y), _grid.u_At(x-1, y));
const float newV = ForwardEuler(dt, _grid.v_Avg_At(x, y), _grid.v_At(x, y), _grid.v_At(x, y+1), _grid.v_At(x, y-1));
Monitoring::MinMaxAvgAt("u-component", newU, x, y);
Monitoring::MinMaxAvgAt("v-component", newV, x, y);
_grid.u_temp_At(x, y) = newU; //u-component (horizontal advection)
_grid.v_temp_At(x, y) = newV; //v-component (vertical advection)*/
}
else if(LiquiefiedParams::integratorChoice == Integrator::BackwardEuler)
{

}
}
}

Expand All @@ -127,13 +174,17 @@ namespace Engine
continue;
}

//u-component (horizontal advection)
const float u_result = ForwardEuler(dt, _grid.deltaX, _grid.u_Avg_At(x, y), _grid.smoke_At(x, y), _grid.smoke_At(x+1, y), _grid.smoke_At(x-1, y));
if(LiquiefiedParams::integratorChoice == Integrator::ForwardEuler)
{
const float uResult = ForwardEuler(dt, _grid.u_Avg_At(x, y), _grid.smoke_At(x, y), _grid.smoke_At(x+1, y), _grid.smoke_At(x-1, y));
const float vResult = ForwardEuler(dt, _grid.v_Avg_At(x, y), _grid.smoke_At(x, y),_grid.smoke_At(x, y+1), _grid.smoke_At(x, y-1));

//v-component (vertical advection)
const float v_result = ForwardEuler(dt, _grid.deltaY, _grid.v_Avg_At(x, y), _grid.smoke_At(x, y),_grid.smoke_At(x, y+1), _grid.smoke_At(x, y-1));
_grid.smoke_temp_At(x, y) = (uResult + vResult) / 2.0f;
}
else if(LiquiefiedParams::integratorChoice == Integrator::BackwardEuler)
{

_grid.smoke_temp_At(x, y) = (u_result + v_result) / 2.0f;
}
}
}

Expand Down Expand Up @@ -168,6 +219,7 @@ namespace Engine
{
float dt = (float)Window::GetDeltaTime_sec();
AddForces(dt);
CorrectForces();
Project(dt);
AdvectVelocity(dt);
AdvectSmoke(dt);
Expand All @@ -176,7 +228,7 @@ namespace Engine
void FluidSimulator::Reset()
{
Init();
Monitor::Reset();
Monitoring::Reset();
}

void FluidSimulator::AddHorizonalTurbine(const uint32 x, const uint32 y, const float power)
Expand Down
Loading

0 comments on commit e89be3c

Please sign in to comment.