From acd6f21a1ce942c7e95252c84e8a950031fd00b5 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Wed, 5 Jul 2017 14:53:41 +0200
Subject: [PATCH] add box filter and filter test
---
bfps/TEST.py | 300 +++++++++++++++++++++++++++
bfps/__main__.py | 7 +-
bfps/cpp/full_code/NSVEparticles.cpp | 2 +
bfps/cpp/full_code/filter_test.cpp | 181 ++++++++++++++++
bfps/cpp/full_code/filter_test.hpp | 69 ++++++
bfps/cpp/kspace.cpp | 42 +++-
bfps/cpp/kspace.hpp | 6 +
setup.py | 3 +-
8 files changed, 605 insertions(+), 5 deletions(-)
create mode 100644 bfps/TEST.py
create mode 100644 bfps/cpp/full_code/filter_test.cpp
create mode 100644 bfps/cpp/full_code/filter_test.hpp
diff --git a/bfps/TEST.py b/bfps/TEST.py
new file mode 100644
index 00000000..90c64d30
--- /dev/null
+++ b/bfps/TEST.py
@@ -0,0 +1,300 @@
+#######################################################################
+# #
+# 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 #
+# #
+#######################################################################
+
+
+
+import os
+import sys
+import shutil
+import subprocess
+import argparse
+import h5py
+import math
+import numpy as np
+import warnings
+
+import bfps
+from ._code import _code
+from bfps import tools
+
+class TEST(_code):
+ """This class is meant to stitch together the C++ code into a final source file,
+ compile it, and handle all job launching.
+ """
+ def __init__(
+ self,
+ work_dir = './',
+ simname = 'test'):
+ _code.__init__(
+ self,
+ work_dir = work_dir,
+ simname = simname)
+ self.host_info = {'type' : 'cluster',
+ 'environment' : None,
+ 'deltanprocs' : 1,
+ 'queue' : '',
+ 'mail_address': '',
+ 'mail_events' : None}
+ self.generate_default_parameters()
+ return None
+ def set_precision(
+ self,
+ fluid_dtype):
+ if fluid_dtype in [np.float32, np.float64]:
+ self.fluid_dtype = fluid_dtype
+ elif fluid_dtype in ['single', 'double']:
+ if fluid_dtype == 'single':
+ self.fluid_dtype = np.dtype(np.float32)
+ elif fluid_dtype == 'double':
+ self.fluid_dtype = np.dtype(np.float64)
+ self.rtype = self.fluid_dtype
+ if self.rtype == np.float32:
+ self.ctype = np.dtype(np.complex64)
+ self.C_field_dtype = 'float'
+ self.fluid_precision = 'single'
+ elif self.rtype == np.float64:
+ self.ctype = np.dtype(np.complex128)
+ self.C_field_dtype = 'double'
+ self.fluid_precision = 'double'
+ return None
+ def write_src(self):
+ self.version_message = (
+ '/***********************************************************************\n' +
+ '* this code automatically generated by bfps\n' +
+ '* version {0}\n'.format(bfps.__version__) +
+ '***********************************************************************/\n\n\n')
+ self.include_list = [
+ '"base.hpp"',
+ '"scope_timer.hpp"',
+ '"fftw_interface.hpp"',
+ '"full_code/main_code.hpp"',
+ '<cmath>',
+ '<iostream>',
+ '<hdf5.h>',
+ '<string>',
+ '<cstring>',
+ '<fftw3-mpi.h>',
+ '<omp.h>',
+ '<cfenv>',
+ '<cstdlib>',
+ '"full_code/{0}.hpp"\n'.format(self.dns_type)]
+ self.main = """
+ int main(int argc, char *argv[])
+ {{
+ bool fpe = (
+ (getenv("BFPS_FPE_OFF") == nullptr) ||
+ (getenv("BFPS_FPE_OFF") != std::string("TRUE")));
+ return main_code< {0} >(argc, argv, fpe);
+ }}
+ """.format(self.dns_type + '<{0}>'.format(self.C_field_dtype))
+ self.includes = '\n'.join(
+ ['#include ' + hh
+ for hh in self.include_list])
+ with open(self.name + '.cpp', 'w') as outfile:
+ outfile.write(self.version_message + '\n\n')
+ outfile.write(self.includes + '\n\n')
+ outfile.write(self.main + '\n')
+ return None
+ def generate_default_parameters(self):
+ # these parameters are relevant for all TEST classes
+ self.parameters['dealias_type'] = int(1)
+ self.parameters['dkx'] = float(1.0)
+ 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 = {}
+ if self.parameters['dealias_type'] == 1:
+ kMx = self.parameters['dkx']*(self.parameters['nx']//2 - 1)
+ kMy = self.parameters['dky']*(self.parameters['ny']//2 - 1)
+ kMz = self.parameters['dkz']*(self.parameters['nz']//2 - 1)
+ else:
+ kMx = self.parameters['dkx']*(self.parameters['nx']//3 - 1)
+ kMy = self.parameters['dky']*(self.parameters['ny']//3 - 1)
+ kMz = self.parameters['dkz']*(self.parameters['nz']//3 - 1)
+ kspace['kM'] = max(kMx, kMy, kMz)
+ kspace['dk'] = min(self.parameters['dkx'],
+ self.parameters['dky'],
+ self.parameters['dkz'])
+ nshells = int(kspace['kM'] / kspace['dk']) + 2
+ kspace['nshell'] = np.zeros(nshells, dtype = np.int64)
+ kspace['kshell'] = np.zeros(nshells, dtype = np.float64)
+ kspace['kx'] = np.arange( 0,
+ self.parameters['nx']//2 + 1).astype(np.float64)*self.parameters['dkx']
+ kspace['ky'] = np.arange(-self.parameters['ny']//2 + 1,
+ self.parameters['ny']//2 + 1).astype(np.float64)*self.parameters['dky']
+ kspace['ky'] = np.roll(kspace['ky'], self.parameters['ny']//2+1)
+ kspace['kz'] = np.arange(-self.parameters['nz']//2 + 1,
+ self.parameters['nz']//2 + 1).astype(np.float64)*self.parameters['dkz']
+ kspace['kz'] = np.roll(kspace['kz'], self.parameters['nz']//2+1)
+ return kspace
+ def get_data_file_name(self):
+ return os.path.join(self.work_dir, self.simname + '.h5')
+ def get_data_file(self):
+ return h5py.File(self.get_data_file_name(), 'r')
+ def write_par(
+ 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:
+ ofile['bfps_info/exec_name'] = self.name
+ kspace = self.get_kspace()
+ for k in kspace.keys():
+ ofile['kspace/' + k] = kspace[k]
+ nshells = kspace['nshell'].shape[0]
+ kspace = self.get_kspace()
+ nshells = kspace['nshell'].shape[0]
+ ofile['checkpoint'] = int(0)
+ return None
+ def job_parser_arguments(
+ self,
+ parser):
+ parser.add_argument(
+ '--ncpu',
+ type = int,
+ dest = 'ncpu',
+ default = -1)
+ parser.add_argument(
+ '--np', '--nprocesses',
+ metavar = 'NPROCESSES',
+ help = 'number of mpi processes to use',
+ type = int,
+ dest = 'nb_processes',
+ default = 4)
+ parser.add_argument(
+ '--ntpp', '--nthreads-per-process',
+ type = int,
+ dest = 'nb_threads_per_process',
+ metavar = 'NTHREADS_PER_PROCESS',
+ help = 'number of threads to use per MPI process',
+ default = 1)
+ parser.add_argument(
+ '--no-submit',
+ action = 'store_true',
+ dest = 'no_submit')
+ parser.add_argument(
+ '--environment',
+ type = str,
+ dest = 'environment',
+ default = None)
+ parser.add_argument(
+ '--minutes',
+ type = int,
+ dest = 'minutes',
+ default = 5,
+ help = 'If environment supports it, this is the requested wall-clock-limit.')
+ parser.add_argument(
+ '--njobs',
+ type = int, dest = 'njobs',
+ default = 1)
+ return None
+ def simulation_parser_arguments(
+ self,
+ parser):
+ parser.add_argument(
+ '--simname',
+ type = str, dest = 'simname',
+ default = 'test')
+ parser.add_argument(
+ '-n', '--grid-size',
+ type = int,
+ dest = 'n',
+ default = 32,
+ metavar = 'N',
+ help = 'code is run by default in a grid of NxNxN')
+ for coord in ['x', 'y', 'z']:
+ parser.add_argument(
+ '--L{0}'.format(coord), '--box-length-{0}'.format(coord),
+ type = float,
+ dest = 'L{0}'.format(coord),
+ default = 2.0,
+ metavar = 'length{0}'.format(coord),
+ help = 'length of the box in the {0} direction will be `length{0} x pi`'.format(coord))
+ parser.add_argument(
+ '--wd',
+ type = str, dest = 'work_dir',
+ default = './')
+ parser.add_argument(
+ '--precision',
+ choices = ['single', 'double'],
+ type = str,
+ default = 'single')
+ return None
+ def add_parser_arguments(
+ self,
+ parser):
+ subparsers = parser.add_subparsers(
+ dest = 'TEST_class',
+ help = 'type of simulation to run')
+ subparsers.required = True
+ parser_filter_test = subparsers.add_parser(
+ 'filter_test',
+ help = 'plain filter test')
+ self.simulation_parser_arguments(parser_filter_test)
+ self.job_parser_arguments(parser_filter_test)
+ self.parameters_to_parser_arguments(parser_filter_test)
+ return None
+ def prepare_launch(
+ self,
+ args = []):
+ opt = _code.prepare_launch(self, args = args)
+ self.set_precision(opt.precision)
+ self.dns_type = opt.TEST_class
+ self.name = self.dns_type + '-' + self.fluid_precision + '-v' + bfps.__version__
+ # merge parameters if needed
+ self.pars_from_namespace(opt)
+ return opt
+ def launch(
+ self,
+ args = [],
+ **kwargs):
+ opt = self.prepare_launch(args = args)
+ self.launch_jobs(opt = opt, **kwargs)
+ return None
+ def launch_jobs(
+ self,
+ opt = None,
+ particle_initial_condition = None):
+ if not os.path.exists(os.path.join(self.work_dir, self.simname + '.h5')):
+ self.write_par(
+ particle_ic = None)
+ self.run(
+ nb_processes = opt.nb_processes,
+ nb_threads_per_process = opt.nb_threads_per_process,
+ njobs = opt.njobs,
+ hours = opt.minutes // 60,
+ minutes = opt.minutes % 60,
+ no_submit = opt.no_submit)
+ return None
+
diff --git a/bfps/__main__.py b/bfps/__main__.py
index fa30289b..c41a6ffb 100644
--- a/bfps/__main__.py
+++ b/bfps/__main__.py
@@ -30,6 +30,7 @@ import argparse
import bfps
from .DNS import DNS
from .PP import PP
+from .TEST import TEST
from .NavierStokes import NavierStokes
from .NSVorticityEquation import NSVorticityEquation
from .FluidResize import FluidResize
@@ -66,7 +67,7 @@ def main():
'NSManyParticles-double']
parser.add_argument(
'base_class',
- choices = ['DNS', 'PP'] +
+ choices = ['DNS', 'PP', 'TEST'] +
NSoptions +
NSVEoptions +
FRoptions +
@@ -86,6 +87,10 @@ def main():
c = PP()
c.launch(args = sys.argv[2:])
return None
+ if opt.base_class == 'TEST':
+ c = TEST()
+ c.launch(args = sys.argv[2:])
+ return None
if 'double' in opt.base_class:
precision = 'double'
else:
diff --git a/bfps/cpp/full_code/NSVEparticles.cpp b/bfps/cpp/full_code/NSVEparticles.cpp
index ba84b394..cb9d17da 100644
--- a/bfps/cpp/full_code/NSVEparticles.cpp
+++ b/bfps/cpp/full_code/NSVEparticles.cpp
@@ -84,6 +84,7 @@ int NSVEparticles<rnumber>::do_stats()
"tracers0", // hdf5 parent group
"velocity" // dataset basename TODO
);
+ DEBUG_MSG("hello after sampling velocity\n");
/// compute acceleration and sample it
this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field);
@@ -93,6 +94,7 @@ int NSVEparticles<rnumber>::do_stats()
(this->simname + "_particles.h5"),
"tracers0",
"acceleration");
+ DEBUG_MSG("hello after sampling acceleration\n");
return EXIT_SUCCESS;
}
diff --git a/bfps/cpp/full_code/filter_test.cpp b/bfps/cpp/full_code/filter_test.cpp
new file mode 100644
index 00000000..6ed16373
--- /dev/null
+++ b/bfps/cpp/full_code/filter_test.cpp
@@ -0,0 +1,181 @@
+#include <string>
+#include <cmath>
+#include "filter_test.hpp"
+#include "scope_timer.hpp"
+
+
+template <typename rnumber>
+int filter_test<rnumber>::initialize(void)
+{
+ this->read_parameters();
+ this->scal_field = new field<rnumber, FFTW, ONE>(
+ nx, ny, nz,
+ this->comm,
+ DEFAULT_FFTW_FLAG);
+ this->kk = new kspace<FFTW, SMOOTH>(
+ this->scal_field->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->read_iteration();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int filter_test<rnumber>::reset_field(
+ int dimension)
+{
+ this->scal_field->real_space_representation = true;
+ *this->scal_field = 0.0;
+ if (this->scal_field->rlayout->starts[0] == 0)
+ {
+ switch(dimension)
+ {
+ case 0:
+ this->scal_field->rval(0) = 1.0 / (
+ (4*acos(0) / (this->nx*this->dkx))*
+ (4*acos(0) / (this->ny*this->dky))*
+ (4*acos(0) / (this->nz*this->dkz)));
+ break;
+ case 1:
+ for (ptrdiff_t xindex = 0; xindex < this->nx; xindex++)
+ this->scal_field->rval(this->scal_field->get_rindex(xindex, 0, 0)) = 1.0 / (
+ (4*acos(0) / (this->ny*this->dky))*
+ (4*acos(0) / (this->nz*this->dkz)));
+ break;
+ case 2:
+ for (ptrdiff_t yindex = 0; yindex < this->ny; yindex++)
+ for (ptrdiff_t xindex = 0; xindex < this->nx; xindex++)
+ {
+ this->scal_field->rval(
+ this->scal_field->get_rindex(xindex, yindex, 0)) = 1.0 / (
+ (4*acos(0) / (this->nz*this->dkz)));
+ }
+ break;
+ default:
+ break;
+ }
+ }
+ this->scal_field->dft();
+ this->scal_field->symmetrize();
+ return EXIT_SUCCESS;
+}
+
+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)
+{
+ std::string filename = this->simname + std::string("_fields.h5");
+ for (int dimension = 0; dimension < 3; dimension++)
+ {
+ this->reset_field(dimension);
+ this->kk->template filter<rnumber, ONE>(
+ this->scal_field->get_cdata(),
+ 4*acos(0) / (this->filter_length),
+ "sharp_Fourier_sphere");
+ this->scal_field->ift();
+ this->scal_field->normalize();
+ this->scal_field->io(
+ filename,
+ "sharp_Fourier_sphere",
+ dimension,
+ false);
+ this->reset_field(dimension);
+ this->kk->template filter<rnumber, ONE>(
+ this->scal_field->get_cdata(),
+ 4*acos(0) / (this->filter_length / sqrt(log(2))),
+ "Gauss");
+ this->scal_field->ift();
+ this->scal_field->normalize();
+ this->scal_field->io(
+ filename,
+ "Gauss",
+ dimension,
+ false);
+ this->reset_field(dimension);
+ this->kk->template filter<rnumber, ONE>(
+ this->scal_field->get_cdata(),
+ 4*acos(0) / (this->filter_length*2),
+ "ball");
+ this->scal_field->ift();
+ this->scal_field->normalize();
+ this->scal_field->io(
+ filename,
+ "ball",
+ dimension,
+ false);
+ }
+ 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>;
+
diff --git a/bfps/cpp/full_code/filter_test.hpp b/bfps/cpp/full_code/filter_test.hpp
new file mode 100644
index 00000000..83ea60e3
--- /dev/null
+++ b/bfps/cpp/full_code/filter_test.hpp
@@ -0,0 +1,69 @@
+/**********************************************************************
+* *
+* 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 FILTER_TEST_HPP
+#define FILTER_TEST_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/direct_numerical_simulation.hpp"
+
+template <typename rnumber>
+class filter_test: public direct_numerical_simulation
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ double filter_length;
+
+ /* other stuff */
+ kspace<FFTW, SMOOTH> *kk;
+ field<rnumber, FFTW, ONE> *scal_field;
+
+ filter_test(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ direct_numerical_simulation(
+ COMMUNICATOR,
+ simulation_name){}
+ ~filter_test(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ int reset_field(int dimension);
+
+ int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void){}
+};
+
+#endif//FILTER_TEST_HPP
+
diff --git a/bfps/cpp/kspace.cpp b/bfps/cpp/kspace.cpp
index 96425fd2..4c296fa2 100644
--- a/bfps/cpp/kspace.cpp
+++ b/bfps/cpp/kspace.cpp
@@ -252,6 +252,37 @@ 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:
+ * \f[
+ * \hat{b}_\ell(\mathbf{k}) = \sin(k \ell / 2) / (k \ell / 2)
+ * \f]
+ */
+template <field_backend be,
+ kspace_dealias_type dt>
+template <typename rnumber,
+ field_components fc>
+void kspace<be, dt>::ball_filter(
+ typename fftw_interface<rnumber>::complex *__restrict__ a,
+ const double sigma)
+{
+ const double prefactor = sigma/2;
+ this->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+ if (k2 > 0)
+ {
+ double argument = sqrt(k2)*prefactor;
+ for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
+ ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= sin(argument) / argument;
+ }
+ });
+}
+
/** \brief Filter a field using a Gaussian kernel.
*
* Filter's mathematical expression in Fourier space is as follows:
@@ -274,7 +305,6 @@ void kspace<be, dt>::Gauss_filter(
ptrdiff_t yindex,
ptrdiff_t zindex,
double k2){
- if (k2 <= this->kM2)
{
for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= exp(prefactor*k2);
@@ -284,8 +314,8 @@ void kspace<be, dt>::Gauss_filter(
/** \brief Filter a field.
*
- * This is a wrapper that can choose between a sharp Fourier spherical filter
- * and a Gaussian 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.
*
@@ -318,6 +348,12 @@ int kspace<be, dt>::filter(
a,
2*acos(0.)/wavenumber);
}
+ else if (filter_type == std::string("ball"))
+ {
+ this->template ball_filter<rnumber, fc>(
+ a,
+ 2*acos(0.)/wavenumber);
+ }
return EXIT_SUCCESS;
}
diff --git a/bfps/cpp/kspace.hpp b/bfps/cpp/kspace.hpp
index 7fa77e8b..eebec838 100644
--- a/bfps/cpp/kspace.hpp
+++ b/bfps/cpp/kspace.hpp
@@ -82,6 +82,12 @@ class kspace
typename fftw_interface<rnumber>::complex *__restrict__ a,
const double sigma);
+ template <typename rnumber,
+ field_components fc>
+ void ball_filter(
+ typename fftw_interface<rnumber>::complex *__restrict__ a,
+ const double sigma);
+
template <typename rnumber,
field_components fc>
int filter(
diff --git a/setup.py b/setup.py
index 1178f547..8562d375 100644
--- a/setup.py
+++ b/setup.py
@@ -88,7 +88,8 @@ print('This is bfps version ' + VERSION)
### lists of files and MANIFEST.in
-src_file_list = ['hdf5_tools',
+src_file_list = ['full_code/filter_test',
+ 'hdf5_tools',
'full_code/get_rfields',
'full_code/NSVE_field_stats',
'full_code/native_binary_to_hdf5',
--
GitLab