Higher order shape function

CMakeLists.txt

@@ -34,6 +34,15 @@ set(DEBUG
 set(precision
     ${default_precision}
     CACHE STRING "Precision")
+
+set(deposit
+    ${default_deposit}
+    CACHE STRING "Deposit")
+
+set(shape_order
+    ${default_shape_order}
+    CACHE STRING "Shape function")
+
 set(pgen
     ${default_pgen}
     CACHE STRING "Problem generator")
@@ -75,6 +84,20 @@ set(precisions
     "single" "double"
     CACHE STRING "Precisions")
 
+set(deposits
+    "zigzag" "esirkepov"
+    CACHE STRING "Deposits")
+
+if(${deposit} STREQUAL "zigzag")
+  set(shape_order
+      ${default_shape_order}
+      CACHE STRING "Shape functions")
+endif()
+
+set(shape_orders
+    "1;2;3;4;5;6;7;8;9;10;11"
+    CACHE STRING "Shape orders")
+
 include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/config.cmake)
 
 # ------------------------- Third-Party Tests ------------------------------ #
@@ -92,6 +115,7 @@ include_directories(${plog_SRC}/include)
 
 # -------------------------------- Main code ------------------------------- #
 set_precision(${precision})
+set_shape_order(${shape_order})
 if("${Kokkos_DEVICES}" MATCHES "CUDA")
   add_compile_options("-D CUDA_ENABLED")
   set(DEVICE_ENABLED ON)

cmake/config.cmake

@@ -16,6 +16,16 @@ function(set_precision precision_name)
   endif()
 endfunction()
 
+# ------------------------------- Shape function --------------------------- #
+function(set_shape_order shape_order)
+  if(${deposit} STREQUAL "esirkepov")
+    if(${shape_order} GREATER 11)
+      message(FATAL_ERROR "Shape order must be between 1 and 11.")
+    endif()
+    add_compile_options("-DSHAPE_ORDER=${shape_order}")
+  endif()
+endfunction()
+
 # ---------------------------- Problem generator --------------------------- #
 function(set_problem_generator pgen_name)
   if(pgen_name STREQUAL ".")

cmake/defaults.cmake

