diff --git a/cmake/parameters.cmake b/cmake/parameters.cmake index ad808295..d0d1a87a 100644 --- a/cmake/parameters.cmake +++ b/cmake/parameters.cmake @@ -19,6 +19,8 @@ set_if_empty(NCHEMISTRY 0) set_if_empty(NTRACER 0) +set_if_empty(NTURBULENCE 0) + set_if_empty(NSTATIC 0) # set_if_empty(TASKLIST ImplicitHydroTasks) diff --git a/cmake/robert.cmake b/cmake/robert.cmake index 70532d60..734cda34 100644 --- a/cmake/robert.cmake +++ b/cmake/robert.cmake @@ -12,3 +12,4 @@ set_if_empty(NUMBER_GHOST_CELLS 3) # canoe configure set(MPI ON) set(PNETCDF ON) +set(TASKLIST ImplicitHydroTasks) diff --git a/configure.hpp.in b/configure.hpp.in index 3dd0135d..d9566e4d 100644 --- a/configure.hpp.in +++ b/configure.hpp.in @@ -6,6 +6,7 @@ enum { NCLOUD = @NCLOUD@, NCHEMISTRY = @NCHEMISTRY@, NTRACER = @NTRACER@, + NTURBULENCE = @NTURBULENCE@, NSTATIC = @NSTATIC@ }; diff --git a/doc/README-parallel-debug.md b/doc/README-parallel-debug.md new file mode 100644 index 00000000..9f14005e --- /dev/null +++ b/doc/README-parallel-debug.md @@ -0,0 +1 @@ +mpirun -n 2 xterm -e gdb robert.debug diff --git a/src/air_parcel.hpp b/src/air_parcel.hpp index bb3a6657..33814526 100644 --- a/src/air_parcel.hpp +++ b/src/air_parcel.hpp @@ -15,7 +15,9 @@ // \brief a collection of all physical data in a computational cell class AirParcel { public: - enum { Size = NHYDRO + NCLOUD + NCHEMISTRY + NTRACER + NSTATIC }; + enum { + Size = NHYDRO + NCLOUD + NCHEMISTRY + NTRACER + NTURBULENCE + NSTATIC + }; enum class Type { MassFrac = 0, MassConc = 1, MoleFrac = 2, MoleConc = 3 }; @@ -24,16 +26,19 @@ class AirParcel { //! data pointers //! hydro data - Real *w; + Real *const w; //! cloud data - Real *c; + Real *const c; //! chemistry data - Real *q; + Real *const q; //! tracer data - Real *x; + Real *const x; + + //! turbulence data + Real *const t; //! static data Real const *s; @@ -42,26 +47,29 @@ class AirParcel { Real const *d; // constructor - explicit AirParcel(Type type = Type::MoleFrac) : mytype_(type) { - w = data_.data(); - c = w + NHYDRO; - q = c + NCLOUD; - x = q + NCHEMISTRY; - s = x + NTRACER; - d = s + NSTATIC; + explicit AirParcel(Type type = Type::MoleFrac) + : mytype_(type), + w(data_.data()), + c(w + NHYDRO), + q(c + NCLOUD), + x(q + NCHEMISTRY), + t(x + NTRACER), + s(t + NTURBULENCE), + d(s + NSTATIC) { std::fill(data_.begin(), data_.end(), 0.0); } // copy constructor - AirParcel(AirParcel const &other) : mytype_(other.mytype_) { - data_ = other.data_; - w = data_.data(); - c = w + NHYDRO; - q = c + NCLOUD; - x = q + NCHEMISTRY; - s = x + NTRACER; - d = s + NSTATIC; - } + AirParcel(AirParcel const &other) + : mytype_(other.mytype_), + data_(other.data_), + w(data_.data()), + c(w + NHYDRO), + q(c + NCLOUD), + x(q + NCHEMISTRY), + t(x + NTRACER), + s(x + NTURBULENCE), + d(s + NSTATIC) {} // Assignment operator AirParcel &operator=(const AirParcel &other) { diff --git a/src/impl.cpp b/src/impl.cpp index ba3e605b..04305590 100644 --- a/src/impl.cpp +++ b/src/impl.cpp @@ -24,6 +24,7 @@ #include "snap/decomposition/decomposition.hpp" #include "snap/implicit/implicit_solver.hpp" #include "snap/thermodynamics/thermodynamics.hpp" +#include "snap/turbulence/turbulence_model.hpp" // microphysics #include "microphysics/microphysics.hpp" @@ -51,7 +52,19 @@ MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) { ptracer = std::make_shared(pmb, pin); // turbulence - // pturb = std::make_shared(pmb, pin); + if (pin->DoesParameterExist("hydro", "turbulence")) { + std::string turbulence_model = pin->GetString("hydro", "turbulence"); + if (turbulence_model == "none") { + pturb = std::make_shared(pmb, pin); + } else if (turbulence_model == "kepsilon") { + pturb = std::make_shared(pmb, pin); + if (NTURBULENCE < 2) + throw NotImplementedError( + "NTURBULENCE must be at least 2 for k-epsilon model"); + } else { + throw NotImplementedError(turbulence_model); + } + } // radiation prad = std::make_shared(pmb, pin); @@ -59,8 +72,8 @@ MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) { // inversion queue fitq = Inversion::NewInversionQueue(pmb, pin); - // reference pressure #ifdef HYDROSTATIC + // reference pressure reference_pressure_ = pin->GetReal("mesh", "ReferencePressure"); // pressure scale height diff --git a/src/snap/turbulence/add_turbulence_flux_divergence.cpp b/src/snap/turbulence/add_turbulence_flux_divergence.cpp deleted file mode 100644 index 173d02e1..00000000 --- a/src/snap/turbulence/add_turbulence_flux_divergence.cpp +++ /dev/null @@ -1,100 +0,0 @@ -// Athena++ headers -#include -#include -#include -#include - -// snap -#include "turbulence_model.hpp" - -// OpenMP header -#ifdef OPENMP_PARALLEL -#include -#endif - -//---------------------------------------------------------------------------------------- -//! \fn void TurbulenceModel::AddFluxDivergence -// \brief Adds flux divergence to weighted average of conservative variables -// from previous step(s) of time integrator algorithm - -// TODO(felker): after the removal of AddCoordTermsDivergence() fn call from -// Hydro::AddFluxDivergence(), the 2x fns could be trivially shared if: -// - flux/s_flux renamed to the same class member name -// - 7x below references of x1face_area_ ... dflx_ private class members (which -// are only ever used in this fn and are members to prevent de/allocating each -// fn call) -// - NHYDRO/NSCALARS is replaced with array_out.GetDim4() - -// ----> Hydro should be derived from Turbulence - -// TODO(felker): remove the following unnecessary private class member? -// field_diffusion.cpp:66: cell_volume_.NewAthenaArray(nc1); - -void TurbulenceModel::addFluxDivergence(const Real wght, - AthenaArray &s_out) { - MeshBlock *pmb = pmy_block; - AthenaArray &x1flux = s_flux[X1DIR]; - AthenaArray &x2flux = s_flux[X2DIR]; - AthenaArray &x3flux = s_flux[X3DIR]; - int is = pmb->is; - int js = pmb->js; - int ks = pmb->ks; - int ie = pmb->ie; - int je = pmb->je; - int ke = pmb->ke; - AthenaArray &x1area = x1face_area_, &x2area = x2face_area_, - &x2area_p1 = x2face_area_p1_, &x3area = x3face_area_, - &x3area_p1 = x3face_area_p1_, &vol = cell_volume_, - &dflx = dflx_; - - int nvar = s_out.GetDim4(); - for (int k = ks; k <= ke; ++k) { - for (int j = js; j <= je; ++j) { - // calculate x1-flux divergence - pmb->pcoord->Face1Area(k, j, is, ie + 1, x1area); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = is; i <= ie; ++i) { - dflx(n, i) = (x1area(i + 1) * x1flux(n, k, j, i + 1) - - x1area(i) * x1flux(n, k, j, i)); - } - } - - // calculate x2-flux divergence - if (pmb->block_size.nx2 > 1) { - pmb->pcoord->Face2Area(k, j, is, ie, x2area); - pmb->pcoord->Face2Area(k, j + 1, is, ie, x2area_p1); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = is; i <= ie; ++i) { - dflx(n, i) += (x2area_p1(i) * x2flux(n, k, j + 1, i) - - x2area(i) * x2flux(n, k, j, i)); - } - } - } - - // calculate x3-flux divergence - if (pmb->block_size.nx3 > 1) { - pmb->pcoord->Face3Area(k, j, is, ie, x3area); - pmb->pcoord->Face3Area(k + 1, j, is, ie, x3area_p1); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = is; i <= ie; ++i) { - dflx(n, i) += (x3area_p1(i) * x3flux(n, k + 1, j, i) - - x3area(i) * x3flux(n, k, j, i)); - } - } - } - - // update conserved variables - pmb->pcoord->CellVolume(k, j, is, ie, vol); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = is; i <= ie; ++i) { - s_out(n, k, j, i) -= wght * dflx(n, i) / vol(i); - } - } - } - } - return; -} diff --git a/src/snap/turbulence/k_epsilon_turbulence.cpp b/src/snap/turbulence/k_epsilon_turbulence.cpp index 6a49907a..df0ce263 100644 --- a/src/snap/turbulence/k_epsilon_turbulence.cpp +++ b/src/snap/turbulence/k_epsilon_turbulence.cpp @@ -6,11 +6,9 @@ // Athena++ headers #include -#include #include -#include // reapply floors to face-centered reconstructed states #include -#include +#include // snap #include "turbulence_model.hpp" @@ -22,6 +20,7 @@ KEpsilonTurbulence::KEpsilonTurbulence(MeshBlock *pmb, ParameterInput *pin) c2_ = pin->GetOrAddReal("turbulence", "kepsilon.c2", 1.92); sigk_ = pin->GetOrAddReal("turbulence", "kepsilon.sigk", 1.0); sige_ = pin->GetOrAddReal("turbulence", "kepsilon.sige", 1.3); + // velocity scale, v ~ 1 m/s, tke ~ v^2 ~ 1 m^2/s^2 Real tke0 = pin->GetOrAddReal("turbulence", "kepsilon.tke0", 1.); @@ -43,23 +42,26 @@ KEpsilonTurbulence::KEpsilonTurbulence(MeshBlock *pmb, ParameterInput *pin) for (int k = kl; k <= ku; ++k) for (int j = jl; j <= ju; ++j) for (int i = il; i <= iu; ++i) { - r(1, k, j, i) = tke0; + w(1, k, j, i) = tke0; // eps ~ tke/T ~ v^2/(dx/v) ~ v^3/dx Real volume = pmb->pcoord->GetCellVolume(k, j, i); Real eps0 = pow(tke0, 1.5) / pow(volume, 1. / 3.); - r(0, k, j, i) = eps0; + w(0, k, j, i) = eps0; } } -void KEpsilonTurbulence::Initialize(AthenaArray const &w) { - MeshBlock *pmb = pmy_block; +void KEpsilonTurbulence::Initialize() { + auto pmb = pmy_block; + auto phydro = pmb->phydro; + for (int k = pmb->ks; k <= pmb->ke; ++k) for (int j = pmb->js; j <= pmb->je; ++j) for (int i = pmb->is; i <= pmb->ie; ++i) { - s(0, k, j, i) = r(0, k, j, i) * w(IDN, k, j, i); - s(1, k, j, i) = r(1, k, j, i) * w(IDN, k, j, i); - mut(k, j, i) = cmu_ * w(IDN, k, j, i) * r(1, k, j, i) * r(1, k, j, i) / - r(0, k, j, i); + Real rhod = phydro->w(IDN, k, j, i); + u(0, k, j, i) = w(0, k, j, i) * rhod; + u(1, k, j, i) = w(1, k, j, i) * rhod; + mut(k, j, i) = + cmu_ * rhod * w(1, k, j, i) * w(1, k, j, i) / w(0, k, j, i); } } @@ -143,22 +145,23 @@ inline Real ShearProduction_(AthenaArray const &w, int k, int j, int i, return result; } -void KEpsilonTurbulence::driveTurbulence(AthenaArray &s, - AthenaArray const &r, - AthenaArray const &w, Real dt) { +void KEpsilonTurbulence::DriveTurbulence(Real dt) { // std::cout << "driving turbulence" << std::endl; - MeshBlock *pmb = pmy_block; + auto pmb = pmy_block; + auto phydro = pmb->phydro; + int ks = pmb->ks, js = pmb->js, is = pmb->is; int ke = pmb->ke, je = pmb->je, ie = pmb->ie; AthenaArray eps, tke; - eps.InitWithShallowSlice(const_cast &>(r), 4, 0, 1); - tke.InitWithShallowSlice(const_cast &>(r), 4, 1, 1); + eps.InitWithShallowSlice(w, 4, 0, 1); + tke.InitWithShallowSlice(w, 4, 1, 1); Real s1, s2, s3; for (int k = ks; k <= ke; ++k) for (int j = js; j <= je; ++j) for (int i = is; i <= ie; ++i) { + Real rhod = phydro->w(IDN, k, j, i); // shear production Real shear = ShearProduction_(w, k, j, i, pmb->pcoord); // std::cout << "shear = " << shear << std::endl; @@ -166,36 +169,35 @@ void KEpsilonTurbulence::driveTurbulence(AthenaArray &s, // turbulent dissipation, de/dt, eq2.2-1 s1 = c1_ * mut(k, j, i) * eps(k, j, i) / tke(k, j, i) * shear; s2 = Laplace_(mut, eps, k, j, i, pmb->pcoord) / sige_; - s3 = - -c2_ * eps(k, j, i) * eps(k, j, i) / tke(k, j, i) * w(IDN, k, j, i); - s(0, k, j, i) += (s1 + s2) * dt; - if (s(0, k, j, i) + s3 * dt > 0) - s(0, k, j, i) += s3 * dt; + s3 = -c2_ * eps(k, j, i) * eps(k, j, i) / tke(k, j, i) * rhod; + u(0, k, j, i) += (s1 + s2) * dt; + if (u(0, k, j, i) + s3 * dt > 0) + u(0, k, j, i) += s3 * dt; else - s(0, k, j, i) /= 2.; - // std::cout << "s = " << s(0,k,j,i) << std::endl; + u(0, k, j, i) /= 2.; + // std::cout << "u = " << u(0,k,j,i) << std::endl; // turbulent kinetic energy, dk/dt, eq 2.2-2 s1 = mut(k, j, i) * shear; s2 = Laplace_(mut, tke, k, j, i, pmb->pcoord) / sigk_; - s3 = -eps(k, j, i) * w(IDN, k, j, i); - s(1, k, j, i) += (s1 + s2) * dt; - if (s(1, k, j, i) + s3 * dt > 0.) - s(1, k, j, i) += s3 * dt; + s3 = -eps(k, j, i) * rhod; + u(1, k, j, i) += (s1 + s2) * dt; + if (u(1, k, j, i) + s3 * dt > 0.) + u(1, k, j, i) += s3 * dt; else - s(1, k, j, i) /= 2.; + u(1, k, j, i) /= 2.; } } -void KEpsilonTurbulence::setDiffusivity(AthenaArray &nu, +void KEpsilonTurbulence::SetDiffusivity(AthenaArray &nu, AthenaArray &kappa, const AthenaArray &prim, const AthenaArray &bcc, int il, int iu, int jl, int ju, int kl, int ku) { AthenaArray eps, tke; - eps.InitWithShallowSlice(r, 4, 0, 1); - tke.InitWithShallowSlice(r, 4, 1, 1); + eps.InitWithShallowSlice(w, 4, 0, 1); + tke.InitWithShallowSlice(w, 4, 1, 1); for (int k = kl; k <= ku; ++k) { for (int j = jl; j <= ju; ++j) { diff --git a/src/snap/turbulence/turbulence_model.cpp b/src/snap/turbulence/turbulence_model.cpp index 21f1a60a..4b746455 100644 --- a/src/snap/turbulence/turbulence_model.cpp +++ b/src/snap/turbulence/turbulence_model.cpp @@ -5,12 +5,14 @@ // Athena++ headers #include -#include -#include -#include -#include #include -#include +#include + +// application +#include + +// canoe +#include // snap #include "turbulence_model.hpp" @@ -18,400 +20,27 @@ // constructor, initializes data structures and parameters TurbulenceModel::TurbulenceModel(MeshBlock *pmb, ParameterInput *pin, int nvar) - : s(nvar, pmb->ncells3, pmb->ncells2, pmb->ncells1), - s1(nvar, pmb->ncells3, pmb->ncells2, pmb->ncells1), - r(nvar, pmb->ncells3, pmb->ncells2, pmb->ncells1), - s_flux{{nvar, pmb->ncells3, pmb->ncells2, pmb->ncells1 + 1}, - {nvar, pmb->ncells3, pmb->ncells2 + 1, pmb->ncells1, - (pmb->pmy_mesh->f2 ? AthenaArray::DataStatus::allocated - : AthenaArray::DataStatus::empty)}, - {nvar, pmb->ncells3 + 1, pmb->ncells2, pmb->ncells1, - (pmb->pmy_mesh->f3 ? AthenaArray::DataStatus::allocated - : AthenaArray::DataStatus::empty)}}, - coarse_s_( - nvar, pmb->ncc3, pmb->ncc2, pmb->ncc1, - (pmb->pmy_mesh->multilevel ? AthenaArray::DataStatus::allocated - : AthenaArray::DataStatus::empty)), - coarse_r_( - nvar, pmb->ncc3, pmb->ncc2, pmb->ncc1, - (pmb->pmy_mesh->multilevel ? AthenaArray::DataStatus::allocated - : AthenaArray::DataStatus::empty)), - sbvar(pmb, &s, &coarse_s_, s_flux, true), - mut(pmb->ncells3, pmb->ncells2, pmb->ncells1), - pmy_block(pmb) { - int nc1 = pmb->ncells1, nc2 = pmb->ncells2, nc3 = pmb->ncells3; - Mesh *pm = pmy_block->pmy_mesh; - - pmb->RegisterMeshBlockData(s); + : pmy_block(pmb) { + if (NTURBULENCE == 0) return; - // If user-requested time integrator is type 3S*, allocate additional memory - // registers - std::string integrator = pin->GetOrAddString("time", "integrator", "vl2"); - if (integrator == "ssprk5_4" || STS_ENABLED) { - // future extension may add "int nregister" to Hydro class - s2.NewAthenaArray(nvar, nc3, nc2, nc1); - } + Application::Logger app("snap"); + app->Log("Initialize Turbulence"); - // "Enroll" in SMR/AMR by adding to vector of pointers in MeshRefinement class - if (pm->multilevel) { - refinement_idx = pmy_block->pmr->AddToRefinement(&s, &coarse_s_); - } + w.InitWithShallowSlice(pmb->pscalars->r, 4, NCLOUD + NCHEMISTRY + NTRACER, + NTURBULENCE); + u.InitWithShallowSlice(pmb->pscalars->s, 4, NCLOUD + NCHEMISTRY + NTRACER, + NTURBULENCE); - // enroll CellCenteredBoundaryVariable object - sbvar.bvar_index = pmb->pbval->bvars.size(); - pmb->pbval->bvars.push_back(&sbvar); - pmb->pbval->bvars_main_int.push_back(&sbvar); + int ncells1 = pmb->ncells1; + int ncells2 = pmb->ncells2; + int ncells3 = pmb->ncells3; - // Allocate memory for scratch arrays - rl_.NewAthenaArray(nvar, nc1); - rr_.NewAthenaArray(nvar, nc1); - rlb_.NewAthenaArray(nvar, nc1); - x1face_area_.NewAthenaArray(nc1 + 1); - if (pm->f2) { - x2face_area_.NewAthenaArray(nc1); - x2face_area_p1_.NewAthenaArray(nc1); - } - if (pm->f3) { - x3face_area_.NewAthenaArray(nc1); - x3face_area_p1_.NewAthenaArray(nc1); - } - cell_volume_.NewAthenaArray(nc1); - dflx_.NewAthenaArray(nvar, nc1); - - mut.ZeroClear(); - s.ZeroClear(); - s1.ZeroClear(); - r.ZeroClear(); + mut.NewAthenaArray(ncells3, ncells2, ncells1); } -//---------------------------------------------------------------------------------------- -//! \fn void Turbulence::calculateFluxes -// \brief Calculate passive scalar fluxes using reconstruction + weighted -// upwinding rule - -void TurbulenceModel::calculateFluxes(AthenaArray &r, const int order) { - MeshBlock *pmb = pmy_block; - - // design decision: do not pass Hydro::flux (for mass flux) via function - // parameters, since 1) it is unlikely that anything else would be passed, 2) - // the current TurbulenceModel class/feature implementation is inherently - // coupled to Hydro class 3) high-order calculation of scalar fluxes will - // require other Hydro flux approximations (flux_fc in calculate_fluxes.cpp is - // currently not saved persistently in Hydro class but each flux dir is temp. - // stored in 4D scratch array scr1_nkji_) - - Hydro &hyd = *(pmb->phydro); - - AthenaArray &x1flux = s_flux[X1DIR]; - AthenaArray mass_flux; - mass_flux.InitWithShallowSlice(hyd.flux[X1DIR], 4, IDN, 1); - int is = pmb->is; - int js = pmb->js; - int ks = pmb->ks; - int ie = pmb->ie; - int je = pmb->je; - int ke = pmb->ke; - int il, iu, jl, ju, kl, ku; - - //-------------------------------------------------------------------------------------- - // i-direction - - // set the loop limits - jl = js, ju = je, kl = ks, ku = ke; - // TODO(felker): fix loop limits for fourth-order hydro - // if (MAGNETIC_FIELDS_ENABLED) { - if (pmb->block_size.nx2 > 1) { - if (pmb->block_size.nx3 == 1) // 2D - jl = js - 1, ju = je + 1, kl = ks, ku = ke; - else // 3D - jl = js - 1, ju = je + 1, kl = ks - 1, ku = ke + 1; - } - // } - - int nvar = s.GetDim4(); - - for (int k = kl; k <= ku; ++k) { - for (int j = jl; j <= ju; ++j) { - // reconstruct L/R states - if (order == 1) { - pmb->precon->DonorCellX1(k, j, is - 1, ie + 1, r, rl_, rr_); - } else if (order == 2) { - pmb->precon->PiecewiseLinearX1(k, j, is - 1, ie + 1, r, rl_, rr_); - } else { - pmb->precon->PiecewiseParabolicX1(k, j, is - 1, ie + 1, r, rl_, rr_); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = is; i <= ie + 1; ++i) { - // apply floor correction - rl_(n, i) = std::max(rl_(n, i), 0.); - rr_(n, i) = std::max(rr_(n, i), 0.); - } - } - } - - computeUpwindFlux(k, j, is, ie + 1, rl_, rr_, mass_flux, x1flux); - } - } - - //------------------------------------------------------------------------------ - // end x1 fourth-order hydro - - //-------------------------------------------------------------------------------------- - // j-direction - - if (pmb->pmy_mesh->f2) { - AthenaArray &x2flux = s_flux[X2DIR]; - mass_flux.InitWithShallowSlice(hyd.flux[X2DIR], 4, IDN, 1); - - // set the loop limits - il = is - 1, iu = ie + 1, kl = ks, ku = ke; - // TODO(felker): fix loop limits for fourth-order hydro - // if (MAGNETIC_FIELDS_ENABLED) { - if (pmb->block_size.nx3 == 1) // 2D - kl = ks, ku = ke; - else // 3D - kl = ks - 1, ku = ke + 1; - // } - - for (int k = kl; k <= ku; ++k) { - // reconstruct the first row - if (order == 1) { - pmb->precon->DonorCellX2(k, js - 1, il, iu, r, rl_, rr_); - } else if (order == 2) { - pmb->precon->PiecewiseLinearX2(k, js - 1, il, iu, r, rl_, rr_); - } else { - pmb->precon->PiecewiseParabolicX2(k, js - 1, il, iu, r, rl_, rr_); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = il; i <= iu; ++i) { - rl_(n, i) = std::max(rl_(n, i), 0.); - // pmb->peos->ApplyPassiveScalarFloors(rr_, n, k, j, i); - } - } - } - for (int j = js; j <= je + 1; ++j) { - // reconstruct L/R states at j - if (order == 1) { - pmb->precon->DonorCellX2(k, j, il, iu, r, rlb_, rr_); - } else if (order == 2) { - pmb->precon->PiecewiseLinearX2(k, j, il, iu, r, rlb_, rr_); - } else { - pmb->precon->PiecewiseParabolicX2(k, j, il, iu, r, rlb_, rr_); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = il; i <= iu; ++i) { - rlb_(n, i) = std::max(rlb_(n, i), 0.); - rr_(n, i) = std::max(rr_(n, i), 0.); - } - } - } - - computeUpwindFlux(k, j, il, iu, rl_, rr_, mass_flux, x2flux); - - // swap the arrays for the next step - rl_.SwapAthenaArray(rlb_); - } - } - } - - //-------------------------------------------------------------------------------------- - // k-direction - - if (pmb->pmy_mesh->f3) { - AthenaArray &x3flux = s_flux[X3DIR]; - mass_flux.InitWithShallowSlice(hyd.flux[X3DIR], 4, IDN, 1); - - // set the loop limits - // TODO(felker): fix loop limits for fourth-order hydro - // if (MAGNETIC_FIELDS_ENABLED) - il = is - 1, iu = ie + 1, jl = js - 1, ju = je + 1; - - for (int j = jl; j <= ju; ++j) { // this loop ordering is intentional - // reconstruct the first row - if (order == 1) { - pmb->precon->DonorCellX3(ks - 1, j, il, iu, r, rl_, rr_); - } else if (order == 2) { - pmb->precon->PiecewiseLinearX3(ks - 1, j, il, iu, r, rl_, rr_); - } else { - pmb->precon->PiecewiseParabolicX3(ks - 1, j, il, iu, r, rl_, rr_); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = il; i <= iu; ++i) { - rl_(n, ks - 1, j, i) = std::max(rl_(n, ks - 1, j, i), 0.); - // pmb->peos->ApplyPassiveScalarFloors(rr_, n, k, j, i); - } - } - } - for (int k = ks; k <= ke + 1; ++k) { - // reconstruct L/R states at k - if (order == 1) { - pmb->precon->DonorCellX3(k, j, il, iu, r, rlb_, rr_); - } else if (order == 2) { - pmb->precon->PiecewiseLinearX3(k, j, il, iu, r, rlb_, rr_); - } else { - pmb->precon->PiecewiseParabolicX3(k, j, il, iu, r, rlb_, rr_); - for (int n = 0; n < nvar; ++n) { -#pragma omp simd - for (int i = il; i <= iu; ++i) { - rlb_(n, i) = std::max(rlb_(n, i), 0.); - rr_(n, i) = std::max(rr_(n, i), 0.); - } - } - } - - computeUpwindFlux(k, j, il, iu, rl_, rr_, mass_flux, x3flux); - - // swap the arrays for the next step - rl_.SwapAthenaArray(rlb_); - } - } - } - - return; -} - -void TurbulenceModel::computeUpwindFlux( - const int k, const int j, const int il, - const int iu, // CoordinateDirection dir, - AthenaArray &rl, AthenaArray &rr, // 2D - AthenaArray &mass_flx, // 3D - AthenaArray &flx_out) { // 4D - for (int n = 0; n < flx_out.GetDim4(); n++) { -#pragma omp simd - for (int i = il; i <= iu; i++) { - Real fluid_flx = mass_flx(k, j, i); - if (fluid_flx >= 0.0) - flx_out(n, k, j, i) = fluid_flx * rl(n, i); - else - flx_out(n, k, j, i) = fluid_flx * rr(n, i); - } - } - return; -} - -void TurbulenceModel::ConservedToPrimitive(AthenaArray &s, - AthenaArray const &w, - AthenaArray &r, - Coordinates *pco, int il, int iu, - int jl, int ju, int kl, int ku) { - int nvar = s.GetDim4(); - for (int n = 0; n < nvar; ++n) - for (int k = kl; k <= ku; ++k) - for (int j = jl; j <= ju; ++j) - for (int i = il; i <= iu; ++i) { - r(n, k, j, i) = s(n, k, j, i) / w(IDN, k, j, i); - } -} - -void TurbulenceModel::PrimitiveToConserved(AthenaArray &r, - AthenaArray const &w, - AthenaArray &s, - Coordinates *pco, int il, int iu, - int jl, int ju, int kl, int ku) { - int nvar = s.GetDim4(); - for (int n = 0; n < nvar; ++n) - for (int k = kl; k <= ku; ++k) - for (int j = jl; j <= ju; ++j) - for (int i = il; i <= iu; ++i) { - s(n, k, j, i) = r(n, k, j, i) * w(IDN, k, j, i); - } -} - -void TurbulenceModel::applyBoundaryCondition(AthenaArray &r, - AthenaArray &s, - AthenaArray const &w, - Coordinates *pco) { - MeshBlock *pmb = pmy_block; - int nvar = r.GetDim4(); - int bis = pmb->is - NGHOST, bie = pmb->ie + NGHOST, bjs = pmb->js, - bje = pmb->je, bks = pmb->ks, bke = pmb->ke; - - // inner x1 - if (pmb->pbval->isPhysicalBoundary(BoundaryFace::inner_x1)) { - for (int n = 0; n < nvar; ++n) - for (int k = bks; k <= bke; ++k) - for (int j = bjs; j <= bje; ++j) - for (int i = pmb->is - NGHOST; i <= pmb->is - 1; ++i) { - r(n, k, j, i) = r(n, k, j, pmb->is); - s(n, k, j, i) = r(n, k, j, i) * w(IDN, k, j, i); - } - } - - // outer x1 - if (pmb->pbval->isPhysicalBoundary(BoundaryFace::outer_x1)) { - for (int n = 0; n < nvar; ++n) - for (int k = bks; k <= bke; ++k) - for (int j = bjs; j <= bje; ++j) - for (int i = pmb->ie + 1; i <= pmb->ie + NGHOST; ++i) { - r(n, k, j, i) = r(n, k, j, pmb->ie); - s(n, k, j, i) = r(n, k, j, i) * w(IDN, k, j, i); - } - } - - if (pmb->pmy_mesh->f2) { - // inner x2 - if (pmb->pbval->isPhysicalBoundary(BoundaryFace::inner_x2)) { - for (int n = 0; n < nvar; ++n) - for (int k = bks; k <= bke; ++k) - for (int j = pmb->js - NGHOST; j <= pmb->js - 1; ++j) - for (int i = bis; i <= bie; ++i) { - r(n, k, j, i) = r(n, k, pmb->js, i); - s(n, k, j, i) = r(n, k, j, i) * w(IDN, k, j, i); - } - } - - // outer x2 - if (pmb->pbval->isPhysicalBoundary(BoundaryFace::outer_x2)) { - for (int n = 0; n < nvar; ++n) - for (int k = bks; k <= bke; ++k) - for (int j = pmb->je + 1; j <= pmb->je + NGHOST; ++j) - for (int i = bis; i <= bie; ++i) { - r(n, k, j, i) = r(n, k, pmb->je, i); - s(n, k, j, i) = r(n, k, j, i) * w(IDN, k, j, i); - } - } - } - - if (pmb->pmy_mesh->f3) { - bjs = pmb->js - NGHOST; - bje = pmb->je + NGHOST; - - // inner x3 - if (pmb->pbval->isPhysicalBoundary(BoundaryFace::inner_x3)) { - for (int n = 0; n < nvar; ++n) - for (int k = pmb->ks - NGHOST; k <= pmb->js - 1; ++k) - for (int j = bjs; j <= bje; ++j) - for (int i = bis; i <= bie; ++i) { - r(n, k, j, i) = r(n, pmb->ks, j, i); - s(n, k, j, i) = r(n, k, j, i) * w(IDN, k, j, i); - } - } - - // outer x3 - if (pmb->pbval->isPhysicalBoundary(BoundaryFace::outer_x3)) { - for (int n = 0; n < nvar; ++n) - for (int k = pmb->ke + 1; k <= pmb->ke + NGHOST; ++k) - for (int j = bjs; j <= bje; ++j) - for (int i = bis; i <= bie; ++i) { - r(n, k, j, i) = r(n, pmb->ke, j, i); - s(n, k, j, i) = r(n, k, j, i) * w(IDN, k, j, i); - } - } - } -} - -size_t TurbulenceModel::getRestartDataSizeInBytes() { - return s.GetSizeInBytes(); -} - -size_t TurbulenceModel::dumpRestartData(char *pdst) { - std::memcpy(pdst, s.data(), s.GetSizeInBytes()); - return s.GetSizeInBytes(); -} +TurbulenceModel::~TurbulenceModel() { + if (NTURBULENCE == 0) return; -size_t TurbulenceModel::loadRestartData(char *psrc) { - std::memcpy(s.data(), psrc, s.GetSizeInBytes()); - // load it into the other memory register(s) too - std::memcpy(s1.data(), psrc, s1.GetSizeInBytes()); - return s.GetSizeInBytes(); + Application::Logger app("snap"); + app->Log("Destroy Turbulence"); } diff --git a/src/snap/turbulence/turbulence_model.hpp b/src/snap/turbulence/turbulence_model.hpp index eed9693b..8e565fb0 100644 --- a/src/snap/turbulence/turbulence_model.hpp +++ b/src/snap/turbulence/turbulence_model.hpp @@ -17,65 +17,25 @@ class ParameterInput; class TurbulenceModel { public: - TurbulenceModel(MeshBlock *pmb, ParameterInput *pin, int nvar); - virtual ~TurbulenceModel() {} + TurbulenceModel(MeshBlock *pmb, ParameterInput *pin, int nvar = 0); + virtual ~TurbulenceModel(); - // public data: - // "conserved vars" = extensive variable - AthenaArray s, s1, s2; // (no more than MAX_NREGISTER allowed) - // "primitive vars" = intensive variable - AthenaArray r; // , r1; - AthenaArray s_flux[3]; // face-averaged flux vector - AthenaArray mut; // dynamic turbulent viscosity + // access members + AthenaArray w, u; - // storage for SMR/AMR - // TODO(KGF): remove trailing underscore or revert to private: - AthenaArray coarse_s_, coarse_r_; - int refinement_idx{-1}; - - CellCenteredBoundaryVariable sbvar; + // public data + AthenaArray mut; // dynamic turbulent viscosity // public functions: - void addFluxDivergence(const Real wght, AthenaArray &s_out); - void calculateFluxes(AthenaArray &s, const int order); - void ConservedToPrimitive(AthenaArray &s, AthenaArray const &w, - AthenaArray &r, Coordinates *pco, int il, - int iu, int jl, int ju, int kl, int ku); - void PrimitiveToConserved(AthenaArray &r, AthenaArray const &w, - AthenaArray &s, Coordinates *pco, int il, - int iu, int jl, int ju, int kl, int ku); - void computeUpwindFlux(const int k, const int j, const int il, - const int iu, // CoordinateDirection dir, - AthenaArray &rl, AthenaArray &rr, - AthenaArray &mass_flx, - AthenaArray &flx_out); - void applyBoundaryCondition(AthenaArray &r, AthenaArray &s, - AthenaArray const &w, Coordinates *pco); - - virtual void driveTurbulence(AthenaArray &s, AthenaArray const &r, - AthenaArray const &w, Real dt) {} - virtual void Initialize(AthenaArray const &w) {} - - virtual void setDiffusivity(AthenaArray &nu, AthenaArray &kappa, + virtual void DriveTurbulence(Real dt) {} + virtual void Initialize() {} + virtual void SetDiffusivity(AthenaArray &nu, AthenaArray &kappa, const AthenaArray &w, const AthenaArray &bc, int il, int iu, int jl, int ju, int kl, int ku) {} - // restart functions - size_t getRestartDataSizeInBytes(); - size_t dumpRestartData(char *pdst); - size_t loadRestartData(char *psrc); - protected: MeshBlock *pmy_block; - // scratch space used to compute fluxes - // 2D scratch arrays - AthenaArray rl_, rr_, rlb_; - // 1D scratch arrays - AthenaArray x1face_area_, x2face_area_, x3face_area_; - AthenaArray x2face_area_p1_, x3face_area_p1_; - AthenaArray cell_volume_; - AthenaArray dflx_; }; using TurbulenceModelPtr = std::shared_ptr; @@ -84,11 +44,10 @@ class KEpsilonTurbulence : public TurbulenceModel { public: KEpsilonTurbulence(MeshBlock *pmb, ParameterInput *pin); ~KEpsilonTurbulence() {} - void driveTurbulence(AthenaArray &s, AthenaArray const &r, - AthenaArray const &w, Real dt); - void Initialize(AthenaArray const &w); - void setDiffusivity(AthenaArray &nu, AthenaArray &kappa, + void DriveTurbulence(Real dt); + void Initialize(); + void SetDiffusivity(AthenaArray &nu, AthenaArray &kappa, const AthenaArray &w, const AthenaArray &bc, int il, int iu, int jl, int ju, int kl, int ku); diff --git a/src/snap/turbulence/turbulence_tasks.cpp b/src/snap/turbulence/turbulence_tasks.cpp deleted file mode 100644 index 4a005e48..00000000 --- a/src/snap/turbulence/turbulence_tasks.cpp +++ /dev/null @@ -1,125 +0,0 @@ -//! \file turbulence_tasks.cpp -//! \brief declared in task_list/task_list.hpp - -// C/C++ headers -#include - -// Athena++ headers -#include -#include - -// canoe -#include - -// tasklist -#include - -// snap -#include "turbulence_model.hpp" - -TaskStatus ImplicitHydroTasks::CalculateTurbulenceFlux(MeshBlock *pmb, - int stage) { - // std::cout << "calculate turbulence flux" << std::endl; - auto pturb = pmb->pimpl->pturb; - if (stage <= nstages) { - pturb->calculateFluxes(pturb->r, 2); - return TaskStatus::next; - } - return TaskStatus::fail; -} - -TaskStatus ImplicitHydroTasks::SendTurbulenceFlux(MeshBlock *pmb, int stage) { - // std::cout << "send turbulence flux" << std::endl; - pmb->pimpl->pturb->sbvar.SendFluxCorrection(); - return TaskStatus::success; -} - -TaskStatus ImplicitHydroTasks::ReceiveTurbulenceFlux(MeshBlock *pmb, - int stage) { - // std::cout << "receive turbulence flux" << std::endl; - if (pmb->pimpl->pturb->sbvar.ReceiveFluxCorrection()) { - return TaskStatus::next; - } else { - return TaskStatus::fail; - } -} - -TaskStatus ImplicitHydroTasks::IntegrateTurbulence(MeshBlock *pmb, int stage) { - // std::cout << "integrate turbulence" << std::endl; - auto pturb = pmb->pimpl->pturb; - if (stage <= nstages) { - // This time-integrator-specific averaging operation logic is identical to - // IntegrateHydro, IntegrateField - Real ave_wghts[5]; - ave_wghts[0] = 1.0; - ave_wghts[1] = stage_wghts[stage - 1].delta; - ave_wghts[2] = 0.0; - ave_wghts[3] = 0.0; - ave_wghts[4] = 0.0; - pmb->WeightedAve(pturb->s1, pturb->s, pturb->s2, pturb->s2, pturb->s2, - ave_wghts); - - ave_wghts[0] = stage_wghts[stage - 1].gamma_1; - ave_wghts[1] = stage_wghts[stage - 1].gamma_2; - ave_wghts[2] = stage_wghts[stage - 1].gamma_3; - if (ave_wghts[0] == 0.0 && ave_wghts[1] == 1.0 && ave_wghts[2] == 0.0) - pturb->s.SwapAthenaArray(pturb->s1); - else - pmb->WeightedAve(pturb->s, pturb->s1, pturb->s2, pturb->s2, pturb->s2, - ave_wghts); - - const Real wght = stage_wghts[stage - 1].beta * pmb->pmy_mesh->dt; - pturb->addFluxDivergence(wght, pturb->s); - - // turbulence forcing - Real dt = (stage_wghts[(stage - 1)].beta) * (pmb->pmy_mesh->dt); - pturb->driveTurbulence(pturb->s, pturb->r, pmb->phydro->w, dt); - return TaskStatus::next; - } - return TaskStatus::fail; -} - -TaskStatus ImplicitHydroTasks::SendTurbulence(MeshBlock *pmb, int stage) { - // std::cout << "send turbulence" << std::endl; - auto pturb = pmb->pimpl->pturb; - if (stage <= nstages) { - // Swap TurbulenceModel quantity in BoundaryVariable interface back to - // conserved var formulation (also needed in SetBoundariesTurbulence() since - // the tasks are independent) - pturb->sbvar.var_cc = &(pturb->s); - pturb->sbvar.SendBoundaryBuffers(); - } else { - return TaskStatus::fail; - } - return TaskStatus::success; -} - -TaskStatus ImplicitHydroTasks::ReceiveTurbulence(MeshBlock *pmb, int stage) { - // std::cout << "recv turbulence" << std::endl; - bool ret; - if (stage <= nstages) { - ret = pmb->pimpl->pturb->sbvar.ReceiveBoundaryBuffers(); - } else { - return TaskStatus::fail; - } - if (ret) { - return TaskStatus::success; - } else { - return TaskStatus::fail; - } - return TaskStatus::success; -} - -TaskStatus ImplicitHydroTasks::SetBoundariesTurbulence(MeshBlock *pmb, - int stage) { - // std::cout << "set turbulence boundary" << std::endl; - auto pturb = pmb->pimpl->pturb; - if (stage <= nstages) { - // Set Turbulence quantity in BoundaryVariable interface to cons var - // formulation - pturb->sbvar.var_cc = &(pturb->s); - pturb->sbvar.SetBoundaries(); - return TaskStatus::success; - } - return TaskStatus::fail; -}