Line data    Source code

       1              : //
       2              : // File BareFieldOperations
       3              : //   Norms and a scalar product for fields
       4              : //
       5              : 
       6              : namespace ippl {
       7              :     /*!
       8              :      * Computes the inner product of two fields
       9              :      * @param f1 first field
      10              :      * @param f2 second field
      11              :      * @return Result of f1^T f2
      12              :      */
      13              :     template <typename BareField>
      14            0 :     typename BareField::value_type innerProduct(const BareField& f1, const BareField& f2) {
      15              :         using T                = typename BareField::value_type;
      16            0 :         constexpr unsigned Dim = BareField::dim;
      17              : 
      18            0 :         T sum                  = 0;
      19            0 :         auto layout            = f1.getLayout();
      20            0 :         auto view1             = f1.getView();
      21            0 :         auto view2             = f2.getView();
      22              :         using exec_space       = typename BareField::execution_space;
      23              :         using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
      24            0 :         ippl::parallel_reduce(
      25            0 :             "Field::innerProduct(Field&, Field&)", f1.getFieldRangePolicy(),
      26            0 :             KOKKOS_LAMBDA(const index_array_type& args, T& val) {
      27            0 :                 val += apply(view1, args) * apply(view2, args);
      28              :             },
      29            0 :             Kokkos::Sum<T>(sum));
      30            0 :         T globalSum = 0;
      31            0 :         layout.comm.allreduce(sum, globalSum, 1, std::plus<T>());
      32            0 :         return globalSum;
      33            0 :     }
      34              : 
      35              :     /*!
      36              :      * Computes the Lp-norm of a field
      37              :      * @param field field
      38              :      * @param p desired norm (default 2)
      39              :      * @return The desired norm of the field
      40              :      */
      41              :     template <typename BareField>
      42           92 :     typename BareField::value_type norm(const BareField& field, int p = 2) {
      43              :         using T                = typename BareField::value_type;
      44           92 :         constexpr unsigned Dim = BareField::dim;
      45              : 
      46           92 :         T local                = 0;
      47           92 :         auto layout            = field.getLayout();
      48           92 :         auto view              = field.getView();
      49              :         using exec_space       = typename BareField::execution_space;
      50              :         using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
      51           92 :         switch (p) {
      52           24 :             case 0: {
      53           24 :                 ippl::parallel_reduce(
      54           48 :                     "Field::norm(0)", field.getFieldRangePolicy(),
      55       477360 :                     KOKKOS_LAMBDA(const index_array_type& args, T& val) {
      56       238656 :                         T myVal = Kokkos::abs(apply(view, args));
      57       238656 :                         if (myVal > val)
      58          844 :                             val = myVal;
      59              :                     },
      60           48 :                     Kokkos::Max<T>(local));
      61           24 :                 T globalMax = 0;
      62           24 :                 layout.comm.allreduce(local, globalMax, 1, std::greater<T>());
      63           24 :                 return globalMax;
      64              :             }
      65           68 :             default: {
      66           68 :                 ippl::parallel_reduce(
      67           68 :                     "Field::norm(int) general", field.getFieldRangePolicy(),
      68       955900 :                     KOKKOS_LAMBDA(const index_array_type& args, T& val) {
      69       477882 :                         val += std::pow(Kokkos::abs(apply(view, args)), p);
      70              :                     },
      71          136 :                     Kokkos::Sum<T>(local));
      72           68 :                 T globalSum = 0;
      73           68 :                 layout.comm.allreduce(local, globalSum, 1, std::plus<T>());
      74           68 :                 return std::pow(globalSum, 1.0 / p);
      75              :             }
      76              :         }
      77           92 :     }
      78              : }  // namespace ippl