testing gaussian field

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

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;
                     }

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);