diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index a29b33a8a07bdf57b593b5cc9bd873c24ae87ee1..4256164b7eb2264e2912c10f95f11ef44c6acab2 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -59,6 +59,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
hid_t group;
hid_t Cspace, Cdset;
int ndims;
+ DEBUG_MSG("about to grow datasets\\n");
// store kspace information
Cdset = H5Dopen(stat_file, "/kspace/kshell", H5P_DEFAULT);
Cspace = H5Dget_space(Cdset);
@@ -77,7 +78,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
H5Dwrite(Cdset, H5T_NATIVE_INT64, H5S_ALL, H5S_ALL, H5P_DEFAULT, fs->nshell);
H5Dclose(Cdset);
Cdset = H5Dopen(stat_file, "/kspace/kM", H5P_DEFAULT);
- H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &fs->kM);
+ H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &fs->kMspec);
H5Dclose(Cdset);
Cdset = H5Dopen(stat_file, "/kspace/dk", H5P_DEFAULT);
H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &fs->dk);
@@ -85,6 +86,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
group = H5Gopen(stat_file, "/statistics", H5P_DEFAULT);
H5Ovisit(group, H5_INDEX_NAME, H5_ITER_NATIVE, grow_statistics_dataset, NULL);
H5Gclose(group);
+ DEBUG_MSG("finished growing datasets\\n");
//endcpp
"""
self.style = {}
@@ -308,6 +310,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.fluid_start += """
//begincpp
char fname[512];
+ DEBUG_MSG("about to allocate fluid_solver\\n");
fs = new fluid_solver<{0}>(
simname,
nx, ny, nz,
@@ -322,6 +325,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
strncpy(fs->forcing_type, forcing_type, 128);
fs->iteration = iteration;
fs->read('v', 'c');
+ DEBUG_MSG("finished reading initial condition\\n");
if (fs->cd->myrank == 0)
{{
H5T_field_complex = H5Tcreate(H5T_COMPOUND, sizeof(tmp_complex_type));
diff --git a/bfps/cpp/fluid_solver_base.cpp b/bfps/cpp/fluid_solver_base.cpp
index a13f8ef29a248daf9ebd6668f1ef6a64a426ba4e..717eb77349d78a88516bca93dc829c844834818b 100644
--- a/bfps/cpp/fluid_solver_base.cpp
+++ b/bfps/cpp/fluid_solver_base.cpp
@@ -24,7 +24,7 @@
-#define NDEBUG
+//#define NDEBUG
#include <cassert>
#include <cmath>
@@ -70,7 +70,7 @@ void fluid_solver_base<rnumber>::cospectrum(cnumber *a, cnumber *b, double *spec
int tmp_int;
CLOOP_K2_NXMODES(
this,
- if (k2 <= this->kM2)
+ if (k2 <= this->kMspec2)
{
tmp_int = int(sqrt(k2)/this->dk)*9;
for (int i=0; i<3; i++)
@@ -99,7 +99,7 @@ void fluid_solver_base<rnumber>::cospectrum(cnumber *a, cnumber *b, double *spec
int tmp_int;
CLOOP_K2_NXMODES(
this,
- if (k2 <= this->kM2)
+ if (k2 <= this->kMspec2)
{
factor = nxmodes*pow(k2, k2exponent);
tmp_int = int(sqrt(k2)/this->dk)*9;
@@ -267,9 +267,9 @@ fluid_solver_base<R>::fluid_solver_base( \
{ \
/* HL07 smooth filter */ \
case 1: \
- this->kMx = this->dkx*(int(2*this->rd->sizes[2] / 5)); \
- this->kMy = this->dky*(int(2*this->rd->sizes[1] / 5)); \
- this->kMz = this->dkz*(int(2*this->rd->sizes[0] / 5)); \
+ this->kMx = this->dkx*(int(this->rd->sizes[2] / 2)-1); \
+ this->kMy = this->dky*(int(this->rd->sizes[1] / 2)-1); \
+ this->kMz = this->dkz*(int(this->rd->sizes[0] / 2)-1); \
break; \
default: \
this->kMx = this->dkx*(int(this->rd->sizes[2] / 3)-1); \
@@ -298,6 +298,17 @@ fluid_solver_base<R>::fluid_solver_base( \
if (this->kM < this->kMy) this->kM = this->kMy; \
if (this->kM < this->kMz) this->kM = this->kMz; \
this->kM2 = this->kM * this->kM; \
+ switch(this->dealias_type) \
+ { \
+ case 1: \
+ this->kMspec = 4*this->kM/5; \
+ this->kMspec2 = this->kMspec*this->kMspec; \
+ break; \
+ default: \
+ this->kMspec = this->kM; \
+ this->kMspec2 = this->kM2; \
+ break; \
+ } \
this->dk = this->dkx; \
if (this->dk > this->dky) this->dk = this->dky; \
if (this->dk > this->dkz) this->dk = this->dkz; \
@@ -306,7 +317,10 @@ fluid_solver_base<R>::fluid_solver_base( \
"kM = %g, kM2 = %g, dk = %g, dk2 = %g\n", \
this->kM, this->kM2, this->dk, this->dk2); \
/* spectra stuff */ \
- this->nshells = int(this->kM / this->dk) + 2; \
+ this->nshells = int(this->kMspec / this->dk) + 2; \
+ DEBUG_MSG( \
+ "kMspec = %g, kMspec2 = %g, nshells = %ld\n", \
+ this->kMspec, this->kMspec2, this->nshells); \
this->kshell = new double[this->nshells]; \
std::fill_n(this->kshell, this->nshells, 0.0); \
this->nshell = new int64_t[this->nshells]; \
@@ -348,7 +362,6 @@ fluid_solver_base<R>::fluid_solver_base( \
template<> \
fluid_solver_base<R>::~fluid_solver_base() \
{ \
- DEBUG_MSG("entered ~fluid_solver_base\n"); \
delete[] this->kshell; \
delete[] this->nshell; \
\
@@ -358,7 +371,6 @@ fluid_solver_base<R>::~fluid_solver_base() \
\
delete this->cd; \
delete this->rd; \
- DEBUG_MSG("exiting ~fluid_solver_base\n"); \
} \
\
template<> \
@@ -483,7 +495,6 @@ void fluid_solver_base<R>::symmetrize(FFTW(complex) *data, const int howmany) \
(*(data + howmany*((yy - this->cd->starts[0])*this->cd->sizes[1] + ii)*this->cd->sizes[2] + cc))[imag_comp]; \
if (ranksrc != rankdst) \
{ \
- DEBUG_MSG("inside fluid_solver_base::symmetrize, about to send/recv data\n"); \
if (this->cd->myrank == ranksrc) \
MPI_Send((void*)buffer, \
howmany*this->cd->sizes[1], MPI_CNUM, rankdst, yy, \
@@ -492,7 +503,6 @@ void fluid_solver_base<R>::symmetrize(FFTW(complex) *data, const int howmany) \
MPI_Recv((void*)buffer, \
howmany*this->cd->sizes[1], MPI_CNUM, ranksrc, yy, \
this->cd->comm, mpistatus); \
- DEBUG_MSG("inside fluid_solver_base::symmetrize, after send/recv data\n"); \
} \
if (this->cd->myrank == rankdst) \
{ \
diff --git a/bfps/cpp/fluid_solver_base.hpp b/bfps/cpp/fluid_solver_base.hpp
index 4cf7c2725ecf887932973ec97322a1a160ac4ee2..f7a3809aa8fd5fd8f964a6c040e6d2d4e42835ea 100644
--- a/bfps/cpp/fluid_solver_base.hpp
+++ b/bfps/cpp/fluid_solver_base.hpp
@@ -62,6 +62,7 @@ class fluid_solver_base
/* mode and dealiasing information */
int dealias_type;
double kMx, kMy, kMz, kM, kM2;
+ double kMspec, kMspec2;
double *kx, *ky, *kz;
std::unordered_map<int, double> Fourier_filter;
double *kshell;
diff --git a/bfps/fluid_base.py b/bfps/fluid_base.py
index a78f30e1ddc86ece23b54d6fcd913ba5911a96e7..facddc7e56c00e2c620b1607a1f7f31c9588b635 100644
--- a/bfps/fluid_base.py
+++ b/bfps/fluid_base.py
@@ -361,9 +361,9 @@ class fluid_particle_base(bfps.code):
def get_kspace(self):
kspace = {}
if self.parameters['dealias_type'] == 1:
- kMx = self.parameters['dkx']*(2*self.parameters['nx']//5)
- kMy = self.parameters['dky']*(2*self.parameters['ny']//5)
- kMz = self.parameters['dkz']*(2*self.parameters['nz']//5)
+ kMx = self.parameters['dkx']*(4*(self.parameters['nx']//2 - 1)//5)
+ kMy = self.parameters['dky']*(4*(self.parameters['ny']//2 - 1)//5)
+ kMz = self.parameters['dkz']*(4*(self.parameters['nz']//2 - 1)//5)
else:
kMx = self.parameters['dkx']*(self.parameters['nx']//3 - 1)
kMy = self.parameters['dky']*(self.parameters['ny']//3 - 1)