Skip to content
GitLab
Commit 07a0d81e authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/filtering' into develop

parents 0523475f 4ba46ebe
Pipeline #15740 passed with stage
in 6 minutes and 13 seconds
Hide whitespace changes
Inline Side-by-side
bfps/TEST.py
View file @ 07a0d81e
......@@ -123,10 +123,6 @@ class TEST(_code):
self.parameters['dky'] = float(1.0)
self.parameters['dkz'] = float(1.0)
self.parameters['filter_length'] = float(1.0)
self.parameters['niter_todo'] = int(1)
self.parameters['niter_stat'] = int(1)
self.parameters['niter_out'] = int(1)
self.parameters['checkpoints_per_file'] = int(1)
return None
def get_kspace(self):
kspace = {}
......@@ -162,9 +158,6 @@ class TEST(_code):
self,
iter0 = 0,
particle_ic = None):
assert (self.parameters['niter_todo'] == 1)
assert (self.parameters['niter_out'] == 1)
assert (self.parameters['niter_stat'] == 1)
assert (iter0 == 0)
_code.write_par(self, iter0 = iter0)
with h5py.File(self.get_data_file_name(), 'r+') as ofile:
......
bfps/cpp/full_code/filter_test.cpp
View file @ 07a0d81e
......@@ -24,7 +24,31 @@ int filter_test<rnumber>::initialize(void)
this->kk->store(stat_file);
H5Fclose(stat_file);
}
this->read_iteration();
return EXIT_SUCCESS;
}
template <typename rnumber>
int filter_test<rnumber>::finalize(void)
{
delete this->scal_field;
delete this->kk;
return EXIT_SUCCESS;
}
template <typename rnumber>
int filter_test<rnumber>::read_parameters()
{
this->test::read_parameters();
hid_t parameter_file;
hid_t dset, memtype, space;
parameter_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
H5P_DEFAULT);
dset = H5Dopen(parameter_file, "/parameters/filter_length", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->filter_length);
H5Dclose(dset);
H5Fclose(parameter_file);
return EXIT_SUCCESS;
}
......@@ -69,15 +93,7 @@ int filter_test<rnumber>::reset_field(
}
template <typename rnumber>
int filter_test<rnumber>::step(void)
{
this->iteration = (this->iteration + this->niter_todo -
(this->iteration % this->niter_todo));
return EXIT_SUCCESS;
}
template <typename rnumber>
int filter_test<rnumber>::write_checkpoint(void)
int filter_test<rnumber>::do_work(void)
{
std::string filename = this->simname + std::string("_fields.h5");
for (int dimension = 0; dimension < 3; dimension++)
......@@ -85,7 +101,7 @@ int filter_test<rnumber>::write_checkpoint(void)
this->reset_field(dimension);
this->kk->template filter<rnumber, ONE>(
this->scal_field->get_cdata(),
4*acos(0) / (this->filter_length),
4*acos(0) / this->filter_length,
"sharp_Fourier_sphere");
this->scal_field->ift();
this->scal_field->normalize();
......@@ -97,7 +113,7 @@ int filter_test<rnumber>::write_checkpoint(void)
this->reset_field(dimension);
this->kk->template filter<rnumber, ONE>(
this->scal_field->get_cdata(),
4*acos(0) / (this->filter_length / sqrt(log(2))),
4*acos(0) / this->filter_length,
"Gauss");
this->scal_field->ift();
this->scal_field->normalize();
......@@ -109,7 +125,7 @@ int filter_test<rnumber>::write_checkpoint(void)
this->reset_field(dimension);
this->kk->template filter<rnumber, ONE>(
this->scal_field->get_cdata(),
4*acos(0) / (this->filter_length*2),
4*acos(0) / this->filter_length,
"ball");
this->scal_field->ift();
this->scal_field->normalize();
......@@ -122,60 +138,6 @@ int filter_test<rnumber>::write_checkpoint(void)
return EXIT_SUCCESS;
}
template <typename rnumber>
int filter_test<rnumber>::finalize(void)
{
delete this->scal_field;
delete this->kk;
return EXIT_SUCCESS;
}
template <typename rnumber>
int filter_test<rnumber>::read_parameters()
{
hid_t parameter_file;
hid_t dset, memtype, space;
char fname[256];
char *string_data;
sprintf(fname, "%s.h5", this->simname.c_str());
parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
dset = H5Dopen(parameter_file, "/parameters/niter_todo", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->niter_todo);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/niter_stat", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->niter_stat);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/niter_out", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->niter_out);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/dealias_type", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dealias_type);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/dkx", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dkx);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/dky", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dky);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/dkz", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dkz);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/nx", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->nx);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/ny", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->ny);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/nz", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->nz);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/filter_length", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->filter_length);
H5Dclose(dset);
H5Fclose(parameter_file);
return EXIT_SUCCESS;
}
template class filter_test<float>;
template class filter_test<double>;
bfps/cpp/full_code/filter_test.hpp
View file @ 07a0d81e
......@@ -31,11 +31,26 @@
#include <cstdlib>
#include "base.hpp"
#include "vorticity_equation.hpp"
#include "full_code/direct_numerical_simulation.hpp"
#include "kspace.hpp"
#include "field.hpp"
#include "full_code/test.hpp"
/** \brief A class for testing filters.
*
* This class applies available filters to three Dirac distributions:
* - nonzero at the origin
* - nonzero on the `x` axis
* - nonzero on the `(x, y)` plane.
* All three distributions are normalized, so simple sanity checks can
* be performed afterwards in a Python script.
*
* While the convolutions can obviously be implemented in Python directly,
* it's better if the functionality is available here directly for easy
* reference.
*/
template <typename rnumber>
class filter_test: public direct_numerical_simulation
class filter_test: public test
{
public:
......@@ -49,20 +64,17 @@ class filter_test: public direct_numerical_simulation
filter_test(
const MPI_Comm COMMUNICATOR,
const std::string &simulation_name):
direct_numerical_simulation(
test(
COMMUNICATOR,
simulation_name){}
~filter_test(){}
int initialize(void);
int step(void);
int do_work(void);
int finalize(void);
int read_parameters(void);
int reset_field(int dimension);
int read_parameters(void);
int write_checkpoint(void);
int do_stats(void){}
};
#endif//FILTER_TEST_HPP
......
bfps/cpp/full_code/test.cpp 0 → 100644
View file @ 07a0d81e
#include <cstdlib>
#include <sys/types.h>
#include <sys/stat.h>
#include "scope_timer.hpp"
#include "hdf5_tools.hpp"
#include "full_code/test.hpp"
int test::main_loop(void)
{
#ifdef USE_TIMINGOUTPUT
TIMEZONE("test::main_loop");
#endif
this->start_simple_timer();
this->do_work();
this->print_simple_timer(
"do_work required " + std::to_string(this->iteration));
return EXIT_SUCCESS;
}
int test::read_parameters()
{
hid_t parameter_file;
hid_t dset, memtype, space;
char fname[256];
char *string_data;
sprintf(fname, "%s.h5", this->simname.c_str());
parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
dset = H5Dopen(parameter_file, "/parameters/dealias_type", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dealias_type);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/dkx", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dkx);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/dky", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dky);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/dkz", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->dkz);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/nx", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->nx);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/ny", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->ny);
H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/nz", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->nz);
H5Dclose(dset);
H5Fclose(parameter_file);
return 0;
}
bfps/cpp/full_code/test.hpp 0 → 100644
View file @ 07a0d81e
/**********************************************************************
* *
* 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 TEST_HPP
#define TEST_HPP
#include <cstdlib>
#include <sys/types.h>
#include <sys/stat.h>
#include <vector>
#include "base.hpp"
#include "full_code/code_base.hpp"
/** \brief base class for miscellaneous tests.
*
* Children of this class can basically do more or less anything inside their
* `do_work` method, which will be executed only once.
*/
class test: public code_base
{
public:
test(
const MPI_Comm COMMUNICATOR,
const std::string &simulation_name):
code_base(
COMMUNICATOR,
simulation_name){}
virtual ~test(){}
virtual int initialize(void) = 0;
virtual int do_work(void) = 0;
virtual int finalize(void) = 0;
int main_loop(void);
virtual int read_parameters(void);
};
#endif//TEST_HPP
bfps/cpp/kspace.cpp
View file @ 07a0d81e
......@@ -254,9 +254,16 @@ void kspace<be, dt>::low_pass(
/** \brief Filter a field using a ball shaped top hat filter.
*
* Filter's mathematical expression in Fourier space is as follows:
* Filter's mathematical expression in real space is as follows:
* \f[
* \phi^b_\ell(r) =
* \frac{1}{\ell^3}\frac{6}{\pi} H(\ell/2 - r)
* \f]
* with the corresponding Fourier space expression:
* \f[
* \hat{b}_\ell(\mathbf{k}) = \sin(k \ell / 2) / (k \ell / 2)
* \hat{\phi^b_\ell}(k) =
* \frac{3}{2(k\ell/2)^3}
* \left(2\sin (k \ell/2) - k \ell \cos (k \ell/2)\right)
* \f]
*/
template <field_backend be,
......@@ -265,9 +272,9 @@ template <typename rnumber,
field_components fc>
void kspace<be, dt>::ball_filter(
typename fftw_interface<rnumber>::complex *__restrict__ a,
const double sigma)
const double ell)
{
const double prefactor = sigma/2;
const double prefactor0 = double(3) / pow(ell/2, 3);
this->CLOOP_K2(
[&](ptrdiff_t cindex,
ptrdiff_t xindex,
......@@ -276,9 +283,12 @@ void kspace<be, dt>::ball_filter(
double k2){
if (k2 > 0)
{
double argument = sqrt(k2)*prefactor;
double argument = sqrt(k2)*ell / 2;
double prefactor = prefactor0 / pow(k2, 1.5);
for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= sin(argument) / argument;
((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= (
prefactor *
(sin(argument) - argument * cos(argument)));
}
});
}
......@@ -316,16 +326,38 @@ void kspace<be, dt>::Gauss_filter(
*
* This is a wrapper that can choose between a sharp Fourier spherical filter,
* a Gaussian filter and a sharp real space spherical filter.
* The cutoff wavenumber \f$k_c\f$ is a parameter, so the low pass Fourier
* operation is straightforward.
*
* For the Gaussian filter, it's mathematical expression in Fourier space is
* as follows:
* Filter expressions in real space are as follows:
* \f{eqnarray*}{
* \phi^b_\ell(r) &=&
* \frac{1}{\ell^3}\frac{6}{\pi} H(\ell/2 - r) \\
* \phi^g_\ell(r) &=&
* \frac{1}{\sigma_\ell^3}\frac{1}{(2\pi)^{3/2}}
* \exp\left(-\frac{1}{2}\left(\frac{r}{\sigma_\ell}\right)^2\right) \\
* \phi^s_\ell(r) &=&
* \frac{1}{2 \pi^2 r^3}
* \left(\sin k_\ell r - k_\ell r \cos k_\ell r\right)
* \f}
* and the corresponding expressions in Fourier space are:
* \f{eqnarray*}{
* \hat{\phi^b_\ell}(k) &=&
* \frac{3}{2(k\ell/2)^3}
* \left(2\sin (k \ell/2) - k \ell \cos (k \ell/2)\right) \\
* \hat{\phi^g_\ell}(k) &=&
* \exp\left(-\frac{1}{2}k^2 \sigma_\ell^2\right) \\
* \hat{\phi^s_\ell}(k) &=& H(k_\ell - k)
* \f}
*
* \f$ k_\ell \f$ is given as a parameter, and then we use
* \f[
* \hat{g}_\ell(\mathbf{k}) = \exp(-k^2 \ell^2 / 2)
* \ell = \pi / k_\ell,
* \sigma_\ell = \pi / k_\ell
* \f]
* And we choose the convention \f$\ell = \frac{\pi}{k_c}\f$.
* This is the same convention used in \cite Buzzicotti2017 .
*
* For the Gaussian filter this is the same convention used in
* \cite Buzzicotti2017 .
*
* See also `filter_calibrated_ell`.
*/
template <field_backend be,
kspace_dealias_type dt>
......@@ -357,6 +389,75 @@ int kspace<be, dt>::filter(
return EXIT_SUCCESS;
}
/** \brief Filter a field.
*
* This is a wrapper that can choose between a sharp Fourier spherical filter,
* a Gaussian filter and a sharp real space spherical filter.
*
* Filter expressions in real space are as follows:
* \f{eqnarray*}{
* \phi^b_\ell(r) &=&
* \frac{1}{\ell^3}\frac{6}{\pi} H(\ell/2 - r) \\
* \phi^g_\ell(r) &=&
* \frac{1}{\sigma_\ell^3}\frac{1}{(2\pi)^{3/2}}
* \exp\left(-\frac{1}{2}\left(\frac{r}{\sigma_\ell}\right)^2\right) \\
* \phi^s_\ell(r) &=&
* \frac{1}{2 \pi^2 r^3}
* \left(\sin k_\ell r - k_\ell r \cos k_\ell r\right)
* \f}
* and the corresponding expressions in Fourier space are:
* \f{eqnarray*}{
* \hat{\phi^b_\ell}(k) &=&
* \frac{3}{2(k\ell/2)^3}
* \left(2\sin (k \ell/2) - k \ell \cos (k \ell/2)\right) \\
* \hat{\phi^g_\ell}(k) &=&
* \exp\left(-\frac{1}{2}k^2 \sigma_\ell^2\right) \\
* \hat{\phi^s_\ell}(k) &=& H(k_\ell - k)
* \f}
*
* \f$\sigma_\ell\f$ and \f$k_\ell\f$ are calibrated such that the energy of
* the large scales is approximately the same (within the inertial range)
* independently of the shape of the filter.
*
* This was done by hand, see [INSERT CITATION HERE] for details, with the
* results:
*
* \f[
* \sigma_\ell = 0.23 \ell,
* k_\ell = 2.8 / \ell
* \f]
*
*/
template <field_backend be,
kspace_dealias_type dt>
template <typename rnumber,
field_components fc>
int kspace<be, dt>::filter_calibrated_ell(
typename fftw_interface<rnumber>::complex *__restrict__ a,
const double ell,
std::string filter_type)
{
if (filter_type == std::string("sharp_Fourier_sphere"))
{
this->template low_pass<rnumber, fc>(
a,
2.8 / ell);
}
else if (filter_type == std::string("Gauss"))
{
this->template Gauss_filter<rnumber, fc>(
a,
0.23*ell);
}
else if (filter_type == std::string("ball"))
{
this->template ball_filter<rnumber, fc>(
a,
ell);
}
return EXIT_SUCCESS;
}
template <field_backend be,
kspace_dealias_type dt>
template <typename rnumber,
......@@ -586,6 +687,32 @@ template int kspace<FFTW, SMOOTH>::filter<double, THREExTHREE>(
const double kmax,
std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<float, ONE>(
typename fftw_interface<float>::complex *__restrict__ a,
const double kmax,
std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<float, THREE>(
typename fftw_interface<float>::complex *__restrict__ a,
const double kmax,
std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<float, THREExTHREE>(
typename fftw_interface<float>::complex *__restrict__ a,
const double kmax,
std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<double, ONE>(
typename fftw_interface<double>::complex *__restrict__ a,
const double kmax,
std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<double, THREE>(
typename fftw_interface<double>::complex *__restrict__ a,
const double kmax,
std::string filter_type);
template int kspace<FFTW, SMOOTH>::filter_calibrated_ell<double, THREExTHREE>(
typename fftw_interface<double>::complex *__restrict__ a,
const double kmax,
std::string filter_type);
template void kspace<FFTW, SMOOTH>::dealias<float, ONE>(
typename fftw_interface<float>::complex *__restrict__ a);
template void kspace<FFTW, SMOOTH>::dealias<float, THREE>(
......
bfps/cpp/kspace.hpp
View file @ 07a0d81e
......@@ -95,6 +95,13 @@ class kspace
const double wavenumber,
std::string filter_type = std::string("Gauss"));
template <typename rnumber,
field_components fc>
int filter_calibrated_ell(
typename fftw_interface<rnumber>::complex *__restrict__ a,
const double wavenumber,
std::string filter_type = std::string("Gauss"));
template <typename rnumber,
field_components fc>
void dealias(typename fftw_interface<rnumber>::complex *__restrict__ a);
......
setup.py
View file @ 07a0d81e
......@@ -88,7 +88,8 @@ print('This is bfps version ' + VERSION)
### lists of files and MANIFEST.in
src_file_list = ['full_code/filter_test',
src_file_list = ['full_code/test',
'full_code/filter_test',
'hdf5_tools',
'full_code/get_rfields',
'full_code/NSVE_field_stats',
......@@ -97,7 +98,6 @@ src_file_list = ['full_code/filter_test',
'full_code/code_base',
'full_code/direct_numerical_simulation',
'full_code/NSVE',
'full_code/NSVEparticles',
'field_binary_IO',
'vorticity_equation',
'field',
......@@ -125,7 +125,8 @@ src_file_list = ['full_code/filter_test',
'spline_n9',
'spline_n10',
'Lagrange_polys',
'scope_timer']
'scope_timer',
'full_code/NSVEparticles']
particle_headers = [
'cpp/particles/particles_distr_mpi.hpp',
......
tests/test_filters.py 0 → 100644
View file @ 07a0d81e
#######################################################################
# #
# Copyright 2015 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 #
# #
#######################################################################
# relevant for results of "bfps TEST filter_test"
import h5py
import numpy as np
import matplotlib.pyplot as plt
def phi_b(
x, ell):
phi = (6. / (np.pi * ell**3)) * np.ones(x.shape, x.dtype)
bindices = np.where(np.abs(x) > ell/2)
phi[bindices] = 0
return phi
def hat_phi_b(
k, ell,
prefactor = 1.0):