Commit 11e0f4d4 authored by Lukas Bentkamp's avatar Lukas Bentkamp
Browse files

Energy considerations

parent a2309ff0
......@@ -742,7 +742,10 @@ class DNS(_code):
opt.nz > opt.n):
opt.n = min(opt.nx, opt.ny, opt.nz)
print("Warning: '-n' parameter changed to minimum of nx, ny, nz. This affects the computation of nu.")
self.parameters['dt'] = (opt.dtfactor / opt.n)
if self.dns_type in ['kraichnan_field']:
self.parameters['dt'] = opt.dtfactor * 0.5 / opt.n**2
else:
self.parameters['dt'] = (opt.dtfactor / opt.n)
self.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
# check value of kMax
kM = opt.n * 0.5
......
......@@ -3,6 +3,7 @@
import numpy as np
from scipy import trapz
from scipy.stats import norm
from scipy.integrate import quad
import h5py
import sys
import time
......@@ -18,11 +19,13 @@ except:
def main():
c = TEST()
# size of grid
n = 128
n = 256
slope = -5./3.
k_cutoff = 16.
coeff = 1.
bin_no = 129
k_cutoff = 30.
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
bin_no = 100
rseed = int(time.time())
c.launch(
......@@ -35,25 +38,29 @@ def main():
'--ntpp', '1',
'--wd', './',
'--histogram_bins', str(bin_no),
'--max_velocity_estimate', '2.',
'--max_velocity_estimate', '8.',
'--spectrum_slope', str(slope),
'--spectrum_k_cutoff', str(k_cutoff),
'--spectrum_coefficient', str(coeff),
'--field_random_seed', str(rseed)] +
sys.argv[1:])
plot_stuff(c.simname)
df = h5py.File(c.simname + '.h5', 'r')
def plot_stuff(simname):
df = h5py.File(simname + '.h5', 'r')
for kk in ['spectrum_slope',
'spectrum_k_cutoff',
'spectrum_coefficient',
'field_random_seed',
'histogram_bins']:
print(kk, df['parameters/' + kk][...])
coeff = df['parameters/spectrum_coefficient'][()]
k_cutoff = df['parameters/spectrum_k_cutoff'][()]
slope = df['parameters/spectrum_slope'][()]
bin_no = df['parameters/histogram_bins'][()]
f = plt.figure()
# EVERYTHING SHOULD BE READ FROM FILE BECAUSE WE CANT BE SURE THE CODE RAN AGAIN
# test spectrum
a = f.add_subplot(121)
kk = df['kspace/kshell'][...]
......@@ -74,8 +81,8 @@ def main():
hist_vel = df['statistics/histograms/velocity'][0, :, :3]
f_vel = hist_vel / np.sum(hist_vel, axis=0, keepdims=True).astype(float) / velbinsize
print('Energy integral: {}'.format(trapz(energy[1:-2], kk[1:-2])))
print('Energy sum: {}'.format(np.sum(energy[1:-2])))
print('Energy analytically: {}'.format(total_energy))
print('Energy sum: {}'.format(np.sum(energy*df['kspace/dk'][()])))
print('Moment sum: {}'.format(df['statistics/moments/velocity'][0,2,3]/2))
print('Velocity variances: {}'.format(trapz(vel[:,None]**2*f_vel, vel[:,None], axis=0)))
......
......@@ -73,27 +73,26 @@ int make_gaussian_random_field(
{
if (k2 > 0)
switch(fc)
{
{
case ONE:
{
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.);
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++)
{ // In the slope, a factor of 1/2 compensates the k^2, a factor of 1/2 goes from the energy to the Fourier coefficient
// and the 2*pi*k^2 divides out the surface of the shell
output_field->cval(cindex,cc,0) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, (slope - 2)/4.) * 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 - 2)/4.) * exp(- sqrt(k2)/k_cutoff/2.);
} // TODO: ADJUST OTHER FIELD TYPES
{
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) = 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.);
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;
}
......@@ -174,7 +173,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(), true);
this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
/// initialize statistics file.
hid_t stat_group, stat_file;
......
......@@ -102,20 +102,14 @@ int kraichnan_field<rnumber>::initialize(void)
"position/0");
this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
this->slope = -5./3.;
this->k_cutoff = 20.; //k_max is approximately N/2 in a N^3 simulation
const double energy = 1.;
this->spectrum_slope = -5./3.;
this->spectrum_k_cutoff = 20.; //k_max is approximately N/2 in a N^3 simulation
// Here we rescale the field to match the desired energy. The divisor
// is the integral of the power law spectrum with exponential cutoff.
// 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
// Projection factor
this->spectrum_coeff = 3./2.;
DEBUG_MSG("Coefficient: %g\n",
this->field_coefficient);
this->spectrum_coeff);
return EXIT_SUCCESS;
}
......@@ -129,32 +123,19 @@ int kraichnan_field<rnumber>::step(void)
this->kk,
this->velocity,
this->iteration,
this->slope,
this->k_cutoff,
this->field_coefficient);
this->spectrum_slope,
this->spectrum_k_cutoff,
this->spectrum_coeff * 3./2.); // incompressibility projection factor
this->velocity->ift();
// PUT THE INCOMPRESSIBILITY PROJECTION FACTOR sqrt(2/3)
// 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: %g\n",
DEBUG_MSG("L2Norm after: %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);
DEBUG_MSG("Field entries: %g\n",
this->velocity->rval(0,0));
DEBUG_MSG("Field entries: %g\n",
this->velocity->rval(0,1));
DEBUG_MSG("Field entries: %g\n",
this->velocity->rval(0,2));
DEBUG_MSG("Field entries: %g\n",
this->velocity->rval(15,2));
// 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->ps->completeLoop(sqrt(this->dt));
this->iteration++;
return EXIT_SUCCESS;
}
......
......@@ -59,9 +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;
double spectrum_slope;
double spectrum_k_cutoff;
double spectrum_coeff;
/* other stuff */
std::unique_ptr<abstract_particles_system<long long int, double>> ps;
......
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