diff --git a/cmake/examples/rtdry.cmake b/cmake/examples/rtdry.cmake index cf4311fc..3c039389 100644 --- a/cmake/examples/rtdry.cmake +++ b/cmake/examples/rtdry.cmake @@ -1,11 +1,5 @@ # configuration for euler equations on cubed sphere -macro(SET_IF_EMPTY _variable) - if("${${_variable}}" STREQUAL "") - set(${_variable} ${ARGN}) - endif() -endmacro() - # athena variables set_if_empty(NUMBER_GHOST_CELLS 3) diff --git a/cmake/parameters.cmake b/cmake/parameters.cmake index f2f43205..a2bfdfb1 100644 --- a/cmake/parameters.cmake +++ b/cmake/parameters.cmake @@ -57,9 +57,3 @@ else() COMPONENTS C REQUIRED) endif() - -if(NOT NBLOCKS OR NOT DEFINED NBLOCKS) - set(USE_NBLOCKS "NOT_USE_NBLOCKS") -else() - set(USE_NBLOCKS "USE_NBLOCKS") -endif() diff --git a/configure.h.in b/configure.h.in index 17d21355..81c93291 100644 --- a/configure.h.in +++ b/configure.h.in @@ -26,6 +26,3 @@ enum { // MPI parallelization (MPI_PARALLEL or NOT_MPI_PARALLEL) #define @MPI_OPTION@ - -// nblocks for exo3 option (USE_NBLOCKS or NOT_USE_NBLOCKS) -#define @USE_NBLOCKS@ diff --git a/examples/2023-Chen-exo3/CMakeLists.txt b/examples/2023-Chen-exo3/CMakeLists.txt index ff3202ca..7e8da1e8 100644 --- a/examples/2023-Chen-exo3/CMakeLists.txt +++ b/examples/2023-Chen-exo3/CMakeLists.txt @@ -15,7 +15,7 @@ else() setup_problem(test_adiabat) setup_problem(hs94) setup_problem(hot_jupiter) - setup_problem(polar_dry) + #setup_problem(polar_dry) endif() # 2. Copy input file to run directory diff --git a/examples/2023-Chen-exo3/earth-dry.yaml b/examples/2023-Chen-exo3/earth-dry.yaml index bf02baf6..b84e499e 100644 --- a/examples/2023-Chen-exo3/earth-dry.yaml +++ b/examples/2023-Chen-exo3/earth-dry.yaml @@ -1,14 +1,15 @@ - -phases: - - name: atm - thermo: ideal-moist - species: [atm] - kinetics: condensation - reactions: none +reference-state: + Tref: 0. + Pref: 1.e5 species: - - name: atm + - name: dry composition: {O: 0.42, N: 1.56, Ar: 0.01} - thermo: - model: constant-cp - cp0: 29.1 J/mol/K + cv_R: 2.5 + +geometry: + type: cartesian + +dynamics: + equation-of-state: + type: ideal-gas diff --git a/examples/2023-Chen-exo3/hot_jupiter.cpp b/examples/2023-Chen-exo3/hot_jupiter.cpp index fcb668b8..7d2b45d4 100644 --- a/examples/2023-Chen-exo3/hot_jupiter.cpp +++ b/examples/2023-Chen-exo3/hot_jupiter.cpp @@ -29,16 +29,13 @@ #include // canoe -#include #include +#include // exo3 #include #include -// snap -#include - #define _sqr(x) ((x) * (x)) #define _qur(x) ((x) * (x) * (x) * (x)) @@ -232,50 +229,24 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) { Rp = pin->GetReal("problem", "Rp"); // construct an isothermal atmosphere - auto pthermo = Thermodynamics::GetInstance(); - auto &w = phydro->w; - std::vector yfrac(1, 1.); + kintera::ThermoX thermo(pimpl->peos->pthermo->options); - for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) { - pthermo->SetMassFractions(yfrac.data()); - pthermo->EquilibrateTP(Ts, p0); - - int i = is; - for (; i <= ie; ++i) { - pthermo->GetPrimitive(w.at(k, j, i)); - pthermo->Extrapolate_inplace(pcoord->dx1f(i), "isothermal", grav, - 0.001); - - // add noise - w(IVY, k, j, i) = 10. * distribution(generator); - w(IVZ, k, j, i) = 10. * distribution(generator); - } - } + auto temp = Ts * torch::ones({ncells3, ncells2}, torch::kDouble); + auto pres = p0 * torch::ones({ncells3, ncells2}, torch::kDouble); + + // half a grid to cell center + pres *= exp(-grav * pcoord->dx1f(is) / (2 * Rd * temp)); + + auto w = get_all(phydro->w); + for (int i = is; i <= ie; ++i) { + w[IPR].select(2, i) = pres; + w[IDN].select(2, i) = pres / (Rd * temp); + pres *= exp(-grav * pcoord->dx1f(i) / (Rd * temp)); + } - // Set up perturbation as the initial condition to break the symmetry of the - // original initial condition. - // Real k3 = 5.; - - // for (int k = ks; k <= ke; ++k) { - // for (int j = js; j <= je; ++j) { - // for (int i = is; i <= ie; ++i) { - // Real theta, phi; - // pexo3->GetLatLon(&theta, &phi, k, j, i); - // Real z = pcoord->x1v(i)-Rp; - // // Set-up isothermal atmosphere - // Real temp = Ts; - // phydro->w(IPR,k,j,i) = p0*exp(-z*grav/(temp*Rd)); - // phydro->w(IDN,k,j,i) = phydro->w(IPR,k,j,i)/ - // (Rd* ( temp + 20.*(distribution(generator) - 0.5) * (1. + - // cos(k3*phi)) ) ); - // //phydro->w(IDN,k,j,i) = phydro->w(IPR,k,j,i)/(Rd*temp); - // phydro->w(IM1,k,j,i) = 0.; - // phydro->w(IM2,k,j,i) = 0.; - // phydro->w(IM3,k,j,i) = 0.; - // } - // } - // } + // add noise + w[IVY] += 10. * torch::randn_like(w[IVY]); + w[IVZ] += 10. * torch::randn_like(w[IVZ]); // transfer to conservative variables // bcc is cell-centered magnetic fields, it is only a place holder here diff --git a/examples/2023-Chen-exo3/hot_jupiter.inp b/examples/2023-Chen-exo3/hot_jupiter.inp index 2bcac0d3..553b3b71 100644 --- a/examples/2023-Chen-exo3/hot_jupiter.inp +++ b/examples/2023-Chen-exo3/hot_jupiter.inp @@ -55,13 +55,13 @@ PressureScaleHeight = 7.56E5 nx1 = 32 -nx2 = 20 -nx3 = 20 +nx2 = 40 +nx3 = 40 grav_acc1 = -8.0 gamma = 1.4 # gamma = C_p/C_v -implicit_flag = 1 +implicit_flag = 9 Rd = 3779. # Gas constant @@ -71,7 +71,7 @@ z_stra = 2.E6 # Stratosphere height, m Gamma_trop = 2.E-4 # Lapse rate, K/m -thermodynamics_config = hot-jupiter.yaml +config_file = earth-dry.yaml Kt = 1.5E5 # Damping timescale, s Ts = 1600. # Surface temperature, K diff --git a/examples/2023-Chen-exo3/hs94.cpp b/examples/2023-Chen-exo3/hs94.cpp index c204b8df..e844d893 100644 --- a/examples/2023-Chen-exo3/hs94.cpp +++ b/examples/2023-Chen-exo3/hs94.cpp @@ -29,16 +29,13 @@ #include // canoe -#include #include +#include // exo3 #include #include -// snap -#include - #define _sqr(x) ((x) * (x)) #define _qur(x) ((x) * (x) * (x) * (x)) @@ -261,39 +258,42 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) { z_iso = pin->GetReal("problem", "z_iso"); // construct an adiabatic atmosphere - auto pthermo = Thermodynamics::GetInstance(); - auto &w = phydro->w; - std::vector yfrac(1, 1.); + kintera::ThermoX thermo(pimpl->peos->pthermo->options); - for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) { - pthermo->SetMassFractions(yfrac.data()); - pthermo->EquilibrateTP(Ts, p0); + auto temp = Ts * torch::ones({ncells3, ncells2}, torch::kDouble); + auto pres = p0 * torch::ones({ncells3, ncells2}, torch::kDouble); + auto xfrac = torch::ones({ncells3, ncells2, 1}, torch::kDouble); - // half a grid to cell center - pthermo->Extrapolate_inplace(pcoord->dx1f(is) / 2., "dry", grav); + // half a grid to cell center + thermo->extrapolate_ad(temp, pres, xfrac, grav, pcoord->dx1f(is) / 2.); - int i = is; - for (; i <= ie; ++i) { - if (pcoord->x1v(i) - Rp > z_iso) break; + auto w = get_all(phydro->w); + int i = is; + for (; i <= ie; ++i) { + auto conc = thermo->compute("TPX->V", {temp, pres, xfrac}); + w[IPR].select(2, i) = pres; + w[IDN].select(2, i) = thermo->compute("V->D", {conc}); + w.slice(0, 1, IVX).select(3, i) = thermo->compute("X->Y", {xfrac}); - pthermo->GetPrimitive(w.at(k, j, i)); - pthermo->Extrapolate_inplace(pcoord->dx1f(i), "dry", grav); + if (pcoord->x1v(i) - Rp > z_iso) break; + thermo->extrapolate_ad(temp, pres, xfrac, grav, pcoord->dx1f(i)); + } - // add noise - w(IVY, k, j, i) = 10. * distribution(generator); - w(IVZ, k, j, i) = 10. * distribution(generator); - } + // isothermal extrapolation + for (; i <= ie; ++i) { + pres *= exp(-grav * pcoord->dx1f(i) / (Rd * temp)); + auto conc = thermo->compute("TPX->V", {temp, pres, xfrac}); + w[IPR].select(2, i) = pres; + w[IDN].select(2, i) = thermo->compute("V->D", {conc}); + w.slice(0, 1, IVX).select(3, i) = thermo->compute("X->Y", {xfrac}); + } - // construct isothermal atmosphere - for (; i <= ie; ++i) { - pthermo->GetPrimitive(w.at(k, j, i)); - pthermo->Extrapolate_inplace(pcoord->dx1f(i), "isothermal", grav); - } - } + // add noise + w[IVY] += 1. * torch::randn_like(w[IVY]); + w[IVZ] += 1. * torch::randn_like(w[IVZ]); - peos->PrimitiveToConserved(w, pfield->bcc, phydro->u, pcoord, is, ie, js, je, - ks, ke); + peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie, + js, je, ks, ke); } void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { diff --git a/examples/2023-Chen-exo3/hs94.inp b/examples/2023-Chen-exo3/hs94.inp index 46094760..b47437f4 100644 --- a/examples/2023-Chen-exo3/hs94.inp +++ b/examples/2023-Chen-exo3/hs94.inp @@ -71,7 +71,7 @@ Rd = 287. cp = 1004. -thermodynamics_config = earth-dry.yaml +config_file = earth-dry.yaml Omega = 7.292E-5 Rp = 6.371E6 diff --git a/examples/2023-Chen-exo3/test_adiabat.cpp b/examples/2023-Chen-exo3/test_adiabat.cpp index 9e349f6b..862b011d 100644 --- a/examples/2023-Chen-exo3/test_adiabat.cpp +++ b/examples/2023-Chen-exo3/test_adiabat.cpp @@ -14,38 +14,32 @@ // canoe #include -#include - -// climath -#include -#include -#include - -// snap -#include -#include +#include Real Ps, Ts, grav; void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { - AllocateUserOutputVariables(3); + AllocateUserOutputVariables(2); SetUserOutputVariableName(0, "temp"); SetUserOutputVariableName(1, "theta"); - SetUserOutputVariableName(2, "mse"); } void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { - auto pthermo = Thermodynamics::GetInstance(); auto &w = phydro->w; + auto temp = get_temp(pimpl->peos, w); + auto pres = get_pres(w); + auto chi = (get_gammad() - 1.) / get_gammad(); + auto theta = temp * pow(Ps / pres, chi); + + auto temp_a = temp.accessor(); + auto theta_a = theta.accessor(); + for (int k = ks; k <= ke; ++k) for (int j = js; j <= je; ++j) for (int i = is; i <= ie; ++i) { - user_out_var(0, k, j, i) = pthermo->GetTemp(w.at(k, j, i)); - user_out_var(1, k, j, i) = potential_temp(pthermo, w.at(k, j, i), Ps); - // msv - user_out_var(2, k, j, i) = - moist_static_energy(pthermo, w.at(k, j, i), grav * pcoord->x1v(i)); + user_out_var(0, k, j, i) = temp_a[k][j][i]; + user_out_var(1, k, j, i) = theta_a[k][j][i]; } } @@ -56,22 +50,25 @@ void Mesh::InitUserMeshData(ParameterInput *pin) { } void MeshBlock::ProblemGenerator(ParameterInput *pin) { - auto pthermo = Thermodynamics::GetInstance(); - auto &w = phydro->w; + kintera::ThermoX thermo(pimpl->peos->pthermo->options); - // construct a reversible adiabat - std::vector yfrac(IVX, 1.); - for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) { - pthermo->SetMassFractions(yfrac.data()); - pthermo->EquilibrateTP(Ts, Ps); + auto temp = Ts * torch::ones({ncells3, ncells2}, torch::kDouble); + auto pres = Ps * torch::ones({ncells3, ncells2}, torch::kDouble); + auto xfrac = torch::ones({ncells3, ncells2, 1}, torch::kDouble); - // half a grid to cell center - pthermo->Extrapolate_inplace(pcoord->dx1f(is) / 2., "reversible", grav); + // half a grid to cell center + thermo->extrapolate_ad(temp, pres, xfrac, grav, pcoord->dx1f(is) / 2.); - for (int i = is; i <= ie; ++i) { - pthermo->GetPrimitive(w.at(k, j, i)); - pthermo->Extrapolate_inplace(pcoord->dx1f(i), "reversible", grav); - } - } + auto w = get_all(phydro->w); + for (int i = is; i <= ie; ++i) { + auto conc = thermo->compute("TPX->V", {temp, pres, xfrac}); + w[IPR].select(2, i) = pres; + w[IDN].select(2, i) = thermo->compute("V->D", {conc}); + w.slice(0, 1, IVX).select(3, i) = thermo->compute("X->Y", {xfrac}); + thermo->extrapolate_ad(temp, pres, xfrac, grav, pcoord->dx1f(i)); + } + + // Change primitive variables to conserved variables + peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie, + js, je, ks, ke); } diff --git a/examples/2023-Chen-exo3/test_adiabat.inp b/examples/2023-Chen-exo3/test_adiabat.inp index 4cdd6ed5..c3f69c83 100644 --- a/examples/2023-Chen-exo3/test_adiabat.inp +++ b/examples/2023-Chen-exo3/test_adiabat.inp @@ -53,10 +53,10 @@ nx3 = 48 grav_acc1 = -9.8 gamma = 1.4 # gamma = C_p/C_v implicit_flag = 1 - - Rd = 287. +config_file = earth-dry.yaml + Ps = 1.E5 Ts = 300. diff --git a/src/exo3/cs_find_block_id.cpp b/src/exo3/cs_find_block_id.cpp index 48ae53bc..123efbbf 100644 --- a/src/exo3/cs_find_block_id.cpp +++ b/src/exo3/cs_find_block_id.cpp @@ -80,7 +80,7 @@ int CubedSphere::FindBlockID(LogicalLocation const &loc) { void CubedSphere::GetLocalIndex(int *lv2_lx2, int *lv2_lx3, int *local_lx2, int *local_lx3, int *bound_lim, LogicalLocation const &loc) { -#ifdef USE_NBLOCKS +#if NBLOCKS > 0 // Updated method, need to manually setup in configure.h, allow 6*n^2 blocks *bound_lim = (int)(sqrt(NBLOCKS / 6) - 0.5); // Find relative location within block diff --git a/src/exo3/cs_transform_ox.cpp b/src/exo3/cs_transform_ox.cpp index 0154f464..af727bfb 100644 --- a/src/exo3/cs_transform_ox.cpp +++ b/src/exo3/cs_transform_ox.cpp @@ -15,8 +15,7 @@ void CubedSphere::TransformOX(int *ox2, int *ox3, int *tox2, int *tox3, int lv2_lx2, lv2_lx3, local_lx2, local_lx3, bound_lim; GetLocalIndex(&lv2_lx2, &lv2_lx3, &local_lx2, &local_lx3, &bound_lim, loc); - - // Hard code the cases... + // Hard code the cases... // The corner cases are not used, abandoned after reading buffers. int target_block = -1; // Block id of target int target_loc_2, target_loc_3; // local x2 and x3 in target block @@ -207,7 +206,7 @@ void CubedSphere::TransformOX(int *ox2, int *ox3, int *tox2, int *tox3, msg << "----------------------------------" << std::endl; ATHENA_ERROR(msg); } -#ifdef USE_NBLOCKS +#if NBLOCKS > 0 lx3_0 = (lx3_0 * (bound_lim + 1)); lx2_0 = (lx2_0 * (bound_lim + 1)); #else diff --git a/src/exo3/cubed_sphere.cpp b/src/exo3/cubed_sphere.cpp index 219f5f9f..e429d977 100644 --- a/src/exo3/cubed_sphere.cpp +++ b/src/exo3/cubed_sphere.cpp @@ -536,7 +536,7 @@ void CubedSphere::sendNeighborBlocks(int ox2, int ox3, int tg_rank, } // Calculate the tag of destination -#ifdef USE_NBLOCKS +#if NBLOCKS > 0 if (tox2 == -1) DirTag = 0 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid; if (tox2 == 1) DirTag = 1 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid; if (tox3 == -1) DirTag = 2 + 4 * pmb->gid + 24 * (bound_lim + 1) * tg_gid; @@ -638,7 +638,7 @@ void CubedSphere::recvNeighborBlocks(int ox2, int ox3, int tg_rank, int dsize = ((kb2 - kb1 + 1) * (jb2 - jb1 + 1) * (ib2 - ib1 + 1) * NWAVE); Real *data = new Real[dsize]; // Calculate the tag for receiving -#ifdef USE_NBLOCKS +#if NBLOCKS > 0 if (ox2 == -1) DirTag = 0 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid; if (ox2 == 1) DirTag = 1 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid; if (ox3 == -1) DirTag = 2 + 4 * tg_gid + 24 * (bound_lim + 1) * pmb->gid; diff --git a/src/exo3/cubed_sphere.hpp b/src/exo3/cubed_sphere.hpp index 6c38f48a..bdc8b2e2 100644 --- a/src/exo3/cubed_sphere.hpp +++ b/src/exo3/cubed_sphere.hpp @@ -7,6 +7,7 @@ // athena #include #include +#include // canoe #include diff --git a/src/impl.cpp b/src/impl.cpp index 1f568df0..b85cb7fe 100644 --- a/src/impl.cpp +++ b/src/impl.cpp @@ -15,6 +15,8 @@ // canoe #include +#include + // snap #include #include @@ -104,6 +106,10 @@ MeshBlock::Impl::Impl(MeshBlock *pmb, ParameterInput *pin) : pmy_block_(pmb) { // planet planet = PlanetFactory::CreateFrom(pmb, pin); + + // radiation data + prad_data = std::make_shared(); + prad_data->cooldown = pin->GetOrAddReal("radiation", "dt", 0.); } MeshBlock::Impl::~Impl() {} diff --git a/src/impl.hpp b/src/impl.hpp index 88b8f935..d6487eda 100644 --- a/src/impl.hpp +++ b/src/impl.hpp @@ -38,6 +38,8 @@ class IdealMoistImpl; class SedHydroImpl; } // namespace snap +struct RadiationData; + //! \class MeshBlock::Impl //! \brief opaque pointer class implements all physics on a MeshBlock class MeshBlock::Impl { @@ -56,6 +58,9 @@ class MeshBlock::Impl { std::shared_ptr planet; std::shared_ptr pexo3; + //// data fields //// + std::shared_ptr prad_data; + public: // constructor and destructor Impl(MeshBlock *pmb, ParameterInput *pin); ~Impl(); diff --git a/src/interface/rad.hpp b/src/interface/rad.hpp new file mode 100644 index 00000000..7a601c1e --- /dev/null +++ b/src/interface/rad.hpp @@ -0,0 +1,16 @@ +// C/C++ +#include +#include + +// torch +#include + +// harp +#include + +struct RadiationData { + double counter = 0.0; + double cooldown = 0.0; + at::Tensor net_flux; + std::map rad_bc; +}; diff --git a/src/snap/riemann/hllc_transform.cpp b/src/snap/riemann/hllc_transform.cpp index d32adb41..e02620d2 100644 --- a/src/snap/riemann/hllc_transform.cpp +++ b/src/snap/riemann/hllc_transform.cpp @@ -12,19 +12,11 @@ #include #include +#include // exo3 #include -// snap -#include - -// microphysics -#include - -// checks -#include - //---------------------------------------------------------------------------------------- //! \fn void Hydro::RiemannSolver //! \brief The HLLC Riemann solver for adiabatic hydrodynamics (use HLLE for @@ -34,9 +26,6 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu, const int ivx, AthenaArray &wl, AthenaArray &wr, AthenaArray &flx, const AthenaArray &dxw) { - auto pthermo = Thermodynamics::GetInstance(); - auto pmicro = pmy_block->pimpl->pmicro; - int ivy = IVX + ((ivx - IVX) + 1) % 3; int ivz = IVX + ((ivx - IVX) + 2) % 3; int dir = ivx - IVX; @@ -64,6 +53,8 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu, } #endif // AFFINE or CUBED_SPHERE + auto peos = pmy_block->pimpl->peos; + #pragma omp simd private(wli, wri, flxi, fl, fr) #pragma distribute_point for (int i = il; i <= iu; ++i) { @@ -87,11 +78,11 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu, // correction for gamma // left - Real gammal = pthermo->GetGamma(wli); + Real gammal = get_gamma(peos, wli); Real kappal = 1. / (gammal - 1.); // right - Real gammar = pthermo->GetGamma(wri); + Real gammar = get_gamma(peos, wri); Real kappar = 1. / (gammar - 1.); //--- Step 2. Compute middle state estimates with PVRS (Toro 10.5.2) @@ -239,21 +230,6 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu, flx(ivy, k, j, i) = flxi[IVY]; flx(ivz, k, j, i) = flxi[IVZ]; flx(IEN, k, j, i) = flxi[IEN]; - - // sedimentation flux - if (ivx == IVX) { - if (i == iu) return; - Real hr = (wri[IPR] * (kappar + 1.) + KE_r) / wri[IDN]; - for (int n = 0; n < NCLOUD + NPRECIP; ++n) { - Real rho = wri[IDN] * wri[1 + NVAPOR + n]; - auto vsed = pmicro->vsedf[0](n, k, j, i); - flx(1 + NVAPOR + n, k, j, i) += vsed * rho; - flx(ivx, k, j, i) += vsed * rho * wri[ivx]; - flx(ivy, k, j, i) += vsed * rho * wri[ivy]; - flx(ivz, k, j, i) += vsed * rho * wri[ivz]; - flx(IEN, k, j, i) += vsed * rho * hr; - } - } } #if defined(AFFINE) || defined(CUBED_SPHERE) // need of deprojection @@ -294,6 +270,4 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu, flx(IVZ, k, j, i) = tz + ty * cth; } #endif // CUBED_SPHERE - - check_hydro_riemann_solver_flux(flx, ivx, k, j, il, iu); } diff --git a/src/tasklist/implicit_hydro_tasks.cpp b/src/tasklist/implicit_hydro_tasks.cpp index 2b0b91c4..c6d6ba38 100644 --- a/src/tasklist/implicit_hydro_tasks.cpp +++ b/src/tasklist/implicit_hydro_tasks.cpp @@ -22,6 +22,7 @@ #include #include +#include // snap #include @@ -219,24 +220,45 @@ TaskStatus ImplicitHydroTasks::AddFluxToConserved(MeshBlock *pmb, int stage) { int ie = pmb->ie, je = pmb->je, ke = pmb->ke; auto prad = pmb->pimpl->prad; + auto peos = pmb->pimpl->peos; + + auto w = get_all(pmb->phydro->w); + auto temp = peos->compute("W->T", {w}); + auto prad_data = pmb->pimpl->prad_data; if (stage <= nstages) { - /*if (stage_wghts[stage - 1].main_stage) { - for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) { - prad->AddRadiativeFlux(pmb->phydro, k, j, is, ie + 1); - } + if (stage_wghts[stage - 1].main_stage) { + if (prad_data->cooldown > 0) { + auto net_flux_a = prad_data->net_flux.accessor(); + // add net-flux + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie + 1; ++i) { + pmb->phydro->flux[X1DIR](IEN, k, j, i) += net_flux_a[k][j][i]; + } + } if (stage == nstages) { // last stage - if (prad->GetCounter() > 0.) { // reuse previous flux - prad->DecrementCounter(pmb->pmy_mesh->dt); - } else { // update radiative flux - prad->ResetCounter(); - for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) prad->CalFlux(pmb, k, j); + if (prad_data->counter > 0.) { // reuse previous flux + prad_data->counter -= pmb->pmy_mesh->dt; + } else if (prad_data->cooldown > 0) { // update radiative flux + std::map atm; + int nlyr = w.size(3); + atm["pres"] = w[IPR].narrow(0, NGHOST, nlyr).view({-1, nlyr}); + atm["temp"] = temp.narrow(0, NGHOST, nlyr).view({-1, nlyr}); + atm["rho"] = w[IDN].narrow(0, NGHOST, nlyr).view({-1, nlyr}); + + auto ivol = + peos->pthermo->compute("DY->V", {w[IDN], w.slice(0, 1, IVX)}); + auto conc = ivol * peos->pthermo->inv_mu; + + prad_data->counter = prad_data->cooldown; + auto [netflux, dnflux, upflux] = + prad->forward(conc, peos->pcoord->dx1f, &prad_data->rad_bc, &atm); + prad_data->net_flux = netflux; } } - }*/ + } return TaskStatus::next; } return TaskStatus::fail;