compute mixed derivatives

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)

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

applications/advection/main.cc

@@ -29,6 +29,7 @@
 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 ///////////////////////////////////////////////////////////////////////////////
 
+#include <algorithm>
 #include <iostream>
 #include <string>
 
@@ -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<T>;
+using T_GRID = GALS::CPU::Grid<T, dim>;
+
+static T compute_dt(const T cfl_max, const T_GRID &grid, const GALS::CPU::Array<T_GRID, TV> &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<T>;
-  using T_GRID = GALS::CPU::Grid<T, dim>;
-
   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<T_GRID, T> 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<T>(0);
 
   if (!is_dt_fixed) {
     dt = grid.dX().min() * static_cast<T>(0.5);
   }
+  std::cout << "here " << ++iii << std::endl;
+
+  // Compute gradient of levelset field.
+  GALS::CPU::Gradient<T, T_GRID, GALS::CPU::ThirdOrder<T, T_GRID>>::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<T, T_GRID, GALS::CPU::ThirdOrder<T, T_GRID>>::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<T, T_GRID, GALS::TEMPORAL_SCHEMES::SEMI_LAGRANGIAN::Euler<T, T_GRID>>::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<T_GRID, TV> &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<T>(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<T>(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;
+}

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

include/gals/cpu/interpolation/hermite.h

@@ -178,7 +178,7 @@ class Hermite<T, GALS::CPU::Grid<T, 1>>
    */
   CPU::InterpolatedFields<CPU::Vec3<T>> interpolate(
       const GALS::CPU::Grid<typename GALS::CPU::Grid<T, 1>::value_type, GALS::CPU::Grid<T, 1>::dim> &grid,
-      const typename GALS::CPU::Grid<T, 1>::position_type &x_interp,
+      const GALS::CPU::Vec3<int> &node_id, const typename GALS::CPU::Grid<T, 1>::position_type &x_interp,
       const CPU::Levelset<GALS::CPU::Grid<T, 1>, T> &levelset, const bool use_gradient_limiting = false);
 
   /*! Interpolate scalar field.
@@ -247,7 +247,7 @@ class Hermite<T, GALS::CPU::Grid<T, 2>>
    */
   CPU::InterpolatedFields<CPU::Vec3<T>> interpolate(
       const GALS::CPU::Grid<typename GALS::CPU::Grid<T, 2>::value_type, GALS::CPU::Grid<T, 2>::dim> &grid,
-      const typename GALS::CPU::Grid<T, 2>::position_type &x_interp,
+      const GALS::CPU::Vec3<int> &node_id, const typename GALS::CPU::Grid<T, 2>::position_type &x_interp,
       const CPU::Levelset<GALS::CPU::Grid<T, 2>, T> &levelset, const bool use_gradient_limiting = false);
 
   /*! Interpolate scalar field.
@@ -316,7 +316,7 @@ class Hermite<T, GALS::CPU::Grid<T, 3>>
    */
   CPU::InterpolatedFields<CPU::Vec3<T>> interpolate(
       const GALS::CPU::Grid<typename GALS::CPU::Grid<T, 3>::value_type, GALS::CPU::Grid<T, 3>::dim> &grid,
-      const typename GALS::CPU::Grid<T, 3>::position_type &x_interp,
+      const GALS::CPU::Vec3<int> &node_id, const typename GALS::CPU::Grid<T, 3>::position_type &x_interp,
       const CPU::Levelset<GALS::CPU::Grid<T, 3>, T> &levelset, const bool use_gradient_limiting = false);
 
   /*! Interpolate scalar field.

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<T_GRID, Vec3<T>>& psi,
+                               Array<T_GRID, VecN<T, T_GRID::num_mixed_derivatives>>& phi_mixed_derivatives);
+
   /*! Return grid.
    *
    * \return grid.

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;
 };
 

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<int> 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;

src/cpu/grid.cc

@@ -215,6 +215,21 @@ void GALS::CPU::Grid<T, DIM>::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 <typename T, int DIM>
+bool GALS::CPU::Grid<T, DIM>::isIndexInDomain(const Vec3<int> 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<double, 1>;
 template class GALS::CPU::Grid<double, 2>;
 template class GALS::CPU::Grid<double, 3>;

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<double>();
   input_fields.m_time->dt = field["dt"].as<double>();
   input_fields.m_time->constant_dt = field["constant_dt"].as<std::string>();
+  input_fields.m_time->cfl_max = field["cfl_max"].as<double>();
+  input_fields.m_time->write_interval = field["write_interval"].as<double>();
 }
 
 void GALS::INPUT_PARSER::Time::operator()(const YAML::Node &field, GALS::INPUT_FIELDS::InputFields *p_input_fields)

src/cpu/interpolation/hermite-1d.cc

@@ -36,7 +36,7 @@
 template <typename T>
 GALS::CPU::InterpolatedFields<GALS::CPU::Vec3<T>> GALS::INTERPOLATION::Hermite<T, GALS::CPU::Grid<T, 1>>::interpolate(
     const GALS::CPU::Grid<typename GALS::CPU::Grid<T, 1>::value_type, GALS::CPU::Grid<T, 1>::dim> &grid,
-    const typename GALS::CPU::Grid<T, 1>::position_type &x_interp,
+    const GALS::CPU::Vec3<int> &node_id, const typename GALS::CPU::Grid<T, 1>::position_type &x_interp,
     const GALS::CPU::Levelset<GALS::CPU::Grid<T, 1>, T> &levelset, const bool use_gradient_limiting)
 {
   GALS::CPU::InterpolatedFields<GALS::CPU::Vec3<T>> hermite_fields;
@@ -96,7 +96,8 @@ void GALS::INTERPOLATION::Hermite<T, GALS::CPU::Grid<T, 1>>::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<int> 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;