Skip to content

Commit

Permalink
Merge branch 'master' of github.com:igfuw/libmpdataxx into outwindow_T
Browse files Browse the repository at this point in the history
  • Loading branch information
@pdziekan
pdziekan committed Sep 10, 2024
2 parents de403f4 + 5dbea9b commit 38f01f9
Show file tree
Hide file tree
Showing 5 changed files with 131 additions and 11 deletions.
1 change: 1 addition & 0 deletions CREDITS
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ The development team consists of (in alphabetic order):
Sylwester Arabas
(core code, library design, OOP-formulae concepts [4],
concurrency, gnuplot output)
Piotr Dziekan
Dorota Jarecka
(shallow-water systems, OOP-formulae concepts [4])
Anna Jaruga
Expand Down
81 changes: 81 additions & 0 deletions libmpdata++/formulae/stress_formulae.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -492,6 +492,37 @@ namespace libmpdataxx
real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1], ijk[2] + 1) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
}

// multiplication of compact vector components by variable anisotropic eddy viscosity
// 3D version
// TODO: when used in tensor multiplication, it is not clear which eddy viscosity components should be used for each term; see Simon and Chow 2021 p. 17
// we opt to use the same approach as therein, i.e. horizontal for all diagonal components; this is controlled by the flag vh
template <int nd, opts_t opts, class arrvec_t, class real_t, class arr_t, class ijk_t>
inline void multiply_vctr_cmpct(const arrvec_t &av,
real_t coeff,
const arrvec_t & k_ma, // two arrays: [0] for horizontal eddy viscosity and [1] for vertical
const arr_t & rho,
const ijk_t &ijk,
const bool vh = false,
typename std::enable_if<nd == 3>::type* = 0)
{
av[0](ijk[0] + h, ijk[1], ijk[2]) *= coeff *
real_t(0.5) * (k_ma[0](ijk[0] + 1, ijk[1], ijk[2]) + k_ma[0](ijk[0], ijk[1], ijk[2])) *
real_t(0.5) * (G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2]) + G<opts, 0>(rho,ijk[0], ijk[1], ijk[2]));

av[1](ijk[0], ijk[1] + h, ijk[2]) *= coeff *
real_t(0.5) * (k_ma[0](ijk[0], ijk[1] + 1, ijk[2]) + k_ma[0](ijk[0], ijk[1], ijk[2])) *
real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2]) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));

if(!vh)
av[2](ijk[0], ijk[1], ijk[2] + h) *= coeff *
real_t(0.5) * (k_ma[1](ijk[0], ijk[1], ijk[2] + 1) + k_ma[1](ijk[0], ijk[1], ijk[2])) *
real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1], ijk[2] + 1) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
else
av[2](ijk[0], ijk[1], ijk[2] + h) *= coeff *
real_t(0.5) * (k_ma[0](ijk[0], ijk[1], ijk[2] + 1) + k_ma[0](ijk[0], ijk[1], ijk[2])) *
real_t(0.5) * (G<opts, 0>(rho, ijk[0], ijk[1], ijk[2] + 1) + G<opts, 0>(rho, ijk[0], ijk[1], ijk[2]));
}

// multiplication of compact tensor components by constant molecular viscosity
// 2D version
template <int nd, class arrvec_t, class real_t, class ijk_t>
Expand Down Expand Up @@ -587,6 +618,56 @@ namespace libmpdataxx
);
}

// multiplication of compact tensor components by variable anisotropic eddy viscosity
// 3D version
// TODO: it is not clear which eddy viscosity components should be used for each term; see Simon and Chow 2021 p. 17
// we opt to use the same approach as therein
template <int nd, opts_t opts, class arrvec_t, class real_t, class arr_t, class ijk_t>
inline void multiply_tnsr_cmpct(const arrvec_t &av,
const real_t coeff,
const arrvec_t &k_ma,
const arr_t &rho,
const ijk_t &ijk,
typename std::enable_if<nd == 3>::type* = 0)
{
multiply_vctr_cmpct<nd, opts>(av, coeff, k_ma, rho, ijk, true);
av[3](ijk[0] + h, ijk[1] + h, ijk[2]) *= coeff *
real_t(0.25) * ( k_ma[0](ijk[0] + 1, ijk[1] , ijk[2])
+ k_ma[0](ijk[0] , ijk[1] , ijk[2])
+ k_ma[0](ijk[0] + 1, ijk[1] + 1, ijk[2])
+ k_ma[0](ijk[0] , ijk[1] + 1, ijk[2])
) *
real_t(0.25) * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1] , ijk[2])
+ G<opts, 0>(rho, ijk[0] , ijk[1] , ijk[2])
+ G<opts, 0>(rho, ijk[0] + 1, ijk[1] + 1, ijk[2])
+ G<opts, 0>(rho, ijk[0] , ijk[1] + 1, ijk[2])
);

