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

[broken] use new sampling class for P2P

Sampled position is different from checkpoint position, I still need to
fix that.
However, there's another issue that needs to be solved first
parent ef790303
Pipeline #20549 passed with stage
in 10 minutes and 16 seconds
......@@ -16,6 +16,7 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
particles_inner_computer<double, long long int> current_particles_inner_computer(inner_v0);
current_particles_inner_computer.setEnable(enable_inner);
this->cutoff = 5.0;
this->ps = particles_system_builder_with_p2p(
this->fs->cvelocity, // (field object)
......@@ -39,6 +40,13 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
"tracers0",
nparticles,
tracers0_integration_steps);
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
long long int, double, 3, 3>(
MPI_COMM_WORLD,
this->ps->getGlobalNbParticles(),
(this->simname + "_particles.h5"),
"tracers0",
"position/0");
return EXIT_SUCCESS;
}
......@@ -76,8 +84,8 @@ int NSVEparticlesP2P<rnumber>::write_checkpoint(void)
template <typename rnumber>
int NSVEparticlesP2P<rnumber>::finalize(void)
{
this->ps.release();
delete this->particles_output_writer_mpi;
delete this->particles_sample_writer_mpi;
this->NSVE<rnumber>::finalize();
return EXIT_SUCCESS;
}
......@@ -95,30 +103,45 @@ int NSVEparticlesP2P<rnumber>::do_stats()
if (!(this->iteration % this->niter_part == 0))
return EXIT_SUCCESS;
// allocate temporary data array
std::unique_ptr<double[]> pdata(new double[3*this->ps->getLocalNbParticles()]);
/// sample position
sample_particles_system_position(
this->ps,
(this->simname + "_particles.h5"), // filename
"tracers0", // hdf5 parent group
"position" // dataset basename TODO
);
/// sample velocity
sample_from_particles_system(*this->tmp_vec_field, // field to save
this->ps,
(this->simname + "_particles.h5"), // filename
"tracers0", // hdf5 parent group
"velocity" // dataset basename TODO
);
/// compute acceleration and sample it
this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field);
this->tmp_vec_field->ift();
sample_from_particles_system(*this->tmp_vec_field,
this->ps,
(this->simname + "_particles.h5"),
"tracers0",
"acceleration");
std::copy(this->ps->getParticlesPositions(),
this->ps->getParticlesPositions()+3*this->ps->getLocalNbParticles(),
pdata.get());
this->particles_sample_writer_mpi->save_dataset(
"tracers0",
"position",
this->ps->getParticlesPositions(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
///// sample velocity
//this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
//this->particles_sample_writer_mpi->save_dataset(
// "tracers0",
// "velocity",
// this->ps->getParticlesPositions(),
// &pdata,
// this->ps->getParticlesIndexes(),
// this->ps->getLocalNbParticles(),
// this->ps->get_step_idx()-1);
///// compute acceleration and sample it
//this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field);
//this->tmp_vec_field->ift();
//this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
//this->particles_sample_writer_mpi->save_dataset(
// "tracers0",
// "acceleration",
// this->ps->getParticlesPositions(),
// &pdata,
// this->ps->getParticlesIndexes(),
// this->ps->getLocalNbParticles(),
// this->ps->get_step_idx()-1);
return EXIT_SUCCESS;
}
......
......@@ -35,6 +35,7 @@
#include "full_code/NSVE.hpp"
#include "particles/particles_system_builder.hpp"
#include "particles/particles_output_hdf5.hpp"
#include "particles/particles_sampling.hpp"
/** \brief Navier-Stokes solver that includes simple Lagrangian tracers.
*
......@@ -65,6 +66,7 @@ class NSVEparticlesP2P: public NSVE<rnumber>
std::unique_ptr<abstract_particles_system<long long int, double>> ps;
// TODO P2P use a reader with particle data
particles_output_hdf5<long long int, double,6,6> *particles_output_writer_mpi;
particles_output_sampling_hdf5<long long int, double,3,3> *particles_sample_writer_mpi;
NSVEparticlesP2P(
......
......@@ -12,45 +12,58 @@ import matplotlib.pyplot as plt
def main():
assert(sys.argv[1] in ['p2p_sampling']);
assert(sys.argv[1] in ['p2p_sampling'])
assert(sys.argv[2] in ['on', 'off'])
niterations = 32
nparticles = 100
nparticles = 1000
njobs = 1
c = DNS()
c.launch(
['NSVEparticlesP2P',
'-n', '32',
'--src-simname', 'B32p1e4',
'--src-wd', bfps.lib_dir + '/test',
'--src-iteration', '0',
'--np', '4',
'--ntpp', '1',
'--niter_todo', '{0}'.format(niterations),
'--niter_out', '{0}'.format(niterations),
'--niter_stat', '1',
'--checkpoints_per_file', '{0}'.format(3),
'--nparticles', '{0}'.format(nparticles),
'--particle-rand-seed', '2',
'--njobs', '{0}'.format(njobs),
'--wd', './'] +
sys.argv[2:])
if sys.argv[2] == 'on':
c = DNS()
c.launch(
['NSVEparticlesP2P',
'-n', '32',
'--src-simname', 'B32p1e4',
'--src-wd', bfps.lib_dir + '/test',
'--src-iteration', '0',
'--np', '4',
'--ntpp', '1',
'--niter_todo', '{0}'.format(niterations),
'--niter_out', '{0}'.format(niterations),
'--niter_stat', '1',
'--checkpoints_per_file', '{0}'.format(3),
'--nparticles', '{0}'.format(nparticles),
'--particle-rand-seed', '2',
'--njobs', '{0}'.format(njobs),
'--wd', './'] +
sys.argv[3:])
if sys.argv[1] == 'p2p_sampling':
cf = h5py.File(
'test_checkpoint_0.h5',
'r')
pf = h5py.File(
os.path.join(
c.work_dir,
c.simname + '_particles.h5'),
'test_particles.h5',
'r')
# show a histogram of the positions
f = plt.figure()
a = f.add_subplot(111)
for iteration in [0, 16, 32, 48, 64]:
for iteration in range(0, niterations*njobs+1, niterations//2):
x = pf['tracers0/position/{0}'.format(iteration)].value
print(x.shape)
#bins, hist = np.histogram(
# a.
hist, bins = np.histogram(
np.sum(x**2, axis = -1).flatten()**.5,
bins = 40)
bb = (bins[:-1] + bins[1:])/2
pp = hist.astype(np.float) / (np.sum(hist) * (bb[1] - bb[0]))
a.plot(bb, pp, label = '{0}'.format(iteration))
a.legend(loc = 'best')
f.tight_layout()
f.savefig('position_histogram.pdf')
plt.close(f)
# compared sampled positions with checkpoint positions
for iteration in range(0, niterations*njobs+1, niterations):
x = pf['tracers0/position/{0}'.format(iteration)].value
s = cf['tracers0/state/{0}'.format(iteration)].value
distance = np.max(np.abs(x - s[..., :3]))
print(iteration, distance)
# show a histogram of the orientations
return None
......
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