Skip to content
Snippets Groups Projects
Commit 29070325 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

broken 

no idea what's wrong...
parent 45b1006a
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
bfps/cpp/fluid_solver.cpp
+ 15
− 8
View file @ 29070325
......@@ -52,7 +52,7 @@ fluid_solver<R>::fluid_solver( \
double DKY, \
double DKZ) : fluid_solver_base<R>(NAME, \
nx , ny , nz, \
DKX, DKY, DKZ) \
DKX, DKY, DKZ, 1) \
{ \
this->cvorticity = FFTW(alloc_complex)(this->cd->local_size);\
this->cvelocity = FFTW(alloc_complex)(this->cd->local_size);\
......@@ -174,7 +174,8 @@ fluid_solver<R>::~fluid_solver() \
template<> \
void fluid_solver<R>::compute_vorticity() \
{ \
this->low_pass_Fourier(this->cu, 3, this->kM); \
/*this->low_pass_Fourier(this->cu, 3, this->kM);*/ \
this->dealias(this->cu, 3); \
CLOOP( \
this->cvorticity[3*cindex+0][0] = -(this->ky[yindex]*this->cu[3*cindex+2][1] - this->kz[zindex]*this->cu[3*cindex+1][1]); \
this->cvorticity[3*cindex+1][0] = -(this->kz[zindex]*this->cu[3*cindex+0][1] - this->kx[xindex]*this->cu[3*cindex+2][1]); \
......@@ -190,7 +191,8 @@ template<> \
void fluid_solver<R>::compute_velocity(C *vorticity) \
{ \
double k2; \
this->low_pass_Fourier(vorticity, 3, this->kM); \
/*this->low_pass_Fourier(vorticity, 3, this->kM);*/ \
this->dealias(vorticity, 3); \
std::fill_n((R*)this->cu, this->cd->local_size*2, 0.0); \
CLOOP( \
k2 = (this->kx[xindex]*this->kx[xindex] + \
......@@ -307,7 +309,8 @@ void fluid_solver<R>::omega_nonlin( \
this->clean_up_real_space(this->ru, 3); \
FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
this->symmetrize(this->cu, 3); \
this->low_pass_Fourier(this->cu, 3, this->kM); \
/*this->low_pass_Fourier(this->cu, 3, this->kM);*/ \
this->dealias(this->cu, 3); \
this->force_divfree(this->cu); \
/* $\imath k \times Fourier(u \times \omega)$ */ \
R tmpx1, tmpy1, tmpz1; \
......@@ -357,7 +360,8 @@ void fluid_solver<R>::step(double dt) \
\
this->force_divfree(this->cv[1]); \
this->omega_nonlin(1); \
this->low_pass_Fourier(this->cv[1], 3, this->kM); \
/*this->low_pass_Fourier(this->cv[1], 3, this->kM);*/ \
this->dealias(this->cv[1], 3); \
CLOOP( \
k2 = (this->kx[xindex]*this->kx[xindex] + \
this->ky[yindex]*this->ky[yindex] + \
......@@ -374,7 +378,8 @@ void fluid_solver<R>::step(double dt) \
\
this->force_divfree(this->cv[2]); \
this->omega_nonlin(2); \
this->low_pass_Fourier(this->cv[2], 3, this->kM); \
/*this->low_pass_Fourier(this->cv[2], 3, this->kM);*/ \
this->dealias(this->cv[2], 3); \
CLOOP( \
k2 = (this->kx[xindex]*this->kx[xindex] + \
this->ky[yindex]*this->ky[yindex] + \
......@@ -389,7 +394,8 @@ void fluid_solver<R>::step(double dt) \
); \
this->symmetrize(this->cu, 3); \
this->force_divfree(this->cv[0]); \
this->low_pass_Fourier(this->cv[0], 3, this->kM); \
/*this->low_pass_Fourier(this->cv[0], 3, this->kM);*/ \
this->dealias(this->cv[0], 3); \
/*this->write_base("cvorticity", this->cvorticity); \
this->read_base("cvorticity", this->cvorticity); \
this->low_pass_Fourier(this->cvorticity, 3, this->kM); \
......@@ -428,7 +434,8 @@ int fluid_solver<R>::read(char field, char representation) \
else \
FFTW(execute)(*((FFTW(plan)*)this->r2c_vorticity )); \
} \
this->low_pass_Fourier(this->cvorticity, 3, this->kM); \
/*this->low_pass_Fourier(this->cvorticity, 3, this->kM);*/ \
this->dealias(this->cvorticity, 3); \
this->force_divfree(this->cvorticity); \
this->symmetrize(this->cvorticity, 3); \
return EXIT_SUCCESS; \
......
bfps/cpp/fluid_solver_base.cpp
+ 59
− 7
View file @ 29070325
......@@ -21,9 +21,9 @@
// code is generally compiled via setuptools, therefore NDEBUG is present
//#ifdef NDEBUG
//#undef NDEBUG
//#endif//NDEBUG
#ifdef NDEBUG
#undef NDEBUG
#endif//NDEBUG
#include <cassert>
#include <cmath>
......@@ -116,7 +116,8 @@ fluid_solver_base<R>::fluid_solver_base( \
int nz, \
double DKX, \
double DKY, \
double DKZ) \
double DKZ, \
int dealias_type) \
{ \
strncpy(this->name, NAME, 256); \
this->name[255] = '\0'; \
......@@ -143,9 +144,20 @@ fluid_solver_base<R>::fluid_solver_base( \
this->kx = new double[this->cd->sizes[2]]; \
this->ky = new double[this->cd->subsizes[0]]; \
this->kz = new double[this->cd->sizes[1]]; \
this->kMx = this->dkx*(int(this->rd->sizes[2] / 3)-1); \
this->kMy = this->dky*(int(this->rd->sizes[1] / 3)-1); \
this->kMz = this->dkz*(int(this->rd->sizes[0] / 3)-1); \
this->dealias_type = dealias_type; \
switch(this->dealias_type) \
{ \
/* HL07 smooth filter */ \
case 1: \
this->kMx = this->dkx*(int(this->rd->sizes[2] / 2)-1); \
this->kMy = this->dky*(int(this->rd->sizes[1] / 2)-1); \
this->kMz = this->dkz*(int(this->rd->sizes[0] / 2)-1); \
break; \
default: \
this->kMx = this->dkx*(int(this->rd->sizes[2] / 3)-1); \
this->kMy = this->dky*(int(this->rd->sizes[1] / 3)-1); \
this->kMz = this->dkz*(int(this->rd->sizes[0] / 3)-1); \
} \
int i, ii; \
for (i = 0; i<this->cd->sizes[2]; i++) \
{ \
......@@ -173,6 +185,14 @@ fluid_solver_base<R>::fluid_solver_base( \
this->dk = this->dkx; \
if (this->dk > this->dky) this->dk = this->dky; \
if (this->dk > this->dkz) this->dk = this->dkz; \
this->dk2 = this->dk*this->dk; \
this->Fourier_filter_size = int(this->kM2 / (this->dk*this->dk))+1; \
this->Fourier_filter = new double[this->Fourier_filter_size]; \
for (int i=0; i < this->Fourier_filter_size; i++) \
{ \
this->Fourier_filter[i] = exp(-36*pow(sqrt(double(i))*this->dk / this->kM, 36)); \
/*DEBUG_MSG("%d %lg\n", i, this->Fourier_filter[i]);*/ \
} \
/* spectra stuff */ \
this->nshells = int(this->kM / this->dk) + 2; \
this->kshell = new double[this->nshells]; \
......@@ -229,6 +249,7 @@ fluid_solver_base<R>::fluid_solver_base( \
template<> \
fluid_solver_base<R>::~fluid_solver_base() \
{ \
delete[] this->Fourier_filter; \
delete[] this->kshell; \
delete[] this->nshell; \
\
......@@ -298,6 +319,37 @@ void fluid_solver_base<R>::low_pass_Fourier(C *a, const int howmany, const doubl
} \
\
template<> \
void fluid_solver_base<R>::dealias(C *a, const int howmany) \
{ \
if (this->dealias_type == 0) \
{ \
this->low_pass_Fourier(a, howmany, this->kM); \
return; \
} \
double k2; \
double tval; \
CLOOP( \
k2 = (this->kx[xindex]*this->kx[xindex] + \
this->ky[yindex]*this->ky[yindex] + \
this->kz[zindex]*this->kz[zindex]); \
tval = this->Fourier_filter[int(k2/this->dk2)]; \
for (int tcounter = 0; tcounter < howmany; tcounter++) \
{ \
a[howmany*cindex+tcounter][0] *= tval; \
a[howmany*cindex+tcounter][1] *= tval; \
} \
if ((xindex == this->rd->sizes[2]/2) || \
(zindex == this->rd->sizes[1]/2) || \
((yindex == this->cd->subsizes[0]) && (this->cd->myrank == this->cd->rank[this->rd->sizes[0]/2]))) \
for (int tcounter = 0; tcounter < howmany; tcounter++) \
{ \
a[howmany*cindex+tcounter][0] = 0; \
a[howmany*cindex+tcounter][1] = 0; \
} \
);\
} \
\
template<> \
void fluid_solver_base<R>::force_divfree(C *a) \
{ \
double k2; \
......
bfps/cpp/fluid_solver_base.hpp
+ 7
− 5
View file @ 29070325
......@@ -50,13 +50,14 @@ class fluid_solver_base
int iteration;
/* physical parameters */
rnumber dkx, dky, dkz, dk;
rnumber dkx, dky, dkz, dk, dk2;
/* mode and dealiasing information */
int dealias_type;
double kMx, kMy, kMz, kM, kM2;
double *kx, *ky, *kz;
//double *Fourier_filter;
//int Fourier_filter_size;
double *Fourier_filter;
int Fourier_filter_size;
double *kshell;
int64_t *nshell;
int nshells;
......@@ -70,11 +71,12 @@ class fluid_solver_base
int nz,
double DKX = 1.0,
double DKY = 1.0,
double DKZ = 1.0);
double DKZ = 1.0,
int dealias_type = 0);
~fluid_solver_base();
void low_pass_Fourier(cnumber *a, int howmany, double kmax);
//void dealias(cnumber *a, int howmany, double kmax);
void dealias(cnumber *a, int howmany);
void force_divfree(cnumber *a);
void symmetrize(cnumber *a, int howmany);
void clean_up_real_space(rnumber *a, int howmany);
......
bfps/resize.py
+ 5
− 1
View file @ 29070325
......@@ -68,7 +68,8 @@ class vorticity_resize(bfps.code):
fs0->iteration = iter0;
fs1->iteration = 0;
fs0->read('v', 'c');
fs0->low_pass_Fourier(fs0->cvorticity, 3, fs0->kM);
//fs0->low_pass_Fourier(fs0->cvorticity, 3, fs0->kM);
fs0->dealias(fs0->cvorticity, 3);
fs0->force_divfree(fs0->cvorticity);
fs0->symmetrize(fs0->cvorticity, 3);
fs0->write('v', 'r');
......@@ -80,6 +81,9 @@ class vorticity_resize(bfps.code):
copy_complex_array(fs0->cd, fs0->cvorticity,
fs1->cd, fs1->cvorticity,
3);
fs1->dealias(fs1->cvorticity, 3);
fs1->force_divfree(fs1->cvorticity);
fs1->symmetrize(fs1->cvorticity, 3);
fs1->write('v', 'c');
fs1->write('v', 'r');
fs1->write('u', 'r');
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment