Commit dd8ef14d authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/public-API' into develop

parents 5f87cc6a 5a0b8278
...@@ -24,22 +24,26 @@ ...@@ -24,22 +24,26 @@
import bfps
import bfps.fluid_base
import bfps.tools
import numpy as np import numpy as np
import pickle import pickle
import os import os
from ._fluid_base import _fluid_particle_base
class fluid_converter(bfps.fluid_base.fluid_particle_base): class FluidConvert(_fluid_particle_base):
"""This class is meant to be used for conversion of native DNS field
representations to real-space representations of velocity/vorticity
fields.
It may be superseeded by streamlined functionality in the future...
"""
def __init__( def __init__(
self, self,
name = 'fluid_converter', name = 'FluidConvert',
work_dir = './', work_dir = './',
simname = 'test', simname = 'test',
fluid_precision = 'single', fluid_precision = 'single',
use_fftw_wisdom = True): use_fftw_wisdom = True):
super(fluid_converter, self).__init__( _fluid_particle_base.__init__(
self,
name = name, name = name,
work_dir = work_dir, work_dir = work_dir,
simname = simname, simname = simname,
...@@ -86,4 +90,24 @@ class fluid_converter(bfps.fluid_base.fluid_particle_base): ...@@ -86,4 +90,24 @@ class fluid_converter(bfps.fluid_base.fluid_particle_base):
""" """
self.fluid_end += 'delete fs;\n' self.fluid_end += 'delete fs;\n'
return None return None
def specific_parser_arguments(
self,
parser):
_fluid_particle_base.specific_parser_arguments(self, parser)
parser.add_argument(
'--src-wd',
type = str,
dest = 'src_work_dir',
default = './')
parser.add_argument(
'--src-simname',
type = str,
dest = 'src_simname',
default = '')
parser.add_argument(
'--src-iteration',
type = int,
dest = 'src_iteration',
default = 0)
return None
...@@ -24,20 +24,27 @@ ...@@ -24,20 +24,27 @@
import bfps import os
import bfps.fluid_base import argparse
import numpy as np import numpy as np
class fluid_resize(bfps.fluid_base.fluid_particle_base): import bfps
from ._fluid_base import _fluid_particle_base
class FluidResize(_fluid_particle_base):
"""This class is meant to resize snapshots of DNS states to new grids.
Typical stuff for DNS of turbulence.
It will become superfluous when HDF5 is used for field I/O.
"""
def __init__( def __init__(
self, self,
name = 'fluid_resize', name = 'FluidResize',
work_dir = './', work_dir = './',
simname = 'test', simname = 'test',
dtype = np.float32, dtype = np.float32,
use_fftw_wisdom = False): use_fftw_wisdom = False):
super(fluid_resize, self).__init__( _fluid_particle_base.__init__(
self,
name = name, name = name,
work_dir = work_dir, work_dir = work_dir,
simname = simname, simname = simname,
...@@ -100,4 +107,50 @@ class fluid_resize(bfps.fluid_base.fluid_particle_base): ...@@ -100,4 +107,50 @@ class fluid_resize(bfps.fluid_base.fluid_particle_base):
//endcpp //endcpp
""" """
return None return None
def specific_parser_arguments(
self,
parser):
_fluid_particle_base.specific_parser_arguments(self, parser)
parser.add_argument(
'-m',
type = int,
dest = 'm',
default = 32,
metavar = 'M',
help = 'resize from N to M')
parser.add_argument(
'--src_wd',
type = str,
dest = 'src_work_dir',
required = True)
parser.add_argument(
'--src_iteration',
type = int,
dest = 'src_iteration',
required = True)
return None
def launch(
self,
args = [],
**kwargs):
parser = argparse.ArgumentParser('bfps ' + type(self).__name__)
self.add_parser_arguments(parser)
opt = parser.parse_args(args)
cmd_line_pars = vars(opt)
for k in ['dst_nx', 'dst_ny', 'dst_nz']:
if type(cmd_line_pars[k]) == type(None):
cmd_line_pars[k] = opt.m
self.pars_from_namespace(opt)
src_file = os.path.join(
opt.src_work_dir,
opt.src_simname + '_cvorticity_i{0:0>5x}'.format(opt.src_iteration))
read_file = os.path.join(
self.work_dir,
opt.src_simname + '_cvorticity_i{0:0>5x}'.format(opt.src_iteration))
if not os.path.exists(read_file):
os.symlink(src_file, read_file)
self.set_host_info(bfps.host_info)
self.write_par(iter0 = opt.src_iteration)
self.run(ncpu = opt.ncpu)
return None
import os
import argparse
import bfps
class Launcher:
def __init__(self):
self.base_class = bfps.NavierStokes
self.parser = argparse.ArgumentParser(prog = 'bfps')
self.parser.add_argument(
'-v', '--version',
action = 'version',
version = '%(prog)s ' + bfps.__version__)
self.parser.add_argument(
'-n',
type = int,
dest = 'n',
default = 32,
metavar = 'N',
help = 'code is run by default in a grid of NxNxN')
self.parser.add_argument(
'--run',
dest = 'run',
action = 'store_true')
self.parser.add_argument(
'--ncpu',
type = int, dest = 'ncpu',
default = 2)
self.parser.add_argument(
'--precision',
type = str, dest = 'precision',
default = 'single')
self.parser.add_argument(
'--simname',
type = str, dest = 'simname',
default = 'test')
self.parser.add_argument(
'--wd',
type = str, dest = 'work_dir',
default = './')
self.parser.add_argument(
'--njobs',
type = int, dest = 'njobs',
default = 1)
self.parser.add_argument(
'--QR-stats',
action = 'store_true',
dest = 'QR_stats',
help = 'add this option if you want to compute velocity gradient and QR stats')
self.parser.add_argument(
'--kMeta',
type = float,
dest = 'kMeta',
default = 2.0)
self.parser.add_argument(
'--dtfactor',
type = float,
dest = 'dtfactor',
default = 0.5,
help = 'dt is computed as DTFACTOR / N')
self.parser.add_argument(
'--environment',
type = str,
dest = 'environment',
default = '')
self.parser.add_argument(
'--src-simname',
type = str,
dest = 'src_simname',
default = '')
self.parser.add_argument(
'--src-iteration',
type = int,
dest = 'src_iteration',
default = 0)
self.parser.add_argument(
'--particle-rand-seed',
type = int,
dest = 'particle_rand_seed',
default = None)
c = self.base_class()
for k in sorted(c.parameters.keys()):
self.parser.add_argument(
'--{0}'.format(k),
type = type(c.parameters[k]),
dest = k,
default = None)
return None
def __call__(
self,
args = None):
opt = self.parser.parse_args(args)
if opt.environment != '':
bfps.host_info['environment'] = opt.environment
opt.work_dir = os.path.join(
os.path.realpath(opt.work_dir),
'N{0:0>4}'.format(opt.n))
c = self.base_class(
work_dir = opt.work_dir,
fluid_precision = opt.precision,
simname = opt.simname,
QR_stats_on = opt.QR_stats)
# with the default Lundgren forcing, I can estimate the dissipation
# with nondefault forcing, figure out the amplitude for this viscosity
# yourself
c.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
c.parameters['dt'] = (opt.dtfactor / opt.n)
if ((c.parameters['niter_todo'] % c.parameters['niter_out']) != 0):
c.parameters['niter_out'] = c.parameters['niter_todo']
if c.QR_stats_on:
# max_Q_estimate and max_R_estimate are just used for the 2D pdf
# therefore I just want them to be small multiples of mean trS2
# I'm already estimating the dissipation with kMeta...
meantrS2 = (opt.n//2 / opt.kMeta)**4 * c.parameters['nu']**2
c.parameters['max_Q_estimate'] = meantrS2
c.parameters['max_R_estimate'] = .4*meantrS2**1.5
# command line parameters will overwrite any defaults
cmd_line_pars = vars(opt)
for k in ['nx', 'ny', 'nz']:
if type(cmd_line_pars[k]) == type(None):
cmd_line_pars[k] = opt.n
for k in c.parameters.keys():
if k in cmd_line_pars.keys():
if not type(cmd_line_pars[k]) == type(None):
c.parameters[k] = cmd_line_pars[k]
c.fill_up_fluid_code()
c.finalize_code()
c.write_src()
c.set_host_info(bfps.host_info)
if opt.run:
if not os.path.exists(os.path.join(c.work_dir, c.simname + '.h5')):
c.write_par()
if c.parameters['nparticles'] > 0:
data = c.generate_tracer_state(species = 0, rseed = opt.particle_rand_seed)
for s in range(1, c.particle_species):
c.generate_tracer_state(species = s, data = data)
init_condition_file = os.path.join(
c.work_dir,
c.simname + '_cvorticity_i{0:0>5x}'.format(0))
if not os.path.exists(init_condition_file):
if len(opt.src_simname) > 0:
src_file = os.path.join(
c.work_dir,
opt.src_simname + '_cvorticity_i{0:0>5x}'.format(opt.src_iteration))
os.symlink(src_file, init_condition_file)
else:
c.generate_vector_field(
write_to_file = True,
spectra_slope = 2.0,
amplitude = 0.25)
c.run(ncpu = opt.ncpu,
njobs = opt.njobs)
return c
...@@ -27,12 +27,17 @@ ...@@ -27,12 +27,17 @@
import os import os
import numpy as np import numpy as np
import h5py import h5py
import argparse
import bfps import bfps
import bfps.fluid_base from ._fluid_base import _fluid_particle_base
import bfps.tools
class NavierStokes(bfps.fluid_base.fluid_particle_base): class NavierStokes(_fluid_particle_base):
"""Objects of this class can be used to generate production DNS codes.
Any functionality that users require should be available through this class,
in the sense that they can implement whatever they need by simply inheriting
this class.
"""
def __init__( def __init__(
self, self,
name = 'NavierStokes', name = 'NavierStokes',
...@@ -46,12 +51,26 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): ...@@ -46,12 +51,26 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.QR_stats_on = QR_stats_on self.QR_stats_on = QR_stats_on
self.frozen_fields = frozen_fields self.frozen_fields = frozen_fields
self.fftw_plan_rigor = fftw_plan_rigor self.fftw_plan_rigor = fftw_plan_rigor
super(NavierStokes, self).__init__( _fluid_particle_base.__init__(
self,
name = name, name = name,
work_dir = work_dir, work_dir = work_dir,
simname = simname, simname = simname,
dtype = fluid_precision, dtype = fluid_precision,
use_fftw_wisdom = use_fftw_wisdom) use_fftw_wisdom = use_fftw_wisdom)
self.parameters['nu'] = 0.1
self.parameters['fmode'] = 1
self.parameters['famplitude'] = 0.5
self.parameters['fk0'] = 1.5
self.parameters['fk1'] = 3.0
self.parameters['forcing_type'] = 'linear'
self.parameters['histogram_bins'] = 256
self.parameters['max_velocity_estimate'] = 1.0
self.parameters['max_vorticity_estimate'] = 1.0
self.parameters['QR2D_histogram_bins'] = 64
self.parameters['max_trS2_estimate'] = 1.0
self.parameters['max_Q_estimate'] = 1.0
self.parameters['max_R_estimate'] = 1.0
self.file_datasets_grow = """ self.file_datasets_grow = """
//begincpp //begincpp
std::string temp_string; std::string temp_string;
...@@ -715,10 +734,87 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): ...@@ -715,10 +734,87 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.parameters['nx']//2+1, self.parameters['nx']//2+1,
3)) 3))
def write_par(self, iter0 = 0): def write_par(self, iter0 = 0):
super(NavierStokes, self).write_par(iter0 = iter0) _fluid_particle_base.write_par(self, iter0 = iter0)
with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r+') as ofile: with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r+') as ofile:
kspace = self.get_kspace() kspace = self.get_kspace()
nshells = kspace['nshell'].shape[0] nshells = kspace['nshell'].shape[0]
for k in ['velocity', 'vorticity']:
time_chunk = 2**20//(8*3*self.parameters['nx'])
time_chunk = max(time_chunk, 1)
ofile.create_dataset('statistics/xlines/' + k,
(1, self.parameters['nx'], 3),
chunks = (time_chunk, self.parameters['nx'], 3),
maxshape = (None, self.parameters['nx'], 3),
dtype = self.dtype,
compression = 'gzip')
time_chunk = 2**20//(8*3*3*nshells)
time_chunk = max(time_chunk, 1)
ofile.create_dataset('statistics/spectra/' + k + '_' + k,
(1, nshells, 3, 3),
chunks = (time_chunk, nshells, 3, 3),
maxshape = (None, nshells, 3, 3),
dtype = np.float64,
compression = 'gzip')
time_chunk = 2**20//(8*4*10)
time_chunk = max(time_chunk, 1)
a = ofile.create_dataset('statistics/moments/' + k,
(1, 10, 4),
chunks = (time_chunk, 10, 4),
maxshape = (None, 10, 4),
dtype = np.float64,
compression = 'gzip')
time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
time_chunk = max(time_chunk, 1)
ofile.create_dataset('statistics/histograms/' + k,
(1,
self.parameters['histogram_bins'],
4),
chunks = (time_chunk,
self.parameters['histogram_bins'],
4),
maxshape = (None,
self.parameters['histogram_bins'],
4),
dtype = np.int64,
compression = 'gzip')
for s in range(self.particle_species):
time_chunk = 2**20 // (8*3*
self.parameters['nparticles']*
self.parameters['tracers{0}_integration_steps'.format(s)])
time_chunk = max(time_chunk, 1)
ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
(1,
self.parameters['tracers{0}_integration_steps'.format(s)],
self.parameters['nparticles'],
3),
maxshape = (None,
self.parameters['tracers{0}_integration_steps'.format(s)],
self.parameters['nparticles'],
3),
chunks = (time_chunk,
self.parameters['tracers{0}_integration_steps'.format(s)],
self.parameters['nparticles'],
3),
dtype = np.float64)
time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
time_chunk = max(time_chunk, 1)
ofile.create_dataset(
'/particles/tracers{0}/velocity'.format(s),
(1,
self.parameters['nparticles'],
3),
chunks = (time_chunk, self.parameters['nparticles'], 3),
maxshape = (None, self.parameters['nparticles'], 3),
dtype = np.float64)
if self.parameters['tracers{0}_acc_on'.format(s)]:
ofile.create_dataset(
'/particles/tracers{0}/acceleration'.format(s),
(1,
self.parameters['nparticles'],
3),
chunks = (time_chunk, self.parameters['nparticles'], 3),
maxshape = (None, self.parameters['nparticles'], 3),
dtype = np.float64)
if self.QR_stats_on: if self.QR_stats_on:
time_chunk = 2**20//(8*3*self.parameters['histogram_bins']) time_chunk = 2**20//(8*3*self.parameters['histogram_bins'])
time_chunk = max(time_chunk, 1) time_chunk = max(time_chunk, 1)
...@@ -827,33 +923,107 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): ...@@ -827,33 +923,107 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.fluid_start += update_fields self.fluid_start += update_fields
self.fluid_loop += update_fields self.fluid_loop += update_fields
return None return None
def specific_parser_arguments(
self,
parser):
_fluid_particle_base.specific_parser_arguments(self, parser)
parser.add_argument(
'--src-wd',
type = str,
dest = 'src_work_dir',
default = './')
parser.add_argument(
'--src-simname',
type = str,
dest = 'src_simname',
default = '')
parser.add_argument(
'--src-iteration',
type = int,
dest = 'src_iteration',
default = 0)
parser.add_argument(
'--precision',
type = str, dest = 'precision',
default = 'single')
parser.add_argument(
'--njobs',
type = int, dest = 'njobs',
default = 1)
parser.add_argument(
'--QR-stats',
action = 'store_true',
dest = 'QR_stats',
help = 'add this option if you want to compute velocity gradient and QR stats')
parser.add_argument(
'--kMeta',
type = float,
dest = 'kMeta',
default = 2.0)
parser.add_argument(
'--dtfactor',
type = float,
dest = 'dtfactor',
default = 0.5,
help = 'dt is computed as DTFACTOR / N')
parser.add_argument(
'--particle-rand-seed',
type = int,
dest = 'particle_rand_seed',
default = None)
return None
def launch(
self,
args = [],
**kwargs):
# with the default Lundgren forcing, I can estimate the dissipation
# with nondefault forcing, figure out the amplitude for this viscosity
# yourself
parser = argparse.ArgumentParser('bfps ' + type(self).__name__)
self.add_parser_arguments(parser)
opt = parser.parse_args(args)
self.QR_stats_on = opt.QR_stats
self.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
self.parameters['dt'] = (opt.dtfactor / opt.n)
if ((self.parameters['niter_todo'] % self.parameters['niter_out']) != 0):
self.parameters['niter_out'] = self.parameters['niter_todo']
if self.QR_stats_on:
# max_Q_estimate and max_R_estimate are just used for the 2D pdf
# therefore I just want them to be small multiples of mean trS2
# I'm already estimating the dissipation with kMeta...
meantrS2 = (opt.n//2 / opt.kMeta)**4 * self.parameters['nu']**2
self.parameters['max_Q_estimate'] = meantrS2
self.parameters['max_R_estimate'] = .4*meantrS2**1.5
def launch( self.pars_from_namespace(opt)
opt, self.fill_up_fluid_code()
nu = None): self.finalize_code()
c = NavierStokes(work_dir = opt.work_dir) self.write_src()
assert((opt.nsteps % 4) == 0) self.set_host_info(bfps.host_info)
c.parameters['nx'] = opt.n if not os.path.exists(os.path.join(self.work_dir, self.simname + '.h5')):
c.parameters['ny'] = opt.n self.write_par()
c.parameters['nz'] = opt.n if self.parameters['nparticles'] > 0:
if type(nu) == type(None): data = self.generate_tracer_state(
c.parameters['nu'] = 5.5*opt.n**(-4./3) species = 0,
else: rseed = opt.particle_rand_seed)
c.parameters['nu'] = nu for s in range(1, self.particle_species):
c.parameters['dt'] = 5e-3 * (64. / opt.n) self.generate_tracer_state(species = s, data = data)
c.parameters['niter_todo'] = opt.nsteps init_condition_file = os.path.join(
c.parameters['niter_part'] = 2 self.work_dir,
c.parameters['famplitude'] = 0.2 self.simname + '_cvorticity_i{0:0>5x}'.format(0))
c.parameters['nparticles'] = 32 if not os.path.exists(init_condition_file):
if opt.particles: if len(opt.src_simname) > 0:
c.add_particles() src_file = os.path.join(
c.add_particles(kcut = 'fs->kM/2') self.work_dir,
c.finalize_code() opt.src_simname + '_cvorticity_i{0:0>5x}'.format(opt.src_iteration))
c.write_src() os.symlink(src_file, init_condition_file)
c.write_par() else:
if opt.run: