Commit 364e016f authored by Chichi Lalescu's avatar Chichi Lalescu
Browse files

NSVEp class compiles

parent 7beb7a4f
......@@ -2,109 +2,6 @@
#include <cmath>
#include "NSVE.hpp"
int grow_single_dataset(hid_t dset, int tincrement)
{
int ndims;
hsize_t space;
space = H5Dget_space(dset);
ndims = H5Sget_simple_extent_ndims(space);
hsize_t *dims = new hsize_t[ndims];
H5Sget_simple_extent_dims(space, dims, NULL);
dims[0] += tincrement;
H5Dset_extent(dset, dims);
H5Sclose(space);
delete[] dims;
return EXIT_SUCCESS;
}
herr_t grow_dataset_visitor(
hid_t o_id,
const char *name,
const H5O_info_t *info,
void *op_data)
{
if (info->type == H5O_TYPE_DATASET)
{
hsize_t dset = H5Dopen(o_id, name, H5P_DEFAULT);
grow_single_dataset(dset, *((int*)(op_data)));
H5Dclose(dset);
}
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVE<rnumber>::read_iteration(void)
{
/* read iteration */
hid_t dset;
hid_t iteration_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
H5P_DEFAULT);
dset = H5Dopen(
iteration_file,
"iteration",
H5P_DEFAULT);
H5Dread(
dset,
H5T_NATIVE_INT,
H5S_ALL,
H5S_ALL,
H5P_DEFAULT,
&this->iteration);
H5Dclose(dset);
dset = H5Dopen(
iteration_file,
"checkpoint",
H5P_DEFAULT);
H5Dread(
dset,
H5T_NATIVE_INT,
H5S_ALL,
H5S_ALL,
H5P_DEFAULT,
&this->checkpoint);
H5Dclose(dset);
H5Fclose(iteration_file);
DEBUG_MSG("simname is %s, iteration is %d and checkpoint is %d\n",
this->simname.c_str(),
this->iteration,
this->checkpoint);
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVE<rnumber>::write_iteration(void)
{
if (this->myrank == 0)
{
hid_t dset = H5Dopen(
this->stat_file,
"iteration",
H5P_DEFAULT);
H5Dwrite(
dset,
H5T_NATIVE_INT,
H5S_ALL,
H5S_ALL,
H5P_DEFAULT,
&this->iteration);
H5Dclose(dset);
dset = H5Dopen(
this->stat_file,
"checkpoint",
H5P_DEFAULT);
H5Dwrite(
dset,
H5T_NATIVE_INT,
H5S_ALL,
H5S_ALL,
H5P_DEFAULT,
&this->checkpoint);
H5Dclose(dset);
}
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVE<rnumber>::initialize(void)
......@@ -246,24 +143,6 @@ int NSVE<rnumber>::finalize(void)
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVE<rnumber>::grow_file_datasets()
{
int file_problems = 0;
hid_t group;
group = H5Gopen(this->stat_file, "/statistics", H5P_DEFAULT);
int tincrement = this->niter_todo / this->niter_stat;
H5Ovisit(
group,
H5_INDEX_NAME,
H5_ITER_NATIVE,
grow_dataset_visitor,
&tincrement);
H5Gclose(group);
return file_problems;
}
template <typename rnumber>
int NSVE<rnumber>::do_stats()
{
......
......@@ -34,25 +34,12 @@
#include "vorticity_equation.hpp"
#include "full_code/direct_numerical_simulation.hpp"
int grow_single_dataset(hid_t dset, int tincrement);
herr_t grow_dataset_visitor(
hid_t o_id,
const char *name,
const H5O_info_t *info,
void *op_data);
template <typename rnumber>
class NSVE: public direct_numerical_simulation
{
public:
/* parameters that are read in read_parameters */
int checkpoints_per_file;
int dealias_type;
double dkx;
double dky;
double dkz;
double dt;
double famplitude;
double fk0;
......@@ -62,20 +49,9 @@ class NSVE: public direct_numerical_simulation
int histogram_bins;
double max_velocity_estimate;
double max_vorticity_estimate;
int niter_out;
int niter_part;
int niter_stat;
int niter_todo;
int nparticles;
double nu;
int nx;
int ny;
int nz;
/* other stuff */
int iteration, checkpoint;
hid_t stat_file;
bool stop_code_now;
vorticity_equation<rnumber, FFTW> *fs;
field<rnumber, FFTW, THREE> *tmp_vec_field;
field<rnumber, FFTW, ONE> *tmp_scal_field;
......@@ -93,10 +69,7 @@ class NSVE: public direct_numerical_simulation
int main_loop(void);
int finalize(void);
int read_iteration(void);
int write_iteration(void);
int read_parameters(void);
int grow_file_datasets(void);
int do_stats(void);
};
......
#include <string>
#include <cmath>
#include "NSVEp.hpp"
template <typename rnumber>
int NSVEp<rnumber>::initialize(void)
{
this->read_iteration();
this->read_parameters();
if (this->myrank == 0)
{
// set caching parameters
hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
DEBUG_MSG("when setting stat_file cache I got %d\n", cache_err);
this->stat_file = H5Fopen(
(this->simname + ".h5").c_str(),
H5F_ACC_RDWR,
fapl);
}
int data_file_problem;
if (this->myrank == 0)
data_file_problem = this->grow_file_datasets();
MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->comm);
if (data_file_problem > 0)
{
std::cerr <<
data_file_problem <<
" problems growing file datasets.\ntrying to exit now." <<
std::endl;
return EXIT_FAILURE;
}
this->fs = new vorticity_equation<rnumber, FFTW>(
simname.c_str(),
nx, ny, nz,
dkx, dky, dkz,
FFTW_MEASURE);
this->tmp_vec_field = new field<rnumber, FFTW, THREE>(
nx, ny, nz,
this->comm,
FFTW_MEASURE);
this->fs->checkpoints_per_file = checkpoints_per_file;
this->fs->nu = nu;
this->fs->fmode = fmode;
this->fs->famplitude = famplitude;
this->fs->fk0 = fk0;
this->fs->fk1 = fk1;
strncpy(this->fs->forcing_type, forcing_type, 128);
this->fs->iteration = this->iteration;
this->fs->checkpoint = this->checkpoint;
this->fs->cvorticity->real_space_representation = false;
this->fs->io_checkpoint();
if (this->myrank == 0 && this->iteration == 0)
this->fs->kk->store(stat_file);
this->ps = particles_system_builder(
fs->cvelocity, // (field object)
fs->kk, // (kspace object, contains dkx, dky, dkz)
tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
nparticles, // to check coherency between parameters and hdf input file
fs->get_current_fname(), // particles input filename
std::string("/tracers0/state/") + std::to_string(fs->iteration), // dataset name for initial input
std::string("/tracers0/rhs/") + std::to_string(fs->iteration), // dataset name for initial input
tracers0_neighbours, // parameter (interpolation no neighbours)
tracers0_smoothness, // parameter
this->comm,
fs->iteration+1);
this->particles_output_writer_mpi = new particles_output_hdf5<double,3,3>(
MPI_COMM_WORLD,
"tracers0",
nparticles,
tracers0_integration_steps);
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVEp<rnumber>::main_loop(void)
{
clock_t time0, time1;
double time_difference, local_time_difference;
time0 = clock();
bool stop_code_now = false;
int max_iter = (this->iteration + this->niter_todo -
(this->iteration % this->niter_todo));
for (; this->iteration < max_iter;)
{
#ifdef USE_TIMINGOUTPUT
const std::string loopLabel = ("code::main_start::loop-" +
std::to_string(this->iteration));
TIMEZONE(loopLabel.c_str());
#endif
if (this->iteration % this->niter_stat == 0)
this->do_stats();
this->fs->compute_velocity(fs->cvorticity);
this->fs->cvelocity->ift();
this->ps->completeLoop(dt);
this->fs->step(dt);
this->iteration = this->fs->iteration;
if (this->fs->iteration % this->niter_out == 0)
{
this->fs->io_checkpoint(false);
this->particles_output_writer_mpi->open_file(fs->get_current_fname());
this->particles_output_writer_mpi->save(
this->ps->getParticlesPositions(),
this->ps->getParticlesRhs(),
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->fs->iteration);
this->particles_output_writer_mpi->close_file();
this->checkpoint = this->fs->checkpoint;
this->write_iteration();
}
if (stop_code_now)
break;
time1 = clock();
local_time_difference = ((
(unsigned int)(time1 - time0)) /
((double)CLOCKS_PER_SEC));
time_difference = 0.0;
MPI_Allreduce(
&local_time_difference,
&time_difference,
1,
MPI_DOUBLE,
MPI_SUM,
MPI_COMM_WORLD);
if (this->myrank == 0)
std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
if (this->myrank == 0)
std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
time0 = time1;
}
if (this->iteration % this->niter_stat == 0)
this->do_stats();
time1 = clock();
local_time_difference = ((
(unsigned int)(time1 - time0)) /
((double)CLOCKS_PER_SEC));
time_difference = 0.0;
MPI_Allreduce(
&local_time_difference,
&time_difference,
1,
MPI_DOUBLE,
MPI_SUM,
MPI_COMM_WORLD);
if (this->myrank == 0)
std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
if (this->myrank == 0)
std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
if (this->iteration % this->niter_out != 0)
{
this->fs->io_checkpoint(false);
this->particles_output_writer_mpi->open_file(fs->get_current_fname());
this->particles_output_writer_mpi->save(
this->ps->getParticlesPositions(),
this->ps->getParticlesRhs(),
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->fs->iteration);
this->particles_output_writer_mpi->close_file();
this->checkpoint = this->fs->checkpoint;
this->write_iteration();
}
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVEp<rnumber>::finalize(void)
{
if (this->myrank == 0)
H5Fclose(this->stat_file);
delete this->fs;
delete this->tmp_vec_field;
this->ps.release();
delete this->particles_output_writer_mpi;
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVEp<rnumber>::do_stats()
{
hid_t stat_group;
if (this->myrank == 0)
stat_group = H5Gopen(
this->stat_file,
"statistics",
H5P_DEFAULT);
else
stat_group = 0;
fs->compute_velocity(fs->cvorticity);
*tmp_vec_field = fs->cvelocity->get_cdata();
tmp_vec_field->compute_stats(
fs->kk,
stat_group,
"velocity",
fs->iteration / niter_stat,
max_velocity_estimate/sqrt(3));
*tmp_vec_field = fs->cvorticity->get_cdata();
tmp_vec_field->compute_stats(
fs->kk,
stat_group,
"vorticity",
fs->iteration / niter_stat,
max_vorticity_estimate/sqrt(3));
if (this->myrank == 0)
H5Gclose(stat_group);
if (myrank == 0)
{
std::string fname = (
std::string("stop_") +
std::string(simname));
{
struct stat file_buffer;
this->stop_code_now = (
stat(fname.c_str(), &file_buffer) == 0);
}
}
MPI_Bcast(
&this->stop_code_now,
1,
MPI_C_BOOL,
0,
MPI_COMM_WORLD);
return EXIT_SUCCESS;
}
template class NSVEp<float>;
template class NSVEp<double>;
/**********************************************************************
* *
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
* This file is part of bfps. *
* *
* bfps is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
* bfps is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with bfps. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
**********************************************************************/
#ifndef NSVEP_HPP
#define NSVEP_HPP
#include <cstdlib>
#include "base.hpp"
#include "vorticity_equation.hpp"
#include "full_code/direct_numerical_simulation.hpp"
#include "particles/particles_system_builder.hpp"
#include "particles/particles_output_hdf5.hpp"
template <typename rnumber>
class NSVEp: public direct_numerical_simulation
{
public:
/* parameters that are read in read_parameters */
double dt;
double famplitude;
double fk0;
double fk1;
double nu;
int fmode;
char forcing_type[512];
int histogram_bins;
double max_velocity_estimate;
double max_vorticity_estimate;
int niter_part;
int nparticles;
int tracers0_integration_steps;
int tracers0_neighbours;
int tracers0_smoothness;
/* other stuff */
vorticity_equation<rnumber, FFTW> *fs;
field<rnumber, FFTW, THREE> *tmp_vec_field;
field<rnumber, FFTW, ONE> *tmp_scal_field;
std::unique_ptr<abstract_particles_system<double>> ps;
particles_output_hdf5<double,3,3> *particles_output_writer_mpi;
NSVEp(
const MPI_Comm COMMUNICATOR,
const std::string &simulation_name):
direct_numerical_simulation(
COMMUNICATOR,
simulation_name){}
~NSVEp(){}
int initialize(void);
int main_loop(void);
int finalize(void);
int read_parameters(void);
int do_stats(void);
};
#endif//NSVEP_HPP
#include "direct_numerical_simulation.hpp"
int grow_single_dataset(hid_t dset, int tincrement)
{
int ndims;
hsize_t space;
space = H5Dget_space(dset);
ndims = H5Sget_simple_extent_ndims(space);
hsize_t *dims = new hsize_t[ndims];
H5Sget_simple_extent_dims(space, dims, NULL);
dims[0] += tincrement;
H5Dset_extent(dset, dims);
H5Sclose(space);
delete[] dims;
return EXIT_SUCCESS;
}
herr_t grow_dataset_visitor(
hid_t o_id,
const char *name,
const H5O_info_t *info,
void *op_data)
{
if (info->type == H5O_TYPE_DATASET)
{
hsize_t dset = H5Dopen(o_id, name, H5P_DEFAULT);
grow_single_dataset(dset, *((int*)(op_data)));
H5Dclose(dset);
}
return EXIT_SUCCESS;
}
direct_numerical_simulation::direct_numerical_simulation(
const MPI_Comm COMMUNICATOR,
const std::string &simulation_name):
......@@ -10,3 +40,93 @@ direct_numerical_simulation::direct_numerical_simulation(
MPI_Comm_size(this->comm, &this->nprocs);
}
int direct_numerical_simulation::grow_file_datasets()
{
int file_problems = 0;
hid_t group;
group = H5Gopen(this->stat_file, "/statistics", H5P_DEFAULT);
int tincrement = this->niter_todo / this->niter_stat;
H5Ovisit(
group,
H5_INDEX_NAME,
H5_ITER_NATIVE,
grow_dataset_visitor,
&tincrement);
H5Gclose(group);
return file_problems;
}
int direct_numerical_simulation::read_iteration(void)
{
/* read iteration */
hid_t dset;
hid_t iteration_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
H5P_DEFAULT);
dset = H5Dopen(
iteration_file,
"iteration",
H5P_DEFAULT);
H5Dread(
dset,
H5T_NATIVE_INT,
H5S_ALL,
H5S_ALL,
H5P_DEFAULT,
&this->iteration);
H5Dclose(dset);
dset = H5Dopen(
iteration_file,
"checkpoint",
H5P_DEFAULT);
H5Dread(
dset,
H5T_NATIVE_INT,
H5S_ALL,
H5S_ALL,
H5P_DEFAULT,
&this->checkpoint);
H5Dclose(dset);
H5Fclose(iteration_file);
DEBUG_MSG("simname is %s, iteration is %d and checkpoint is %d\n",
this->simname.c_str(),
this->iteration,
this->checkpoint);
return EXIT_SUCCESS;