some fixes

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

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>

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

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;