@@ -19,6 +19,12 @@ set(default_engine
 set(default_precision
     "single"
     CACHE INTERNAL "Default precision")
+set(default_deposit
+    "zigzag"
+    CACHE INTERNAL "Default deposit")
+set(default_shape_order
+    1
+    CACHE INTERNAL "Default shape function order")
 set(default_pgen
     "."
     CACHE INTERNAL "Default problem generator")

cmake/report.cmake

@@ -37,6 +37,24 @@ printchoices(
   "${Blue}"
   PRECISION_REPORT
   46)
+printchoices(
+  "Deposit"
+  "deposit"
+  "${deposits}"
+  ${deposit}
+  ${default_deposit}
+  "${Blue}"
+  DEPOSIT_REPORT
+  46)
+printchoices(
+  "Shape order"
+  "shape_order"
+  "${shape_orders}"
+  ${shape_order}
+  ${default_shape_order}
+  "${Blue}"
+  SHAPEFUNCTION_REPORT
+  46)
 printchoices(
   "Output"
   "output"
@@ -113,6 +131,12 @@ string(
   ${PRECISION_REPORT}
   "\n"
   "  "
+  ${DEPOSIT_REPORT}
+  "\n"
+  "  "
+  ${SHAPEFUNCTION_REPORT}
+  "\n"
+  "  "
   ${OUTPUT_REPORT}
   "\n")
 

cmake/styling.cmake

@@ -140,27 +140,30 @@ function(
   else()
     padto("${rstring}" " " ${Padding} rstring)
 
-    set(new_choices ${Choices})
-    foreach(ch IN LISTS new_choices)
-      string(REPLACE ${ch} "${Dim}${ch}${ColorReset}" new_choices
-                     "${new_choices}")
-    endforeach()
-    set(Choices ${new_choices})
     if(${Value} STREQUAL "ON")
       set(col ${Green})
     elseif(${Value} STREQUAL "OFF")
       set(col ${Red})
     else()
       set(col ${Color})
     endif()
-    if(NOT "${Value}" STREQUAL "")
-      string(REPLACE ${Value} "${col}${Value}${ColorReset}" Choices
-                     "${Choices}")
-    endif()
-    if(NOT "${Default}" STREQUAL "")
-      string(REPLACE ${Default} "${Underline}${Default}${ColorReset}" Choices
-                     "${Choices}")
-    endif()
+    set(new_choices "")
+    foreach(ch IN LISTS Choices)
+      set(elem "${ch}")
+      if((NOT "${Value}" STREQUAL "") AND (${ch} STREQUAL ${Value}))
+        set(elem "${col}${ch}${ColorReset}")
+      else()
+        set(elem "${Dim}${ch}${ColorReset}")
+      endif()
+      if((NOT "${Default}" STREQUAL "") AND (${ch} STREQUAL ${Default}))
+        set(elem "${Underline}${elem}${ColorReset}")
+      endif()
+      string(APPEND new_choices "${elem};")
+    endforeach()
+    string(LENGTH "${new_choices}" nlen)
+    math(EXPR nlen "${nlen} - 1")
+    string(SUBSTRING "${new_choices}" 0 ${nlen} new_choices)
+    set(Choices ${new_choices})
     string(REPLACE ";" "/" Choices "${Choices}")
     string(APPEND rstring "${Choices}")
   endif()

input.example.toml

@@ -201,15 +201,35 @@
   #   @default: 0
   current_filters = ""
 
-  [algorithms.toggles]
-    # Toggle for the field solver
+  [algorithms.deposit]
+    # Enable the current deposition
     #   @type: bool
     #   @default: true
-    fieldsolver = ""
-    # Toggle for the current deposition
+    enable = ""
+    # Order of the particle shape function
+    #   @type: int
+    #   @default: 1
+    order = ""
+
+  # @TODO: fix fieldsolver params below
+  [algorithms.fieldsolver]
+    # Enable the EM fieldsolver
     #   @type: bool
     #   @default: true
-    deposit = ""
+    enable = ""
+    # Yee - all 0.0 - default
+    # 1D
+    deltax = -0.065
+    # 2D
+    deltay = -0.065
+    betaxy = -0.065
+    betayx = -0.065
+    # 3D - not yet tested
+    deltaz = 0.0
+    betaxz = 0.0
+    betazx = 0.0
+    betayz = 0.0
+    betazy = 0.0
 
   [algorithms.timestep]
     # Courant-Friedrichs-Lewy number

src/engines/srpic.hpp

@@ -76,9 +76,9 @@ namespace ntt {
 
     void step_forward(timer::Timers& timers, domain_t& dom) override {
       const auto fieldsolver_enabled = m_params.template get<bool>(
-        "algorithms.toggles.fieldsolver");
+        "algorithms.fieldsolver.enable");
       const auto deposit_enabled = m_params.template get<bool>(
-        "algorithms.toggles.deposit");
+        "algorithms.deposit.enable");
       const auto clear_interval = m_params.template get<std::size_t>(
         "particles.clear_interval");
 
@@ -505,9 +505,30 @@ namespace ntt {
       }
     }
 
+    template <unsigned short O>
+    void deposit_with(const Particles<M::Dim, M::CoordType>& species,
+                      const M&                               metric,
+                      const scatter_ndfield_t<M::Dim, 3>&    scatter_cur,
+                      real_t                                 dt) {
+      // clang-format off
+      Kokkos::parallel_for("CurrentsDeposit",
+                          species.rangeActiveParticles(),
+                          kernel::DepositCurrents_kernel<SimEngine::SRPIC, M, O>(
+                            scatter_cur,
+                            species.i1, species.i2, species.i3,
+                            species.i1_prev, species.i2_prev, species.i3_prev,
+                            species.dx1, species.dx2, species.dx3,
+                            species.dx1_prev, species.dx2_prev, species.dx3_prev,
+                            species.ux1, species.ux2, species.ux3,
+                            species.phi, species.weight, species.tag,
+                            metric, (real_t)(species.charge()), dt));
+      // clang-format on
+    }
+
     void CurrentsDeposit(domain_t& domain) {
       auto scatter_cur = Kokkos::Experimental::create_scatter_view(
         domain.fields.cur);
+      auto shape_order = m_params.template get<int>("algorithms.deposit.order");
       for (auto& species : domain.species) {
         if ((species.pusher() == PrtlPusher::NONE) or (species.npart() == 0) or
             cmp::AlmostZero_host(species.charge())) {
@@ -520,20 +541,8 @@ namespace ntt {
                       species.npart(),
                       (double)species.charge()),
           HERE);
-        // clang-format off
-        Kokkos::parallel_for("CurrentsDeposit",
-                             species.rangeActiveParticles(),
-                             kernel::DepositCurrents_kernel<SimEngine::SRPIC, M>(
-                               scatter_cur,
-                               species.i1, species.i2, species.i3,
-                               species.i1_prev, species.i2_prev, species.i3_prev,
-                               species.dx1, species.dx2, species.dx3,
-                               species.dx1_prev, species.dx2_prev, species.dx3_prev,
-                               species.ux1, species.ux2, species.ux3,
-                               species.phi, species.weight, species.tag,
-                               domain.mesh.metric,
-                               (real_t)(species.charge()), dt));
-        // clang-format on
+
+        deposit_with<SHAPE_ORDER>(species, domain.mesh.metric, scatter_cur, dt);
       }
       Kokkos::Experimental::contribute(domain.fields.cur, scatter_cur);
     }

src/framework/parameters.cpp

@@ -409,13 +409,16 @@ namespace ntt {
                       "current_filters",
                       defaults::current_filters));
 
-    /* [algorithms.toggles] ------------------------------------------------- */
-    set("algorithms.toggles.fieldsolver",
-        toml::find_or(toml_data, "algorithms", "toggles", "fieldsolver", true));
-    set("algorithms.toggles.deposit",
-        toml::find_or(toml_data, "algorithms", "toggles", "deposit", true));
+    /* [algorithms.deposit] ------------------------------------------------- */
+    set("algorithms.deposit.enable",
+        toml::find_or(toml_data, "algorithms", "deposit", "enable", true));
+    set("algorithms.deposit.order",
+        toml::find_or(toml_data, "algorithms", "deposit", "order", 1));
 
     /* [algorithms.fieldsolver] --------------------------------------------- */
+    set("algorithms.fieldsolver.enable",
+        toml::find_or(toml_data, "algorithms", "fieldsolver", "enable", true));
+
     set("algorithms.fieldsolver.delta_x",
         toml::find_or(toml_data,
                       "algorithms",

src/global/global.h

@@ -113,7 +113,15 @@ namespace files {
 
 namespace ntt {
 
+#if !defined(SHAPE_ORDER)
+  #define SHAPE_ORDER 0
   inline constexpr std::size_t N_GHOSTS = 2;
+#else  // SHAPE_ORDER
+  inline constexpr std::size_t N_GHOSTS = static_cast<std::size_t>(
+                                            (SHAPE_ORDER + 1) / 2) +
+                                          1;
+#endif // SHAPE_ORDER
+
 // Coordinate shift to account for ghost cells
 #define COORD(I)                                                               \
   (static_cast<real_t>(static_cast<int>((I)) - static_cast<int>(N_GHOSTS)))

src/global/utils/numeric.h

@@ -36,9 +36,13 @@ inline constexpr float TWO    = 2.0f;
 inline constexpr float THREE  = 3.0f;
 inline constexpr float FOUR   = 4.0f;
 inline constexpr float FIVE   = 5.0f;
+inline constexpr float SIX    = 6.0f;
 inline constexpr float TWELVE = 12.0f;
 inline constexpr float ZERO   = 0.0f;
 inline constexpr float HALF   = 0.5f;
+inline constexpr float THIRD  = 0.333333f;
+inline constexpr float THREE_FOURTHS = 0.75f;
+inline constexpr float THREE_HALFS   = 1.5f;
 inline constexpr float INV_2  = 0.5f;
 inline constexpr float INV_4  = 0.25f;
 inline constexpr float INV_8  = 0.125f;
@@ -51,9 +55,13 @@ inline constexpr double TWO    = 2.0;
 inline constexpr double THREE  = 3.0;
 inline constexpr double FOUR   = 4.0;
 inline constexpr double FIVE   = 5.0;
+inline constexpr double SIX    = 6.0;
 inline constexpr double TWELVE = 12.0;
 inline constexpr double ZERO   = 0.0;
 inline constexpr double HALF   = 0.5;
+inline constexpr double THIRD  = 0.3333333333333333;
+inline constexpr double THREE_FOURTHS = 0.75;
+inline constexpr float  THREE_HALFS   = 1.5;
 inline constexpr double INV_2  = 0.5;
 inline constexpr double INV_4  = 0.25;
 inline constexpr double INV_8  = 0.125;