diff --git a/CMakeLists.txt b/CMakeLists.txt index c195f74..c59e135 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,7 +31,7 @@ ENDIF() MARK_AS_ADVANCED(BUILD_COVERAGE) -SET(CMAKE_CXX_STANDARD 14) +SET(CMAKE_CXX_STANDARD 17) SET(CMAKE_CXX_STANDARD_REQUIRED ON) IF (BUILD_TESTS) diff --git a/README.md b/README.md index 5383399..c4aa92c 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ #### Dependencies * [CMake] -* Compiler that supports C++14 +* Compiler that supports C++17 ```sh git clone https://github.com/acrlakshman/gradient-augmented-levelset-cuda --recursive diff --git a/applications/advection/main.cc b/applications/advection/main.cc index ddbe361..98e7df4 100644 --- a/applications/advection/main.cc +++ b/applications/advection/main.cc @@ -29,6 +29,7 @@ // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. /////////////////////////////////////////////////////////////////////////////// +#include #include #include @@ -42,6 +43,13 @@ #include "gals/utilities/grid.h" #include "gals/utilities/vec3.h" +const int dim = 2; +using T = double; +using TV = GALS::CPU::Vec3; +using T_GRID = GALS::CPU::Grid; + +static T compute_dt(const T cfl_max, const T_GRID &grid, const GALS::CPU::Array &velocity); + int main(int argc, char **argv) { std::cout << "Inside applications/advection" << std::endl; @@ -60,11 +68,6 @@ int main(int argc, char **argv) } } - const int dim = 2; - using T = double; - using TV = GALS::CPU::Vec3; - using T_GRID = GALS::CPU::Grid; - GALS::INPUT_FIELDS::InputFields input_fields; GALS::INPUT_PARSER::InputParser input_parser; @@ -118,8 +121,11 @@ int main(int argc, char **argv) auto &psi_prev = levelset.psiPrev(); auto &phi_mixed_derivatives_prev = levelset.phiMixedDerivativesPrev(); GALS::ANALYTICAL_FIELDS::Levelset analytical_levelset(grid, levelset_inputs); + auto iii = 0; + std::cout << "here " << ++iii << std::endl; analytical_levelset.compute(positions, levelset); + std::cout << "here " << ++iii << std::endl; // Build time data. const auto &time_inputs = *(input_fields.m_time); @@ -128,39 +134,107 @@ int main(int argc, char **argv) T t_end = time_inputs.end; T dt = time_inputs.dt; bool is_dt_fixed = std::strcmp(time_inputs.constant_dt.c_str(), "NO"); + T cfl_max = time_inputs.cfl_max; + T next_write_time = time_inputs.write_interval; T sim_time = static_cast(0); if (!is_dt_fixed) { dt = grid.dX().min() * static_cast(0.5); } + std::cout << "here " << ++iii << std::endl; + + // Compute gradient of levelset field. + GALS::CPU::Gradient>::compute(phi, psi); // Time loop + int write_count = 1; bool run_sim = true; while (run_sim) { - sim_time += dt; - // Compute velocity and its gradient at current time. analytical_velocity.compute(positions, sim_time, levelset_velocity); + if (!is_dt_fixed) dt = compute_dt(cfl_max, grid, velocity_field); + std::cout << "=======================================\n" << std::endl; + std::cout << "Time step = " << dt << std::endl; + sim_time += dt; + + iii = 1000; + std::cout << "here " << ++iii << std::endl; + // Compute gradient of levelset field. GALS::CPU::Gradient>::compute(phi, psi); phi_prev = phi; psi_prev = psi; + // std::cout << "phi_prev = " << phi_prev << std::endl; + // std::cout << "psi_prev = " << psi_prev << std::endl; + std::cout << "here " << ++iii << std::endl; // Compute levelset mixed derivatives for `_prev` fields. - // TODO (lakshman) - // levelset.computeMixedDerivatives(psi_prev, phi_mixed_derivatives_prev); + levelset.computeMixedDerivatives(psi_prev, phi_mixed_derivatives_prev); + // std::cout << "phi_mixed_derivatives_prev = " << phi_mixed_derivatives_prev << std::endl; + std::cout << "here " << ++iii << std::endl; // Advect levelset. GALS::CPU::Temporal>::compute( dt, levelset_velocity, levelset); + std::cout << "here " << ++iii << std::endl; + + if (sim_time >= next_write_time) { + auto write_dir = file_utils.getRootDirectory() + std::to_string(write_count); + file_utils.createDirectory(write_dir); + + // Write velocity to a file. + file_utils.write(std::string(write_dir + "/velocity"), velocity_field); + + // Write levelset to a file. + file_utils.write(std::string(write_dir + "/phi"), levelset.phi()); + + ++write_count; + next_write_time += time_inputs.write_interval; + } if (GALS::is_equal(sim_time, t_end) || sim_time > t_end) run_sim = false; } - // Write velocity to a file. - file_utils.createDirectory(file_utils.getRootDirectory()); - file_utils.write(std::string(file_utils.getRootDirectory() + "velocity"), velocity_field); - return 0; } + +T compute_dt(const T cfl_max, const T_GRID &grid, const GALS::CPU::Array &velocity) +{ + T dt = 0.; + + const auto dx = grid.dX(); + const auto dx_min = dx.min(); + + // accessing grid details + const auto mask = grid.getMask(); + const int pad = grid.getPadding(); + const auto num_cells = grid.numCells(); + + // defining working+ghost domain extent + int i_min = 0; + int j_min = 0; + int k_min = 0; + int i_max = num_cells[0]; + int j_max = num_cells[1]; + int k_max = num_cells[2]; + + T max_u_mag = static_cast(0.); + + for (int i = i_min; i < i_max; ++i) + for (int j = j_min; j < j_max; ++j) + for (int k = k_min; k < k_max; ++k) { + // std::cout << "velocity(" << i << "," << j << "," << k << "): " << velocity(i, j, k) << std::endl; + auto vel = velocity(i, j, k).mag(); + max_u_mag = std::max(max_u_mag, vel); + } + T one_by_max_u_mag = static_cast(1.) / max_u_mag; + + std::cout << "one_by_max_u_mag = " << one_by_max_u_mag << std::endl; + std::cout << "\tcfl_max = " << cfl_max << std::endl; + std::cout << "\tdx = " << dx << std::endl; + std::cout << "\tdx_min = " << dx_min << std::endl; + dt = cfl_max * dx_min * one_by_max_u_mag; + + return dt; +} diff --git a/docs/pages/README.md b/docs/pages/README.md index f6887fd..d4ace37 100644 --- a/docs/pages/README.md +++ b/docs/pages/README.md @@ -16,7 +16,7 @@ Implementation of Gradient Augmented Levelset method in CPU and GPU. #### Dependencies * [CMake] -* Compiler that supports C++14 +* Compiler that supports C++17 ```sh git clone https://github.com/acrlakshman/gradient-augmented-levelset-cuda --recursive diff --git a/include/gals/cpu/interpolation/hermite.h b/include/gals/cpu/interpolation/hermite.h index cf637c2..eda88ce 100644 --- a/include/gals/cpu/interpolation/hermite.h +++ b/include/gals/cpu/interpolation/hermite.h @@ -178,7 +178,7 @@ class Hermite> */ CPU::InterpolatedFields> interpolate( const GALS::CPU::Grid::value_type, GALS::CPU::Grid::dim> &grid, - const typename GALS::CPU::Grid::position_type &x_interp, + const GALS::CPU::Vec3 &node_id, const typename GALS::CPU::Grid::position_type &x_interp, const CPU::Levelset, T> &levelset, const bool use_gradient_limiting = false); /*! Interpolate scalar field. @@ -247,7 +247,7 @@ class Hermite> */ CPU::InterpolatedFields> interpolate( const GALS::CPU::Grid::value_type, GALS::CPU::Grid::dim> &grid, - const typename GALS::CPU::Grid::position_type &x_interp, + const GALS::CPU::Vec3 &node_id, const typename GALS::CPU::Grid::position_type &x_interp, const CPU::Levelset, T> &levelset, const bool use_gradient_limiting = false); /*! Interpolate scalar field. @@ -316,7 +316,7 @@ class Hermite> */ CPU::InterpolatedFields> interpolate( const GALS::CPU::Grid::value_type, GALS::CPU::Grid::dim> &grid, - const typename GALS::CPU::Grid::position_type &x_interp, + const GALS::CPU::Vec3 &node_id, const typename GALS::CPU::Grid::position_type &x_interp, const CPU::Levelset, T> &levelset, const bool use_gradient_limiting = false); /*! Interpolate scalar field. diff --git a/include/gals/cpu/levelset.h b/include/gals/cpu/levelset.h index 795ce7c..15f2582 100644 --- a/include/gals/cpu/levelset.h +++ b/include/gals/cpu/levelset.h @@ -84,6 +84,19 @@ class Levelset */ ~Levelset(); + /*! Compute mixed derivatives. + * + * Mixed derivatives of the levelset is computed using second order gradient approximation on + * the levelset gradient. + * + * \param psi gradient of levelset. + * \param phi_mixed_derivatives mixed derivatives of the levelset. + * + * \return Void. + */ + void computeMixedDerivatives(const Array>& psi, + Array>& phi_mixed_derivatives); + /*! Return grid. * * \return grid. diff --git a/include/gals/input-fields/time.h b/include/gals/input-fields/time.h index 931fd1a..01f71ce 100644 --- a/include/gals/input-fields/time.h +++ b/include/gals/input-fields/time.h @@ -38,7 +38,7 @@ namespace GALS namespace INPUT_FIELDS { struct Time { - double start, end, dt; + double start, end, dt, cfl_max, write_interval; std::string constant_dt; }; diff --git a/include/gals/utilities/grid.h b/include/gals/utilities/grid.h index 7972fd8..c87a2f8 100644 --- a/include/gals/utilities/grid.h +++ b/include/gals/utilities/grid.h @@ -238,6 +238,12 @@ class Grid */ void generate(T x_min, T x_max, T y_min, T y_max, T z_min, T z_max); + /*! Returns true if index in the domain. + * + * \param node_id node index of type NodeId. + */ + bool isIndexInDomain(const Vec3 node_id, const bool use_padding = false) const; + private: int m_dimension, m_nx, m_ny, m_nz, m_pad; size_t m_total_cells; diff --git a/src/cpu/grid.cc b/src/cpu/grid.cc index 00a7f95..dba6470 100644 --- a/src/cpu/grid.cc +++ b/src/cpu/grid.cc @@ -215,6 +215,21 @@ void GALS::CPU::Grid::generate(T x_min, T x_max, T y_min, T y_max, T z_m for (int d = 0; d < DIM; ++d) m_total_cells *= this->numCells()[d]; } +template +bool GALS::CPU::Grid::isIndexInDomain(const Vec3 node_id, const bool use_padding) const +{ + int pad = use_padding ? m_pad : 0; + + int i_min = -pad * m_mask[0], j_min = -pad * m_mask[1], k_min = -pad * m_mask[2]; + int i_max = m_nx + pad * m_mask[0], j_max = m_ny + pad * m_mask[1], k_max = m_nz + pad * m_mask[2]; + + if (node_id[0] < i_min || node_id[1] >= i_max) return false; + if (m_dimension > 1 && (node_id[1] < j_min || node_id[1] >= j_max)) return false; + if (m_dimension > 2 && (node_id[2] < k_min || node_id[2] >= k_max)) return false; + + return true; +} + template class GALS::CPU::Grid; template class GALS::CPU::Grid; template class GALS::CPU::Grid; diff --git a/src/cpu/input-parser/time.cc b/src/cpu/input-parser/time.cc index 458a58b..b0e764d 100644 --- a/src/cpu/input-parser/time.cc +++ b/src/cpu/input-parser/time.cc @@ -45,6 +45,8 @@ void GALS::INPUT_PARSER::Time::parse(const YAML::Node &field, GALS::INPUT_FIELDS input_fields.m_time->end = field["end"].as(); input_fields.m_time->dt = field["dt"].as(); input_fields.m_time->constant_dt = field["constant_dt"].as(); + input_fields.m_time->cfl_max = field["cfl_max"].as(); + input_fields.m_time->write_interval = field["write_interval"].as(); } void GALS::INPUT_PARSER::Time::operator()(const YAML::Node &field, GALS::INPUT_FIELDS::InputFields *p_input_fields) diff --git a/src/cpu/interpolation/hermite-1d.cc b/src/cpu/interpolation/hermite-1d.cc index 056c91c..ffcbb3d 100644 --- a/src/cpu/interpolation/hermite-1d.cc +++ b/src/cpu/interpolation/hermite-1d.cc @@ -36,7 +36,7 @@ template GALS::CPU::InterpolatedFields> GALS::INTERPOLATION::Hermite>::interpolate( const GALS::CPU::Grid::value_type, GALS::CPU::Grid::dim> &grid, - const typename GALS::CPU::Grid::position_type &x_interp, + const GALS::CPU::Vec3 &node_id, const typename GALS::CPU::Grid::position_type &x_interp, const GALS::CPU::Levelset, T> &levelset, const bool use_gradient_limiting) { GALS::CPU::InterpolatedFields> hermite_fields; @@ -96,7 +96,8 @@ void GALS::INTERPOLATION::Hermite>::compute( for (int i = 0; i < num_cells_interp[0]; ++i) for (int j = 0; j < num_cells_interp[1]; ++j) for (int k = 0; k < num_cells_interp[2]; ++k) { - const auto &hermite_fields = this->interpolate(grid, x_interp(i, j, k), levelset); + GALS::CPU::Vec3 node_id(i, j, k); + const auto &hermite_fields = this->interpolate(grid, node_id, x_interp(i, j, k), levelset); levelset.phiInterpPrev()(i, j, k) = hermite_fields.phi_interpolated; levelset.psiInterpPrev()(i, j, k) = hermite_fields.psi_interpolated; diff --git a/src/cpu/interpolation/hermite-2d.cc b/src/cpu/interpolation/hermite-2d.cc index 7555f92..ec9c716 100644 --- a/src/cpu/interpolation/hermite-2d.cc +++ b/src/cpu/interpolation/hermite-2d.cc @@ -36,7 +36,7 @@ template GALS::CPU::InterpolatedFields> GALS::INTERPOLATION::Hermite>::interpolate( const GALS::CPU::Grid::value_type, GALS::CPU::Grid::dim> &grid, - const typename GALS::CPU::Grid::position_type &x_interp, + const GALS::CPU::Vec3 &node_id, const typename GALS::CPU::Grid::position_type &x_interp, const GALS::CPU::Levelset, T> &levelset, const bool use_gradient_limiting) { GALS::CPU::InterpolatedFields> hermite_fields; @@ -50,12 +50,23 @@ GALS::CPU::InterpolatedFields> GALS::INTERPOLATION::Hermite base_i_jp1 = GALS::CPU::Vec3(base_i_j[0], base_i_j[1] + axis_vectors(1, 1), 0); const GALS::CPU::Vec3 base_ip1_jp1 = GALS::CPU::Vec3(base_i_j[0] + axis_vectors(0, 0), base_i_j[1] + axis_vectors(1, 1), 0); + // std::cout << "x_interp = " << x_interp << std::endl; + // std::cout << "\t\t\tbase i j = " << base_i_j << std::endl; + // std::cout << "here" << std::endl; + + hermite_fields.phi_interpolated = levelset.phiPrev()(node_id); + hermite_fields.psi_interpolated = levelset.psiPrev()(node_id); + + if (!grid.isIndexInDomain(base_i_j) || !grid.isIndexInDomain(base_ip1_j) || !grid.isIndexInDomain(base_i_jp1) || + !grid.isIndexInDomain(base_ip1_jp1)) + return hermite_fields; const typename T_GRID::position_type x_base = grid(base_i_j); const auto &dx = grid.dX(); const auto &one_over_dx = grid.oneOverDX(); GALS::CPU::Vec3 eta = (x_interp - x_base) * one_over_dx; + // std::cout << "here" << std::endl; const ControlPoints &control_points_bottom = GALS::INTERPOLATION::get_control_points( levelset.phiPrev()(base_i_j), levelset.psiPrev()(base_i_j)[axis_x], levelset.phiPrev()(base_ip1_j), @@ -69,20 +80,27 @@ GALS::CPU::InterpolatedFields> GALS::INTERPOLATION::Hermite &control_points_right = GALS::INTERPOLATION::get_control_points( levelset.phiPrev()(base_ip1_j), levelset.psiPrev()(base_ip1_j)[axis_y], levelset.phiPrev()(base_ip1_jp1), levelset.psiPrev()(base_ip1_jp1)[axis_y], dx[axis_y], use_gradient_limiting); + // std::cout << "here" << std::endl; T c_1 = levelset.psiPrev()(base_i_j)[axis_y] + dx[axis_x] * one_third * levelset.phiMixedDerivativesPrev()(base_i_j)[0]; + // std::cout << "c_1" << std::endl; T c_2 = levelset.psiPrev()(base_ip1_j)[axis_x] + dx[axis_y] * one_third * levelset.phiMixedDerivativesPrev()(base_ip1_j)[0]; + // std::cout << "c_2" << std::endl; T c_2_prime = levelset.psiPrev()(base_ip1_jp1)[axis_x] - dx[axis_y] * one_third * levelset.phiMixedDerivativesPrev()(base_ip1_jp1)[0]; + // std::cout << "c_2p" << std::endl; T c_3 = levelset.psiPrev()(base_i_jp1)[axis_y] + dx[axis_x] * one_third * levelset.phiMixedDerivativesPrev()(base_i_jp1)[0]; + // std::cout << "c_3" << std::endl; + // std::cout << "here" << std::endl; T bx[] = {B0(eta[0]), B1(eta[0]), B2(eta[0]), B3(eta[0])}; T bx_prime[] = {B0_Prime(eta[0]), B1_Prime(eta[0]), B2_Prime(eta[0]), B3_Prime(eta[0])}; T by[] = {B0(eta[1]), B1(eta[1]), B2(eta[1]), B3(eta[1])}; T by_prime[] = {B0_Prime(eta[1]), B1_Prime(eta[1]), B2_Prime(eta[1]), B3_Prime(eta[1])}; + // std::cout << "here" << std::endl; T control_points_all[] = {control_points_left.c_30, control_points_left.c_21, @@ -100,12 +118,14 @@ GALS::CPU::InterpolatedFields> GALS::INTERPOLATION::Hermite> GALS::INTERPOLATION::Hermite> GALS::INTERPOLATION::Hermite>::compute( const GALS::CPU::Vec3 dx = grid.dX(); const auto &axis_vectors = GALS::CPU::Grid::axis_vectors; + std::cout << "Num cells: " << num_cells_interp << std::endl; for (int i = 0; i < num_cells_interp[0]; ++i) for (int j = 0; j < num_cells_interp[1]; ++j) for (int k = 0; k < num_cells_interp[2]; ++k) { - const auto &hermite_fields = this->interpolate(grid, x_interp(i, j, k), levelset); + GALS::CPU::Vec3 node_id(i, j, k); + const auto &hermite_fields = this->interpolate(grid, node_id, x_interp(i, j, k), levelset); levelset.phiInterpPrev()(i, j, k) = hermite_fields.phi_interpolated; levelset.psiInterpPrev()(i, j, k) = hermite_fields.psi_interpolated; + + // std::cout << "\tgrid.X(" << i << "," << j << "," << k << "): " << grid(i,j,k) << std::endl; + // std::cout << "\t\tlevelset.phi(" << i << "," << j << "," << k << "): " << levelset.phiPrev()(i,j,k) << + // std::endl; std::cout << "\t\tlevelset.psi(" << i << "," << j << "," << k << "): " << levelset.psiPrev()(i,j,k) + // << std::endl; std::cout << "\t\tx_interp(" << i << "," << j << "," << k << "): " << x_interp(i,j,k) << + // std::endl; std::cout << "\t\tphi_interp(" << i << "," << j << "," << k << "): " << + // levelset.phiInterpPrev()(i,j,k) << std::endl; std::cout << "\t\tpsi_interp(" << i << "," << j << "," << k << + // "): " << levelset.psiInterpPrev()(i,j,k) << std::endl; } } diff --git a/src/cpu/interpolation/hermite-3d.cc b/src/cpu/interpolation/hermite-3d.cc index 2acd806..3963597 100644 --- a/src/cpu/interpolation/hermite-3d.cc +++ b/src/cpu/interpolation/hermite-3d.cc @@ -36,7 +36,7 @@ template GALS::CPU::InterpolatedFields> GALS::INTERPOLATION::Hermite>::interpolate( const GALS::CPU::Grid::value_type, GALS::CPU::Grid::dim>& grid, - const typename GALS::CPU::Grid::position_type& x_interp, + const GALS::CPU::Vec3& node_id, const typename GALS::CPU::Grid::position_type& x_interp, const GALS::CPU::Levelset, T>& levelset, const bool use_gradient_limiting) { GALS::CPU::InterpolatedFields> hermite_fields; @@ -393,7 +393,8 @@ void GALS::INTERPOLATION::Hermite>::compute( for (int i = 0; i < num_cells_interp[0]; ++i) for (int j = 0; j < num_cells_interp[1]; ++j) for (int k = 0; k < num_cells_interp[2]; ++k) { - const auto& hermite_fields = this->interpolate(grid, x_interp(i, j, k), levelset); + GALS::CPU::Vec3 node_id(i, j, k); + const auto& hermite_fields = this->interpolate(grid, node_id, x_interp(i, j, k), levelset); levelset.phiInterpPrev()(i, j, k) = hermite_fields.phi_interpolated; levelset.psiInterpPrev()(i, j, k) = hermite_fields.psi_interpolated; diff --git a/src/cpu/levelset.cc b/src/cpu/levelset.cc index dba0755..32d6c8e 100644 --- a/src/cpu/levelset.cc +++ b/src/cpu/levelset.cc @@ -30,6 +30,7 @@ /////////////////////////////////////////////////////////////////////////////// #include "gals/cpu/levelset.h" +#include "gals/cpu/gradient.h" #include @@ -52,6 +53,74 @@ GALS::CPU::Levelset::~Levelset() { } +template +void GALS::CPU::Levelset::computeMixedDerivatives( + const Array>& psi, Array>& phi_mixed_derivatives) +{ + // For now gradient computations are hardcoded in this function. + const Vec3 num_cells = psi.numCells(); + const T_GRID& grid = psi.grid(); + const Vec3 dx = grid.dX(); + const auto one_over_dx = grid.oneOverDX(); + const auto& axis_vectors = GALS::CPU::Grid::axis_vectors; + + if constexpr (T_GRID::dim == 2) { + int axis_x = 0, axis_y = 1; + + // Computing \frac{\partial}{\partial y} \left( \frac{\partial}{\partial x} \right). + for (int i = 0; i < num_cells[0]; ++i) + for (int j = 0; j < num_cells[1]; ++j) + for (int k = 0; k < num_cells[2]; ++k) { + phi_mixed_derivatives(i, j, k)[0] = + (psi(i + axis_vectors(axis_y, 0), j + axis_vectors(axis_y, 1), k + axis_vectors(axis_y, 2))[axis_x] - + psi(i - axis_vectors(axis_y, 0), j - axis_vectors(axis_y, 1), k - axis_vectors(axis_y, 2))[axis_x]) * + one_over_dx[axis_y] * static_cast(0.5); + } + } else if constexpr (T_GRID::dim == 3) { + int axis_x = 0, axis_y = 1, axis_z = 2; + + // Computing \frac{\partial}{\partial x} \left( \frac{\partial}{\partial y} \right). + for (int i = 0; i < num_cells[0]; ++i) + for (int j = 0; j < num_cells[1]; ++j) + for (int k = 0; k < num_cells[2]; ++k) { + phi_mixed_derivatives(i, j, k)[0] = + (psi(i + axis_vectors(axis_x, 0), j + axis_vectors(axis_x, 1), k + axis_vectors(axis_x, 2))[axis_y] - + psi(i - axis_vectors(axis_x, 0), j - axis_vectors(axis_x, 1), k - axis_vectors(axis_x, 2))[axis_y]) * + one_over_dx[axis_x] * static_cast(0.5); + } + + // Computing \frac{\partial}{\partial y} \left( \frac{\partial}{\partial z} \right). + for (int i = 0; i < num_cells[0]; ++i) + for (int j = 0; j < num_cells[1]; ++j) + for (int k = 0; k < num_cells[2]; ++k) { + phi_mixed_derivatives(i, j, k)[1] = + (psi(i + axis_vectors(axis_y, 0), j + axis_vectors(axis_y, 1), k + axis_vectors(axis_y, 2))[axis_z] - + psi(i - axis_vectors(axis_y, 0), j - axis_vectors(axis_y, 1), k - axis_vectors(axis_y, 2))[axis_z]) * + one_over_dx[axis_y] * static_cast(0.5); + } + + // Computing \frac{\partial}{\partial z} \left( \frac{\partial}{\partial x} \right). + for (int i = 0; i < num_cells[0]; ++i) + for (int j = 0; j < num_cells[1]; ++j) + for (int k = 0; k < num_cells[2]; ++k) { + phi_mixed_derivatives(i, j, k)[2] = + (psi(i + axis_vectors(axis_z, 0), j + axis_vectors(axis_z, 1), k + axis_vectors(axis_z, 2))[axis_x] - + psi(i - axis_vectors(axis_z, 0), j - axis_vectors(axis_z, 1), k - axis_vectors(axis_z, 2))[axis_x]) * + one_over_dx[axis_z] * static_cast(0.5); + } + + // Computing \frac{\partial}{\partial y} \left( \frac{\partial}{\partial z} \right). + for (int i = 0; i < num_cells[0]; ++i) + for (int j = 0; j < num_cells[1]; ++j) + for (int k = 0; k < num_cells[2]; ++k) { + // phi_mixed_derivatives(i, j, k)[1] = + //(psi(i + axis_vectors(axis_y, 0), j + axis_vectors(axis_y, 1), k + axis_vectors(axis_y, 2))[axis_z] - + // psi(i - axis_vectors(axis_y, 0), j - axis_vectors(axis_y, 1), k - axis_vectors(axis_y, 2))[axis_z]) * + // one_over_dx[axis_y] * static_cast(0.5); + } + } +} + template void GALS::CPU::Levelset::print() { diff --git a/src/cpu/temporal-schemes/semi-lagrangian/euler.cc b/src/cpu/temporal-schemes/semi-lagrangian/euler.cc index e17ad39..a763c09 100644 --- a/src/cpu/temporal-schemes/semi-lagrangian/euler.cc +++ b/src/cpu/temporal-schemes/semi-lagrangian/euler.cc @@ -56,16 +56,25 @@ void GALS::TEMPORAL_SCHEMES::SEMI_LAGRANGIAN::Euler::compute( const GALS::CPU::Vec3 num_cells = grid.numCells(); auto &phi = levelset.phi(); auto &psi = levelset.psi(); + // std::cout << "phi = " << phi << std::endl; // Compute x_root `root of the characteristic`. GALS::CPU::Array> x_root(grid); for (int i = 0; i < num_cells[0]; ++i) for (int j = 0; j < num_cells[1]; ++j) - for (int k = 0; k < num_cells[2]; ++k) x_root(i, j, k) = grid(i, j, k) - velocity(i, j, k) * dt; + for (int k = 0; k < num_cells[2]; ++k) { + x_root(i, j, k) = grid(i, j, k) - velocity(i, j, k) * dt; + // std::cout << "grid(" << i << "," << j << "," << k << "): " << grid(i, j, k) << std::endl; + // std::cout << "\tvelocity(" << i << "," << j << "," << k << "): " << velocity(i,j,k) << std::endl; + // std::cout << "\tdt = " << dt << std::endl; + // std::cout << "\tx_root(" << i << "," << j << "," << k << "): " << x_root(i,j,k) << std::endl; + } + // std::cout << "x_root = " << x_root << std::endl; // Compute phi_interp_prev, psi_interp_prev at x_root. GALS::CPU::Interpolate>::compute(x_root, levelset); + // std::cout << "after interp" << std::endl; // Compute x_root_grad. GALS::CPU::Array> x_root_grad(grid); @@ -75,9 +84,12 @@ void GALS::TEMPORAL_SCHEMES::SEMI_LAGRANGIAN::Euler::compute( for (int i = 0; i < num_cells[0]; ++i) for (int j = 0; j < num_cells[1]; ++j) for (int k = 0; k < num_cells[2]; ++k) x_root_grad(i, j, k) = identity_mat - velocity_grad(i, j, k) * dt; + // std::cout << "x_root_grad = " << x_root_grad << std::endl; const auto &phi_interp_prev = levelset.phiInterpPrev(); const auto &psi_interp_prev = levelset.psiInterpPrev(); + // std::cout << "phi_interp_prev = " << phi_interp_prev << std::endl; + // std::cout << "psi_interp_prev = " << psi_interp_prev << std::endl; // Update phi and psi. for (int i = 0; i < num_cells[0]; ++i) @@ -86,6 +98,8 @@ void GALS::TEMPORAL_SCHEMES::SEMI_LAGRANGIAN::Euler::compute( phi(i, j, k) = phi_interp_prev(i, j, k); psi(i, j, k) = x_root_grad(i, j, k).dot(psi_interp_prev(i, j, k)); } + // std::cout << "phi = " << phi << std::endl; + // std::cout << "psi = " << psi << std::endl; } template class GALS::TEMPORAL_SCHEMES::SEMI_LAGRANGIAN::Euler>; diff --git a/src/cpu/vec3.cc b/src/cpu/vec3.cc index b9a67b5..00626e1 100644 --- a/src/cpu/vec3.cc +++ b/src/cpu/vec3.cc @@ -76,8 +76,8 @@ const double GALS::CPU::Vec3::mag() const template const T GALS::CPU::Vec3::min() const { - if (m_data[0] < m_data[1] && m_data[0] < m_data[2]) return m_data[0]; - if (m_data[1] < m_data[0] && m_data[1] < m_data[2]) return m_data[1]; + if (m_data[0] <= m_data[1] && m_data[0] <= m_data[2]) return m_data[0]; + if (m_data[1] <= m_data[0] && m_data[1] <= m_data[2]) return m_data[1]; return m_data[2]; } diff --git a/tests/inputs b/tests/inputs index 3af485b..ad3b772 100644 --- a/tests/inputs +++ b/tests/inputs @@ -17,6 +17,8 @@ time: constant_dt: NO # NO/YES # NO: time step will be dynamically calculated. # YES: `dt` will be used for entire simulation. + cfl_max: 0.5 # Max CFL number for simulation. + write_interval: 1e-4 # Write interval for output files. velocity: name: CIRCULAR # Available: diff --git a/tools/octave/Visualize.m b/tools/octave/Visualize.m index dd1a10b..96dd628 100644 --- a/tools/octave/Visualize.m +++ b/tools/octave/Visualize.m @@ -1,24 +1,28 @@ -% Visualize fields. -% -% Fields that will be read are: -% phi -% psi -% velocity -% -% Data formats: x y z (or) x y z ... - -function Visualize(path, dim, nx, ny, nz) - - phi_str = strcat(path, '/phi'); - psi_str = strcat(path, '/psi'); - velocity_str = strcat(path, '/velocity'); - - if exist(phi_str, 'file') == 2 - ShowScalar(phi_str, dim, nx, ny, nz); +% Visualize fields.% % Fields that will be read are : % phi % psi % velocity % % + Data formats + : x y z(or) x y z... + + function Visualize(path, dim, nx, ny, nz) + + num_dirs = length(ls(path)); + +fprintf('num_dirs = %d\n', num_dirs); + for + d = 1 : num_dirs + + phi_str = strcat(path, '/', num2str(d), '/phi'); + psi_str = strcat(path, '/', num2str(d), '/psi'); + velocity_str = strcat(path, '/', num2str(d), '/velocity'); + + if + exist(phi_str, 'file') == 2 ShowScalar(phi_str, dim, nx, ny, nz); end - - if exist(velocity_str, 'file') == 2 - ShowVector(velocity_str, dim); + + if exist(velocity_str, 'file') == 2 ShowVector(velocity_str, dim); + end + + fprintf('d = %d\n', d); + end end @@ -36,8 +40,9 @@ function ShowScalar(path, dim, nx, ny, nz) a = reshape(ts{ 4 }, nx, ny); if (dim == 2) - figure() - surfc(x, y, a); + figure(1) contour(x, y, a, [0 0]); + axis square; + % surfc(x, y, a); end end @@ -57,8 +62,9 @@ function ShowVector(path, dim) w = ts{6}; if (dim == 2) - figure(); + figure(2); quiver(x, y, u, v); - end + axis square; + end -end + end