add L2norm methods

bfps/cpp/field.cpp

@@ -1152,6 +1152,42 @@ void field<rnumber, be, fc>::compute_stats(
     }
 }
 
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+template <kspace_dealias_type dt>
+double field<rnumber, be, fc>::L2norm(
+        kspace<be, dt> *kk)
+{
+    TIMEZONE("field::L2norm");
+    if (!this->real_space_representation)
+        return kk->template L2norm<rnumber, fc>(this->get_cdata());
+    else
+    {
+        shared_array<double> local_m2_threaded(1, [&](double* local_moment){
+            std::fill_n(local_moment, 1, 0);});
+
+        this->RLOOP(
+                [&](ptrdiff_t rindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex){
+                double *local_m2 = local_m2_threaded.getMine();
+                for (unsigned int i=0; i<ncomp(fc); i++)
+                    local_m2[0] += this->data[rindex*ncomp(fc)+i]*this->data[rindex*ncomp(fc)+i];
+            });
+
+        local_m2_threaded.mergeParallel();
+        double m2;
+        MPI_Allreduce(
+                (void*)local_m2_threaded.getMasterData(),
+                &m2,
+                1,
+                MPI_DOUBLE, MPI_SUM, this->comm);
+        return m2 / this->npoints;
+    }
+}
+
 template <typename rnumber,
           field_backend be,
           field_components fc1,
@@ -1655,6 +1691,34 @@ template void field<double, FFTW, THREExTHREE>::compute_stats<SMOOTH>(
         kspace<FFTW, SMOOTH> *,
         const hid_t, const std::string, const hsize_t, const double);
 
+template double field<float, FFTW, ONE>::L2norm<TWO_THIRDS>(
+        kspace<FFTW, TWO_THIRDS> *);
+template double field<float, FFTW, THREE>::L2norm<TWO_THIRDS>(
+        kspace<FFTW, TWO_THIRDS> *);
+template double field<float, FFTW, THREExTHREE>::L2norm<TWO_THIRDS>(
+        kspace<FFTW, TWO_THIRDS> *);
+
+template double field<double, FFTW, ONE>::L2norm<TWO_THIRDS>(
+        kspace<FFTW, TWO_THIRDS> *);
+template double field<double, FFTW, THREE>::L2norm<TWO_THIRDS>(
+        kspace<FFTW, TWO_THIRDS> *);
+template double field<double, FFTW, THREExTHREE>::L2norm<TWO_THIRDS>(
+        kspace<FFTW, TWO_THIRDS> *);
+
+template double field<float, FFTW, ONE>::L2norm<SMOOTH>(
+        kspace<FFTW, SMOOTH> *);
+template double field<float, FFTW, THREE>::L2norm<SMOOTH>(
+        kspace<FFTW, SMOOTH> *);
+template double field<float, FFTW, THREExTHREE>::L2norm<SMOOTH>(
+        kspace<FFTW, SMOOTH> *);
+
+template double field<double, FFTW, ONE>::L2norm<SMOOTH>(
+        kspace<FFTW, SMOOTH> *);
+template double field<double, FFTW, THREE>::L2norm<SMOOTH>(
+        kspace<FFTW, SMOOTH> *);
+template double field<double, FFTW, THREExTHREE>::L2norm<SMOOTH>(
+        kspace<FFTW, SMOOTH> *);
+
 template int compute_gradient<float, FFTW, THREE, THREExTHREE, SMOOTH>(
         kspace<FFTW, SMOOTH> *,
         field<float, FFTW, THREE> *,

bfps/cpp/field.hpp

@@ -249,6 +249,9 @@ class field
                 const std::string dset_name,
                 const hsize_t toffset,
                 const double max_estimate);