av[4](ijk[0] + h, ijk[1], ijk[2] + h) *= coeff *
real_t(0.25) * ( k_ma[1](ijk[0] + 1, ijk[1], ijk[2] )
+ k_ma[1](ijk[0] , ijk[1], ijk[2] )
+ k_ma[1](ijk[0] + 1, ijk[1], ijk[2] + 1)
+ k_ma[1](ijk[0] , ijk[1], ijk[2] + 1)
) *
real_t(0.25) * ( G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2] )
+ G<opts, 0>(rho, ijk[0] , ijk[1], ijk[2] )
+ G<opts, 0>(rho, ijk[0] + 1, ijk[1], ijk[2] + 1)
+ G<opts, 0>(rho, ijk[0] , ijk[1], ijk[2] + 1)
);

av[5](ijk[0], ijk[1] + h, ijk[2] + h) *= coeff *
real_t(0.25) * ( k_ma[1](ijk[0], ijk[1] , ijk[2] + 1)
+ k_ma[1](ijk[0], ijk[1] , ijk[2] )
+ k_ma[1](ijk[0], ijk[1] + 1, ijk[2] + 1)
+ k_ma[1](ijk[0], ijk[1] + 1, ijk[2] )
) *
real_t(0.25) * ( G<opts, 0>(rho, ijk[0], ijk[1] , ijk[2] + 1)
+ G<opts, 0>(rho, ijk[0], ijk[1] , ijk[2] )
+ G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2] + 1)
+ G<opts, 0>(rho, ijk[0], ijk[1] + 1, ijk[2] )
);
}

// flux divergence
// 2D version
template <int nd, opts_t opts, class arrvec_t, class arr_t, class ijk_t, class dijk_t>
Expand Down
11 changes: 8 additions & 3 deletions libmpdata++/output/detail/output_common.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ namespace libmpdataxx

int do_record_cnt = 0;
typename parent_t::real_t record_time;
const typename parent_t::advance_arg_t outfreq;
const typename parent_t::advance_arg_t outfreq, outstart;
const int outwindow;
const std::string outdir;

Expand Down Expand Up @@ -158,7 +158,7 @@ namespace libmpdataxx
record_time = this->time;
for (int t = 0; t < outwindow; ++t)
{
if ((this->timestep - t) % static_cast<int>(outfreq) == 0) record_all();
if ((this->timestep - t) % static_cast<int>(outfreq) == 0 && (this->timestep - t) >= static_cast<int>(outstart)) record_all();
}
}
}
Expand All @@ -170,7 +170,8 @@ namespace libmpdataxx

