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

simplify io in final code

parent 068de907
Branches
Tags
No related merge requests found
......@@ -51,12 +51,14 @@ void fluid_solver<rnumber>::impose_zero_modes()
\
template<> \
fluid_solver<R>::fluid_solver( \
const char *NAME, \
int nx, \
int ny, \
int nz, \
double DKX, \
double DKY, \
double DKZ) : fluid_solver_base<R>(nx , ny , nz, \
double DKZ) : fluid_solver_base<R>(NAME, \
nx , ny , nz, \
DKX, DKY, DKZ) \
{ \
this->cvorticity = FFTW(alloc_complex)(this->cd->local_size);\
......@@ -306,6 +308,52 @@ void fluid_solver<R>::step(double dt) \
this->force_divfree(this->cv[0]); \
\
this->iteration++; \
} \
\
template<> \
int fluid_solver<R>::read(char field, char representation) \
{ \
if ((field == 'v') && (representation == 'c')) \
return this->read_base("cvorticity", this->cvorticity); \
if ((field == 'v') && (representation == 'r')) \
{ \
int read_result = this->read_base("rvorticity", this->rvorticity); \
if (!(read_result == EXIT_SUCCESS)) \
return read_result; \
else \
{ \
FFTW(execute)(*((FFTW(plan)*)this->r2c_vorticity )); \
return EXIT_SUCCESS; \
} \
} \
if ((field == 'u') && (representation == 'c')) \
return this->read_base("cvelocity", this->cvelocity); \
if ((field == 'u') && (representation == 'r')) \
return this->read_base("rvelocity", this->rvelocity); \
return EXIT_FAILURE; \
} \
\
template<> \
int fluid_solver<R>::write(char field, char representation) \
{ \
if ((field == 'v') && (representation == 'c')) \
return this->write_base("cvorticity", this->cvorticity); \
if ((field == 'v') && (representation == 'r')) \
{ \
FFTW(execute)(*((FFTW(plan)*)this->c2r_vorticity )); \
clip_zero_padding<R>(this->rd, this->rvorticity, 3); \
return this->write_base("rvorticity", this->rvorticity); \
} \
this->compute_velocity(this->cvorticity); \
if ((field == 'u') && (representation == 'c')) \
return this->write_base("cvelocity", this->cvelocity); \
if ((field == 'u') && (representation == 'r')) \
{ \
FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \
clip_zero_padding<R>(this->rd, this->rvelocity, 3); \
return this->write_base("rvelocity", this->rvelocity); \
} \
return EXIT_FAILURE; \
}
/*****************************************************************************/
......@@ -319,12 +367,12 @@ FLUID_SOLVER_DEFINITIONS(
fftwf_complex,
MPI_REAL4,
MPI_COMPLEX8)
FLUID_SOLVER_DEFINITIONS(
FFTW_MANGLE_DOUBLE,
double,
fftw_complex,
MPI_REAL8,
MPI_COMPLEX16)
//FLUID_SOLVER_DEFINITIONS(
// FFTW_MANGLE_DOUBLE,
// double,
// fftw_complex,
// MPI_REAL8,
// MPI_COMPLEX16)
/*****************************************************************************/
......
......@@ -64,6 +64,7 @@ class fluid_solver:public fluid_solver_base<rnumber>
/* methods */
fluid_solver(
const char *NAME,
int nx,
int ny,
int nz,
......@@ -77,6 +78,9 @@ class fluid_solver:public fluid_solver_base<rnumber>
void omega_nonlin(int src);
void step(double dt);
void impose_zero_modes(void);
int read(char field, char representation);
int write(char field, char representation);
};
#endif//FLUID_SOLVER
......
......@@ -21,6 +21,7 @@
#include <cassert>
#include <cmath>
#include <cstring>
#include "fluid_solver_base.hpp"
#include "fftw_tools.hpp"
......@@ -33,6 +34,7 @@
\
template<> \
fluid_solver_base<R>::fluid_solver_base( \
const char *NAME, \
int nx, \
int ny, \
int nz, \
......@@ -40,6 +42,8 @@ fluid_solver_base<R>::fluid_solver_base( \
double DKY, \
double DKZ) \
{ \
strncpy(this->name, NAME, 256); \
this->name[255] = '\0'; \
this->iteration = 0; \
\
int ntmp[4]; \
......@@ -259,6 +263,38 @@ void fluid_solver_base<R>::symmetrize(C *data, const int howmany) \
FFTW(free)(buffer); \
delete mpistatus; \
} \
\
template<> \
int fluid_solver_base<R>::read_base(const char *fname, R *data) \
{ \
char full_name[512]; \
sprintf(full_name, "%s_%s_i%.5x", this->name, fname, this->iteration); \
return this->rd->read(full_name, (void*)data); \
} \
\
template<> \
int fluid_solver_base<R>::read_base(const char *fname, C *data) \
{ \
char full_name[512]; \
sprintf(full_name, "%s_%s_i%.5x", this->name, fname, this->iteration); \
return this->cd->read(full_name, (void*)data); \
} \
\
template<> \
int fluid_solver_base<R>::write_base(const char *fname, R *data) \
{ \
char full_name[512]; \
sprintf(full_name, "%s_%s_i%.5x", this->name, fname, this->iteration); \
return this->rd->write(full_name, (void*)data); \
} \
\
template<> \
int fluid_solver_base<R>::write_base(const char *fname, C *data) \
{ \
char full_name[512]; \
sprintf(full_name, "%s_%s_i%.5x", this->name, fname, this->iteration); \
return this->cd->write(full_name, (void*)data); \
}
/*****************************************************************************/
......@@ -272,12 +308,12 @@ FLUID_SOLVER_BASE_DEFINITIONS(
fftwf_complex,
MPI_REAL4,
MPI_COMPLEX8)
FLUID_SOLVER_BASE_DEFINITIONS(
FFTW_MANGLE_DOUBLE,
double,
fftw_complex,
MPI_REAL8,
MPI_COMPLEX16)
//FLUID_SOLVER_BASE_DEFINITIONS(
// FFTW_MANGLE_DOUBLE,
// double,
// fftw_complex,
// MPI_REAL8,
// MPI_COMPLEX16)
/*****************************************************************************/
......
......@@ -46,6 +46,7 @@ class fluid_solver_base
ptrdiff_t normalization_factor;
/* simulation parameters */
char name[256];
int iteration;
/* physical parameters */
......@@ -63,6 +64,7 @@ class fluid_solver_base
/* methods */
fluid_solver_base(
const char *NAME,
int nx,
int ny,
int nz,
......@@ -75,6 +77,10 @@ class fluid_solver_base
void force_divfree(cnumber *a);
void symmetrize(cnumber *a, int howmany);
rnumber correl_vec(cnumber *a, cnumber *b);
int read_base(const char *fname, rnumber *data);
int read_base(const char *fname, cnumber *data);
int write_base(const char *fname, rnumber *data);
int write_base(const char *fname, cnumber *data);
};
......
......@@ -108,7 +108,7 @@ class stat_test(code):
//begincpp
fluid_solver<float> *fs;
char fname[512];
fs = new fluid_solver<float>(nx, ny, nz);
fs = new fluid_solver<float>(simname, nx, ny, nz);
fs->nu = nu;
fs->iteration = iter0;
if (myrank == fs->cd->io_myrank)
......@@ -116,28 +116,14 @@ class stat_test(code):
sprintf(fname, "%s_stats.bin", simname);
stat_file = fopen(fname, "wb");
}
sprintf(fname, "%s_kvorticity_i%.5x", simname, fs->iteration);
fs->cd->read(
fname,
(void*)fs->cvorticity);
fs->read('v', 'c');
fs->low_pass_Fourier(fs->cvorticity, 3, fs->kM);
fs->force_divfree(fs->cvorticity);
fs->symmetrize(fs->cvorticity, 3);
t = 0.0;
do_stats(fs);
fftwf_execute(*((fftwf_plan*)fs->c2r_velocity));
sprintf(fname, "%s_rvelocity_i%.5x", simname, fs->iteration);
clip_zero_padding<float>(fs->rd, fs->rvelocity, 3);
fs->rd->write(
fname,
(void*)fs->rvelocity);
fftwf_execute(*((fftwf_plan*)fs->c2r_vorticity));
sprintf(fname, "%s_rvorticity_i%.5x", simname, fs->iteration);
clip_zero_padding<float>(fs->rd, fs->rvorticity, 3);
fs->rd->write(
fname,
(void*)fs->rvorticity);
fs->write('u', 'r');
fs->write('v', 'r');
for (; fs->iteration < iter0 + niter_todo;)
{
fs->step(dt);
......@@ -145,22 +131,9 @@ class stat_test(code):
do_stats(fs);
}
fclose(stat_file);
sprintf(fname, "%s_kvorticity_i%.5x", simname, fs->iteration);
fs->cd->write(
fname,
(void*)fs->cvorticity);
fftwf_execute(*((fftwf_plan*)fs->c2r_vorticity));
sprintf(fname, "%s_rvorticity_i%.5x", simname, fs->iteration);
clip_zero_padding<float>(fs->rd, fs->rvorticity, 3);
fs->rd->write(
fname,
(void*)fs->rvorticity);
fftwf_execute(*((fftwf_plan*)fs->c2r_velocity));
sprintf(fname, "%s_rvelocity_i%.5x", simname, fs->iteration);
clip_zero_padding<float>(fs->rd, fs->rvelocity, 3);
fs->rd->write(
fname,
(void*)fs->rvelocity);
fs->write('v', 'c');
fs->write('v', 'r');
fs->write('u', 'r');
delete fs;
//endcpp
"""
......@@ -195,12 +168,12 @@ if __name__ == '__main__':
c.parameters['dt'] = 0.004
c.parameters['niter_todo'] = opt.nsteps
c.write_par(simname = 'test1')
Kdata0.tofile("test1_kvorticity_i00000")
Kdata0.tofile("test1_cvorticity_i00000")
c.run(ncpu = opt.ncpu, simname = 'test1')
c.parameters['dt'] /= 2
c.parameters['niter_todo'] *= 2
c.write_par(simname = 'test2')
Kdata0.tofile("test2_kvorticity_i00000")
Kdata0.tofile("test2_cvorticity_i00000")
c.run(ncpu = opt.ncpu, simname = 'test2')
Rdata = np.fromfile(
'test1_rvorticity_i00000',
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment