diff --git a/CREDITS b/CREDITS index 448267ee..1489656a 100644 --- a/CREDITS +++ b/CREDITS @@ -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 diff --git a/libmpdata++/formulae/stress_formulae.hpp b/libmpdata++/formulae/stress_formulae.hpp index b009dccf..7724a73c 100644 --- a/libmpdata++/formulae/stress_formulae.hpp +++ b/libmpdata++/formulae/stress_formulae.hpp @@ -492,6 +492,37 @@ namespace libmpdataxx real_t(0.5) * (G(rho, ijk[0], ijk[1], ijk[2] + 1) + G(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 + 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::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(rho, ijk[0] + 1, ijk[1], ijk[2]) + G(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(rho, ijk[0], ijk[1] + 1, ijk[2]) + G(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(rho, ijk[0], ijk[1], ijk[2] + 1) + G(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(rho, ijk[0], ijk[1], ijk[2] + 1) + G(rho, ijk[0], ijk[1], ijk[2])); + } + // multiplication of compact tensor components by constant molecular viscosity // 2D version template @@ -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 + 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::type* = 0) + { + multiply_vctr_cmpct(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(rho, ijk[0] + 1, ijk[1] , ijk[2]) + + G(rho, ijk[0] , ijk[1] , ijk[2]) + + G(rho, ijk[0] + 1, ijk[1] + 1, ijk[2]) + + G(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(rho, ijk[0] + 1, ijk[1], ijk[2] ) + + G(rho, ijk[0] , ijk[1], ijk[2] ) + + G(rho, ijk[0] + 1, ijk[1], ijk[2] + 1) + + G(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(rho, ijk[0], ijk[1] , ijk[2] + 1) + + G(rho, ijk[0], ijk[1] , ijk[2] ) + + G(rho, ijk[0], ijk[1] + 1, ijk[2] + 1) + + G(rho, ijk[0], ijk[1] + 1, ijk[2] ) + ); + } + // flux divergence // 2D version template diff --git a/libmpdata++/output/detail/output_common.hpp b/libmpdata++/output/detail/output_common.hpp index 47920259..1a7d968d 100644 --- a/libmpdata++/output/detail/output_common.hpp +++ b/libmpdata++/output/detail/output_common.hpp @@ -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; @@ -158,7 +158,7 @@ namespace libmpdataxx record_time = this->time; for (int t = 0; t < outwindow; ++t) { - if ((this->timestep - t) % static_cast(outfreq) == 0) record_all(); + if ((this->timestep - t) % static_cast(outfreq) == 0 && (this->timestep - t) >= static_cast(outstart)) record_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 outvars; std::string outdir; @@ -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, {"", ""}}}; diff --git a/libmpdata++/output/gnuplot.hpp b/libmpdata++/output/gnuplot.hpp index 9ffd8f4d..432e1163 100644 --- a/libmpdata++/output/gnuplot.hpp +++ b/libmpdata++/output/gnuplot.hpp @@ -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 << ", '-'"; diff --git a/libmpdata++/output/hdf5.hpp b/libmpdata++/output/hdf5.hpp index cbf3a81c..56766fdc 100644 --- a/libmpdata++/output/hdf5.hpp +++ b/libmpdata++/output/hdf5.hpp @@ -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 coord(nt_out); - - for(int i=0; ioutstart == 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)); @@ -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) { @@ -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