Commit ef3d23f4 authored by Chichi Lalescu's avatar Chichi Lalescu
Browse files

add back autocorrelation functionality

parent d319b5e6
......@@ -46,6 +46,26 @@ void fluid_solver_base<rnumber>::clean_up_real_space(rnumber *a, int howmany)
std::fill_n(a+rindex+this->rd->subsizes[2]*howmany, 2*howmany, 0.0);
}
template <class rnumber>
rnumber fluid_solver_base<rnumber>::autocorrel(cnumber *a)
{
double *spec = fftw_alloc_real(this->nshells*9);
double sum_local, sum;
this->cospectrum(a, a, spec);
sum_local = 0.0;
for (int n = 0; n < this->nshells; n++)
sum_local += spec[n*9] + spec[n*9 + 4] + spec[n*9 + 8];
fftw_free(spec);
MPI_Allreduce(
&sum_local,
&sum,
1,
MPI_DOUBLE,
MPI_SUM,
this->cd->comm);
return sum;
}
template <class rnumber>
void fluid_solver_base<rnumber>::cospectrum(cnumber *a, cnumber *b, double *spec)
{
......
......@@ -83,6 +83,7 @@ class fluid_solver_base
void clean_up_real_space(rnumber *a, int howmany);
void cospectrum(cnumber *a, cnumber *b, double *spec);
void cospectrum(cnumber *a, cnumber *b, double *spec, const double k2exponent);
rnumber autocorrel(cnumber *a);
void compute_rspace_stats(rnumber *a,
double *moments,
ptrdiff_t *hist,
......
......@@ -126,7 +126,7 @@ class fluid_particle_base(bfps.code):
'local_time_difference = ((unsigned int)(time1 - time0))/((double)CLOCKS_PER_SEC);\n' +
'time_difference = 0.0;\n' +
'MPI_Allreduce(&local_time_difference, &time_difference, ' +
'1, MPI_DOUBLE, MPI_SUM, fs->rd->comm);\n' +
'1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);\n' +
'if (myrank == 0) std::cout << "iteration " ' +
'<< iteration << " took " ' +
'<< time_difference/nprocs << " seconds" << std::endl;\n' +
......
......@@ -70,15 +70,15 @@ class fluid_resize(bfps.fluid_base.fluid_particle_base):
fs1->iteration = 0;
fs0->read('v', 'c');
double a, b;
a = 0.5*fs0->correl_vec(fs0->cvelocity, fs0->cvelocity);
b = 0.5*fs0->correl_vec(fs0->cvorticity, fs0->cvorticity);
a = 0.5*fs0->autocorrel(fs0->cvelocity);
b = 0.5*fs0->autocorrel(fs0->cvorticity);
DEBUG_MSG("old field %d %g %g\\n", fs0->iteration, a, b);
copy_complex_array<{0}>(fs0->cd, fs0->cvorticity,
fs1->cd, fs1->cvorticity,
3);
fs1->write('v', 'c');
a = 0.5*fs1->correl_vec(fs1->cvelocity, fs1->cvelocity);
b = 0.5*fs1->correl_vec(fs1->cvorticity, fs1->cvorticity);
a = 0.5*fs1->autocorrel(fs1->cvelocity);
b = 0.5*fs1->autocorrel(fs1->cvorticity);
DEBUG_MSG("new field %d %g %g\\n", fs1->iteration, a, b);
niter_todo = 0;
//endcpp
......
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