diff --git a/bfps/cpp/field.cpp b/bfps/cpp/field.cpp index e4e08f54eaef12459f80df44ea99565e9a1c6461..489903b715b14f47a0f7deba6884b3a899bff664 100644 --- a/bfps/cpp/field.cpp +++ b/bfps/cpp/field.cpp @@ -1077,17 +1077,17 @@ void field<rnumber, be, fc>::symmetrize() fftw_interface<rnumber>::free(buffer); delete mpistatus; /* put asymmetric data to 0 */ - /*if (this->clayout->myrank == this->clayout->rank[0][this->clayout->sizes[0]/2]) + if (this->clayout->myrank == this->clayout->rank[0][this->clayout->sizes[0]/2]) { - tindex = ncomp(fc)*(this->clayout->sizes[0]/2 - this->clayout->starts[0])*this->clayout->sizes[1]*this->clayout->sizes[2]; + ptrdiff_t tindex = ncomp(fc)*(this->clayout->sizes[0]/2 - this->clayout->starts[0])*this->clayout->sizes[1]*this->clayout->sizes[2]; for (ii = 0; ii < this->clayout->sizes[1]; ii++) { std::fill_n((rnumber*)(data + tindex), ncomp(fc)*2*this->clayout->sizes[2], 0.0); tindex += ncomp(fc)*this->clayout->sizes[2]; } } - tindex = ncomp(fc)*(); - std::fill_n((rnumber*)(data + tindex), ncomp(fc)*2, 0.0);*/ + //tindex = ncomp(fc)*(); + //std::fill_n((rnumber*)(data + tindex), ncomp(fc)*2, 0.0); } template <typename rnumber, @@ -1184,7 +1184,7 @@ double field<rnumber, be, fc>::L2norm( &m2, 1, MPI_DOUBLE, MPI_SUM, this->comm); - return m2 / this->npoints; + return sqrt(m2 / this->npoints); } } diff --git a/bfps/cpp/full_code/field_test.cpp b/bfps/cpp/full_code/field_test.cpp index b07f9b393a6aa9d10ea86d2b8e4e09d51cab48e2..b6c65a4a036f9a9fe31ee213d691937a6af4db33 100644 --- a/bfps/cpp/full_code/field_test.cpp +++ b/bfps/cpp/full_code/field_test.cpp @@ -76,9 +76,13 @@ int field_test<rnumber>::do_work(void) }); *scal_field_alt = scal_field->get_rdata(); + double L2r = scal_field->L2norm(kk); scal_field->dft(); + double L2c = scal_field->L2norm(kk); scal_field->ift(); scal_field->normalize(); + DEBUG_MSG("L2r = %g, L2c = %g\n", + L2r, L2c / scal_field->npoints); double max_error = 0; scal_field->RLOOP( @@ -93,6 +97,16 @@ int field_test<rnumber>::do_work(void) DEBUG_MSG("maximum error is %g\n", max_error); + scal_field->dft(); + kk->template dealias<rnumber, ONE>(scal_field->get_cdata()); + scal_field->symmetrize(); + scal_field->normalize(); + L2c = scal_field->L2norm(kk); + scal_field->ift(); + L2r = scal_field->L2norm(kk); + DEBUG_MSG("L2r = %g, L2c = %g\n", + L2r, L2c); + // deallocate delete kk; delete scal_field; diff --git a/bfps/cpp/kspace.cpp b/bfps/cpp/kspace.cpp index e3dabfe98de8d25983a3da3c5cb3b70c86675d93..72f8d70d0bebed55bf1513473e13f5100c216b39 100644 --- a/bfps/cpp/kspace.cpp +++ b/bfps/cpp/kspace.cpp @@ -616,14 +616,12 @@ double kspace<be, dt>::L2norm( 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++){ + for (hsize_t i=0; i<ncomp(fc); i++){ 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])); + (a[ncomp(fc)*cindex + i][0] * a[ncomp(fc)*cindex + i][0]) + + (a[ncomp(fc)*cindex + i][1] * a[ncomp(fc)*cindex + i][1])); } } }); @@ -636,7 +634,7 @@ double kspace<be, dt>::L2norm( &L2, 1, MPI_DOUBLE, MPI_SUM, this->layout->comm); - return L2 * this->dkx * this->dky * this->dkz; + return sqrt(L2 * this->dkx * this->dky * this->dkz); }