From 5caa1e6a50d0f1613d5ced8ba5f9536b54dd37b3 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Thu, 20 Jun 2024 08:33:00 +0200
Subject: [PATCH] [wip] adds working folder for unit test of forcing
---
examples/forcing_unit_test/README.rst | 40 +++
.../forcing_unit_test/fluid_solver_check.cpp | 311 ++++++++++++++++++
.../forcing_unit_test/fluid_solver_check.hpp | 100 ++++++
examples/forcing_unit_test/sanity_check.py | 100 ++++++
4 files changed, 551 insertions(+)
create mode 100644 examples/forcing_unit_test/README.rst
create mode 100644 examples/forcing_unit_test/fluid_solver_check.cpp
create mode 100644 examples/forcing_unit_test/fluid_solver_check.hpp
create mode 100644 examples/forcing_unit_test/sanity_check.py
diff --git a/examples/forcing_unit_test/README.rst b/examples/forcing_unit_test/README.rst
new file mode 100644
index 00000000..0e85aaf6
--- /dev/null
+++ b/examples/forcing_unit_test/README.rst
@@ -0,0 +1,40 @@
+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)
diff --git a/examples/forcing_unit_test/fluid_solver_check.cpp b/examples/forcing_unit_test/fluid_solver_check.cpp
new file mode 100644
index 00000000..6993d920
--- /dev/null
+++ b/examples/forcing_unit_test/fluid_solver_check.cpp
@@ -0,0 +1,311 @@
+/******************************************************************************
+* *
+* 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>;
+
diff --git a/examples/forcing_unit_test/fluid_solver_check.hpp b/examples/forcing_unit_test/fluid_solver_check.hpp
new file mode 100644
index 00000000..dc58515a
--- /dev/null
+++ b/examples/forcing_unit_test/fluid_solver_check.hpp
@@ -0,0 +1,100 @@
+/******************************************************************************
+* *
+* 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
+
diff --git a/examples/forcing_unit_test/sanity_check.py b/examples/forcing_unit_test/sanity_check.py
new file mode 100644
index 00000000..4f3cc120
--- /dev/null
+++ b/examples/forcing_unit_test/sanity_check.py
@@ -0,0 +1,100 @@
+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()
--
GitLab