diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py index c2d33afbd3853ba2760432bd1a47c0718ca29c07..f28abf74bf4557815f4805a8b769de9adb52ff7a 100644 --- a/TurTLE/test/test_Gaussian_field.py +++ b/TurTLE/test/test_Gaussian_field.py @@ -17,10 +17,9 @@ def main(): c = TEST() # size of grid n = 64 - - slope = 5./3. - k_cutoff = 10. - coeff = np.sqrt(1./(4*np.pi*k_cutoff**(3+slope)*gamma(3+slope))) + slope = 1. + k_cutoff = 20. + coeff = np.sqrt(1./0.93058) bin_no = 33 c.launch( @@ -49,7 +48,9 @@ def main(): kk = df['kspace/kshell'][...] phi_ij = df['statistics/spectra/velocity_velocity'][0] energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2 - a.plot(kk, energy) # plot analytic (read from parameter file) + + a.plot(kk[1:-2], energy[1:-2]) + a.plot(kk[1:-2], coeff*kk[1:-2]**slope*np.exp(-kk[1:-2]/k_cutoff), ls='--', c='C0') a.set_xscale('log') a.set_yscale('log') # test isotropy diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp index ae0a1df784671192c1b537968ea14cc35dfa3676..7e26954d41034a94a01e49cd1da7eebccf9f25fb 100644 --- a/cpp/full_code/Gauss_field_test.cpp +++ b/cpp/full_code/Gauss_field_test.cpp @@ -47,7 +47,7 @@ int make_gaussian_random_field( TIMEZONE("make_gaussian_random_field"); // initialize a separate random number generator for each thread std::vector<std::mt19937_64> rgen; - std::uniform_real_distribution<rnumber> rdist; + std::normal_distribution<rnumber> rdist; rgen.resize(omp_get_max_threads()); // seed random number generators such that no seed is ever repeated if we change the value of rseed. // basically use a multi-dimensional array indexing technique to come up with actual seed. @@ -62,6 +62,7 @@ int make_gaussian_random_field( } output_field->real_space_representation = false; *output_field = 0.0; + DEBUG_MSG("slope: %g\n", slope); // inside loop use only thread-local random number generator kk->CLOOP_K2([&]( ptrdiff_t cindex, @@ -75,26 +76,23 @@ 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/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.); + output_field->cval(cindex,0) = coefficient * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); + output_field->cval(cindex,1) = coefficient * rdist(rgen[omp_get_thread_num()]) * 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/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.); - } + output_field->cval(cindex,cc,0) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4. - 1) * exp(- sqrt(k2)/k_cutoff/2.); + output_field->cval(cindex,cc,1) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4. - 1) * exp(- sqrt(k2)/k_cutoff/2.); + } // TODO: ADJUST OTHER FIELD TYPES break; case THREExTHREE: for (int cc = 0; cc<3; cc++) 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/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.); + output_field->cval(cindex,cc,ccc,0) = coefficient * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); + output_field->cval(cindex,cc,ccc,1) = coefficient * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.); } break; } diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp index 2884c9b7fdb761b6d01e938caa324a21d60de8a5..61192d7ed35b7f6ece154b019cd3814178e9374e 100644 --- a/cpp/full_code/kraichnan_field.cpp +++ b/cpp/full_code/kraichnan_field.cpp @@ -103,14 +103,16 @@ int kraichnan_field<rnumber>::initialize(void) 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 + this->k_cutoff = 20.; //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))); + // The prefactor sqrt(3/2) compensates the projection to incompressibity. + this->field_coefficient = sqrt(3./2.* energy / + (1.10369)); + // 0.93058 for k_cutoff = 10 + // 1.10369 for k_cutoff = 20 DEBUG_MSG("Coefficient: %g\n", this->field_coefficient); @@ -132,11 +134,13 @@ int kraichnan_field<rnumber>::step(void) this->field_coefficient); this->velocity->ift(); - this->kk->template project_divfree<rnumber>(this->velocity->get_cdata(), true); + // PUT THE INCOMPRESSIBILITY PROJECTION FACTOR sqrt(2/3) + this->kk->template project_divfree<rnumber>(this->velocity->get_cdata()); DEBUG_MSG("L2Norm: %g\n", this->velocity->L2norm(this->kk)); // What does the ->template do? // ADJUST TOTAL ENERGY, SLOPE AND CUTOFF + // we actually want to propagate with sqrt(dt) instead of dt DEBUG_MSG("Iteration: %d\n", this->iteration);