Branch data     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                 :          36 :     typename BareField::value_type norm(const BareField& field, int p = 2) {
      43                 :             :         using T                = typename BareField::value_type;
      44                 :          36 :         constexpr unsigned Dim = BareField::dim;
      45                 :             : 
      46                 :          36 :         T local                = 0;
      47   [ +  -  +  - ]:          36 :         auto layout            = field.getLayout();
      48         [ +  - ]:          36 :         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         [ +  + ]:          36 :         switch (p) {
      52                 :          12 :             case 0: {
      53   [ +  -  +  - ]:          12 :                 ippl::parallel_reduce(
      54         [ +  - ]:          24 :                     "Field::norm(0)", field.getFieldRangePolicy(),
      55         [ +  - ]:      477336 :                     KOKKOS_LAMBDA(const index_array_type& args, T& val) {
      56                 :      238656 :                         T myVal = std::abs(apply(view, args));
      57   [ +  +  +  -  :      238656 :                         if (myVal > val)
          +  +  +  -  +  
          +  +  +  +  +  
          +  +  +  +  +  
             +  +  +  +  
                      + ]
           [ #  #  #  # ]
                 [ #  # ]
      58                 :         614 :                             val = myVal;
      59                 :             :                     },
      60         [ +  - ]:          24 :                     Kokkos::Max<T>(local));
      61                 :          12 :                 T globalMax = 0;
      62         [ +  - ]:          12 :                 layout.comm.allreduce(local, globalMax, 1, std::greater<T>());
      63                 :          12 :                 return globalMax;
      64                 :             :             }
      65                 :          24 :             default: {
      66   [ +  -  +  - ]:          24 :                 ippl::parallel_reduce(
      67                 :          24 :                     "Field::norm(int) general", field.getFieldRangePolicy(),
      68   [ +  -  +  -  :      954672 :                     KOKKOS_LAMBDA(const index_array_type& args, T& val) {
                   -  - ]
      69                 :      477312 :                         val += std::pow(std::abs(apply(view, args)), p);
      70                 :             :                     },
      71         [ +  - ]:          48 :                     Kokkos::Sum<T>(local));
      72                 :          24 :                 T globalSum = 0;
      73         [ +  - ]:          24 :                 layout.comm.allreduce(local, globalSum, 1, std::plus<T>());
      74                 :          24 :                 return std::pow(globalSum, 1.0 / p);
      75                 :             :             }
      76                 :             :         }
      77                 :          36 :     }
      78                 :             : }  // namespace ippl