LevelsetVelocity is used to compute analytical velocities

* Velocity gradient can be computed analytically for solid body rotation field.

Signed-off-by: Lakshman Anumolu <[email protected]>

Author: acrlakshman

applications/advection/main.cc

@@ -101,9 +101,10 @@ int main(int argc, char **argv)
 
   // Construct velocity.
   const auto &velocity_inputs = *(input_fields.m_velocity);
-  GALS::CPU::Array<T_GRID, GALS::CPU::Vec3<T>> velocity_field(grid);
+  GALS::CPU::LevelsetVelocity<T_GRID, T> levelset_velocity(grid);
+  const auto &velocity_field = levelset_velocity.velocity();
   GALS::ANALYTICAL_FIELDS::Velocity<T_GRID, T> velocity(grid, velocity_inputs);
-  velocity.compute(positions, velocity_field);
+  velocity.compute(positions, levelset_velocity);
 
   const auto &time_inputs = *(input_fields.m_time);
 

include/gals/analytical-fields/velocity.h

@@ -35,6 +35,7 @@
 #include <string>
 #include <vector>
 
+#include "gals/cpu/levelset-velocity.h"
 #include "gals/input-fields/velocity.h"
 #include "gals/utilities/array.h"
 #include "gals/utilities/grid.h"
@@ -84,10 +85,10 @@ class Velocity
   /*! Compute velocity field.
    *
    * \param positions array of positions where velocity needs to be computed.
-   * \param velocity velocity field which will be updated by this function.
+   * \param levelset_velocity velocity field which will be updated by this function.
    */
   void compute(const GALS::CPU::Array<T_GRID, GALS::CPU::Vec3<T>>& positions,
-               GALS::CPU::Array<T_GRID, GALS::CPU::Vec3<T>>& velocity);
+               GALS::CPU::LevelsetVelocity<T_GRID, T>& levelset_velocity);
 
  private:
   const T_GRID& m_grid;

src/cpu/analytical-fields/velocity.cc

@@ -46,13 +46,15 @@ GALS::ANALYTICAL_FIELDS::Velocity<T_GRID, T>::~Velocity()
 template <typename T_GRID, typename T>
 void GALS::ANALYTICAL_FIELDS::Velocity<T_GRID, T>::compute(
     const GALS::CPU::Array<T_GRID, GALS::CPU::Vec3<T>>& positions,
-    GALS::CPU::Array<T_GRID, GALS::CPU::Vec3<T>>& velocity)
+    GALS::CPU::LevelsetVelocity<T_GRID, T>& levelset_velocity)
 {
   const GALS::CPU::Vec3<int> num_cells = positions.numCells();
   const auto& velocity_name = m_inputs.name;
 
   switch (velocity_name_map[velocity_name]) {
     case VelocityFieldNames::CIRCULAR: {
+      auto& velocity = levelset_velocity.velocity();
+
       if (T_GRID::dim == 1)
         GALS_FUNCTION_NOT_IMPLEMENTED(
             "GALS::ANALYTICAL_FIELDS::Velocity::compute: CIRCULAR velocity field for 1D grid.");
@@ -71,6 +73,24 @@ void GALS::ANALYTICAL_FIELDS::Velocity<T_GRID, T>::compute(
             velocity(i, j, k) = velocity_node;
           }
 
+      // Compute velocity gradient using analytical expressions.
+      if (!std::strcmp(m_inputs.gradient_scheme, "ANALYTICAL")) {
+        auto& velocity_gradient = levelset_velocity.velocityGradient();
+
+        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) {
+              for (int cmpt = 0; cmpt < T_GRID::dim; ++cmp) {
+                for (int derv = 0; derv < T_GRID::dim; ++derv) {
+                  if (cmpt == 0 && derv == 0) velocity_gradient[cmpt][derv] = static_cast<T>(0);
+                  if (cmpt == 0 && derv == 1) velocity_gradient[cmpt][derv] = -coeff;
+                  if (cmpt == 1 && derv == 0) velocity_gradient[cmpt][derv] = coeff;
+                  if (cmpt == 1 && derv == 1) velocity_gradient[cmpt][derv] = static_cast<T>(0);
+                }
+              }
+            }
+      }
+
       break;
     }
 

tests/analytical-fields-velocity.cc

@@ -103,9 +103,11 @@ TEST(GALS, ANALYTICAL_FIELDS_VELOCITY_2D)
         positions(i, j, k) = grid(i, j, k);
       }
 
-  GALS::CPU::Array<T_GRID, GALS::CPU::Vec3<T>> velocity_field(grid);
+  GALS::CPU::LevelsetVelocity<T_GRID, T> levelset_velocity(grid);
+  const auto& velocity_field = levelset_velocity.velocity();
+
   GALS::ANALYTICAL_FIELDS::Velocity<T_GRID, T> velocity(grid, velocity_inputs);
-  velocity.compute(positions, velocity_field);
+  velocity.compute(positions, levelset_velocity);
 
   // TODO (lakshman): Compare these values with results from matlab.
   // for (int i = 0; i < num_cells[0]; ++i)