+        template <kspace_dealias_type dt>
+        double L2norm(
+                kspace<be, dt> *kk);
         inline void impose_zero_mode()
         {
             if (this->clayout->myrank == this->clayout->rank[0][0] &&

bfps/cpp/kspace.cpp

@@ -597,6 +597,48 @@ void kspace<be, dt>::cospectrum(
     }
 }
 
+template <field_backend be,
+          kspace_dealias_type dt>
+template <typename rnumber,
+          field_components fc>
+double kspace<be, dt>::L2norm(
+        const rnumber(* __restrict a)[2])
+{
+    TIMEZONE("field::L2norm");
+    shared_array<double> L2_local_thread(1,[&](double* spec_local){
+        std::fill_n(spec_local, 1, 0);
+    });
+
+    this->CLOOP_K2_NXMODES(
+            [&](ptrdiff_t cindex,
+                ptrdiff_t xindex,
+                ptrdiff_t yindex,
+                ptrdiff_t zindex,
+                double k2,
+                int nxmodes){
+            if (k2 <= this->kM2)
+            {
+                double* L2_local = L2_local_thread.getMine();
+                for (hsize_t i=0; i<ncomp(fc); i++)
+                for (hsize_t j=0; j<ncomp(fc); j++){
+                    L2_local[0] += nxmodes * (
+                        (a[ncomp(fc)*cindex + i][0] * a[ncomp(fc)*cindex + j][0]) +
+                        (a[ncomp(fc)*cindex + i][1] * a[ncomp(fc)*cindex + j][1]));
+                }
+            }
+            });
+
+    L2_local_thread.mergeParallel();
+
+    double L2;
+    MPI_Allreduce(
+            L2_local_thread.getMasterData(),
+            &L2,
+            1,
+            MPI_DOUBLE, MPI_SUM, this->layout->comm);
+    return L2 * this->dkx * this->dky * this->dkz;
+}
+
 
 template class kspace<FFTW, TWO_THIRDS>;
 template class kspace<FFTW, SMOOTH>;
@@ -801,6 +843,32 @@ template void kspace<FFTW, SMOOTH>::cospectrum<double, THREExTHREE>(
         const std::string dset_name,
         const hsize_t toffset);
 
+template double kspace<FFTW, TWO_THIRDS>::L2norm<float, ONE>(
+        const typename fftw_interface<float>::complex *__restrict__ a);
+template double kspace<FFTW, TWO_THIRDS>::L2norm<float, THREE>(
+        const typename fftw_interface<float>::complex *__restrict__ a);
+template double kspace<FFTW, TWO_THIRDS>::L2norm<float, THREExTHREE>(
+        const typename fftw_interface<float>::complex *__restrict__ a);
+template double kspace<FFTW, TWO_THIRDS>::L2norm<double, ONE>(
+        const typename fftw_interface<double>::complex *__restrict__ a);
+template double kspace<FFTW, TWO_THIRDS>::L2norm<double, THREE>(
+        const typename fftw_interface<double>::complex *__restrict__ a);
+template double kspace<FFTW, TWO_THIRDS>::L2norm<double, THREExTHREE>(
+        const typename fftw_interface<double>::complex *__restrict__ a);
+
+template double kspace<FFTW, SMOOTH>::L2norm<float, ONE>(
+        const typename fftw_interface<float>::complex *__restrict__ a);
+template double kspace<FFTW, SMOOTH>::L2norm<float, THREE>(
+        const typename fftw_interface<float>::complex *__restrict__ a);
+template double kspace<FFTW, SMOOTH>::L2norm<float, THREExTHREE>(
+        const typename fftw_interface<float>::complex *__restrict__ a);
+template double kspace<FFTW, SMOOTH>::L2norm<double, ONE>(
+        const typename fftw_interface<double>::complex *__restrict__ a);
+template double kspace<FFTW, SMOOTH>::L2norm<double, THREE>(
+        const typename fftw_interface<double>::complex *__restrict__ a);
+template double kspace<FFTW, SMOOTH>::L2norm<double, THREExTHREE>(
+        const typename fftw_interface<double>::complex *__restrict__ a);
+
 template void kspace<FFTW, SMOOTH>::force_divfree<float>(
        typename fftw_interface<float>::complex *__restrict__ a);
 template void kspace<FFTW, SMOOTH>::force_divfree<double>(

bfps/cpp/kspace.hpp

@@ -114,6 +114,12 @@ class kspace
                 const hid_t group,
                 const std::string dset_name,
                 const hsize_t toffset);
+
+        template <typename rnumber,
+                  field_components fc>
+        double L2norm(
+                const rnumber(* __restrict__ a)[2]);
+
         template <class func_type>
         void CLOOP(func_type expression)
         {