Skip to content
Snippets Groups Projects
Commit ac6ef6ca authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

make code with different tracer species

parent 21a33811
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
bfps/NavierStokes.py
+ 96
− 92
View file @ ac6ef6ca
...@@ -27,10 +27,9 @@ import pickle ...@@ -27,10 +27,9 @@ import pickle
class NavierStokes(bfps.code): class NavierStokes(bfps.code):
def __init__( def __init__(
self, self,
name = 'NavierStokes', name = 'NavierStokes'):
particles_on = False):
super(NavierStokes, self).__init__() super(NavierStokes, self).__init__()
self.particles_on = particles_on self.particle_species = 0
self.name = name self.name = name
self.parameters['dkx'] = 1.0 self.parameters['dkx'] = 1.0
self.parameters['dky'] = 1.0 self.parameters['dky'] = 1.0
...@@ -40,35 +39,23 @@ class NavierStokes(bfps.code): ...@@ -40,35 +39,23 @@ class NavierStokes(bfps.code):
self.parameters['nu'] = 0.1 self.parameters['nu'] = 0.1
self.parameters['famplitude'] = 1.0 self.parameters['famplitude'] = 1.0
self.parameters['fmode'] = 1 self.parameters['fmode'] = 1
self.fluid_includes = '' self.parameters['nparticles'] = 0
self.fluid_includes = '#include "fluid_solver.hpp"\n'
self.fluid_variables = '' self.fluid_variables = ''
self.fluid_definitions = '' self.fluid_definitions = ''
self.fluid_start = '' self.fluid_start = ''
self.fluid_loop = '' self.fluid_loop = ''
self.fluid_end = '' self.fluid_end = ''
self.particle_includes = '' self.particle_includes = '#include "tracers.hpp"\n'
self.particle_variables = '' self.particle_variables = ''
self.particle_definitions = '' self.particle_definitions = ''
self.particle_start = '' self.particle_start = ''
self.particle_loop = '' self.particle_loop = ''
self.particle_end = '' self.particle_end = ''
self.stats_dtype = np.dtype([('iteration', np.int32),
('t', np.float64),
('energy', np.float64),
('enstrophy', np.float64)])
pickle.dump(
self.stats_dtype,
open(self.name + '_dtype.pickle', 'w'))
self.fill_up_fluid_code() self.fill_up_fluid_code()
return None return None
def fill_up_fluid_code(self): def write_fluid_stats(self):
return None self.fluid_definitions += """
self.includes += '#include <cstring>\n'
self.includes += '#include "fftw_tools.hpp"\n'
self.variables += ('double t;\n' +
'FILE *stat_file;\n'
'double stats[2];\n')
self.definitions += """
//begincpp //begincpp
void do_stats(fluid_solver<float> *fsolver) void do_stats(fluid_solver<float> *fsolver)
{ {
...@@ -81,30 +68,27 @@ class NavierStokes(bfps.code): ...@@ -81,30 +68,27 @@ class NavierStokes(bfps.code):
fwrite((void*)&t, sizeof(double), 1, stat_file); fwrite((void*)&t, sizeof(double), 1, stat_file);
fwrite((void*)stats, sizeof(double), 2, stat_file); fwrite((void*)stats, sizeof(double), 2, stat_file);
} }
fs->write_spectrum("velocity", fs->cvelocity);
} }
//endcpp //endcpp
""" """
if self.particles_on: self.stats_dtype = np.dtype([('iteration', np.int32),
self.parameters['nparticles'] = 1 ('t', np.float64),
self.includes += '#include "tracers.hpp"\n' ('energy', np.float64),
self.variables += 'FILE *traj_file;\n' ('enstrophy', np.float64)])
self.definitions += """ pickle.dump(
//begincpp self.stats_dtype,
void out_traj(tracers<float> *tracers) open(self.name + '_dtype.pickle', 'w'))
{ return None
if (myrank == 0) def fill_up_fluid_code(self):
{ self.fluid_includes += '#include <cstring>\n'
fwrite((void*)tracers->state, sizeof(double), tracers->array_size, traj_file); self.fluid_variables += ('double t;\n' +
} 'fluid_solver<float> *fs;\n' +
} 'FILE *stat_file;\n' +
//endcpp 'double stats[2];\n')
""" self.write_fluid_stats()
self.variables += self.cdef_pars() self.fluid_start += """
self.definitions += self.cread_pars()
self.main = """
//begincpp //begincpp
fluid_solver<float> *fs;
tracers<float> *ps;
char fname[512]; char fname[512];
fs = new fluid_solver<float>( fs = new fluid_solver<float>(
simname, simname,
...@@ -124,69 +108,86 @@ class NavierStokes(bfps.code): ...@@ -124,69 +108,86 @@ class NavierStokes(bfps.code):
stat_file = fopen(fname, "wb"); stat_file = fopen(fname, "wb");
} }
t = 0.0; t = 0.0;
do_stats(fs);
//endcpp //endcpp
""" """
if self.particles_on: self.fluid_loop += """
self.main += """
//begincpp
sprintf(fname, "%s_tracers", simname);
ps = new tracers<float>(
fname, fs,
nparticles, 1, 1,
fs->ru);
ps->dt = dt;
ps->iteration = iter0;
ps->read();
fs->compute_velocity(fs->cvorticity);
fftwf_execute(*((fftwf_plan*)fs->c2r_velocity));
ps->update_field();
if (myrank == 0)
{
sprintf(fname, "%s_traj.bin", ps->name);
traj_file = fopen(fname, "wb");
}
out_traj(ps);
//endcpp
"""
self.main += """
//begincpp //begincpp
do_stats(fs);
fs->write_spectrum("velocity", fs->cvelocity);
for (; fs->iteration < iter0 + niter_todo;)
{
fs->step(dt); fs->step(dt);
t += dt; t += dt;
do_stats(fs); do_stats(fs);
//endcpp //endcpp
""" """
if self.particles_on: self.fluid_end += """
self.main += """
//begincpp //begincpp
fs->compute_velocity(fs->cvorticity);
fftwf_execute(*((fftwf_plan*)fs->c2r_velocity));
ps->update_field();
ps->Euler();
ps->iteration++;
ps->synchronize();
//endcpp
"""
self.main += """
//begincpp
}
if (myrank == 0) if (myrank == 0)
{ {
fclose(stat_file); fclose(stat_file);
} }
fs->write('v', 'c'); fs->write('v', 'c');
fs->write_spectrum("velocity", fs->cvelocity); delete fs;
//endcpp //endcpp
""" """
if self.particles_on: return None
self.main += """ def add_particles(
ps->write(); self,
delete ps; neighbours = 1,
""" smoothness = 1,
self.main += 'delete fs;\n' kcut = 'fs->kM'):
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.particle_variables += 'tracers<float> *ps{0};'.format(self.particle_species)
self.particle_variables += 'FILE *traj_file{0};\n'.format(self.particle_species)
#self.particle_definitions
update_field = ('fs->compute_velocity(fs->cvorticity);\n' +
'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, "%s_tracers{0}", simname);\n' +
'ps{0} = new tracers<float>(\n' +
'fname, fs,\n' +
'nparticles, neighbours{0}, smoothness{0},\n' +
'fs->ru);\n' +
'ps{0}->dt = dt;\n' +
'ps{0}->iteration = iter0;\n' +
update_field +
'ps{0}->read();\n' +
'if (myrank == 0)\n' +
'{{\n' +
' sprintf(fname, "%s_traj{0}.bin", ps{0}->name);\n' +
' traj_file{0} = fopen(fname, "wb");\n' +
' fwrite((void*)ps{0}->state, sizeof(double), ps{0}->array_size, traj_file{0});\n' +
'}}\n').format(self.particle_species)
self.particle_loop += (update_field +
'ps{0}->Euler();\n' +
'ps{0}->iteration++;\n' +
'ps{0}->synchronize();\n').format(self.particle_species)
self.particle_end += ('ps{0}->write();\n' +
'delete ps{0};\n').format(self.particle_species)
self.particle_species += 1
return None
def finalize_code(self):
self.variables += self.cdef_pars()
self.definitions+= self.cread_pars()
self.includes += self.fluid_includes
self.variables += self.fluid_variables
self.definitions+= self.fluid_definitions
if self.particle_species > 0:
self.includes += self.particle_includes
self.variables += self.particle_variables
self.definitions += self.particle_definitions
self.main = self.fluid_start
if self.particle_species > 0:
self.main += self.particle_start
self.main += 'for (; fs->iteration < iter0 + niter_todo;)\n{\n'
self.main += self.fluid_loop
if self.particle_species > 0:
self.main += self.particle_loop
self.main += '\n}\n'
if self.particle_species > 0:
self.main += self.particle_end
self.main += self.fluid_end
return None return None
def plot_vel_cut( def plot_vel_cut(
self, self,
...@@ -255,10 +256,11 @@ class NavierStokes(bfps.code): ...@@ -255,10 +256,11 @@ class NavierStokes(bfps.code):
self, self,
rseed = 34982, rseed = 34982,
simname = None, simname = None,
iteration = 0): iteration = 0,
species = 0):
data = np.random.random(self.parameters['nparticles']*3)*2*np.pi data = np.random.random(self.parameters['nparticles']*3)*2*np.pi
if not (type(simname) == type(None)): if not (type(simname) == type(None)):
data.tofile(simname + "_tracers_state_i{0:0>5x}".format(iteration)) data.tofile(simname + "_tracers{0}_state_i{1:0>5x}".format(species, iteration))
return data return data
import subprocess import subprocess
...@@ -270,7 +272,7 @@ def test(opt): ...@@ -270,7 +272,7 @@ def test(opt):
if opt.clean: if opt.clean:
subprocess.call(['make', 'clean']) subprocess.call(['make', 'clean'])
return None return None
c = NavierStokes(particles_on = True) c = NavierStokes()
c.parameters['nx'] = opt.n c.parameters['nx'] = opt.n
c.parameters['ny'] = opt.n c.parameters['ny'] = opt.n
c.parameters['nz'] = opt.n c.parameters['nz'] = opt.n
...@@ -279,6 +281,8 @@ def test(opt): ...@@ -279,6 +281,8 @@ def test(opt):
c.parameters['niter_todo'] = opt.nsteps c.parameters['niter_todo'] = opt.nsteps
c.parameters['famplitude'] = 0.0 c.parameters['famplitude'] = 0.0
c.parameters['nparticles'] = 32 c.parameters['nparticles'] = 32
c.add_particles()
c.finalize_code()
c.write_src() c.write_src()
c.write_par(simname = 'test') c.write_par(simname = 'test')
c.generate_vector_field(simname = 'test') c.generate_vector_field(simname = 'test')
......
src/fluid_solver.cpp
+ 6
− 0
View file @ ac6ef6ca
...@@ -240,6 +240,12 @@ void fluid_solver<R>::compute_velocity(C *vorticity) \ ...@@ -240,6 +240,12 @@ void fluid_solver<R>::compute_velocity(C *vorticity) \
} \ } \
\ \
template<> \ template<> \
void fluid_solver<R>::ift_velocity() \
{ \
FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \
} \
\
template<> \
void fluid_solver<R>::add_forcing(\ void fluid_solver<R>::add_forcing(\
C *field, R factor) \ C *field, R factor) \
{ \ { \
......
src/fluid_solver.hpp
+ 1
− 1
View file @ ac6ef6ca
...@@ -22,7 +22,6 @@ ...@@ -22,7 +22,6 @@
#include <stdlib.h> #include <stdlib.h>
#include <iostream> #include <iostream>
#include "field_descriptor.hpp" #include "field_descriptor.hpp"
#include "vector_field.hpp"
#include "fluid_solver_base.hpp" #include "fluid_solver_base.hpp"
#ifndef FLUID_SOLVER #ifndef FLUID_SOLVER
...@@ -77,6 +76,7 @@ class fluid_solver:public fluid_solver_base<rnumber> ...@@ -77,6 +76,7 @@ class fluid_solver:public fluid_solver_base<rnumber>
void compute_vorticity(void); void compute_vorticity(void);
void compute_velocity(rnumber (*vorticity)[2]); void compute_velocity(rnumber (*vorticity)[2]);
void ift_velocity();
void omega_nonlin(int src); void omega_nonlin(int src);
void step(double dt); void step(double dt);
void impose_zero_modes(void); void impose_zero_modes(void);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment