From 8d0440926f17567e07948fbf8f29b1a0d82f7411 Mon Sep 17 00:00:00 2001 From: lisyarus Date: Wed, 1 Dec 2021 21:23:25 +0300 Subject: [PATCH] Shallow water steady state solver (wip) --- .../include/psemek/phys/shallow_water.hpp | 141 ++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 libs/phys/include/psemek/phys/shallow_water.hpp diff --git a/libs/phys/include/psemek/phys/shallow_water.hpp b/libs/phys/include/psemek/phys/shallow_water.hpp new file mode 100644 index 00000000..693b7dd7 --- /dev/null +++ b/libs/phys/include/psemek/phys/shallow_water.hpp @@ -0,0 +1,141 @@ +#pragma once + +#include +#include +#include + +namespace psemek::phys +{ + + template + std::pair shallow_water_steady_state_single_step(util::array const & depth, util::array & water_height, + util::array & x_velocity, util::array & y_velocity, + T gravity, std::size_t max_iterations, T max_error) + { + std::size_t const n = depth.width() - 2; + std::size_t const m = depth.height() - 2; + + util::array delta_height({n + 2, m + 2}, 0.f); + util::array delta_x_velocity({n + 1, m + 2}, 0.f); + util::array delta_y_velocity({n + 2, m + 1}, 0.f); + + auto h_dual = [&](auto i, auto j) { return geom::dual{water_height(i, j), {delta_height(i, j), 0.f}}; }; + auto ux_dual = [&](auto i, auto j) { return geom::dual{x_velocity(i, j + 1), {delta_x_velocity(i, j + 1), 0.f}}; }; + auto uy_dual = [&](auto i, auto j) { return geom::dual{y_velocity(i + 1, j), {delta_y_velocity(i + 1, j), 0.f}}; }; + + auto h_diag = [&](auto i, auto j) { return geom::dual{water_height(i, j), {0.f, 1.f}}; }; + auto ux_diag = [&](auto i, auto j) { return geom::dual{x_velocity(i, j + 1), {0.f, 1.f}}; }; + auto uy_diag = [&](auto i, auto j) { return geom::dual{y_velocity(i + 1, j), {0.f, 1.f}}; }; + + T error = 0; + + T omega = T{0.1}; + + std::size_t iteration = 0; + for (; iteration < max_iterations; ++iteration) + { + // Single iteration of Gauss-Seidel + error = 0; + + // Mass conservation equations + for (std::size_t j = 1; j < m + 1; ++j) + { + for (std::size_t i = 1; i < n + 1; ++i) + { + geom::dual value = + + (h_dual(i + 1, j ) + h_diag(i, j)) * ux_dual(i , j - 1) / 2.f + - (h_dual(i - 1, j ) + h_diag(i, j)) * ux_dual(i - 1, j - 1) / 2.f + + (h_dual(i , j + 1) + h_diag(i, j)) * uy_dual(i - 1, j ) / 2.f + - (h_dual(i , j - 1) + h_diag(i, j)) * uy_dual(i - 1, j - 1) / 2.f + ; + +// error += geom::sqr(value.scalar - value.delta[0] - value.delta[1] * delta_height(i, j)); + error += geom::sqr(value.scalar); + + if (value.delta[1] != 0.f) + delta_height(i, j) = geom::lerp(delta_height(i, j), (-value.scalar - value.delta[0]) / value.delta[1], omega); + } + } + + // X-momentum conservation equations + for (std::size_t j = 0; j < m; ++j) + { + for (std::size_t i = 1; i < n; ++i) + { + geom::dual value = + + h_dual(i + 1, j + 1) * (geom::sqr(ux_dual(i + 1, j ) + geom::sqr(ux_diag(i, j)))) / 2.f + - h_dual(i , j + 1) * (geom::sqr(ux_dual(i - 1, j ) + geom::sqr(ux_diag(i, j)))) / 2.f + + + (h_dual(i, j + 1) + h_dual(i + 1, j + 1) + h_dual(i, j + 2) + h_dual(i + 1, j + 2)) + * (ux_diag(i, j) + ux_dual(i, j + 1)) * (uy_dual(i - 1, j + 1) + uy_dual(i, j + 1)) / 16.f + + - (h_dual(i, j + 1) + h_dual(i + 1, j + 1) + h_dual(i, j) + h_dual(i - 1, j)) + * (ux_diag(i, j) + ux_dual(i, j - 1)) * (uy_dual(i - 1, j) + uy_dual(i, j)) / 16.f + + + gravity * (geom::sqr(h_dual(i + 1, j + 1)) - geom::sqr(h_dual(i, j + 1))) / 2.f + + - gravity * (h_dual(i + 1, j + 1) + h_dual(i, j + 1)) * (depth(i + 1, j + 1) - depth(i, j + 1)) / 2.f + ; + + value = 0; + +// error += geom::sqr(value.scalar - value.delta[0] - value.delta[1] * delta_x_velocity(i, j + 1)); + error += geom::sqr(value.scalar); + + if (value.delta[1] != 0.f) + delta_x_velocity(i, j + 1) = geom::lerp(delta_x_velocity(i, j + 1), (value.scalar - value.delta[0]) / value.delta[1], omega); + } + } + + // Y-momentum conservation equations + for (std::size_t j = 1; j < m; ++j) + { + for (std::size_t i = 0; i < n; ++i) + { + geom::dual value = + + h_dual(i + 1, j + 1) * (geom::sqr(uy_dual(i , j + 1) + geom::sqr(uy_diag(i, j)))) / 2.f + - h_dual(i + 1, j ) * (geom::sqr(uy_dual(i , j - 1) + geom::sqr(uy_diag(i, j)))) / 2.f + + + (h_dual(i + 1, j) + h_dual(i + 1, j + 1) + h_dual(i + 2, j) + h_dual(i + 2, j + 1)) + * (uy_diag(i, j) + uy_dual(i + 1, j)) * (ux_dual(i + 1, j - 1) + ux_dual(i + 1, j)) / 16.f + + - (h_dual(i + 1, j) + h_dual(i + 1, j + 1) + h_dual(i, j) + h_dual(i, j - 1)) + * (uy_diag(i, j) + uy_dual(i - 1, j)) * (ux_dual(i, j - 1) + ux_dual(i, j)) / 16.f + + + gravity * (geom::sqr(h_dual(i + 1, j + 1)) - geom::sqr(h_dual(i + 1, j))) / 2.f + + - gravity * (h_dual(i + 1, j + 1) + h_dual(i + 1, j)) * (depth(i + 1, j + 1) - depth(i + 1, j)) / 2.f + ; + + value = 0; + +// error += geom::sqr(value.scalar - value.delta[0] - value.delta[1] * delta_y_velocity(i + 1, j)); + error += geom::sqr(value.scalar); + + if (value.delta[1] != 0.f) + delta_y_velocity(i + 1, j) = geom::lerp(delta_y_velocity(i + 1, j), (-value.scalar - value.delta[0]) / value.delta[1], omega); + } + } + + int temp = 42; + (void)temp; + + (void)max_error; + } + + for (std::size_t j = 0; j < m + 2; ++j) + for (std::size_t i = 0; i < n + 2; ++i) + water_height(i, j) += delta_height(i, j); + + for (std::size_t j = 0; j < m + 2; ++j) + for (std::size_t i = 0; i < n + 1; ++i) + x_velocity(i, j) += delta_x_velocity(i, j); + + for (std::size_t j = 0; j < m + 1; ++j) + for (std::size_t i = 0; i < n + 2; ++i) + y_velocity(i, j) += delta_y_velocity(i, j); + + return {iteration, error}; + } + +}