struct rt_params_t : parent_t::rt_params_t
{
typename parent_t::advance_arg_t outfreq = 1;
typename parent_t::advance_arg_t outfreq = 1,
outstart = 0; // TODO: make it work for var_dt
int outwindow = 1;
std::map<int, info_t> outvars;
std::string outdir;
Expand All @@ -184,11 +185,15 @@ namespace libmpdataxx
) :
parent_t(args, p),
outfreq(p.outfreq),
outstart(p.outstart),
outwindow(p.outwindow),
outvars(p.outvars),
outdir(p.outdir),
intrp_vars(args.mem->tmp[__FILE__][0])
{
assert(int(outstart) % int(outfreq) == 0);
assert(int(outstart)==0 || this->var_dt==0);

// default value for outvars
if (this->outvars.size() == 0 && parent_t::n_eqns == 1)
outvars = {{0, {"", ""}}};
Expand Down
1 change: 1 addition & 0 deletions libmpdata++/output/gnuplot.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,7 @@ namespace libmpdataxx

for (int t = 0; t <= nt; t+=p.outfreq)
{
if(t > 0 && t < p.outstart) continue;
for (const auto &v : this->outvars)
{
*gp << ", '-'";
Expand Down
48 changes: 40 additions & 8 deletions libmpdata++/output/hdf5.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,17 +198,19 @@ namespace libmpdataxx

// T
{
const hsize_t
nt_out = (nt / this->outfreq + 1) * this->outwindow // incl. t=0
- (this->outwindow - (nt % this->outfreq) - 1); // timsteps from outwindow that go beyond nt
// incl. t=0 and t=outstart
const hsize_t nt_out = ((nt - this->outstart) / this->outfreq + (this->outstart == 0 ? 1 : 2)) * this->outwindow // for outstart>0 we still want to store t=0
- (this->outwindow - (nt % this->outfreq) - 1); // timsteps from outwindow that go beyond nt

float dt = this->dt;

blitz::Array<typename solver_t::real_t, 1> coord(nt_out);

for(int i=0; i<nt_out; ++i)
if(this->outstart == 0)
coord(blitz::Range(0,nt_out-1)) = (this->var_dt ? this->outfreq : this->outfreq * this->dt) * (floor(blitz::firstIndex() / this->outwindow) * this->outfreq + fmod(blitz::firstIndex(), this->outwindow) + this->outstart/this->outfreq);
else
{
int timestep = int(i / this->outwindow) * this->outfreq + (i % this->outwindow);
coord(i) = timestep * (this->var_dt ? this->outfreq : this->outfreq * this->dt);
coord(blitz::Range(1,nt_out-1)) = (this->var_dt ? this->outfreq : this->outfreq * this->dt) * (floor(blitz::firstIndex() / this->outwindow) * this->outfreq + fmod(blitz::firstIndex(), this->outwindow) + this->outstart/this->outfreq);
coord(0)=0;
}

auto curr_dim = (*hdfp).createDataSet("T", flttype_output, H5::DataSpace(1, &nt_out));
Expand Down Expand Up @@ -371,6 +373,22 @@ namespace libmpdataxx
aux.write(data, flttype_solver, H5::DataSpace(parent_t::n_dims, shape.data()), space, dxpl_id);
}

// for 1-D arrays
void record_aux_hlpr(const std::string &name, typename solver_t::real_t *data, hsize_t size, H5::H5File hdf)
{
assert(this->rank == 0);

auto aux = hdf.createDataSet(
name,
flttype_output,
H5::DataSpace(1, &size)
);

auto space = aux.getSpace();
space.selectHyperslab(H5S_SELECT_SET, &size, &zero);
aux.write(data, flttype_solver, H5::DataSpace(1, &size), space, dxpl_id);
}

// for discontiguous array with halos
void record_aux_dsc_hlpr(const std::string &name, const typename solver_t::arr_t &arr, H5::H5File hdf, bool srfc = false)
{
Expand Down Expand Up @@ -522,10 +540,24 @@ namespace libmpdataxx
// has to be called after const file was created (i.e. after start())
void record_aux_const(const std::string &name, typename solver_t::real_t *data)
{
H5::H5File hdfcp(const_file, H5F_ACC_RDWR); // reopen the const file
H5::H5File hdfcp(const_file, H5F_ACC_RDWR
#if defined(USE_MPI)
, H5P_DEFAULT, fapl_id
#endif
); // reopen the const file
record_aux_hlpr(name, data, hdfcp);
}

void record_aux_const(const std::string &name, typename solver_t::real_t *data, const int &size)
{
H5::H5File hdfcp(const_file, H5F_ACC_RDWR
#if defined(USE_MPI)
, H5P_DEFAULT, fapl_id
#endif
); // reopen the const file
record_aux_hlpr(name, data, size, hdfcp);
}

void record_aux_const(const std::string &name, const std::string &group_name, typename solver_t::real_t data)
{
H5::H5File hdfcp(const_file, H5F_ACC_RDWR
Expand Down

0 comments on commit 38f01f9

Please sign in to comment.