[wip] adds working folder for unit test of forcing

examples/forcing_unit_test/README.rst

+Forcing unit test
+=================
+
+For this test, velocity field is initialized as a single mode.
+Numerical solution is compared with analytic solution for the different
+available forcing methods.
+
+
+**Linear forcing**
+
+initial condition
+
+.. math::
+    u(x,0) = U (e^{i\bar{k}\cdot x} + e^{-i\bar{k}\cdot x})
+    U \perp \bar{k}
+
+solution
+
+.. math::
+    u(x,t) = exp((-\nu\bar{k}^2 + \alpha)t) u(x,0)
+
+where :math:`\alpha` is the linear coefficient.
+
+**Fixed energy injection rate**
+
+initial condition
+
+.. math::
+    u(x,0) = U (e^{i\bar{k}\cdot x} + e^{-i\bar{k}\cdot x})
+    U \perp \bar{k}
+
+solution has same form as for linear forcing, but
+
+.. math::
+    U^2(t) = V^2 + (U^2(0) - V^2) exp(-2\nu\bar{k}^2 t)
+
+with
+
+.. math::
+    V^2 = \epsilon / (2\nu\bar{k}^2)

examples/forcing_unit_test/fluid_solver_check.cpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2024 TurTLE team                                                 *
+*                                                                             *
+*  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 "fluid_solver_check.hpp"
+#include "scope_timer.hpp"
+#include "shared_array.hpp"
+#include "fftw_tools.hpp"
+#include "full_code/NSVE.hpp"
+#include "full_code/NSE.hpp"
+
+#include <string>
+#include <cmath>
+#include <filesystem>
+
+template <typename rnumber>
+int fluid_solver_check<rnumber>::initialize(void)
+{
+    TIMEZONE("fluid_solver_check::intialize");
+    this->nx = 32;
+    this->ny = 32;
+    this->nz = 32;
+    this->dkx = 1.0;
+    this->dky = 1.0;
+    this->dkz = 1.0;
+    this->dealias_type = SMOOTH;
+
+    this->spectrum_dissipation = 0.4;
+    this->spectrum_Lint = 1.3;
+    this->spectrum_etaK = 0.3;
+    this->spectrum_large_scale_const = 6.78;
+    this->spectrum_small_scale_const = 0.40;
+    this->random_seed = 2;
+
+    this->dt                   = 0.1;
+    this->famplitude           = 0.5;
+    this->friction_coefficient = 0.1;
+    this->fk0                  = 1.4;
+    this->fk1                  = 2.3;
+    this->energy               = 0.5;
+    this->injection_rate       = 0.4;
+    this->variation_strength   = 0.25;
+    this->variation_time_scale = 1.0;
+    this->fmode                = 2;
+    this->forcing_type         = "linear";
+    this->nu                   = 0.1;
+
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int fluid_solver_check<rnumber>::finalize(void)
+{
+    TIMEZONE("fluid_solver_check::finalize");
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber,
+          field_components fc>
+int get_error(
+        field<rnumber, FFTW, fc>* a,
+        field<rnumber, FFTW, fc>* b,
+        double& diff_L2,
+        double& diff_max,
+        double& aunit_L2,
+        double& bunit_L2)
+{
+    diff_L2  = 0;
+    diff_max = 0;
+    aunit_L2 = 0;
+    bunit_L2 = 0;
+    shared_array<double> aunit_L2_threaded(
+            1,
+            [&](double *val_tmp){
+                *val_tmp = double(0.0);
+            });
+    shared_array<double> bunit_L2_threaded(
+            1,
+            [&](double *val_tmp){
+                *val_tmp = double(0.0);
+            });
+    shared_array<double> diff_L2_threaded(
+            1,
+            [&](double *val_tmp){
+                *val_tmp = double(0.0);
+            });
+    shared_array<double> diff_max_threaded(
+            1,
+            [&](double *val_tmp){
+                *val_tmp = double(0.0);
+            });
+    switch (fc) {
+        case THREE:
+            {
+            a->RLOOP(
+                [&](ptrdiff_t rindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex){
+                for (auto cc = 0; cc < 3; cc++) {
+                    const double difference = std::abs(
+                            a->rval(rindex, cc) - b->rval(rindex, cc));
+                    double *diff_max = diff_max_threaded.getMine();
+                    double *diff_L2  = diff_L2_threaded.getMine();
+                    double *aunit_L2 = aunit_L2_threaded.getMine();
+                    double *bunit_L2 = bunit_L2_threaded.getMine();
+                    *diff_max  = std::max(*diff_max, difference);
+                    *diff_L2  += difference * difference;
+                    *aunit_L2 += a->rval(rindex, cc) * a->rval(rindex, cc);
+                    *bunit_L2 += b->rval(rindex, cc) * b->rval(rindex, cc);
+            }});
+            }
+            break;
+        case ONE:
+            {
+            a->RLOOP(
+                [&](ptrdiff_t rindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex){
+                const double difference =
+                    std::abs(a->rval(rindex) - b->rval(rindex));
+                double *diff_max = diff_max_threaded.getMine();
+                double *diff_L2  = diff_L2_threaded.getMine();
+                double *aunit_L2 = aunit_L2_threaded.getMine();
+                double *bunit_L2 = bunit_L2_threaded.getMine();
+                *diff_max  = std::max(*diff_max, difference);
+                *diff_L2  += difference * difference;
+                *aunit_L2 += (a->rval(rindex) * a->rval(rindex));
+                *bunit_L2 += (b->rval(rindex) * b->rval(rindex));
+                });
+            }
+            break;
+        }
+
+    aunit_L2_threaded.mergeParallel();
+    bunit_L2_threaded.mergeParallel();
+    diff_L2_threaded.mergeParallel();
+    diff_max_threaded.mergeParallel(
+            [&](const int idx, const double& v1, const double& v2) -> double {
+                return std::max(v1, v2);
+            });
+
+    MPI_Allreduce(
+            aunit_L2_threaded.getMasterData(),
+            &aunit_L2,
+            1,
+            MPI_DOUBLE,
+            MPI_SUM,
+            a->comm);
+    MPI_Allreduce(
+            bunit_L2_threaded.getMasterData(),
+            &bunit_L2,
+            1,
+            MPI_DOUBLE,
+            MPI_SUM,
+            a->comm);
+    MPI_Allreduce(
+            diff_L2_threaded.getMasterData(),
+            &diff_L2,
+            1,
+            MPI_DOUBLE,
+            MPI_SUM,
+            a->comm);
+    MPI_Allreduce(
+            diff_max_threaded.getMasterData(),
+            &diff_max,
+            1,
+            MPI_DOUBLE,
+            MPI_MAX,
+            a->comm);
+    diff_L2 = std::sqrt(diff_L2 / (ncomp(fc)*a->npoints));
+    aunit_L2 = std::sqrt(aunit_L2 / (ncomp(fc)*a->npoints));
+    bunit_L2 = std::sqrt(bunit_L2 / (ncomp(fc)*a->npoints));
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int compare_solvers(
+        NSE<rnumber>*  nse,
+        NSVE<rnumber>* nsve)
+{
+    assert(nse->get_nx() == nsve->get_nx());
+    assert(nse->get_ny() == nsve->get_ny());
+    assert(nse->get_nz() == nsve->get_nz());
+
+    nse->print_debug_info();
+    nsve->print_debug_info();
+
+    // fix initial conditions
+    compute_curl(
+        nse->kk,
+        nse->velocity,
+        nsve->fs->cvorticity);
+    *nsve->fs->cvelocity = 0.0;
+
+    // step
+    //nse->Euler_step();
+    //nsve->fs->Euler_step(nsve->fs->dt);
+    nse->step();
+    nsve->fs->step(nsve->fs->dt);
+
+    // compare results
+    double diff_L2, diff_max;
+    double nse_uL2, nsve_uL2;
+    nsve->fs->compute_velocity(nsve->fs->cvorticity);
+
+    nse->velocity->ift();
+    nsve->fs->cvelocity->ift();
+    get_error(
+        nse->velocity,
+        nsve->fs->cvelocity,
+        diff_L2,
+        diff_max,
+        nse_uL2,
+        nsve_uL2);
+    DEBUG_MSG("#### compare solvers information\n");
+    DEBUG_MSG("nse->forcing_type = %s, nsve->forcing_type = %s\n",
+            nse->forcing_type.c_str(), nsve->fs->forcing_type.c_str());
+    DEBUG_MSG("diff_L2 = %g, diff_max = %g, nse_uL2 = %g, nsve_uL2 = %g\n",
+            diff_L2, diff_max, nse_uL2, nsve_uL2);
+    DEBUG_MSG("diff_L2 / nse_uL2 = %g, diff_max / nse_uL2 = %g\n",
+            diff_L2 / nse_uL2, diff_max / nse_uL2);
+    const bool difference_small =
+        diff_max / std::sqrt(nse_uL2*nsve_uL2) < 1e-5;
+    if (difference_small)
+        return EXIT_SUCCESS;
+    else
+        return EXIT_FAILURE;
+}
+
+template <typename rnumber>
+int fluid_solver_check<rnumber>::do_work(void)
+{
+    TIMEZONE("fluid_solver_check::do_work");
+
+    NSE<rnumber>* nse = new NSE<rnumber>(this->get_communicator(), this->simname);
+    clone_parameters(nse);
+    nse->allocate();
+
+    NSVE<rnumber>* nsve = new NSVE<rnumber>(this->get_communicator(), this->simname);
+    clone_parameters(nsve);
+    nsve->allocate();
+
+    // temporary field
+    field<rnumber, FFTW, THREE>* temp_vec_field = new field<rnumber, FFTW, THREE>(
+            nse->get_nx(), nse->get_ny(), nse->get_nz(),
+            nse->get_communicator(),
+            FFTW_ESTIMATE);
+    make_gaussian_random_field(
+            nse->kk,
+            temp_vec_field,
+            this->random_seed,
+            this->spectrum_dissipation,
+            this->spectrum_Lint,
+            this->spectrum_etaK,
+            this->spectrum_large_scale_const,
+            this->spectrum_small_scale_const,
+            3./2.);
+    nse->kk->template project_divfree<rnumber>(temp_vec_field->get_cdata());
+
+    *(nse->velocity) = temp_vec_field->get_cdata();
+    nse->forcing_type = "fixed_energy_injection_rate";
+    nsve->fs->forcing_type = "fixed_energy_injection_rate";
+    const int return_result_feir = compare_solvers(nse, nsve);
+    *(nse->velocity) = temp_vec_field->get_cdata();
+    nse->forcing_type = "linear";
+    nsve->fs->forcing_type = "linear";
+    const int return_result_linear = compare_solvers(nse, nsve);
+    *(nse->velocity) = temp_vec_field->get_cdata();
+    nse->injection_rate = 0;
+    nsve->fs->injection_rate = 0;
+    nse->famplitude = 0;
+    nsve->fs->famplitude = 0;
+    const int return_result_decay = compare_solvers(nse, nsve);
+
+    delete temp_vec_field;
+    delete nse;
+    delete nsve;
+
+    if ((return_result_linear == EXIT_SUCCESS)
+     && (return_result_feir   == EXIT_SUCCESS)
+     && (return_result_decay  == EXIT_SUCCESS))
+        return EXIT_SUCCESS;
+    else
+        return EXIT_FAILURE;
+}
+
+template class fluid_solver_check<float>;
+template class fluid_solver_check<double>;
+

examples/forcing_unit_test/fluid_solver_check.hpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2024 TurTLE team                                                 *
+*                                                                             *
+*  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 FLUID_SOLVER_CHECK_HPP
+#define FLUID_SOLVER_CHECK_HPP
+
+/** \brief Error analysis for fluid solver.
+ */
+
+#include "full_code/test.hpp"
+
+template <typename rnumber>
+class fluid_solver_check: public test
+{
+    private:
+        double max_velocity_estimate;
+        double spectrum_dissipation;
+        double spectrum_Lint;
+        double spectrum_etaK;
+        double spectrum_large_scale_const;
+        double spectrum_small_scale_const;
+        int output_incompressible_field;
+        int random_seed;
+
+        double dt;
+        double famplitude;
+        double friction_coefficient;
+        double fk0;
+        double fk1;
+        double energy;
+        double injection_rate;
+        double variation_strength;
+        double variation_time_scale;
+        int fmode;
+        std::string forcing_type;
+        double nu;
+
+    public:
+        fluid_solver_check(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simname):
+            test(
+                    COMMUNICATOR,
+                    simname)
+        {}
+        ~fluid_solver_check(){}
+
+        int initialize(void);
+        int do_work(void);
+        int finalize(void);
+
+        template <class dtype>
+        inline int clone_parameters(dtype *dst) {
+            dst->nx                   = this->nx;
+            dst->ny                   = this->ny;
+            dst->nz                   = this->nz;
+            dst->dkx                  = this->dkx;
+            dst->dky                  = this->dky;
+            dst->dkz                  = this->dkz;
+            dst->dealias_type         = this->dealias_type;
+
+            dst->dt                   = this->dt;
+            dst->famplitude           = this->famplitude;
+            dst->friction_coefficient = this->friction_coefficient;
+            dst->fk0                  = this->fk0;
+            dst->fk1                  = this->fk1;
+            dst->energy               = this->energy;
+            dst->injection_rate       = this->injection_rate;
+            dst->fmode                = this->fmode;
+            dst->forcing_type         = this->forcing_type;
+            dst->nu                   = this->nu;
+            return EXIT_SUCCESS;
+        }
+
+};
+
+#endif//FLUID_SOLVER_CHECK_HPP
+

examples/forcing_unit_test/sanity_check.py

+import os
+import sys
+import getpass
+import numpy as np
+import h5py
+
+import TurTLE
+import TurTLE.PP
+import TurTLE.DNS
+import TurTLE.TEST
+
+username = getpass.getuser()
+homefolder = os.path.expanduser('~')
+cpp_location = None
+
+import pathlib
+cpp_location = pathlib.Path(__file__).parent.resolve()
+
+
+class sanity_check(TurTLE.TEST):
+    def __init__(
+            self,
+            **kwargs):
+        TurTLE.PP.__init__(self, **kwargs)
+        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 = [
+                '"full_code/main_code.hpp"',
+                '<cmath>']
+        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 = 'fluid_solver_check'
+        self.dns_type = 'fluid_solver_check'
+        self.definitions = ''
+        self.definitions += open(
+            os.path.join(
+                cpp_location, 'fluid_solver_check.hpp'), 'r').read()
+        self.definitions += open(
+            os.path.join(
+                cpp_location, 'fluid_solver_check.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')
+        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_ac = subparsers.add_parser(
+                'fluid_solver_check')
+        self.simulation_parser_arguments(parser_ac)
+        self.job_parser_arguments(parser_ac)
+        self.parameters_to_parser_arguments(parser_ac)
+        self.parameters_to_parser_arguments(
+            parser_ac,
+            parameters = self.extra_postprocessing_parameters(
+                'fluid_solver_check'))
+        return None
+    def extra_postprocessing_parameters(
+            self,
+            dns_type = None):
+        assert(dns_type == 'fluid_solver_check')
+        pars = {}
+        return pars
+
+import matplotlib.pyplot as plt
+
+def main():
+    bla = sanity_check()
+    bla.launch(
+            ['fluid_solver_check',
+             '--np', '1',
+             '--ntpp', '1',
+             '--simname', 'fluid_solver_onestep_sanity_check']
+            + sys.argv[1:])
+    return None
+
+if __name__ == '__main__':
+    main()