diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9dae34e5e18ae4d9fb7106aea800412b0ca7bc5f..bfc65f8a34e21538c41ff111b8325936f142f17e 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -43,9 +43,8 @@ build-gcc-openmpi:
- *export_GCC_compilers
- *build
variables:
- COMPILER: "gcc"
+ COMPILER: "gcc/11"
MPI: "openmpi"
- MPICXX: "mpicxx"
tags:
- docker
artifacts:
@@ -63,9 +62,8 @@ build-gcc-impi:
export MPI_HOME=${I_MPI_ROOT}
- *build
variables:
- COMPILER: "gcc"
- MPI: "impi"
- MPICXX: "mpigxx"
+ COMPILER: "gcc/11"
+ MPI: "impi/2021.5"
tags:
- docker
artifacts:
@@ -83,9 +81,8 @@ build-intel:
export MPI_HOME=${I_MPI_ROOT}
- *build
variables:
- COMPILER: "intel"
- MPI: "impi"
- MPICXX: "mpiicpc"
+ COMPILER: "intel/21.5.0"
+ MPI: "impi/2021.5"
tags:
- docker
artifacts:
@@ -105,9 +102,8 @@ test-gcc-impi:
tags:
- docker
variables:
- COMPILER: "gcc"
- MPI: "impi"
- MPICXX: "mpigxx"
+ COMPILER: "gcc/11"
+ MPI: "impi/2021.5"
needs:
- job: build-gcc-impi
artifacts: true
@@ -124,7 +120,6 @@ test-gcc-impi:
# variables:
# COMPILER: "gcc"
# MPI: "openmpi"
-# MPICXX: "mpicxx"
# needs:
# - job: build-gcc-openmpi
# artifacts: true
@@ -140,9 +135,8 @@ test-intel:
tags:
- docker
variables:
- COMPILER: "intel"
- MPI: "impi"
- MPICXX: "mpiicpc"
+ COMPILER: "intel/21.5.0"
+ MPI: "impi/2021.5"
needs:
- job: build-intel
artifacts: true
@@ -153,7 +147,7 @@ build-doc:
- *load_modules
- mkdir build-doc
- cd build-doc
- - module load git gcc graphviz doxygen
+ - module load git gcc/11 graphviz doxygen
- *export_GCC_compilers
- pip install breathe
- yum -y install gd
@@ -163,9 +157,8 @@ build-doc:
tags:
- docker
variables:
- COMPILER: "gcc"
- MPI: "impi"
- MPICXX: "mpigxx"
+ COMPILER: "gcc/11"
+ MPI: "impi/2021.5"
artifacts:
paths:
- build-doc/
diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 4b26dd81590097d2f9b40a4fb677ecfbcbcaf745..d69af1ead1f0006a71c4b520a8ffa5a2df89a1f7 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -149,6 +149,10 @@ class DNS(_code):
self.parameter_description['famplitude'] = 'Forcing parameter: amplitude of Kolmogorov forcing, in code units.'
self.parameters['friction_coefficient'] = float(0.5)
self.parameter_description['friction_coefficient'] = 'Forcing parameter: drag coefficient, in code units.'
+ self.parameters['variation_strength'] = float(0.25)
+ self.parameter_description['variation_strength'] = 'Forcing parameter: amplitude of time variation of injection rate, in code units.'
+ self.parameters['variation_time_scale'] = float(1.)
+ self.parameter_description['variation_time_scale'] = 'Forcing parameter: time scale of variation of injection rate, in code units.'
self.parameters['energy'] = float(0.5)
self.parameter_description['energy'] = 'Forcing parameter: if fluid is forced by enforcing a constant energy, this is the value (in code units).'
self.parameters['injection_rate'] = float(0.4)
diff --git a/TurTLE/test/particle_set/NSVE_tracers_test.cpp b/TurTLE/test/particle_set/NSVE_tracers_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0f542e6221fbd518477ddc6b3d99a37bedde6fe0
--- /dev/null
+++ b/TurTLE/test/particle_set/NSVE_tracers_test.cpp
@@ -0,0 +1,292 @@
+/**********************************************************************
+* *
+* Copyright 2022 Max Planck Computing and Data Facility *
+* *
+* This file is part of TurTLE. *
+* *
+* TurTLE 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. *
+* *
+* TurTLE 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 TurTLE. If not, see <http://www.gnu.org/licenses/> *
+* *
+* Contact: Cristian.Lalescu@mpcdf.mpg.de *
+* *
+**********************************************************************/
+
+
+
+//#include "NSVE_tracers_test.hpp"
+
+#include <string>
+#include <cmath>
+
+using std::make_pair;
+
+template <typename rnumber>
+template <int neighbours, int smoothness>
+int NSVE_tracers_test<rnumber>::create_tracers()
+{
+ TIMEZONE("NSVE_tracers_test::create_solver");
+ assert(this->pset[make_pair(0, 0)] != nullptr);
+ this->pset[make_pair(neighbours, smoothness)] = new particle_set<3, neighbours, smoothness>(
+ this->fs->cvelocity->rlayout,
+ this->dkx,
+ this->dky,
+ this->dkz);
+ *(this->pset[make_pair(neighbours, smoothness)]) = this->pset[make_pair(0, 0)];
+ this->psolver[make_pair(neighbours, smoothness)] = new particle_solver(
+ *(this->pset[make_pair(neighbours, smoothness)]),
+ 0);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_tracers_test<rnumber>::initialize(void)
+{
+ TIMEZONE("NSVE_tracers_test::intialize");
+ this->NSVE<rnumber>::initialize();
+
+ // allocate field for acceleration
+ this->acceleration = new field<rnumber, FFTW, THREE>(
+ this->fs->cvelocity->rlayout->sizes[2],
+ this->fs->cvelocity->rlayout->sizes[1],
+ this->fs->cvelocity->rlayout->sizes[0],
+ this->fs->cvelocity->rlayout->comm,
+ this->fs->cvelocity->fftw_plan_rigor);
+
+ // allocate previous velocity
+ this->previous_velocity = new field<rnumber, FFTW, THREE>(
+ this->fs->cvelocity->rlayout->sizes[2],
+ this->fs->cvelocity->rlayout->sizes[1],
+ this->fs->cvelocity->rlayout->sizes[0],
+ this->fs->cvelocity->rlayout->comm,
+ this->fs->cvelocity->fftw_plan_rigor);
+ this->fs->compute_velocity(this->fs->cvorticity);
+ *this->previous_velocity = *this->fs->cvelocity;
+ this->previous_velocity->ift();
+
+ this->fti = new field_tinterpolator<rnumber, FFTW, THREE, LINEAR>();
+ this->trhs = new tracer_with_collision_counter_rhs<rnumber, FFTW, LINEAR>();
+
+ // set two fields for the temporal interpolation
+ this->fti->set_field(this->previous_velocity, 0);
+ this->fti->set_field(this->fs->cvelocity, 1);
+ this->trhs->setVelocity(this->fti);
+
+ // create linear interpolation particle set
+ this->pset[make_pair(0, 0)] = new particle_set<3, 0, 0>(
+ this->fs->cvelocity->rlayout,
+ this->dkx,
+ this->dky,
+ this->dkz);
+
+ // set up initial condition reader
+ particles_input_hdf5<long long int, double, 3, 3> particle_reader(
+ this->comm,
+ this->fs->get_current_fname(),
+ std::string("/tracers0/state/") + std::to_string(this->fs->iteration), // dataset name for initial input
+ std::string("/tracers0/rhs/") + std::to_string(this->fs->iteration), // dataset name for initial input
+ this->pset[make_pair(0, 0)]->getSpatialLowLimitZ(),
+ this->pset[make_pair(0, 0)]->getSpatialUpLimitZ());
+
+ // read initial condition with linear interpolation particle set
+ this->pset[make_pair(0, 0)]->init(particle_reader);
+
+ // create solver for linear interpolation particle set
+ this->psolver[make_pair(0, 0)] = new particle_solver(
+ *(this->pset[make_pair(0, 0)]),
+ 0);
+
+ // create Lagrange interpolation particles
+ this->template create_tracers<1, 0>();
+ this->template create_tracers<2, 0>();
+ this->template create_tracers<3, 0>();
+ this->template create_tracers<1, 1>();
+ this->template create_tracers<1, 2>();
+ this->template create_tracers<1, 3>();
+ this->template create_tracers<2, 3>();
+ this->template create_tracers<1, 4>();
+ this->template create_tracers<2, 4>();
+ this->template create_tracers<3, 4>();
+
+ this->particles_output_writer_mpi = new particles_output_hdf5<
+ long long int, double, 3>(
+ MPI_COMM_WORLD,
+ "tracers0",
+ nparticles,
+ 0);
+ this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+ long long int, double, double, 3>(
+ MPI_COMM_WORLD,
+ this->pset[make_pair(0, 0)]->getTotalNumberOfParticles(),
+ (this->simname + "_particles.h5"),
+ "tracers_n0_m0",
+ "position/0");
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_tracers_test<rnumber>::step(void)
+{
+ TIMEZONE("NSVE_tracers_test::step");
+ // compute next step Navier-Stokes
+ this->NSVE<rnumber>::step();
+ // update velocity 1
+ this->fs->compute_velocity(this->fs->cvorticity);
+ this->fs->cvelocity->ift();
+ // compute particle step using velocity 0 and velocity 1
+ for (auto it=this->psolver.begin(); it != this->psolver.end(); ++it)
+ {
+ it->second->cRK4(this->dt, *(this->trhs));
+ it->second->setIteration(this->fs->iteration);
+ }
+ // update velocity 0
+ *this->previous_velocity = *this->fs->cvelocity;
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_tracers_test<rnumber>::write_checkpoint(void)
+{
+ TIMEZONE("NSVE_tracers_test::write_checkpoint");
+ this->NSVE<rnumber>::write_checkpoint();
+
+ for (auto it=this->psolver.begin(); it != this->psolver.end(); ++it)
+ {
+ this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
+ this->particles_output_writer_mpi->update_particle_species_name(
+ "tracers_n" +
+ std::to_string(it->first.first) +
+ "_m" +
+ std::to_string(it->first.second));
+ it->second->setIteration(this->fs->iteration);
+ it->second->template writeCheckpoint<3>(this->particles_output_writer_mpi);
+ }
+ this->particles_output_writer_mpi->close_file();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_tracers_test<rnumber>::finalize(void)
+{
+ TIMEZONE("NSVE_tracers_test::finalize");
+ delete this->particles_sample_writer_mpi;
+ delete this->particles_output_writer_mpi;
+ for (auto it=this->psolver.begin(); it != this->psolver.end(); ++it)
+ {
+ delete it->second;
+ }
+ for (auto it=this->pset.begin(); it != this->pset.end(); ++it)
+ {
+ delete it->second;
+ }
+ delete this->trhs;
+ delete this->fti;
+ delete this->previous_velocity;
+ delete this->acceleration;
+ this->NSVE<rnumber>::finalize();
+ return EXIT_SUCCESS;
+}
+
+/** \brief Compute fluid stats and sample fields at particle locations.
+ */
+
+template <typename rnumber>
+int NSVE_tracers_test<rnumber>::do_stats()
+{
+ TIMEZONE("NSVE_tracers_test::do_stats");
+ /// fluid stats go here
+ this->NSVE<rnumber>::do_stats();
+
+
+ /// either one of two conditions suffices to compute statistics:
+ /// 1) current iteration is a multiple of niter_part
+ /// 2) we are within niter_part_fine_duration/2 of a multiple of niter_part_fine_period
+ if (!(this->iteration % this->niter_part == 0 ||
+ ((this->iteration + this->niter_part_fine_duration/2) % this->niter_part_fine_period <=
+ this->niter_part_fine_duration)))
+ return EXIT_SUCCESS;
+
+ // ensure velocity field is computed, and is in real space
+ if (!(this->iteration % this->niter_stat == 0))
+ {
+ // we need to compute velocity field manually, because it didn't happen in NSVE::do_stats()
+ this->fs->compute_velocity(this->fs->cvorticity);
+ *this->tmp_vec_field = *(this->fs->cvelocity);
+ this->tmp_vec_field->ift();
+ }
+
+ // compute acceleration
+ this->fs->compute_Lagrangian_acceleration(this->acceleration);
+ this->acceleration->ift();
+ DEBUG_MSG("hello, finished with acceleration\n");
+
+ for (auto it=this->pset.begin(); it != this->pset.end(); ++it)
+ {
+ const int neighbours = it->first.first;
+ const int smoothness = it->first.second;
+ abstract_particle_set *ps = it->second;
+ const std::string species_name = (
+ "tracers_n" +
+ std::to_string(neighbours) +
+ "_m" +
+ std::to_string(smoothness));
+ DEBUG_MSG("hello, doing %s\n", species_name.c_str());
+
+ // sample position
+ ps->writeStateTriplet(
+ 0,
+ this->particles_sample_writer_mpi,
+ species_name,
+ "position",
+ this->iteration);
+ MPI_Barrier(this->comm);
+
+ // sample velocity
+ ps->writeSample(
+ this->tmp_vec_field,
+ this->particles_sample_writer_mpi,
+ species_name,
+ "velocity",
+ this->iteration);
+ MPI_Barrier(this->comm);
+
+ // sample velocity
+ ps->writeSample(
+ this->acceleration,
+ this->particles_sample_writer_mpi,
+ species_name,
+ "acceleration",
+ this->iteration);
+ MPI_Barrier(this->comm);
+ }
+
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_tracers_test<rnumber>::read_parameters(void)
+{
+ TIMEZONE("NSVE_tracers_test::read_parameters");
+ this->NSVE<rnumber>::read_parameters();
+ hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+ this->niter_part = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part");
+ this->niter_part_fine_period = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_period");
+ this->niter_part_fine_duration = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_duration");
+ this->nparticles = hdf5_tools::read_value<int>(parameter_file, "parameters/nparticles");
+ H5Fclose(parameter_file);
+ MPI_Barrier(this->comm);
+ return EXIT_SUCCESS;
+}
+
+template class NSVE_tracers_test<float>;
+template class NSVE_tracers_test<double>;
+
diff --git a/TurTLE/test/particle_set/NSVE_tracers_test.hpp b/TurTLE/test/particle_set/NSVE_tracers_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..01ed7d983ac7b26149fa5e5412678e107c6a5f85
--- /dev/null
+++ b/TurTLE/test/particle_set/NSVE_tracers_test.hpp
@@ -0,0 +1,99 @@
+/**********************************************************************
+* *
+* Copyright 2022 Max Planck Computing and Data Facility *
+* *
+* This file is part of TurTLE. *
+* *
+* TurTLE 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. *
+* *
+* TurTLE 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 TurTLE. If not, see <http://www.gnu.org/licenses/> *
+* *
+* Contact: Cristian.Lalescu@mpcdf.mpg.de *
+* *
+**********************************************************************/
+
+
+
+#ifndef NSVE_TRACERS_TEST_HPP
+#define NSVE_TRACERS_TEST_HPP
+
+
+
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+#include "particles/interpolation/field_tinterpolator.hpp"
+#include "particles/interpolation/particle_set.hpp"
+#include "particles/particle_solver.hpp"
+#include "particles/rhs/tracer_with_collision_counter_rhs.hpp"
+
+#include <cstdlib>
+#include <map>
+
+/** \brief Test of tracer integration.
+ *
+ * Child of Navier Stokes vorticity equation solver, this class calls all the
+ * methods from `NSVE`, and in addition integrates simple Lagrangian tracers
+ * in the resulting velocity field.
+ *
+ * The code is meant to compare results of using different interpolation methods.
+ */
+
+template <typename rnumber>
+class NSVE_tracers_test: public NSVE<rnumber>
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ int niter_part;
+ int niter_part_fine_period;
+ int niter_part_fine_duration;
+ long long int nparticles;
+
+ /* other stuff */
+ field_tinterpolator<rnumber, FFTW, THREE, LINEAR> *fti;
+ tracer_rhs<rnumber, FFTW, LINEAR> *trhs;
+
+ std::map<std::pair<int, int>, abstract_particle_set *> pset;
+ std::map<std::pair<int, int>, particle_solver *> psolver;
+
+ field<rnumber, FFTW, THREE> *previous_velocity;
+ field<rnumber, FFTW, THREE> *acceleration;
+
+ particles_output_hdf5<long long int, double, 3> *particles_output_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
+
+ NSVE_tracers_test(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ NSVE<rnumber>(
+ COMMUNICATOR,
+ simulation_name){}
+ ~NSVE_tracers_test(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+
+ template <int neighbours, int smoothness>
+ int create_tracers();
+};
+
+#endif//NSVE_TRACERS_TEST_HPP
+
diff --git a/TurTLE/test/test_tracer_integration.py b/TurTLE/test/test_tracer_integration.py
new file mode 100644
index 0000000000000000000000000000000000000000..d9087e0d4a52aca9217862f8fd169f21a174b868
--- /dev/null
+++ b/TurTLE/test/test_tracer_integration.py
@@ -0,0 +1,182 @@
+import os
+import sys
+import argparse
+import getpass
+
+import numpy as np
+import h5py
+import matplotlib.pyplot as plt
+
+import TurTLE
+import TurTLE._base
+import TurTLE.DNS
+from TurTLE._code import _code
+
+cpp_location = os.path.join(
+ TurTLE.data_dir, 'particle_set')
+
+
+class ADNS(TurTLE.DNS):
+ def __init__(
+ self,
+ **kwargs):
+ TurTLE.DNS.__init__(self, **kwargs)
+ self.list_of_particle_codes.append('NSVE_tracers_test')
+ self.extra_parameters['NSVE_tracers_test'] = {}
+ return None
+ def write_src(
+ self):
+ self.version_message = (
+ '/***********************************************************************\n' +
+ '* this code automatically generated by TurTLE\n' +
+ '* version {0}\n'.format(TurTLE.__version__) +
+ '***********************************************************************/\n\n\n')
+ self.include_list = [
+ '"base.hpp"',
+ '"scope_timer.hpp"',
+ '"fftw_interface.hpp"',
+ '"full_code/main_code.hpp"',
+ '"full_code/NSVEparticles.hpp"',
+ '<cmath>',
+ '<iostream>',
+ '<hdf5.h>',
+ '<string>',
+ '<cstring>',
+ '<fftw3-mpi.h>',
+ '<omp.h>',
+ '<cfenv>',
+ '<cstdlib>']
+ 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])
+ self.name = 'NSVE_tracers_test'
+ self.dns_type = 'NSVE_tracers_test'
+ self.definitions += open(
+ os.path.join(
+ cpp_location, 'NSVE_tracers_test.hpp'), 'r').read()
+ self.definitions += open(
+ os.path.join(
+ cpp_location, 'NSVE_tracers_test.cpp'), 'r').read()
+ with open(self.name + '.cpp', 'w') as outfile:
+ outfile.write(self.version_message + '\n\n')
+ outfile.write(self.includes + '\n\n')
+ outfile.write(self.definitions + '\n\n')
+ outfile.write(self.main + '\n')
+ self.check_current_vorticity_exists = True
+ return None
+
+ def write_par(
+ self,
+ iter0 = 0):
+ TurTLE.DNS.write_par(self, iter0 = iter0)
+ with h5py.File(self.get_data_file_name(), 'r+') as ofile:
+ #TODO move collision stats to particle stats
+ for i in range(self.parameters['niter_out']):
+ ofile.create_dataset('statistics/collisions/'+str(i),
+ shape=(1,),
+ dtype = np.int64)
+ return None
+
+ def launch(
+ self,
+ args = [],
+ **kwargs):
+ opt = self.prepare_launch(args = args)
+ if not os.path.exists(self.get_data_file_name()):
+ self.write_par()
+ self.generate_initial_condition(opt = opt)
+ with h5py.File(self.get_particle_file_name(), 'w') as particle_file:
+ for (n, m) in [
+ (0, 0),
+ (1, 0), (2, 0), (3, 0),
+ (1, 1),
+ (1, 2),
+ (1, 3), (2, 3),
+ (1, 4), (2, 4), (3, 4)]:
+ particle_file.create_group('tracers_n{0}_m{1}/position'.format(n, m))
+ particle_file.create_group('tracers_n{0}_m{1}/velocity'.format(n, m))
+ particle_file.create_group('tracers_n{0}_m{1}/acceleration'.format(n, m))
+ self.generate_tracer_state(integration_steps = 0,
+ rseed = opt.particle_rand_seed,
+ ncomponents = 3)
+ 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,
+ no_debug = opt.no_debug)
+ return None
+
+def main():
+ bla = ADNS()
+ if not os.path.exists('test.h5'):
+ bla.launch([
+ 'NSVE_tracers_test',
+ '--niter_todo', '20',
+ '--niter_stat', '10',
+ '--niter_out', '20',
+ '--nparticles', '1000',
+ '--kMeta', '3.0',
+ ] + sys.argv[1:])
+ else:
+ bla.read_parameters()
+ df = bla.get_particle_file()
+ f = plt.figure(figsize=(16,4))
+ ax1 = f.add_subplot(141)
+ ax2 = f.add_subplot(142)
+ ax3 = f.add_subplot(143)
+ ax4 = f.add_subplot(144)
+ for species in df.keys():
+ iterations = range(bla.parameters['niter_out']//2, bla.parameters['niter_out'])
+ pos = np.array(
+ [df[species + '/position/{0}'.format(ii)][()]
+ for ii in iterations])
+ vel = np.array(
+ [df[species + '/velocity/{0}'.format(ii)][()]
+ for ii in iterations])
+ acc = np.array(
+ [df[species + '/acceleration/{0}'.format(ii)][()]
+ for ii in iterations])
+ d1x = (pos[2:] - pos[:-2]) / (2*bla.parameters['dt'])
+ d2x = (pos[2:] - 2*pos[1:-1] + pos[:-2]) / (bla.parameters['dt']**2)
+ d1v = (vel[2:] - vel[:-2]) / (2*bla.parameters['dt'])
+ vv = vel[1:-1]
+ aa = acc[1:-1]
+ errv = np.abs(np.sum((vv - d1x)**2, axis = 2) / np.sum(vv**2, axis = 2))**0.5
+ erra1 = np.abs(np.sum((aa - d1v)**2, axis = 2) / np.sum(aa**2, axis = 2))**0.5
+ erra2 = np.abs(np.sum((aa - d2x)**2, axis = 2) / np.sum(aa**2, axis = 2))**0.5
+ erra3 = np.abs(np.sum((d1v - d2x)**2, axis = 2) / np.sum(aa**2, axis = 2))**0.5
+ ax1.hist(errv.flatten(), label = species, bins = 101, histtype = 'step', density = True)
+ ax2.hist(erra1.flatten(), label = species, bins = 101, histtype = 'step', density = True)
+ ax3.hist(erra2.flatten(), label = species, bins = 101, histtype = 'step', density = True)
+ ax4.hist(erra3.flatten(), label = species, bins = 101, histtype = 'step', density = True)
+ for aa in [ax1, ax2, ax3, ax4]:
+ aa.set_xscale('log')
+ aa.set_yscale('log')
+ aa.legend(loc = 'best', fontsize = 6)
+ aa.tick_params(axis = 'both', which = 'both', labelsize = 6, pad = 0)
+ ax2.set_xlim(1e-4, 2)
+ ax3.set_xlim(1e-4, 2)
+ ax4.set_xlim(1e-4, 2)
+ ax1.set_title('vel vs d1x', fontsize = 6)
+ ax2.set_title('acc vs d1v', fontsize = 6)
+ ax3.set_title('acc vs d2x', fontsize = 6)
+ ax4.set_title('d1v vs d2x', fontsize = 6)
+ f.tight_layout()
+ f.savefig('tmp.pdf')
+ plt.close(f)
+ return None
+
+if __name__ == '__main__':
+ main()
diff --git a/cpp/full_code/NSVE.cpp b/cpp/full_code/NSVE.cpp
index 09ca7f286be9d1d17ba3bce5125187bf88fb8024..5c0d360c0cdced78c17b72b11b30d813239d75c8 100644
--- a/cpp/full_code/NSVE.cpp
+++ b/cpp/full_code/NSVE.cpp
@@ -77,10 +77,13 @@ int NSVE<rnumber>::initialize(void)
this->fs->fmode = fmode;
this->fs->famplitude = famplitude;
this->fs->friction_coefficient = friction_coefficient;
+ this->fs->variation_strength = variation_strength;
+ this->fs->variation_time_scale = variation_time_scale;
this->fs->energy = energy;
this->fs->injection_rate = injection_rate;
this->fs->fk0 = fk0;
this->fs->fk1 = fk1;
+ this->fs->dt = dt;
strncpy(this->fs->forcing_type, forcing_type, 128);
this->fs->iteration = this->iteration;
this->fs->checkpoint = this->checkpoint;
@@ -181,6 +184,8 @@ int NSVE<rnumber>::read_parameters(void)
this->fmode = hdf5_tools::read_value<int>(parameter_file, "parameters/fmode");
this->famplitude = hdf5_tools::read_value<double>(parameter_file, "parameters/famplitude");
this->friction_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/friction_coefficient");
+ this->variation_strength = hdf5_tools::read_value<double>(parameter_file, "parameters/variation_strength");
+ this->variation_time_scale = hdf5_tools::read_value<double>(parameter_file, "parameters/variation_time_scale");
this->injection_rate = hdf5_tools::read_value<double>(parameter_file, "parameters/injection_rate");
this->fk0 = hdf5_tools::read_value<double>(parameter_file, "parameters/fk0");
this->fk1 = hdf5_tools::read_value<double>(parameter_file, "parameters/fk1");
diff --git a/cpp/full_code/NSVE.hpp b/cpp/full_code/NSVE.hpp
index 2f29f91fef5c990f1676dc547e30458c0ff97a0b..15a388dc4f8e7fdd1db8e11b3cff1806094e18d3 100644
--- a/cpp/full_code/NSVE.hpp
+++ b/cpp/full_code/NSVE.hpp
@@ -51,6 +51,8 @@ class NSVE: public direct_numerical_simulation
double dt;
double famplitude;
double friction_coefficient;
+ double variation_strength;
+ double variation_time_scale;
double fk0;
double fk1;
double energy;
diff --git a/cpp/full_code/joint_acc_vel_stats.cpp b/cpp/full_code/joint_acc_vel_stats.cpp
index 8fe55d2025e1cf7c0c1b91144db2aab56ba6b4bc..1daf7ed0f3d6b1cc26be4854ca0bb33b93bd0d9a 100644
--- a/cpp/full_code/joint_acc_vel_stats.cpp
+++ b/cpp/full_code/joint_acc_vel_stats.cpp
@@ -51,14 +51,10 @@ int joint_acc_vel_stats<rnumber>::initialize(void)
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
H5P_DEFAULT);
- hid_t 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);
+ this->niter_out = hdf5_tools::read_value<int>(parameter_file, "/parameters/niter_out");
if (H5Lexists(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT))
{
- dset = H5Dopen(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT);
- H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->checkpoints_per_file);
- H5Dclose(dset);
+ this->checkpoints_per_file = hdf5_tools::read_value<int>(parameter_file, "/parameters/checkpoints_per_file");
}
else
this->checkpoints_per_file = 1;
@@ -72,12 +68,8 @@ int joint_acc_vel_stats<rnumber>::initialize(void)
this->iteration_list = hdf5_tools::read_vector<int>(
parameter_file,
"/joint_acc_vel_stats/parameters/iteration_list");
- dset = H5Dopen(parameter_file, "joint_acc_vel_stats/parameters/max_acceleration_estimate", H5P_DEFAULT);
- H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->max_acceleration_estimate);
- H5Dclose(dset);
- dset = H5Dopen(parameter_file, "joint_acc_vel_stats/parameters/max_velocity_estimate", H5P_DEFAULT);
- H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->max_velocity_estimate);
- H5Dclose(dset);
+ this->max_acceleration_estimate = hdf5_tools::read_value<int>(parameter_file, "joint_acc_vel_stats/parameters/max_acceleration_estimate");
+ this->max_velocity_estimate = hdf5_tools::read_value<int>(parameter_file, "joint_acc_vel_stats/parameters/max_velocity_estimate");
H5Fclose(parameter_file);
// the following ensures the file is free for rank 0 to open in read/write mode
MPI_Barrier(this->comm);
@@ -196,7 +188,7 @@ int joint_acc_vel_stats<rnumber>::work_on_current_iteration(void)
template <typename rnumber>
int joint_acc_vel_stats<rnumber>::finalize(void)
{
- DEBUG_MSG("entered joint_acc_vel_stats::finalize\n");
+ TIMEZONE("joint_acc_vel_stats::finalize");
delete this->ve;
delete this->kk;
if (this->myrank == 0)
diff --git a/cpp/full_code/postprocess.cpp b/cpp/full_code/postprocess.cpp
index 0f80890c144ba6a64ffcfbb72e9b24ede476aea2..2fe20c34ef524419e6e591b6d381421f12944d3b 100644
--- a/cpp/full_code/postprocess.cpp
+++ b/cpp/full_code/postprocess.cpp
@@ -69,6 +69,8 @@ int postprocess::read_parameters()
this->friction_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/friction_coefficient");
this->fk0 = hdf5_tools::read_value<double>(parameter_file, "parameters/fk0");
this->fk1 = hdf5_tools::read_value<double>(parameter_file, "parameters/fk1");
+ this->variation_strength = hdf5_tools::read_value<double>(parameter_file, "parameters/variation_strength");
+ this->variation_time_scale = hdf5_tools::read_value<double>(parameter_file, "parameters/variation_time_scale");
this->energy = hdf5_tools::read_value<double>(parameter_file, "parameters/energy");
this->injection_rate = hdf5_tools::read_value<double>(parameter_file, "parameters/injection_rate");
std::string tmp = hdf5_tools::read_string(parameter_file, "parameters/forcing_type");
@@ -83,6 +85,7 @@ template <typename rnumber,
int postprocess::copy_parameters_into(
vorticity_equation<rnumber, be> *dst)
{
+ TIMEZONE("postprocess::copy_parameters_into");
dst->nu = this->nu;
dst->fmode = this->fmode;
dst->famplitude = this->famplitude;
@@ -91,7 +94,9 @@ int postprocess::copy_parameters_into(
dst->injection_rate = this->injection_rate;
dst->energy = this->energy;
dst->friction_coefficient = this->friction_coefficient;
- strncpy(this->forcing_type, dst->forcing_type, 128);
+ dst->variation_strength = this->variation_strength;
+ dst->variation_time_scale = this->variation_time_scale;
+ strncpy(dst->forcing_type, this->forcing_type, 128);
return EXIT_SUCCESS;
}
diff --git a/cpp/full_code/postprocess.hpp b/cpp/full_code/postprocess.hpp
index da22c91e9f7ca187cf7b3038994c0852e255924b..4a8fab396149ccddc008d15dedbc15ce968b8a70 100644
--- a/cpp/full_code/postprocess.hpp
+++ b/cpp/full_code/postprocess.hpp
@@ -53,6 +53,8 @@ class postprocess: public code_base
double fk1;
double energy;
double injection_rate;
+ double variation_strength;
+ double variation_time_scale;
int fmode;
char forcing_type[512];
double nu;
diff --git a/cpp/full_code/test_particle_integration.cpp b/cpp/full_code/test_particle_integration.cpp
index a1d3cbaf957ffb6dab9070746ac70d3d22f57695..d914aa921cb15fd34b21e326a77ac0d85ffbc19e 100644
--- a/cpp/full_code/test_particle_integration.cpp
+++ b/cpp/full_code/test_particle_integration.cpp
@@ -174,7 +174,7 @@ int test_particle_integration<rnumber>::do_work()
0.40, // default
3*3./2); // 3/2 to account for divfree call below
// 3 because I want a larger amplitude
- this->kk->force_divfree<rnumber>(this->velocity_back->get_cdata());
+ this->kk->template force_divfree<rnumber>(this->velocity_back->get_cdata());
// take field to real space
this->velocity_back->ift();
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index 95d94fa1798a806cc97fd92197e96757ba174564..af907584ddd55dedc7d8a6d5a6b59f6f01262928 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -26,14 +26,18 @@
#ifndef ABSTRACT_PARTICLE_SET_HPP
#define ABSTRACT_PARTICLE_SET_HPP
+
+
+#include "field.hpp"
+#include "particles/p2p/p2p_distr_mpi.hpp"
+#include "particles/particles_sampling.hpp"
+#include "particles/abstract_particles_input.hpp"
+
#include <array>
#include <cassert>
#include <vector>
#include <memory>
#include <string>
-#include "field.hpp"
-#include "particles/p2p/p2p_distr_mpi.hpp"
-#include "particles/particles_sampling.hpp"
@@ -77,6 +81,12 @@ class abstract_particle_set
virtual std::vector<hsize_t> getParticleFileLayout() = 0;
+ virtual particle_rnumber getSpatialLowLimitZ() const = 0;
+ virtual particle_rnumber getSpatialUpLimitZ() const = 0;
+ virtual int init(abstract_particles_input<partsize_t, particle_rnumber>& particles_input) = 0;
+ virtual int operator=(abstract_particle_set *src) = 0;
+
+
// get p2p computer
virtual p2p_distr_mpi<partsize_t, particle_rnumber>* getP2PComputer() = 0;
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 071614ad52369d596565a50fa2e8ac34d5c22e99..79335800ee1de0cbc77666fbf459eb3295c827b4 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -709,12 +709,12 @@ class particle_set: public abstract_particle_set
return EXIT_SUCCESS;
}
- particle_rnumber getSpatialLowLimitZ()
+ particle_rnumber getSpatialLowLimitZ() const
{
return this->pInterpolator.getSpatialLowLimitZ();
}
- particle_rnumber getSpatialUpLimitZ()
+ particle_rnumber getSpatialUpLimitZ() const
{
return this->pInterpolator.getSpatialUpLimitZ();
}
diff --git a/cpp/vorticity_equation.cpp b/cpp/vorticity_equation.cpp
index e6c500b90277f2d586bf92ef5aa3c645f2949c3b..c32c03319d1e35796a847a30d1b9590999d3678a 100644
--- a/cpp/vorticity_equation.cpp
+++ b/cpp/vorticity_equation.cpp
@@ -159,6 +159,7 @@ vorticity_equation<rnumber, be>::vorticity_equation(
this->friction_coefficient = 1.0;
this->fk0 = 2.0;
this->fk1 = 4.0;
+ this->dt = 0.1;
}
template <class rnumber,
@@ -284,7 +285,8 @@ template <class rnumber,
field_backend be>
void vorticity_equation<rnumber, be>::add_forcing(
field<rnumber, be, THREE> *dst,
- field<rnumber, be, THREE> *vort_field)
+ field<rnumber, be, THREE> *vort_field,
+ double t)
{
TIMEZONE("vorticity_equation::add_forcing");
if (strcmp(this->forcing_type, "Kolmogorov") == 0)
@@ -334,7 +336,8 @@ void vorticity_equation<rnumber, be>::add_forcing(
return;
}
if ((strcmp(this->forcing_type, "fixed_energy_injection_rate") == 0) ||
- (strcmp(this->forcing_type, "fixed_energy_injection_rate_and_drag") == 0))
+ (strcmp(this->forcing_type, "fixed_energy_injection_rate_and_drag") == 0) ||
+ (strcmp(this->forcing_type, "sinusoidal_energy_injection_rate") == 0))
{
// first, compute energy in shell
shared_array<double> local_energy_in_shell(1);
@@ -373,7 +376,14 @@ void vorticity_equation<rnumber, be>::add_forcing(
// now, modify amplitudes
if (energy_in_shell < 10*std::numeric_limits<rnumber>::epsilon())
energy_in_shell = 1;
- double temp_famplitude = this->injection_rate / energy_in_shell;
+ double temp_famplitude = 0;
+ if (strcmp(this->forcing_type, "sinusoidal_energy_injection_rate") == 0){
+ temp_famplitude = (this->injection_rate + this->variation_strength
+ *sin(t/this->variation_time_scale))
+ / energy_in_shell;
+ }else{
+ temp_famplitude = this->injection_rate / energy_in_shell;
+ }
this->add_field_band(
dst, vort_field,
this->fk0, this->fk1,
@@ -472,7 +482,7 @@ void vorticity_equation<rnumber, be>::impose_forcing(
template <class rnumber,
field_backend be>
void vorticity_equation<rnumber, be>::omega_nonlin(
- int src)
+ int src, double t)
{
//DEBUG_MSG("vorticity_equation::omega_nonlin(%d)\n", src);
assert(src >= 0 && src < 3);
@@ -529,7 +539,7 @@ void vorticity_equation<rnumber, be>::omega_nonlin(
this->u->cval(cindex, cc, i) = tmp[cc][i];
}
);
- this->add_forcing(this->u, this->v[src]);
+ this->add_forcing(this->u, this->v[src], t);
this->kk->template force_divfree<rnumber>(this->u->get_cdata());
this->u->symmetrize();
}
@@ -538,10 +548,15 @@ template <class rnumber,
field_backend be>
void vorticity_equation<rnumber, be>::step(double dt)
{
- //DEBUG_MSG("vorticity_equation::step\n");
+ // This is an implementation of a memory-saving, third-order
+ // Runge-Kutta method. See also
+ // https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Explicit_methods
+ // where it is called "Third-order Strong Stability Preserving Runge-Kutta (SSPRK3)".
+ // The intermediate values of t passed to `omega_nonlin` can be read from
+ // the corresponding Butcher tableau.
TIMEZONE("vorticity_equation::step");
*this->v[1] = 0.0;
- this->omega_nonlin(0);
+ this->omega_nonlin(0, this->iteration*dt);
this->kk->CLOOP_K2(
[&](const ptrdiff_t cindex,
const ptrdiff_t xindex,
@@ -560,7 +575,7 @@ void vorticity_equation<rnumber, be>::step(double dt)
);
this->impose_forcing(this->v[1], this->v[0]);
- this->omega_nonlin(1);
+ this->omega_nonlin(1, (this->iteration + 1)*dt);
this->kk->CLOOP_K2(
[&](const ptrdiff_t cindex,
const ptrdiff_t xindex,
@@ -581,7 +596,7 @@ void vorticity_equation<rnumber, be>::step(double dt)
);
this->impose_forcing(this->v[2], this->v[0]);
- this->omega_nonlin(2);
+ this->omega_nonlin(2, (this->iteration + 1./2.)*dt);
// store old vorticity
*this->v[1] = *this->v[0];
this->kk->CLOOP_K2(
@@ -719,6 +734,9 @@ void vorticity_equation<rnumber, be>::compute_Lagrangian_acceleration(
this->compute_velocity(this->cvorticity);
acceleration->real_space_representation = false;
*acceleration = 0.0;
+ // add curl of forcing
+ this->add_forcing(acceleration, this->cvorticity, this->iteration*this->dt);
+ // loop over all modes, inverting curl and then adding the rest of the terms
this->kk->CLOOP_K2(
[&](const ptrdiff_t cindex,
const ptrdiff_t xindex,
@@ -727,27 +745,36 @@ void vorticity_equation<rnumber, be>::compute_Lagrangian_acceleration(
const double k2){
if (k2 <= this->kk->kM2)
{
- ptrdiff_t tindex = 3*cindex;
- for (int cc=0; cc<3; cc++)
- for (int i=0; i<2; i++)
- acceleration->get_cdata()[tindex+cc][i] = \
- - this->nu*k2*this->cvelocity->get_cdata()[tindex+cc][i];
- if (strcmp(this->forcing_type, "linear") == 0)
+ // first invert curl
+ if (k2 > 0)
{
- double knorm = sqrt(k2);
- if ((this->fk0 <= knorm) &&
- (this->fk1 >= knorm))
- for (int c=0; c<3; c++)
- for (int i=0; i<2; i++)
- acceleration->get_cdata()[tindex+c][i] += \
- this->famplitude*this->cvelocity->get_cdata()[tindex+c][i];
+ const rnumber rax = acceleration->cval(cindex, 0, 0);
+ const rnumber cax = acceleration->cval(cindex, 0, 1);
+ const rnumber ray = acceleration->cval(cindex, 1, 0);
+ const rnumber cay = acceleration->cval(cindex, 1, 1);
+ const rnumber raz = acceleration->cval(cindex, 2, 0);
+ const rnumber caz = acceleration->cval(cindex, 2, 1);
+ acceleration->cval(cindex,0,0) = -(kk->ky[yindex]*caz - kk->kz[zindex]*cay) / k2;
+ acceleration->cval(cindex,0,1) = (kk->ky[yindex]*raz - kk->kz[zindex]*ray) / k2;
+ acceleration->cval(cindex,1,0) = -(kk->kz[zindex]*cax - kk->kx[xindex]*caz) / k2;
+ acceleration->cval(cindex,1,1) = (kk->kz[zindex]*rax - kk->kx[xindex]*raz) / k2;
+ acceleration->cval(cindex,2,0) = -(kk->kx[xindex]*cay - kk->ky[yindex]*cax) / k2;
+ acceleration->cval(cindex,2,1) = (kk->kx[xindex]*ray - kk->ky[yindex]*rax) / k2;
}
- acceleration->get_cdata()[tindex+0][0] += this->kk->kx[xindex]*pressure->get_cdata()[cindex][1];
- acceleration->get_cdata()[tindex+1][0] += this->kk->ky[yindex]*pressure->get_cdata()[cindex][1];
- acceleration->get_cdata()[tindex+2][0] += this->kk->kz[zindex]*pressure->get_cdata()[cindex][1];
- acceleration->get_cdata()[tindex+0][1] -= this->kk->kx[xindex]*pressure->get_cdata()[cindex][0];
- acceleration->get_cdata()[tindex+1][1] -= this->kk->ky[yindex]*pressure->get_cdata()[cindex][0];
- acceleration->get_cdata()[tindex+2][1] -= this->kk->kz[zindex]*pressure->get_cdata()[cindex][0];
+
+ // add dissipative term
+ for (int cc=0; cc<3; cc++)
+ for (int i=0; i<2; i++)
+ acceleration->cval(cindex,cc, i) -= \
+ this->nu*k2*this->cvelocity->cval(cindex,cc,i);
+
+ // add nonlinear term
+ acceleration->cval(cindex,0,0) += this->kk->kx[xindex]*pressure->cval(cindex,1);
+ acceleration->cval(cindex,0,1) -= this->kk->kx[xindex]*pressure->cval(cindex,0);
+ acceleration->cval(cindex,1,0) += this->kk->ky[yindex]*pressure->cval(cindex,1);
+ acceleration->cval(cindex,1,1) -= this->kk->ky[yindex]*pressure->cval(cindex,0);
+ acceleration->cval(cindex,2,0) += this->kk->kz[zindex]*pressure->cval(cindex,1);
+ acceleration->cval(cindex,2,1) -= this->kk->kz[zindex]*pressure->cval(cindex,0);
}
});
if (own_pressure)
diff --git a/cpp/vorticity_equation.hpp b/cpp/vorticity_equation.hpp
index 32c6547b9781d23bf180d2b28a13d5cfdcd4e69a..90803b4b29a65951cac5afe93bfdd49cccf722da 100644
--- a/cpp/vorticity_equation.hpp
+++ b/cpp/vorticity_equation.hpp
@@ -70,12 +70,15 @@ class vorticity_equation
/* physical parameters */
double nu;
+ double dt;
int fmode; // for Kolmogorov flow
double famplitude; // both for Kflow and band forcing
double fk0, fk1; // for band forcing
double injection_rate; // for fixed energy injection rate
double energy; // for fixed energy
double friction_coefficient; // for Kolmogorov_and_drag
+ double variation_strength; // for time-varying forcing
+ double variation_time_scale; // for time-varying forcing
char forcing_type[128];
/* constructor, destructor */
@@ -91,7 +94,7 @@ class vorticity_equation
virtual ~vorticity_equation(void) noexcept(false);
/* solver essential methods */
- virtual void omega_nonlin(int src);
+ virtual void omega_nonlin(int src, double t = 0.);
virtual void step(double dt);
void impose_zero_modes(void);
@@ -103,7 +106,8 @@ class vorticity_equation
*
*/
void add_forcing(field<rnumber, be, THREE> *dst,
- field<rnumber, be, THREE> *src_vorticity);
+ field<rnumber, be, THREE> *src_vorticity,
+ double t = 0.);
void add_Kolmogorov_forcing(field<rnumber, be, THREE> *dst,
const int fmode,
diff --git a/documentation/index.html.in b/documentation/index.html.in
index 8614030531341f1e033c295771b5e1ccbc7793d7..e565b34a2913fbf66ddb4aa961e4047fd7da3567 100644
--- a/documentation/index.html.in
+++ b/documentation/index.html.in
@@ -17,13 +17,14 @@
and the
<a href="https://www.wilczek.physik.uni-bayreuth.de/en/team/index.html">
University of Bayreuth</a>
- in collaboration with the Application Support Group of the Max Planck
- Computing and Data Facility.
+ in collaboration with the Application Support Group of the <a href=https://www.mpcdf.mpg.de/>Max Planck
+ Computing and Data Facility</a>.
</p>
<p>
TurTLE is a framework for
pseudo-spectral simulations of turbulence, with a highly optimized particle tracking method.
It consists of a C++ core library, as well as a Python "wrapper" library.
+ TurTLE is open sourced (GPLv3), you may access it <a href="https://gitlab.mpcdf.mpg.de/TurTLE/turtle">here</a>.
</p>
<p>
The documentation consists of two parts: