From e0bf5553ad6ccca8602e2d61d6cb880cf8e884cc Mon Sep 17 00:00:00 2001 From: Lukas Bentkamp <bentkamp@lashgarak.lmp.ds.mpg.de> Date: Wed, 30 Oct 2019 17:03:34 +0100 Subject: [PATCH] some fixes --- cpp/full_code/Gauss_field_test.cpp | 12 +++---- cpp/full_code/Gauss_field_test.hpp | 2 +- cpp/full_code/kraichnan_field.cpp | 58 +++++++++++++++--------------- cpp/full_code/kraichnan_field.hpp | 3 ++ 4 files changed, 38 insertions(+), 37 deletions(-) diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp index 5005000b..3c8405e1 100644 --- a/cpp/full_code/Gauss_field_test.cpp +++ b/cpp/full_code/Gauss_field_test.cpp @@ -76,16 +76,16 @@ int make_gaussian_random_field( case ONE: { double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]); - output_field->cval(cindex,0) = coefficient * cos(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff); - output_field->cval(cindex,1) = coefficient * sin(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff); + output_field->cval(cindex,0) = coefficient * cos(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); + output_field->cval(cindex,1) = coefficient * sin(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); break; } case THREE: for (int cc = 0; cc<3; cc++) { double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]); - output_field->cval(cindex,cc,0) = coefficient * cos(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff); - output_field->cval(cindex,cc,1) = coefficient * sin(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff); + output_field->cval(cindex,cc,0) = coefficient * cos(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); + output_field->cval(cindex,cc,1) = coefficient * sin(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); } break; case THREExTHREE: @@ -93,8 +93,8 @@ int make_gaussian_random_field( for (int ccc = 0; ccc<3; ccc++) { double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]); - output_field->cval(cindex,cc,ccc,0) = coefficient * cos(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff); - output_field->cval(cindex,cc,ccc,1) = coefficient * sin(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff); + output_field->cval(cindex,cc,ccc,0) = coefficient * cos(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); + output_field->cval(cindex,cc,ccc,1) = coefficient * sin(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); } break; } diff --git a/cpp/full_code/Gauss_field_test.hpp b/cpp/full_code/Gauss_field_test.hpp index d14c3cfb..185351cc 100644 --- a/cpp/full_code/Gauss_field_test.hpp +++ b/cpp/full_code/Gauss_field_test.hpp @@ -50,7 +50,7 @@ int make_gaussian_random_field( field<rnumber, be, fc> *output_field, const int rseed = 0, const double slope = -5./3., - const double k_cutoff = 10., + const double k2_cutoff = 10., const double coefficient = 1.); template <typename rnumber> diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp index a0f81f52..2884c9b7 100644 --- a/cpp/full_code/kraichnan_field.cpp +++ b/cpp/full_code/kraichnan_field.cpp @@ -62,37 +62,15 @@ int kraichnan_field<rnumber>::initialize(void) return EXIT_FAILURE; } - /* - this->vorticity = new field<rnumber, FFTW, THREE>( - nx, ny, nz, - this->comm, - fftw_planner_string_to_flag[this->fftw_plan_rigor]); - this->vorticity->real_space_representation = false; -*/ this->velocity = new field<rnumber, FFTW, THREE>( nx, ny, nz, this->comm, fftw_planner_string_to_flag[this->fftw_plan_rigor]); this->velocity->real_space_representation = false; - //read vorticity field - /*this->vorticity->io( - this->get_current_fname(), - "vorticity", - this->iteration, - true);*/ + this->kk = new kspace<FFTW, SMOOTH>( this->velocity->clayout, this->dkx, this->dky, this->dkz); // IS THE VELOCITY LAYOUT THE SAME AS THE VORTICITY ONE? -/* - // compute the velocity field and store - invert_curl( - this->kk, - this->vorticity, - velocity); - */ - // transform velocity and vorticity fields to real space - //this->velocity->ift(); - //this->vorticity->ift(); // initialize particles this->ps = particles_system_builder( @@ -124,6 +102,19 @@ int kraichnan_field<rnumber>::initialize(void) "position/0"); this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout()); + this->slope = -5./3.; + this->k_cutoff = 10.; //k_max is approximately N/2 in a N^3 simulation + const double energy = 1.; + assert(this->slope > -3); + + // Here we rescale the field to match the desired energy. The divisor + // is the integral of the power law spectrum with exponential cutoff. + this->field_coefficient = sqrt(energy / + (4*M_PI*pow(this->k_cutoff, 3+this->slope)*tgamma(3+this->slope))); + + DEBUG_MSG("Coefficient: %g\n", + this->field_coefficient); + return EXIT_SUCCESS; } @@ -132,19 +123,24 @@ int kraichnan_field<rnumber>::step(void) { TIMEZONE("kraichnan_file::step"); - make_gaussian_random_field(this->kk, this->velocity, this->iteration, -5./3., 10., 0.01); + make_gaussian_random_field( + this->kk, + this->velocity, + this->iteration, + this->slope, + this->k_cutoff, + this->field_coefficient); this->velocity->ift(); - //this->kk->template project_divfree<rnumber>(this->velocity->get_cdata(), true); + this->kk->template project_divfree<rnumber>(this->velocity->get_cdata(), true); + DEBUG_MSG("L2Norm: %g\n", + this->velocity->L2norm(this->kk)); // What does the ->template do? - // What is cdata and rdata? // ADJUST TOTAL ENERGY, SLOPE AND CUTOFF DEBUG_MSG("Iteration: %d\n", this->iteration); - DEBUG_MSG("Real space: %d\n", - this->velocity->real_space_representation); - DEBUG_MSG("Field entries: %g\n", + DEBUG_MSG("Field entries: %g\n", this->velocity->rval(0,0)); DEBUG_MSG("Field entries: %g\n", this->velocity->rval(0,1)); @@ -152,7 +148,9 @@ int kraichnan_field<rnumber>::step(void) this->velocity->rval(0,2)); DEBUG_MSG("Field entries: %g\n", this->velocity->rval(15,2)); - this->ps->completeLoop(sqrt(this->dt)); + // In principle this should be sqrt(dt), but to match the NS energy, + // we absorb this into the velocity field. + this->ps->completeLoop(this->dt); this->iteration++; return EXIT_SUCCESS; } diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp index 93b044c5..a6be715f 100644 --- a/cpp/full_code/kraichnan_field.hpp +++ b/cpp/full_code/kraichnan_field.hpp @@ -59,6 +59,9 @@ class kraichnan_field: public direct_numerical_simulation /* other stuff */ kspace<FFTW, SMOOTH> *kk; field<rnumber, FFTW, THREE> *velocity; + double slope; + double k_cutoff; + double field_coefficient; /* other stuff */ std::unique_ptr<abstract_particles_system<long long int, double>> ps; -- GitLab