Skip to content
GitLab
Commit a78310d6 authored by Chichi Lalescu's avatar Chichi Lalescu
Browse files

plain test works with double precision

parent 49cd5bf1
Hide whitespace changes
Inline Side-by-side
bfps/NavierStokes.py
View file @ a78310d6
......@@ -32,11 +32,13 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self,
name = 'NavierStokes',
work_dir = './',
simname = 'test'):
simname = 'test',
fluid_precision = 'single'):
super(NavierStokes, self).__init__(
name = name,
work_dir = work_dir,
simname = simname)
simname = simname,
dtype = fluid_precision)
self.file_datasets_create = """
//begincpp
std::string temp_string;
......@@ -209,13 +211,17 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
"""
return None
def fill_up_fluid_code(self):
if self.dtype == np.float32:
C_dtype = 'float'
else:
C_dtype = 'double'
self.fluid_includes += '#include <cstring>\n'
self.fluid_variables += 'fluid_solver<float> *fs;\n'
self.fluid_variables += 'fluid_solver<{0}> *fs;\n'.format(C_dtype)
self.write_fluid_stats()
self.fluid_start += """
//begincpp
char fname[512];
fs = new fluid_solver<float>(
fs = new fluid_solver<{0}>(
simname,
nx, ny, nz,
dkx, dky, dkz);
......@@ -228,7 +234,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
fs->iteration = iteration;
fs->read('v', 'c');
//endcpp
"""
""".format(C_dtype)
self.fluid_loop += """
//begincpp
fs->step(dt);
......@@ -250,11 +256,15 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
smoothness = 1,
integration_steps = 2,
kcut = 'fs->kM'):
if self.dtype == np.float32:
C_dtype = 'float'
else:
C_dtype = 'double'
self.parameters['neighbours{0}'.format(self.particle_species)] = neighbours
self.parameters['smoothness{0}'.format(self.particle_species)] = smoothness
self.parameters['kcut{0}'.format(self.particle_species)] = kcut
self.parameters['integration_steps{0}'.format(self.particle_species)] = integration_steps
self.particle_variables += 'tracers<float> *ps{0};'.format(self.particle_species)
self.particle_variables += 'tracers<{0}> *ps{1};'.format(C_dtype, self.particle_species)
self.file_datasets_create += ('tmp_dspace = H5::DataSpace(2, dims, maxdims);\n' +
'temp_string = std::string("/particles/") + std::string(ps{0}->name);\n' +
'group = data_file.openGroup(temp_string);\n' +
......@@ -286,16 +296,16 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut) +
'fs->ift_velocity();\n' +
'ps{0}->update_field();\n').format(self.particle_species)
self.particle_start += ('sprintf(fname, "tracers{0}");\n' +
'ps{0} = new tracers<float>(\n' +
self.particle_start += ('sprintf(fname, "tracers{1}");\n' +
'ps{1} = new tracers<{0}>(\n' +
'fname, fs,\n' +
'nparticles,\n' +
'neighbours{0}, smoothness{0}, niter_part, integration_steps{0},\n' +
'neighbours{1}, smoothness{1}, niter_part, integration_steps{1},\n' +
'fs->ru);\n' +
'ps{0}->dt = dt;\n' +
'ps{0}->iteration = iteration;\n' +
'ps{1}->dt = dt;\n' +
'ps{1}->iteration = iteration;\n' +
update_field +
'ps{0}->read(&data_file);\n').format(self.particle_species)
'ps{1}->read(&data_file);\n').format(C_dtype, self.particle_species)
self.particle_loop += (update_field +
'ps{0}->step();\n' +
'if (ps{0}->iteration % niter_part == 0)\n' +
......@@ -347,16 +357,15 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
# plot particles
if particles_on and self.particle_species > 0:
traj = self.read_traj()
fig = plt.figure(figsize = (12, 12))
a = fig.add_subplot(111, projection = '3d')
for t in range(self.parameters['nparticles']):
a.plot(traj['state'][0, :, t, 0],
traj['state'][0, :, t, 1],
traj['state'][0, :, t, 2], color = 'blue')
a.plot(traj['state'][1, :, t, 0],
traj['state'][1, :, t, 1],
traj['state'][1, :, t, 2], color = 'red', dashes = (1, 1))
a.plot(self.trajectories[0][:, t, 0],
self.trajectories[0][:, t, 1],
self.trajectories[0][:, t, 2], color = 'blue')
a.plot(self.trajectories[1][:, t, 0],
self.trajectories[1][:, t, 1],
self.trajectories[1][:, t, 2], color = 'red', dashes = (1, 1))
fig.savefig('traj.pdf', format = 'pdf')
return None
def compute_statistics(self, iter0 = 0):
......@@ -385,7 +394,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.statistics['diss' + suffix])**.5
self.statistics['kshell'] = data_file['kspace/kshell'].value
self.statistics['kM'] = data_file['kspace/kM'].value
self.trajectories = self.read_traj()
self.trajectories = [data_file['particles/' + key + '/state'].value
for key in data_file['particles'].keys()]
return None
def plot_spectrum(
self,
......
bfps/cpp/fftw_tools.cpp
View file @ a78310d6
......@@ -41,7 +41,7 @@ int copy_complex_array(
min_fast_dim =
(fi->sizes[2] > fo->sizes[2]) ?
fo->sizes[2] : fi->sizes[2];
MPI_Datatype MPI_CNUM = (sizeof(rnumber) == 4) ? MPI_COMPLEX8 : MPI_COMPLEX16;
MPI_Datatype MPI_CNUM = (sizeof(rnumber) == 4) ? MPI_COMPLEX : MPI_COMPLEX16;
// clean up destination, in case we're padding with zeros
// (even if only for one dimension)
......@@ -164,7 +164,7 @@ int get_descriptors_3D(
field_descriptor<rnumber> **fc)
{
MPI_Datatype MPI_RNUM = (sizeof(rnumber) == 4) ? MPI_FLOAT : MPI_DOUBLE;
MPI_Datatype MPI_CNUM = (sizeof(rnumber) == 4) ? MPI_COMPLEX8 : MPI_COMPLEX16;
MPI_Datatype MPI_CNUM = (sizeof(rnumber) == 4) ? MPI_COMPLEX : MPI_COMPLEX16;
int ntmp[3];
ntmp[0] = n0;
ntmp[1] = n1;
......
bfps/cpp/field_descriptor.cpp
View file @ a78310d6
......@@ -19,9 +19,9 @@
************************************************************************/
// code is generally compiled via setuptools, therefore NDEBUG is present
//#ifdef NDEBUG
//#undef NDEBUG
//#endif//NDEBUG
#ifdef NDEBUG
#undef NDEBUG
#endif//NDEBUG
#include <stdlib.h>
#include <algorithm>
......@@ -95,7 +95,14 @@ field_descriptor<rnumber>::field_descriptor(
tsubsizes[ndims-1] *= sizeof(rnumber);
tstarts[ndims-1] *= sizeof(rnumber);
if (this->mpi_dtype == MPI_COMPLEX ||
this->mpi_dtype == MPI_COMPLEX16)
this->mpi_dtype == MPI_C_DOUBLE_COMPLEX)
{
tsizes[ndims-1] *= 2;
tsubsizes[ndims-1] *= 2;
tstarts[ndims-1] *= 2;
}
if (this->mpi_dtype == MPI_DOUBLE ||
this->mpi_dtype == MPI_C_DOUBLE_COMPLEX)
{
tsizes[ndims-1] *= 2;
tsubsizes[ndims-1] *= 2;
......@@ -236,8 +243,8 @@ int field_descriptor<rnumber>::read(
if (sizeof(rnumber)==4)
ttype = MPI_COMPLEX;
else if (sizeof(rnumber)==8)
ttype = MPI_DOUBLE_COMPLEX;
DEBUG_MSG("aloha 00\n");
ttype = MPI_C_DOUBLE_COMPLEX;
DEBUG_MSG("aloha from field_descriptor read\n");
char representation[] = "native";
if (this->subsizes[0] > 0)
{
......@@ -245,11 +252,11 @@ int field_descriptor<rnumber>::read(
MPI_Info_create(&info);
MPI_File f;
ptrdiff_t read_size = this->local_size*sizeof(rnumber);
DEBUG_MSG("aloha %ld\n", read_size);
DEBUG_MSG("read size is %ld\n", read_size);
char ffname[200];
if (this->mpi_dtype == ttype)
read_size *= 2;
DEBUG_MSG("aloha %ld\n", read_size);
DEBUG_MSG("read size is %ld\n", read_size);
sprintf(ffname, "%s", fname);
MPI_File_open(
......@@ -285,7 +292,7 @@ int field_descriptor<rnumber>::write(
if (sizeof(rnumber)==4)
ttype = MPI_COMPLEX;
else if (sizeof(rnumber)==8)
ttype = MPI_DOUBLE_COMPLEX;
ttype = MPI_C_DOUBLE_COMPLEX;
char representation[] = "native";
if (this->subsizes[0] > 0)
{
......
bfps/cpp/fluid_solver.cpp
View file @ a78310d6
......@@ -454,12 +454,12 @@ FLUID_SOLVER_DEFINITIONS(
fftwf_complex,
MPI_FLOAT,
MPI_COMPLEX)
//FLUID_SOLVER_DEFINITIONS(
// FFTW_MANGLE_DOUBLE,
// double,
// fftw_complex,
// MPI_DOUBLE,
// MPI_C_DOUBLE_COMPLEX)
FLUID_SOLVER_DEFINITIONS(
FFTW_MANGLE_DOUBLE,
double,
fftw_complex,
MPI_DOUBLE,
MPI_C_DOUBLE_COMPLEX)
/*****************************************************************************/
......@@ -467,5 +467,6 @@ FLUID_SOLVER_DEFINITIONS(
/*****************************************************************************/
/* finally, force generation of code for single precision */
template class fluid_solver<float>;
template class fluid_solver<double>;
/*****************************************************************************/
bfps/cpp/fluid_solver_base.cpp
View file @ a78310d6
......@@ -472,12 +472,12 @@ FLUID_SOLVER_BASE_DEFINITIONS(
fftwf_complex,
MPI_FLOAT,
MPI_COMPLEX)
//FLUID_SOLVER_BASE_DEFINITIONS(
// FFTW_MANGLE_DOUBLE,
// double,
// fftw_complex,
// MPI_DOUBLE,
// MPI_C_DOUBLE_COMPLEX)
FLUID_SOLVER_BASE_DEFINITIONS(
FFTW_MANGLE_DOUBLE,
double,
fftw_complex,
MPI_DOUBLE,
MPI_C_DOUBLE_COMPLEX)
/*****************************************************************************/
......@@ -485,5 +485,6 @@ FLUID_SOLVER_BASE_DEFINITIONS(
/*****************************************************************************/
/* finally, force generation of code for single precision */
template class fluid_solver_base<float>;
template class fluid_solver_base<double>;
/*****************************************************************************/
bfps/cpp/slab_field_particles.cpp
View file @ a78310d6
......@@ -820,4 +820,5 @@ void slab_field_particles<rnumber>::rFFTW_to_buffered(rnumber *src, rnumber *dst
/*****************************************************************************/
/* finally, force generation of code for single precision */
template class slab_field_particles<float>;
template class slab_field_particles<double>;
/*****************************************************************************/
bfps/cpp/tracers.cpp
View file @ a78310d6
......@@ -39,7 +39,7 @@ void tracers<rnumber>::jump_estimate(double *jump)
double *tjump = new double[this->nparticles];
int *xg = new int[this->array_size];
double *xx = new double[this->array_size];
float *vel = this->data + this->buffer_size;
rnumber *vel = this->data + this->buffer_size;
double tmp[3];
/* get grid coordinates */
this->get_grid_coordinates(this->state, xg, xx);
......@@ -75,7 +75,7 @@ void tracers<rnumber>::get_rhs(double *x, double *y)
/* get grid coordinates */
int *xg = new int[this->array_size];
double *xx = new double[this->array_size];
float *vel = this->data + this->buffer_size;
rnumber *vel = this->data + this->buffer_size;
this->get_grid_coordinates(x, xg, xx);
/* perform interpolation */
for (int p=0; p<this->nparticles; p++) if (this->fs->rd->myrank == this->computing[p])
......@@ -146,12 +146,19 @@ TRACERS_DEFINITIONS(
fftwf_complex,
MPI_FLOAT,
MPI_COMPLEX)
TRACERS_DEFINITIONS(
FFTW_MANGLE_DOUBLE,
double,
fftw_complex,
MPI_DOUBLE,
MPI_C_DOUBLE_COMPLEX)
/*****************************************************************************/
/*****************************************************************************/
/* finally, force generation of code for single precision */
/* finally, force generation of code */
template class tracers<float>;
template class tracers<double>;
/*****************************************************************************/
bfps/fluid_base.py
View file @ a78310d6
......@@ -31,12 +31,20 @@ class fluid_particle_base(bfps.code):
self,
name = 'solver',
work_dir = './',
simname = 'test'):
simname = 'test',
dtype = np.float32):
super(fluid_particle_base, self).__init__(
work_dir = work_dir,
simname = simname)
self.name = name
self.particle_species = 0
if dtype in [np.float32, np.float64]:
self.dtype = dtype
elif dtype in ['single', 'double']:
if dtype == 'single':
self.dtype = np.float32
elif dtype == 'double':
self.dtype = np.float64
self.parameters['dkx'] = 1.0
self.parameters['dky'] = 1.0
self.parameters['dkz'] = 1.0
......@@ -128,7 +136,7 @@ class fluid_particle_base(bfps.code):
self.simname + '_r' + field + '_i{0:0>5x}'.format(iteration))
return np.memmap(
filename,
dtype = np.float32,
dtype = self.dtype,
shape = (self.parameters['nz'],
self.parameters['ny'],
self.parameters['nx'], 3))
......@@ -147,7 +155,7 @@ class fluid_particle_base(bfps.code):
self.parameters['nz'],
self.parameters['ny'],
self.parameters['nx']),
dtype = Rdata.dtype)
dtype = self.dtype)
for i in range(3):
new_data[i] = Rdata[..., i]
if type(ofile) == type(None):
......@@ -201,9 +209,9 @@ class fluid_particle_base(bfps.code):
iteration = 0,
field_name = 'vorticity',
write_to_file = False):
if precision == 'single':
if self.dtype == np.float32:
dtype = np.complex64
elif precision == 'double':
elif self.dtype == np.float64:
dtype = np.complex128
np.random.seed(rseed)
Kdata00 = bfps.tools.generate_data_3D(
......@@ -271,30 +279,6 @@ class fluid_particle_base(bfps.code):
self.work_dir,
"tracers{0}_state_i{1:0>5x}".format(species, iteration)))
return data
def read_spec(
self,
field = 'velocity'):
k = np.fromfile(
os.path.join(
self.work_dir,
self.simname + '_kshell'),
dtype = np.float64)
spec_dtype = np.dtype([('iteration', np.int32),
('val', np.float64, k.shape[0])])
spec = np.fromfile(
os.path.join(
self.work_dir,
self.simname + '_' + field + '_spec'),
dtype = spec_dtype)
return k, spec
def read_traj(self):
if self.particle_species == 0:
return None
with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r') as data_file:
traj_list = []
for t in range(self.particle_species):
traj_list.append(data_file['particles/tracers{0}/state'.format(t)].value)
return traj_list
def generate_initial_condition(self):
self.generate_vector_field(write_to_file = True)
for species in range(self.particle_species):
......
bfps/fluid_converter.py
View file @ a78310d6
......@@ -31,11 +31,13 @@ class fluid_converter(bfps.fluid_base.fluid_particle_base):
self,
name = 'fluid_converter',
work_dir = './',
simname = 'test'):
simname = 'test',
dtype = np.float32):
super(fluid_converter, self).__init__(
name = name,
work_dir = work_dir,
simname = simname)
simname = simname,
dtype = dtype)
self.parameters['write_rvelocity'] = 1
self.parameters['write_rvorticity'] = 1
self.parameters['fluid_name'] = 'test'
......@@ -43,12 +45,16 @@ class fluid_converter(bfps.fluid_base.fluid_particle_base):
self.finalize_code()
return None
def fill_up_fluid_code(self):
if self.dtype == np.float32:
C_dtype = 'float'
else:
C_dtype = 'double'
self.fluid_includes += '#include <cstring>\n'
self.fluid_variables += ('double t;\n' +
'fluid_solver<float> *fs;\n')
'fluid_solver<{0}> *fs;\n').format(C_dtype)
self.fluid_definitions += """
//begincpp
void do_conversion(fluid_solver<float> *bla)
void do_conversion(fluid_solver<{0}> *bla)
{
bla->read('v', 'c');
if (write_rvelocity)
......@@ -57,17 +63,17 @@ class fluid_converter(bfps.fluid_base.fluid_particle_base):
bla->write('v', 'r');
}
//endcpp
"""
""".format(C_dtype)
self.fluid_start += """
//begincpp
fs = new fluid_solver<float>(
fs = new fluid_solver<{0}>(
fluid_name,
nx, ny, nz,
dkx, dky, dkz);
fs->iteration = iteration;
do_conversion(fs);
//endcpp
"""
""".format(C_dtype)
self.fluid_loop += """
//begincpp
fs->iteration++;
......
bfps/fluid_resize.py
View file @ a78310d6
......@@ -21,16 +21,20 @@
import bfps
import bfps.fluid_base
import numpy as np
class fluid_resize(bfps.fluid_base.fluid_particle_base):
def __init__(
self,
name = 'fluid_converter',
work_dir = './',
simname = 'test'):
simname = 'test',
dtype = np.float32):
super(fluid_resize, self).__init__(
name = name,
work_dir = work_dir,
simname = simname)
simname = simname,
dtype = dtype)
self.parameters['dst_iter'] = 0
self.parameters['dst_nx'] = 32
self.parameters['dst_ny'] = 32
......@@ -45,16 +49,20 @@ class fluid_resize(bfps.fluid_base.fluid_particle_base):
def fill_up_fluid_code(self):
self.fluid_includes += '#include <cstring>\n'
self.fluid_includes += '#include "fftw_tools.hpp"\n'
if self.dtype == np.float32:
C_dtype = 'float'
else:
C_dtype = 'double'
self.fluid_variables += ('double t;\n' +
'fluid_solver<float> *fs0, *fs1;\n')
'fluid_solver<' + C_dtype + '> *fs0, *fs1;\n')
self.fluid_start += """
//begincpp
char fname[512];
fs0 = new fluid_solver<float>(
fs0 = new fluid_solver<{0}>(
simname,
nx, ny, nz,
dkx, dky, dkz);
fs1 = new fluid_solver<float>(
fs1 = new fluid_solver<{0}>(
dst_simname,
dst_nx, dst_ny, dst_nz,
dst_dkx, dst_dky, dst_dkz);
......@@ -74,7 +82,7 @@ class fluid_resize(bfps.fluid_base.fluid_particle_base):
DEBUG_MSG("new field %d %g %g\\n", fs1->iteration, a, b);
niter_todo = 0;
//endcpp
"""
""".format(C_dtype)
self.fluid_end += """
//begincpp
delete fs0;
......
test.py
View file @ a78310d6
......@@ -135,7 +135,8 @@ def double(opt):
new_simname = 'N{0:0>3x}'.format(opt.n*2)
c = fluid_resize(
work_dir = opt.work_dir + '/resize',
simname = old_simname)
simname = old_simname,
dtype = opt.precision)
c.parameters['nx'] = opt.n
c.parameters['ny'] = opt.n
c.parameters['nz'] = opt.n
......@@ -165,7 +166,9 @@ def NSlaunch(
opt,
nu = None,
tracer_state_file = None):
c = bfps.NavierStokes(work_dir = opt.work_dir)
c = bfps.NavierStokes(
work_dir = opt.work_dir,
fluid_precision = opt.precision)
assert((opt.nsteps % 4) == 0)
c.parameters['nx'] = opt.n
c.parameters['ny'] = opt.n
......@@ -473,17 +476,24 @@ def plain(opt):
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('--run', dest = 'run', action = 'store_true')
parser.add_argument('--double', dest = 'double', action = 'store_true')
parser.add_argument('--initialize', dest = 'initialize', action = 'store_true')
parser.add_argument('--convergence', dest = 'convergence', action = 'store_true')
parser.add_argument('--plain', dest = 'plain', action = 'store_true')
parser.add_argument('--io', dest = 'io', action = 'store_true')
parser.add_argument('--iteration', type = int, dest = 'iteration', default = 0)
parser.add_argument('--ncpu', type = int, dest = 'ncpu', default = 2)
parser.add_argument('--nsteps', type = int, dest = 'nsteps', default = 16)
parser.add_argument('--njobs', type = int, dest = 'njobs', default = 1)
parser.add_argument('-n', type = int, dest = 'n', default = 64)
parser.add_argument('--wd', type = str, dest = 'work_dir', default = 'data')
parser.add_argument('-n',
type = int, dest = 'n', default = 64)
parser.add_argument('--iteration',
type = int, dest = 'iteration', default = 0)
parser.add_argument('--ncpu',
type = int, dest = 'ncpu', default = 2)
parser.add_argument('--nsteps',
type = int, dest = 'nsteps', default = 16)
parser.add_argument('--njobs',
type = int, dest = 'njobs', default = 1)
parser.add_argument('--wd',
type = str, dest = 'work_dir', default = 'data')
parser.add_argument('--precision',
type = str, dest = 'precision', default = 'single')
opt = parser.parse_args()
if opt.convergence:
convergence_test(opt)
......
tox_io.ini
View file @ a78310d6
[tox]
envlist = py27
envlist = py27, py34
[testenv]
whitelist_externals =
echo
......@@ -9,5 +9,6 @@ changedir =
{envtmpdir}
commands =
cp {toxinidir}/test.py {envtmpdir}
python test.py -n 32 --run --initialize --ncpu 2 --io
python test.py -n 32 --run --initialize --ncpu 2 --io --precision single --wd "data/single"
python test.py -n 32 --run --initialize --ncpu 2 --io --precision double --wd "data/double"
tox_plain.ini
View file @ a78310d6
......@@ -9,5 +9,6 @@ changedir =
{envtmpdir}
commands =
cp {toxinidir}/test.py {envtmpdir}
python test.py -n 32 --run --initialize --ncpu 2 --plain --nsteps 64
python test.py -n 32 --run --initialize --ncpu 2 --plain --nsteps 64 --precision single --wd "data/single"
python test.py -n 32 --run --initialize --ncpu 2 --plain --nsteps 64 --precision double --wd "data/double"
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment