Commit 2952398f authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/ksqu_energy_spectrum' into develop

parents d4f77830 70d2f4ec
Pipeline #68269 passed with stage
in 6 minutes and 49 seconds
......@@ -188,6 +188,7 @@ class TEST(_code):
scal_rspace_stats = []
if self.dns_type in ['Gauss_field_test']:
vec_spectra_stats.append('velocity')
vec_spectra_stats.append('k*velocity')
vec4_rspace_stats.append('velocity')
tens_rspace_stats.append('velocity_gradient')
scal_rspace_stats.append('velocity_divergence')
......
......@@ -5,7 +5,7 @@ from scipy import trapz
from scipy.stats import norm
from scipy.integrate import quad
import h5py
import sys
import sys, os
import time
import TurTLE
......@@ -17,24 +17,28 @@ except:
plt = None
def main():
c = TEST()
# size of grid
n = 256
n = 1024
slope = -5./3.
k_cutoff = 30.
k_cutoff = 64.
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]
total_energy = quad(func, 0.6, k_cutoff*8)[0]
coeff = 1./total_energy
bin_no = 100
rseed = int(time.time())
simname = 'Gaussianity_test'
c.launch(
if not os.path.exists(simname+'.h5'):
c = TEST()
opt = c.launch(
['Gauss_field_test',
'--nx', str(n),
'--ny', str(n),
'--nz', str(n),
'--simname', 'Gaussianity_test',
'--simname', simname,
'--np', '4',
'--environment', 'short',
'--minutes', '60',
'--ntpp', '1',
'--wd', './',
'--histogram_bins', str(bin_no),
......@@ -44,7 +48,7 @@ def main():
'--spectrum_coefficient', str(coeff),
'--field_random_seed', str(rseed)] +
sys.argv[1:])
plot_stuff(c.simname, total_energy = total_energy)
plot_stuff(simname, total_energy = total_energy)
return None
def plot_stuff(simname, total_energy = 1.):
......@@ -82,6 +86,7 @@ def plot_stuff(simname, total_energy = 1.):
f_vel = hist_vel / np.sum(hist_vel, axis=0, keepdims=True).astype(float) / velbinsize
print('Energy analytically: {}'.format(total_energy))
print(np.sum(energy*np.arange(len(energy))**2))
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)))
......@@ -101,6 +106,17 @@ def plot_stuff(simname, total_energy = 1.):
df['statistics/moments/velocity_divergence'][0, 2]))
print('Gradient second moment is: {0}'.format(
df['statistics/moments/velocity_gradient'][0, 2].mean()))
print('----------- k2-premultiplied spectrum -----------')
k2func = lambda k, k_c=k_cutoff, s=slope : k**(2+s)*np.exp(-k/k_c)
k2sum_analytic = quad(k2func, 0, k_cutoff*20)[0]
print('Analytically: {}'.format(k2sum_analytic))
k2spec_trace = (
df['statistics/spectra/k*velocity_k*velocity'][..., 0, 0]
+ df['statistics/spectra/k*velocity_k*velocity'][..., 1, 1]
+ df['statistics/spectra/k*velocity_k*velocity'][..., 2, 2])
print('Energy sum: {}'.format(np.sum(k2spec_trace*df['kspace/dk'][()])/2./coeff))
df.close()
return None
......
......@@ -144,6 +144,14 @@ int Gauss_field_test<rnumber>::do_work(void)
this->vector_field,
this->scalar_field);
/// compute integral of premultiplied energy spectrum
this->kk->cospectrum<rnumber, THREE>(
this->vector_field->get_cdata(),
stat_group,
"k*velocity_k*velocity",
0,
2.);
/// compute basic statistics of Gaussian field
this->vector_field->compute_stats(
this->kk,
......
......@@ -684,7 +684,8 @@ void kspace<be, dt>::cospectrum(
const rnumber(* __restrict b)[2],
const hid_t group,
const std::string dset_name,
const hsize_t toffset)
const hsize_t toffset,
const double wavenumber_exp)
{
TIMEZONE("field::cospectrum2");
shared_array<double> spec_local_thread(this->nshells*ncomp(fc)*ncomp(fc),[&](double* spec_local){
......@@ -704,7 +705,8 @@ void kspace<be, dt>::cospectrum(
int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
for (hsize_t i=0; i<ncomp(fc); i++)
for (hsize_t j=0; j<ncomp(fc); j++){
spec_local[tmp_int + i*ncomp(fc)+j] += nxmodes * (
spec_local[tmp_int + i*ncomp(fc)+j] += pow(k2,wavenumber_exp/2.)*
nxmodes * (
(a[ncomp(fc)*cindex + i][0] * b[ncomp(fc)*cindex + j][0]) +
(a[ncomp(fc)*cindex + i][1] * b[ncomp(fc)*cindex + j][1]));
}
......@@ -762,7 +764,8 @@ void kspace<be, dt>::cospectrum(
const rnumber(* __restrict a)[2],
const hid_t group,
const std::string dset_name,
const hsize_t toffset)
const hsize_t toffset,
const double wavenumber_exp)
{
TIMEZONE("field::cospectrum1");
shared_array<double> spec_local_thread(this->nshells*ncomp(fc)*ncomp(fc),[&](double* spec_local){
......@@ -780,9 +783,11 @@ void kspace<be, dt>::cospectrum(
{
double* spec_local = spec_local_thread.getMine();
int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
for (hsize_t i=0; i<ncomp(fc); i++)
for (hsize_t j=0; j<ncomp(fc); j++){
spec_local[tmp_int + i*ncomp(fc)+j] += nxmodes * (
spec_local[tmp_int + i*ncomp(fc)+j] += pow(k2, wavenumber_exp/2.)*
nxmodes * (
(a[ncomp(fc)*cindex + i][0] * a[ncomp(fc)*cindex + j][0]) +
(a[ncomp(fc)*cindex + i][1] * a[ncomp(fc)*cindex + j][1]));
}
......@@ -1007,136 +1012,160 @@ template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREExTHREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, ONE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREExTHREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, ONE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREExTHREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, ONE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREExTHREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, ONE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp);
template double kspace<FFTW, TWO_THIRDS>::L2norm<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ a);
......
......@@ -129,7 +129,8 @@ class kspace
const rnumber(* __restrict__ b)[2],
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp = 0);
template <typename rnumber,
field_components fc>
......@@ -137,7 +138,8 @@ class kspace
const rnumber(* __restrict__ a)[2],
const hid_t group,
const std::string dset_name,
const hsize_t toffset);
const hsize_t toffset,
const double wavenumber_exp = 0);
template <typename rnumber,
field_components fc>
......
Markdown is supported
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