-
Cristian Lalescu authoredCristian Lalescu authored
test_interpolation.cpp 9.38 KiB
/******************************************************************************
* *
* Copyright 2019 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 *
* *
******************************************************************************/
#include "full_code/test_interpolation.hpp"
template <typename rnumber>
int test_interpolation<rnumber>::read_parameters(void)
{
TIMEZONE("test_interpolation::read_parameters");
this->test::read_parameters();
hid_t parameter_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
H5P_DEFAULT);
this->nparticles = hdf5_tools::read_value<int>(
parameter_file, "/parameters/nparticles");
this->tracers0_integration_steps = hdf5_tools::read_value<int>(
parameter_file, "/parameters/tracers0_integration_steps");
this->tracers0_neighbours = hdf5_tools::read_value<int>(
parameter_file, "/parameters/tracers0_neighbours");
this->tracers0_smoothness = hdf5_tools::read_value<int>(
parameter_file, "/parameters/tracers0_smoothness");
H5Fclose(parameter_file);
return EXIT_SUCCESS;
}
template <typename rnumber>
int test_interpolation<rnumber>::initialize(void)
{
TIMEZONE("test_interpolation::initialize");
this->read_parameters();
this->vorticity = new field<rnumber, FFTW, THREE>(
this->nx, this->ny, this->nz,
this->comm,
FFTW_ESTIMATE);
this->vorticity->real_space_representation = false;
this->velocity = new field<rnumber, FFTW, THREE>(
this->nx, this->ny, this->nz,
this->comm,
FFTW_ESTIMATE);
this->nabla_u = new field<rnumber, FFTW, THREExTHREE>(
this->nx, this->ny, this->nz,
this->comm,
FFTW_ESTIMATE);
this->kk = new kspace<FFTW, SMOOTH>(
this->vorticity->clayout, this->dkx, this->dky, this->dkz);
if (this->myrank == 0)
{
hid_t stat_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDWR,
H5P_DEFAULT);
this->kk->store(stat_file);
H5Fclose(stat_file);
}
this->ps = particles_system_builder(
this->velocity, // (field object)
this->kk, // (kspace object, contains dkx, dky, dkz)
this->tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
(long long int)nparticles, // to check coherency between parameters and hdf input file
this->simname + "_input.h5", // particles input filename
std::string("/tracers0/state/0"), // dataset name for initial input
std::string("/tracers0/rhs/0") , // dataset name for initial input
this->tracers0_neighbours, // parameter (interpolation no neighbours)
this->tracers0_smoothness, // parameter
this->comm,
1);
this->particles_output_writer_mpi = new particles_output_hdf5<
long long int, double, 3>(
MPI_COMM_WORLD,
"tracers0",
nparticles,
this->tracers0_integration_steps);
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
long long int, double, 3>(
MPI_COMM_WORLD,
this->ps->getGlobalNbParticles(),
(this->simname + "_output.h5"),
"tracers0",
"position/0");
return EXIT_SUCCESS;
}
template <typename rnumber>
int test_interpolation<rnumber>::finalize(void)
{
TIMEZONE("test_interpolation::finalize");
delete this->nabla_u;
delete this->velocity;
delete this->vorticity;
delete this->ps.release();
delete this->kk;
delete particles_output_writer_mpi;
delete particles_sample_writer_mpi;
return EXIT_SUCCESS;
}
template <typename rnumber>
int test_interpolation<rnumber>::do_work()
{
TIMEZONE("test_interpolation::do_work");
*this->nabla_u = 0.0;
this->velocity->real_space_representation = false;
this->vorticity->real_space_representation = false;
this->nabla_u->real_space_representation = false;
// read vorticity field
this->vorticity->io(
this->simname + std::string("_input.h5"),
"vorticity",
0, true);
this->kk->template force_divfree<rnumber>(this->vorticity->get_cdata());
// compute velocity
invert_curl(this->kk, this->vorticity, this->velocity);
// compute velocity gradient
compute_gradient(this->kk, this->velocity, this->nabla_u);
// go to real space
this->vorticity->ift();
this->velocity->ift();
this->nabla_u->ift();
DEBUG_MSG("some vorticity values: %g %g %g\n",
this->vorticity->rval(20, 1),
this->vorticity->rval(200, 2),
this->vorticity->rval(741, 0));
DEBUG_MSG("corresponding velocity gradient to vorticity values: %g %g %g\n",
this->nabla_u->rval( 20, 2, 0) - this->nabla_u->rval( 20, 0, 2),
this->nabla_u->rval(200, 1, 0) - this->nabla_u->rval(200, 0, 1),
this->nabla_u->rval(741, 1, 2) - this->nabla_u->rval(741, 2, 1));
// allocate interpolation arrays
std::unique_ptr<double[]> p3data;
std::unique_ptr<double[]> p9data;
if(this->ps->getLocalNbParticles()){
p3data.reset(new double[3*this->ps->getLocalNbParticles()]);
p9data.reset(new double[9*this->ps->getLocalNbParticles()]);
}
/// sample position
std::copy(this->ps->getParticlesState(),
this->ps->getParticlesState()+3*this->ps->getLocalNbParticles(),
p3data.get());
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"position",
this->ps->getParticlesState(),
&p3data,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
/// sample velocity at particles' position
std::fill_n(p3data.get(), 3*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->velocity, p3data.get());
if(p3data){
DEBUG_MSG("first vel value is %g\n", p3data.get()[0]);
}
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"velocity",
this->ps->getParticlesState(),
&p3data,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
/// sample vorticity at particles' position
std::fill_n(p3data.get(), 3*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->vorticity, p3data.get());
if(p3data){
DEBUG_MSG("first vort value is %g\n", p3data.get()[0]);
}
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"vorticity",
this->ps->getParticlesState(),
&p3data,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1); /// sample velocity gradient at particles' position
std::fill_n(p9data.get(), 9*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->nabla_u, p9data.get());
if(p9data){
DEBUG_MSG("first vel gradient value is %g\n", p9data.get()[0]);
}
this->particles_sample_writer_mpi->template save_dataset<9>(
"tracers0",
"velocity_gradient",
this->ps->getParticlesState(),
&p9data,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
// deallocate temporary arrays
delete[] p3data.release();
delete[] p9data.release();
return EXIT_SUCCESS;
}
template class test_interpolation<float>;
template class test_interpolation<double>;