Commit 481916bc authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

moves `make_gaussian_random_field` to field.?pp

parent 11e0f4d4
......@@ -22,9 +22,9 @@ def main():
n = 256
slope = -5./3.
k_cutoff = 30.
func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
total_energy = quad(func, 1, k_cutoff*4)[0]
coeff = 3./2./total_energy
coeff = 1./total_energy
bin_no = 100
rseed = int(time.time())
......@@ -44,10 +44,10 @@ def main():
'--spectrum_coefficient', str(coeff),
'--field_random_seed', str(rseed)] +
sys.argv[1:])
plot_stuff(c.simname, total_energy = total_energy)
return None
plot_stuff(c.simname)
def plot_stuff(simname):
def plot_stuff(simname, total_energy = 1.):
df = h5py.File(simname + '.h5', 'r')
for kk in ['spectrum_slope',
'spectrum_k_cutoff',
......@@ -59,7 +59,7 @@ def plot_stuff(simname):
k_cutoff = df['parameters/spectrum_k_cutoff'][()]
slope = df['parameters/spectrum_slope'][()]
bin_no = df['parameters/histogram_bins'][()]
f = plt.figure()
# test spectrum
a = f.add_subplot(121)
......@@ -67,7 +67,7 @@ def plot_stuff(simname):
print('dk: {}'.format(df['kspace/dk'][()]))
phi_ij = df['statistics/spectra/velocity_velocity'][0] / df['kspace/dk'][()]
energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
a.scatter(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')
......
......@@ -29,6 +29,7 @@
#include <cstdlib>
#include <algorithm>
#include <cassert>
#include <random>
#include "field.hpp"
#include "scope_timer.hpp"
#include "shared_array.hpp"
......@@ -2032,6 +2033,76 @@ field<rnumber, be, fc> &field<rnumber, be, fc>::operator=(
return *this;
}
template <typename rnumber,
field_backend be,
field_components fc,
kspace_dealias_type dt>
int make_gaussian_random_field(
kspace<be, dt> *kk,
field<rnumber, be, fc> *output_field,
const int rseed,
const double slope,
const double k_cutoff,
const double coefficient)
{
TIMEZONE("make_gaussian_random_field");
// initialize a separate random number generator for each thread
std::vector<std::mt19937_64> rgen;
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.
for (int thread_id=0; thread_id < omp_get_max_threads(); thread_id++)
{
int current_seed = (
rseed*omp_get_max_threads()*output_field->clayout->nprocs +
output_field->clayout->myrank*omp_get_max_threads() +
thread_id);
DEBUG_MSG("in make_gaussian_random_field, thread_id = %d, current_seed = %d\n", thread_id, current_seed);
rgen[thread_id].seed(current_seed);
}
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,
ptrdiff_t xindex,
ptrdiff_t yindex,
ptrdiff_t zindex,
double k2)
{
if (k2 > 0)
switch(fc)
{
case ONE:
{
output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
break;
}
case THREE:
for (int cc = 0; cc<3; cc++)
{
output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
}
break;
case THREExTHREE:
for (int cc = 0; cc<3; cc++)
for (int ccc = 0; ccc<3; ccc++)
{
output_field->cval(cindex,cc,ccc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
}
break;
}
});
output_field->symmetrize();
return EXIT_SUCCESS;
}
template class field<float, FFTW, ONE>;
template class field<float, FFTW, THREE>;
template class field<float, FFTW, THREExTHREE>;
......@@ -2189,3 +2260,49 @@ template int joint_rspace_3PDF<double, FFTW>(
const std::vector<double>,
const std::vector<double>);
template int make_gaussian_random_field<
float,
FFTW,
ONE,
SMOOTH>(
kspace<FFTW, SMOOTH> *kk,
field<float, FFTW, ONE> *output_field,
const int rseed,
const double slope,
const double k_cutoff,
const double coefficient);
template int make_gaussian_random_field<
double,
FFTW,
ONE,
SMOOTH>(
kspace<FFTW, SMOOTH> *kk,
field<double, FFTW, ONE> *output_field,
const int rseed,
const double slope,
const double k_cutoff,
const double coefficient);
template int make_gaussian_random_field<
float,
FFTW,
THREE,
SMOOTH>(
kspace<FFTW, SMOOTH> *kk,
field<float, FFTW, THREE> *output_field,
const int rseed,
const double slope,
const double k_cutoff,
const double coefficient);
template int make_gaussian_random_field<
double,
FFTW,
THREE,
SMOOTH>(
kspace<FFTW, SMOOTH> *kk,
field<double, FFTW, THREE> *output_field,
const int rseed,
const double slope,
const double k_cutoff,
const double coefficient);
......@@ -374,5 +374,20 @@ int joint_rspace_3PDF(
const std::vector<double> max_f2_estimate,
const std::vector<double> max_f3_estimate);
/** \brief Generate a Gaussian random field
*
* */
template <typename rnumber,
field_backend be,
field_components fc,
kspace_dealias_type dt>
int make_gaussian_random_field(
kspace<be, dt> *kk,
field<rnumber, be, fc> *output_field,
const int rseed = 0,
const double slope = -5./3.,
const double k2_cutoff = 10.,
const double coefficient = 1.);
#endif//FIELD_HPP
......@@ -32,75 +32,7 @@
template <typename rnumber,
field_backend be,
field_components fc,
kspace_dealias_type dt>
int make_gaussian_random_field(
kspace<be, dt> *kk,
field<rnumber, be, fc> *output_field,
const int rseed,
const double slope,
const double k_cutoff,
const double coefficient)
{
TIMEZONE("make_gaussian_random_field");
// initialize a separate random number generator for each thread
std::vector<std::mt19937_64> rgen;
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.
for (int thread_id=0; thread_id < omp_get_max_threads(); thread_id++)
{
int current_seed = (
rseed*omp_get_max_threads()*output_field->clayout->nprocs +
output_field->clayout->myrank*omp_get_max_threads() +
thread_id);
DEBUG_MSG("in make_gaussian_random_field, thread_id = %d, current_seed = %d\n", thread_id, current_seed);
rgen[thread_id].seed(current_seed);
}
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,
ptrdiff_t xindex,
ptrdiff_t yindex,
ptrdiff_t zindex,
double k2)
{
if (k2 > 0)
switch(fc)
{
case ONE:
{
output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
break;
}
case THREE:
for (int cc = 0; cc<3; cc++)
{
output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
}
break;
case THREExTHREE:
for (int cc = 0; cc<3; cc++)
for (int ccc = 0; ccc<3; ccc++)
{
output_field->cval(cindex,cc,ccc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
}
break;
}
});
output_field->symmetrize();
return EXIT_SUCCESS;
}
template <typename rnumber>
int Gauss_field_test<rnumber>::initialize(void)
......@@ -173,7 +105,7 @@ int Gauss_field_test<rnumber>::do_work(void)
this->spectrum_coefficient);
/// impose divergence free condition while maintaining the energy of the field
this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata(), true);
/// initialize statistics file.
hid_t stat_group, stat_file;
......
......@@ -39,20 +39,6 @@
*
*/
template <typename rnumber,
field_backend be,
field_components fc,
kspace_dealias_type dt>
int make_gaussian_random_field(
kspace<be, dt> *kk,
field<rnumber, be, fc> *output_field,
const int rseed = 0,
const double slope = -5./3.,
const double k2_cutoff = 10.,
const double coefficient = 1.);
template <typename rnumber>
class Gauss_field_test: public test
{
......
......@@ -28,7 +28,6 @@
#include "kraichnan_field.hpp"
#include "scope_timer.hpp"
#include "fftw_tools.hpp"
#include "Gauss_field_test.hpp"
template <typename rnumber>
......@@ -65,12 +64,12 @@ int kraichnan_field<rnumber>::initialize(void)
this->velocity = new field<rnumber, FFTW, THREE>(
nx, ny, nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->velocity->real_space_representation = false;
// the velocity layout should be the same as the vorticity one
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?
// initialize particles
this->ps = particles_system_builder(
......@@ -128,13 +127,13 @@ int kraichnan_field<rnumber>::step(void)
this->spectrum_coeff * 3./2.); // incompressibility projection factor
this->velocity->ift();
// PUT THE INCOMPRESSIBILITY PROJECTION FACTOR sqrt(2/3)
DEBUG_MSG("L2Norm before: %g\n",
// PUT THE INCOMPRESSIBILITY PROJECTION FACTOR sqrt(2/3)
DEBUG_MSG("L2Norm before: %g\n",
this->velocity->L2norm(this->kk));
this->kk->template project_divfree<rnumber>(this->velocity->get_cdata());
DEBUG_MSG("L2Norm after: %g\n",
DEBUG_MSG("L2Norm after: %g\n",
this->velocity->L2norm(this->kk));
this->ps->completeLoop(sqrt(this->dt));
this->iteration++;
return EXIT_SUCCESS;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment