diff --git a/README b/README index f05caadb..d59c7683 100644 --- a/README +++ b/README @@ -4,6 +4,13 @@ when vapor pressure are too low (~1e-10 Pa) at low temperatuers, the saturationa when vapor pressures are too high (1e5 Pa) eg. for CO2, the saturation adjustment fails immediately (all nan) -exo can use cores 24/96/384/1536 -also permit 216/600/864/1176/ if specify NBLOCKS in the .cmake file +exo can use cores 24/96/384/1536: 6*4^n +also permit 216/600/864/1176/1536/1944/2400... (24*n^2) if specify NBLOCKS in the .cmake file +That is, if meshblock is x, the mesh should be x*4*n, and x*6*n, for example, if x is 10, mesh could be 40/60, or 240/360 etc. +meshblock = 5 is not very efficient, maybe meshblock is 8 or higher + +gammal in svp is delta[n] = (cp_ratio_mass[n] - cp_ratio_mass[i]) * mu_ratio[i] / (1. - 1. / gammad); + +cp_ratio_mass is rcp +mu_ratio is eps diff --git a/cmake/examples/ice_h2s.cmake b/cmake/examples/ice_h2s.cmake new file mode 100644 index 00000000..b0fd9320 --- /dev/null +++ b/cmake/examples/ice_h2s.cmake @@ -0,0 +1,20 @@ +# configure file for test jupiter crm + +macro(SET_IF_EMPTY _variable) + if("${${_variable}}" STREQUAL "") + set(${_variable} ${ARGN}) + endif() +endmacro() + +# athena variables +set_if_empty(NUMBER_GHOST_CELLS 3) + +# canoe configure +set(NVAPOR 1) +set(NCLOUD 2) +set(NPHASE_LEGACY 3) +set(PNETCDF ON) +set(MPI ON) +set(TASKLIST ImplicitHydroTasks) +#set_if_empty(RSOLVER hllc_transform) +set(RSOLVER lmars) diff --git a/cmake/examples/minichem.cmake b/cmake/examples/minichem.cmake index 93a3e5a8..f4895240 100644 --- a/cmake/examples/minichem.cmake +++ b/cmake/examples/minichem.cmake @@ -10,10 +10,8 @@ endmacro() set_if_empty(NUMBER_GHOST_CELLS 3) set(MINICHEM ON) -set(NCHEMISTRY 13) -#set(NVAPOR 1) -#set(NCLOUD 2) -#set(NPHASE_LEGACY 3) +set(NCHEM 13) +# set(NVAPOR 1) set(NCLOUD 2) set(NPHASE_LEGACY 3) set(NETCDF ON) set(PNETCDF ON) set(MPI ON) diff --git a/cmake/examples/rtcldexo.cmake b/cmake/examples/rtcldexo.cmake index 5f9adc11..61423a5c 100644 --- a/cmake/examples/rtcldexo.cmake +++ b/cmake/examples/rtcldexo.cmake @@ -9,14 +9,14 @@ endmacro() # athena variables set_if_empty(NUMBER_GHOST_CELLS 3) -set(CUBED_SPHERE ON) -#set(NBLOCKS 600) +#set(CUBED_SPHERE ON) +#set(NBLOCKS 2400) #set(COORDINATE_SYSTEM "spherical_polar") -set(NVAPOR 1) -set(NCLOUD 2) +set(NVAPOR 2) +set(NCLOUD 4) set(NPHASE_LEGACY 3) set(NETCDF ON) -#set(PNETCDF OFF) +set(PNETCDF ON) set(MPI ON) set(DISORT ON) set(PYTHON_BINDINGS ON) diff --git a/cmake/examples/rtdry.cmake b/cmake/examples/rtdry.cmake index 3c039389..f29191fd 100644 --- a/cmake/examples/rtdry.cmake +++ b/cmake/examples/rtdry.cmake @@ -3,11 +3,8 @@ # athena variables set_if_empty(NUMBER_GHOST_CELLS 3) -set(CUBED_SPHERE ON) -set(NBLOCKS 864) -# set(COORDINATE_SYSTEM "spherical_polar") -set(NETCDF OFF) -set(PNETCDF ON) +set(NETCDF ON) +set(PNETCDF OFF) set(MPI ON) set(DISORT ON) set(PYTHON_BINDINGS ON) diff --git a/cmake/minichem.cmake b/cmake/minichem.cmake index 9850ecc2..c22f3fe4 100644 --- a/cmake/minichem.cmake +++ b/cmake/minichem.cmake @@ -7,4 +7,4 @@ set(patch_command set(PACKAGE_NAME minichem) set(REPO_URL "https://github.com/chengcli/mini_chem") set(REPO_TAG "3372a500f038f83228c9f8f944b3fb6b2dedc572") -add_package(${PACKAGE_NAME} ${REPO_URL} ${REPO_TAG} "${patch_commend}" OFF) +add_package(${PACKAGE_NAME} ${REPO_URL} ${REPO_TAG} "${patch_command}" OFF) diff --git a/cmake/modules/FindKintera.cmake b/cmake/modules/FindKintera.cmake index 18efb724..5abc6ef1 100644 --- a/cmake/modules/FindKintera.cmake +++ b/cmake/modules/FindKintera.cmake @@ -4,8 +4,6 @@ # # KINTERA_LIBRARY: Link these libraries when using KINTERA # -# VAPORS_LIBRARY: Link these libraries when using KINTERA -# # KINTERA_FOUND: True if kintera found # # Normal usage would be: @@ -14,7 +12,7 @@ # # include_directories(${KINTERA_INCLUDE_DIR}) # -# target_link_libraries(${KINTERA_LIBRARY} ${VAPORS_LIBRARY}) +# target_link_libraries(${KINTERA_LIBRARY}) include(FindPackageHandleStandardArgs) @@ -102,22 +100,11 @@ find_library( $ENV{KINTERA_DIR}/lib $ENV{KINTERA_ROOT}/lib) -find_library( - VAPORS_LIBRARY vapors_release - HINTS ${kintera_lib_dir} - /opt/homebrew/lib - /usr/lib/x86_64-linux-gnu/ - KINTERA_DIR/lib - KINTERA_LIB - $ENV{KINTERA_LIB} - $ENV{KINTERA_DIR}/lib - $ENV{KINTERA_ROOT}/lib) - if(${CUDAToolKit_FOUND}) - set(kintera_required_vars VAPORS_LIBRARY KINTERA_LIBRARY KINTERA_CUDA_LIBRARY + set(kintera_required_vars KINTERA_LIBRARY KINTERA_CUDA_LIBRARY KINTERA_INCLUDE_DIR) else() - set(kintera_required_vars VAPORS_LIBRARY KINTERA_LIBRARY KINTERA_INCLUDE_DIR) + set(kintera_required_vars KINTERA_LIBRARY KINTERA_INCLUDE_DIR) endif() mark_as_advanced(${kintera_required_vars}) diff --git a/cmake/modules/FindSnap.cmake b/cmake/modules/FindSnap.cmake index e153d073..2afd8961 100644 --- a/cmake/modules/FindSnap.cmake +++ b/cmake/modules/FindSnap.cmake @@ -4,8 +4,6 @@ # # SNAP_LIBRARY: Link these libraries when using SNAP # -# BC_LIBRARY: Link these libraries when using SNAP -# # SNAP_FOUND: True if snap found # # Normal usage would be: @@ -14,7 +12,7 @@ # # include_directories(${SNAP_INCLUDE_DIR}) # -# target_link_libraries(${SNAP_LIBRARY} ${BC_LIBRARY}) +# target_link_libraries(${SNAP_LIBRARY}) include(FindPackageHandleStandardArgs) @@ -100,22 +98,10 @@ find_library( $ENV{SNAP_DIR}/lib $ENV{SNAP_ROOT}/lib) -find_library( - BC_LIBRARY bc_release - HINTS ${snap_lib_dir} - /opt/homebrew/lib - /usr/lib/x86_64-linux-gnu/ - SNAP_DIR/lib - SNAP_LIB - $ENV{SNAP_LIB} - $ENV{SNAP_DIR}/lib - $ENV{SNAP_ROOT}/lib) - if(${CUDAToolKit_FOUND}) - set(snap_required_vars BC_LIBRARY SNAP_LIBRARY SNAP_CUDA_LIBRARY - SNAP_INCLUDE_DIR) + set(snap_required_vars SNAP_LIBRARY SNAP_CUDA_LIBRARY SNAP_INCLUDE_DIR) else() - set(snap_required_vars BC_LIBRARY SNAP_LIBRARY SNAP_INCLUDE_DIR) + set(snap_required_vars SNAP_LIBRARY SNAP_INCLUDE_DIR) endif() mark_as_advanced(${snap_required_vars}) diff --git a/examples/2024-XZhang-cloud-rt/CMakeLists.txt b/examples/2024-XZhang-cloud-rt/CMakeLists.txt index 2142c0c7..9844e3d5 100644 --- a/examples/2024-XZhang-cloud-rt/CMakeLists.txt +++ b/examples/2024-XZhang-cloud-rt/CMakeLists.txt @@ -1,11 +1,13 @@ #if (${NVAPOR} EQUAL 1 AND ${NCLOUD} EQUAL 2 AND ${NPHASE_LEGACY} EQUAL 3) if (${NVAPOR} EQUAL 1) #setup_problem(rt_silicate) - setup_problem(beta_rtsi) + setup_problem(ice_h2s) + #setup_problem(beta_rtsi) #setup_problem(rt_h2o) elseif (${NVAPOR} EQUAL 2) #setup_problem(rt_h2o_nh3) - setup_problem(jupiter_crm) + #setup_problem(jupiter_crm) + setup_problem(rtcldu3d) elseif (${NVAPOR} EQUAL 3) setup_problem(beta_rtjup3cld) else() diff --git a/examples/2024-XZhang-cloud-rt/ice_h2s.cpp b/examples/2024-XZhang-cloud-rt/ice_h2s.cpp new file mode 100644 index 00000000..8514ac13 --- /dev/null +++ b/examples/2024-XZhang-cloud-rt/ice_h2s.cpp @@ -0,0 +1,195 @@ +/* ------------------------------------------------------------------------------------- + * SNAP Example Program + * + * Contributer: + * Cheng Li, University of Michigan + * + * Year: 2023 + * Contact: chengcli@umich.edu + * Reference: Test Jupiter CRM + * ------------------------------------------------------------------------------------- + */ + +// athena +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// application +#include +#include + +// canoe +#include +#include + +// climath +#include +#include + +// snap +#include + +// special includes +#include "svp_functions.hpp" + +Real grav, P0, T0, Tmin, prad, hrate; +int iH2S; + +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { + AllocateUserOutputVariables(4 + NVAPOR); + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + SetUserOutputVariableName(2, "thetav"); + SetUserOutputVariableName(3, "mse"); + SetUserOutputVariableName(4, "rh_H2S"); +} + +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + 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(this, k, j, i); + user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); + // theta_v + user_out_var(2, k, j, i) = + user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); + // mse + user_out_var(3, k, j, i) = + pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i); + // relative humidity + user_out_var(4, k, j, i) = pthermo->RelativeHumidity(this, 1, k, j, i); + } +} + +void Forcing(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, AthenaArray const &r, + AthenaArray const &bcc, AthenaArray &du, + AthenaArray &s) { + int is = pmb->is, js = pmb->js, ks = pmb->ks; + int ie = pmb->ie, je = pmb->je, ke = pmb->ke; + auto pthermo = Thermodynamics::GetInstance(); + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); + + air.ToMoleFraction(); + Real cv = + get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass()); + +// if (w(IPR, k, j, i) < prad) { +// du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv * +// (1. + 1.E-4 * sin(2. * M_PI * rand() / RAND_MAX)); +// } + } +} + +void Mesh::InitUserMeshData(ParameterInput *pin) { + grav = -pin->GetReal("hydro", "grav_acc1"); + + P0 = pin->GetReal("problem", "P0"); + T0 = pin->GetReal("problem", "T0"); + + Tmin = pin->GetReal("problem", "Tmin"); +// prad = pin->GetReal("problem", "prad"); +// hrate = pin->GetReal("problem", "hrate") / 86400.; + + // index + auto pindex = IndexMap::GetInstance(); + iH2S = pindex->GetVaporId("H2S"); + EnrollUserExplicitSourceFunction(Forcing); +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + srand(Globals::my_rank + time(0)); + + Application::Logger app("main"); + app->Log("ProblemGenerator: jupiter_crm"); + + auto pthermo = Thermodynamics::GetInstance(); + + // mesh limits + Real x1min = pmy_mesh->mesh_size.x1min; + Real x1max = pmy_mesh->mesh_size.x1max; + + // request temperature and pressure + app->Log("request T", T0); + app->Log("request P", P0); + + // thermodynamic constants + Real gamma = pin->GetReal("hydro", "gamma"); + Real Rd = pthermo->GetRd(); + Real cp = gamma / (gamma - 1.) * Rd; + + // set up an adiabatic atmosphere + int max_iter = 400, iter = 0; + Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); + + AirParcel air(AirParcel::Type::MoleFrac); + + // estimate surface temperature and pressure + Real Ts = T0 - grav / cp * x1min; + Real Ps = P0 * pow(Ts / T0, cp / Rd); + Real xH2S = pin->GetReal("problem", "qH2S.ppmv") / 1.E6; + + while (iter++ < max_iter) { + // read in vapors + air.w[iH2S] = xH2S; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // stop at just above P0 + for (int i = is; i <= ie; ++i) { + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + if (air.w[IPR] < P0) break; + } + + // make up for the difference + Ts += T0 - air.w[IDN]; + if (std::abs(T0 - air.w[IDN]) < Ttol) break; + + app->Log("Iteration #", iter); + app->Log("T", air.w[IDN]); + } + + if (iter > max_iter) { + throw RuntimeError("ProblemGenerator", "maximum iteration reached"); + } + + // construct atmosphere from bottom up + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + air.SetZero(); + air.w[iH2S] = xH2S; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); + + int i = is; + for (; i <= ie; ++i) { + if (air.w[IDN] < Tmin) break; + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + //pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + } + + // Replace adiabatic atmosphere with isothermal atmosphere if temperature + // is too low + for (; i <= ie; ++i) { + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + } + } +} diff --git a/examples/2024-XZhang-cloud-rt/jupiter_crm.cpp b/examples/2024-XZhang-cloud-rt/jupiter_crm.cpp new file mode 100644 index 00000000..020a077e --- /dev/null +++ b/examples/2024-XZhang-cloud-rt/jupiter_crm.cpp @@ -0,0 +1,279 @@ +//======================================================================================== +// Athena++ astrophysical MHD code +// Copyright(C) 2014 James M. Stone and other code +// contributors Licensed under the 3-clause BSD License, see LICENSE file for +// details +//======================================================================================== +//! \file hs94.cpp +// \brief Problem generator for Held-Suarez-94 GCM bench mark. +// +// REFERENCE: I.M Held & M.J Suarez, "A Proposal for the Intercomparison of the +// Dynamical Cores of Atmospheric General Circulation Models" + +// C++ headers +#include +#include +#include +#include +#include + +// athena +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// canoe +#include +#include +#include +#include + +// exo3 +#include +#include + +// snap +#include + +// application +#include +#include + +// climath +#include +#include + +// special includes +#include + + +#define _sqr(x) ((x) * (x)) +#define _qur(x) ((x) * (x) * (x) * (x)) + +using namespace std; + +Real grav, P0, T0, Tmin, prad, hrate, Omega, Rp; +int iH2O, iNH3; + +std::default_random_engine generator; +std::normal_distribution distribution(0.0, 1.0); + +namespace cs = CubedSphereUtility; + +void Forcing(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, AthenaArray const &r, + AthenaArray const &bcc, AthenaArray &du, + AthenaArray &s) { + int is = pmb->is, js = pmb->js, ks = pmb->ks; + int ie = pmb->ie, je = pmb->je, ke = pmb->ke; + auto pthermo = Thermodynamics::GetInstance(); + + auto pexo3 = pmb->pimpl->pexo3; + Real om_earth = Omega; + + 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) { + } + + + /* Sponge Layer + 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) { + Real pres = w(IPR, k, j, i); + if (pres < piso) { // sponge layer at top + Real tau = sponge_tau * pow(pres / piso, 2); + u(IVX, k, j, i) = u(IVX, k, j, i) / (1 + dt / tau); + u(IVY, k, j, i) = u(IVY, k, j, i) / (1 + dt / tau); + u(IVZ, k, j, i) = u(IVZ, k, j, i) / (1 + dt / tau); + } + } + } + }*/ + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + // coriolis force + Real lat, lon; + pexo3->GetLatLon(&lat, &lon, k, j, i); + // coriolis force + Real f = 2. * om_earth * sin(lat); + Real f2 = 2. * om_earth * cos(lat); + Real U, V; + + pexo3->GetUV(&U, &V, w(IVY, k, j, i), w(IVZ, k, j, i), k, j, i); + + Real m1 = w(IDN, k, j, i) * w(IVX, k, j, i); + Real m2 = w(IDN, k, j, i) * U; + Real m3 = w(IDN, k, j, i) * V; + + // du(IM1, k, j, i) += dt * f * m2; + Real ll_acc_U = f * m3; //- f2 * m1; + Real ll_acc_V = -f * m2; + Real acc1, acc2, acc3; + pexo3->GetVyVz(&acc2, &acc3, ll_acc_U, ll_acc_V, k, j, i); + pexo3->ContravariantVectorToCovariant(j, k, acc2, acc3, &acc2, &acc3); + du(IM2, k, j, i) += dt * acc2; + du(IM3, k, j, i) += dt * acc3; + } +} + +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { + AllocateUserOutputVariables(8 + NVAPOR); + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + SetUserOutputVariableName(2, "thetav"); + SetUserOutputVariableName(3, "mse"); + SetUserOutputVariableName(4, "lat"); + SetUserOutputVariableName(5, "lon"); + SetUserOutputVariableName(6, "U"); + SetUserOutputVariableName(7, "V"); + for (int n = 1; n <= NVAPOR; ++n) { + std::string name = "rh" + std::to_string(n); + SetUserOutputVariableName(7 + n, name.c_str()); + } +} + +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + 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(this, k, j, i); + user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); + // theta_v + user_out_var(2, k, j, i) = + user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); + // mse + user_out_var(3, k, j, i) = + pthermo->MoistStaticEnergy(this, grav * (pcoord->x1v(i)-Rp), k, j, i); + for (int n = 1; n <= NVAPOR; ++n) + user_out_var(7 + n, k, j, i) = + pthermo->RelativeHumidity(this, n, k, j, i); + Real x = tan(pcoord->x2v(j)); + + Real Vy = phydro->w(IVY, k, j, i); + Real Vz = phydro->w(IVZ, k, j, i); + Real lat, lon; + pimpl->pexo3->GetLatLon(&lat, &lon, k, j, i); + user_out_var(4, k, j, i) = lat / PI * 180.0; + user_out_var(5, k, j, i) = lon / PI * 180.0; + Real U, V; + pimpl->pexo3->GetUV(&U, &V, Vy, Vz, k, j, i); + user_out_var(6, k, j, i) = U; + user_out_var(7, k, j, i) = V; + } +} + +void Mesh::InitUserMeshData(ParameterInput *pin) { + grav = -pin->GetReal("hydro", "grav_acc1"); + + P0 = pin->GetReal("problem", "P0"); + T0 = pin->GetReal("problem", "T0"); + Omega = pin->GetReal("problem", "Omega"); + Rp = pin->GetReal("problem", "Rp"); + + Tmin = pin->GetReal("problem", "Tmin"); + prad = pin->GetReal("problem", "prad"); + hrate = pin->GetReal("problem", "hrate") / 86400.; + + // index + auto pindex = IndexMap::GetInstance(); + iH2O = pindex->GetVaporId("H2O"); + iNH3 = pindex->GetVaporId("NH3"); + EnrollUserExplicitSourceFunction(Forcing); +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + srand(Globals::my_rank + time(0)); + + Application::Logger app("main"); + app->Log("ProblemGenerator: jupiter_crm"); + + auto pthermo = Thermodynamics::GetInstance(); + + // mesh limits + Real x1min = pmy_mesh->mesh_size.x1min; + Real x1max = pmy_mesh->mesh_size.x1max; + + // request temperature and pressure + app->Log("request T", T0); + app->Log("request P", P0); + + // thermodynamic constants + Real gamma = pin->GetReal("hydro", "gamma"); + Real Rd = pthermo->GetRd(); + Real cp = gamma / (gamma - 1.) * Rd; + + // set up an adiabatic atmosphere + int max_iter = 400, iter = 0; + Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); + + AirParcel air(AirParcel::Type::MoleFrac); + + // estimate surface temperature and pressure + Real Ts = T0 - grav / cp * (x1min - Rp); + Real Ps = P0 * pow(Ts / T0, cp / Rd); + Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6; + Real xNH3 = pin->GetReal("problem", "qNH3.ppmv") / 1.E6; + + while (iter++ < max_iter) { + // read in vapors + air.w[iH2O] = xH2O; + air.w[iNH3] = xNH3; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // stop at just above P0 + for (int i = is; i <= ie; ++i) { + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + if (air.w[IPR] < P0) break; + } + + // make up for the difference + Ts += T0 - air.w[IDN]; + if (std::abs(T0 - air.w[IDN]) < Ttol) break; + + app->Log("Iteration #", iter); + app->Log("T", air.w[IDN]); + } + + if (iter > max_iter) { + throw RuntimeError("ProblemGenerator", "maximum iteration reached"); + } + + // construct atmosphere from bottom up + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + air.SetZero(); + air.w[iH2O] = xH2O; + air.w[iNH3] = xNH3; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); + + int i = is; + for (; i <= ie; ++i) { + if (air.w[IDN] < Tmin) break; + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5); + } + + // Replace adiabatic atmosphere with isothermal atmosphere if temperature + // is too low + for (; i <= ie; ++i) { + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + } + } +} diff --git a/examples/2024-XZhang-cloud-rt/rtcldu3d.cpp b/examples/2024-XZhang-cloud-rt/rtcldu3d.cpp new file mode 100644 index 00000000..2c102735 --- /dev/null +++ b/examples/2024-XZhang-cloud-rt/rtcldu3d.cpp @@ -0,0 +1,288 @@ +// C++ headers +#include +#include +#include +#include +#include + +// athena +#include +#include +#include +#include +#include +#include +#include +#include + +// canoe +#include +#include +#include +#include + +// exo3 +#include +#include + +// snap +#include + +// application +#include +#include + +// climath +#include +#include + +// special includes +#include "svp_functions.hpp" + +// astro +#include + +// harp +#include "harp/radiation.hpp" + +#define _sqr(x) ((x) * (x)) +#define _qur(x) ((x) * (x) * (x) * (x)) + +using namespace std; + +Real grav, P0, T0, Tmin, Omega, radius; +Real Omega_f, center_deg, beta_f; +int iH2S, iCH4; + +std::default_random_engine generator; +std::normal_distribution distribution(0.0, 1.0); + +namespace cs = CubedSphereUtility; + +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { +#ifdef CUBED_SPHERE + AllocateUserOutputVariables(11 + NVAPOR); +#else + AllocateUserOutputVariables(5 + NVAPOR); +#endif + + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + SetUserOutputVariableName(2, "thetav"); + SetUserOutputVariableName(3, "mse"); + SetUserOutputVariableName(4, "pres"); + for (int n = 1; n <= NVAPOR; ++n) { + std::string name = "rh" + std::to_string(n); + SetUserOutputVariableName(4 + n, name.c_str()); + } + +#ifdef CUBED_SPHERE + SetUserOutputVariableName(5 + NVAPOR, "lat"); + SetUserOutputVariableName(6 + NVAPOR, "lon"); + SetUserOutputVariableName(7 + NVAPOR, "vlat"); + SetUserOutputVariableName(8 + NVAPOR, "vlon"); + SetUserOutputVariableName(9 + NVAPOR, "w"); + SetUserOutputVariableName(10 + NVAPOR, "zenith"); +#endif +} + +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + 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(this, k, j, i); + user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); + // theta_v + user_out_var(2, k, j, i) = + user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); + // mse + user_out_var(3, k, j, i) = pthermo->MoistStaticEnergy( + this, grav * (pcoord->x1v(i) - radius), k, j, i); + user_out_var(4, k, j, i) = phydro->w(IPR, k, j, i); + for (int n = 1; n <= NVAPOR; ++n) + user_out_var(4 + n, k, j, i) = + pthermo->RelativeHumidity(this, n, k, j, i); + } + +#ifdef CUBED_SPHERE + auto pexo3 = pimpl->pexo3; + Real lat, lon; + Real U, V; + Direction ray = pimpl->prad->GetRayInput(0); + Real zenith; + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + pexo3->GetLatLon(&lat, &lon, k, j, i); + pexo3->GetUV(&U, &V, phydro->w(IVY, k, j, i), phydro->w(IVZ, k, j, i), + k, j, i); + user_out_var(5 + NVAPOR, k, j, i) = lat; + user_out_var(6 + NVAPOR, k, j, i) = lon; + user_out_var(7 + NVAPOR, k, j, i) = U; + user_out_var(8 + NVAPOR, k, j, i) = V; + user_out_var(9 + NVAPOR, k, j, i) = phydro->w(IVX, k, j, i); + + ray = pimpl->planet->ParentZenithAngle(pmy_mesh->time, M_PI / 2. - lat, + lon); + zenith = std::acos(ray.mu) / M_PI * 180.0; + user_out_var(10 + NVAPOR, k, j, i) = zenith; + } +#endif +} + +void Forcing(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, AthenaArray const &r, + AthenaArray const &bcc, AthenaArray &du, + AthenaArray &s) { + int is = pmb->is, js = pmb->js, ks = pmb->ks; + int ie = pmb->ie, je = pmb->je, ke = pmb->ke; + auto pthermo = Thermodynamics::GetInstance(); + +#ifndef CUBED_SPHERE + 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) { + Real dx2 = pmb->pcoord->dx2v(j); + Real dx3 = pmb->pcoord->dx3v(k); + // Coriolis force + Real x2 = pmb->pcoord->x2v(j); + // Load Latitude (constant f here) + Real theta0 = center_deg * M_PI / 180.; + // Coriolis force (shallow approximation) + Real local_f = 2. * sin(theta0) * Omega_f + x2 * beta_f; + + Real m2 = w(IDN, k, j, i) * w(IVY, k, j, i); + Real m3 = w(IDN, k, j, i) * w(IVZ, k, j, i); + du(IM2, k, j, i) -= dt * local_f * m3; + du(IM3, k, j, i) += dt * local_f * m2; + } +#endif + +#ifdef CUBED_SPHERE + auto pexo3 = pmb->pimpl->pexo3; + 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) { + Real lat, lon; + pexo3->GetLatLon(&lat, &lon, k, j, i); + + // coriolis force + Real f = 2. * Omega * sin(lat); + Real f2 = 2. * Omega * cos(lat); + Real U, V; + + pexo3->GetUV(&U, &V, w(IVY, k, j, i), w(IVZ, k, j, i), k, j, i); + + Real m1 = w(IDN, k, j, i) * w(IVX, k, j, i); + Real m2 = w(IDN, k, j, i) * U; + Real m3 = w(IDN, k, j, i) * V; + + Real ll_acc_U = f * m3; + Real ll_acc_V = -f * m2; + Real acc2, acc3; + pexo3->GetVyVz(&acc2, &acc3, ll_acc_U, ll_acc_V, k, j, i); + pexo3->ContravariantVectorToCovariant(j, k, acc2, acc3, &acc2, &acc3); + du(IM2, k, j, i) += dt * acc2; + du(IM3, k, j, i) += dt * acc3; + } +#endif +} + +void Mesh::InitUserMeshData(ParameterInput *pin) { + grav = -pin->GetReal("hydro", "grav_acc1"); + P0 = pin->GetReal("problem", "P0"); + T0 = pin->GetReal("problem", "T0"); + Omega = pin->GetReal("problem", "Omega"); + radius = pin->GetReal("problem", "radius"); + Tmin = pin->GetReal("problem", "Tmin"); + center_deg = pin->GetReal("problem", "center_deg"); + beta_f = pin->GetReal("problem", "beta_f"); + Omega_f = pin->GetReal("problem", "Omega_f"); + + // index + auto pindex = IndexMap::GetInstance(); + iH2S = pindex->GetVaporId("H2S"); + iCH4 = pindex->GetVaporId("CH4"); + EnrollUserExplicitSourceFunction(Forcing); +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + srand(Globals::my_rank + time(0)); + auto pthermo = Thermodynamics::GetInstance(); + AirParcel air(AirParcel::Type::MoleFrac); + + // mesh limits + Real x1min = pmy_mesh->mesh_size.x1min; + Real x1max = pmy_mesh->mesh_size.x1max; + + // thermodynamic constants + Real gammad = pin->GetReal("hydro", "gamma"); + Real Rd = pthermo->GetRd(); + Real cp = gammad / (gammad - 1.) * Rd; + + // set up an adiabatic atmosphere + int max_iter = 400, iter = 0; + Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); + + // estimate surface temperature and pressure + Real Ts = T0 - grav / cp * (x1min - radius); + Real Ps = P0 * pow(Ts / T0, cp / Rd); + Real xH2S = pin->GetReal("problem", "qH2S.ppmv") / 1.E6; + Real xCH4 = pin->GetReal("problem", "qCH4.ppmv") / 1.E6; + + while (iter++ < max_iter) { + // read in vapors + air.w[iH2S] = xH2S; + air.w[iCH4] = xCH4; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // stop at just above P0 + for (int i = is; i <= ie; ++i) { + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + //pthermo->Extrapolate(&air, pcoord->dx1f(i), "reversible", grav); + if (air.w[IPR] < P0) break; + } + + // make up for the difference + Ts += T0 - air.w[IDN]; + if (std::abs(T0 - air.w[IDN]) < Ttol) break; + } + + if (iter > max_iter) { + throw RuntimeError("ProblemGenerator", "maximum iteration reached"); + } + + // construct atmosphere from bottom up + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + air.SetZero(); + air.w[iH2S] = xH2S; + air.w[iCH4] = xCH4; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); + + int i = is; + for (; i <= ie; ++i) { + air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); + if (air.w[IDN] < Tmin) break; + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + //pthermo->Extrapolate(&air, pcoord->dx1f(i), "reversible", grav); + } + + // Replace adiabatic atmosphere with isothermal atmosphere if temperature + // is too low + for (; i <= ie; ++i) { + air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + } + } +} diff --git a/examples/2024-XZhang-cloud-rt/rtdry.cpp b/examples/2024-XZhang-cloud-rt/rtdry.cpp index 2b29b303..495d8af1 100644 --- a/examples/2024-XZhang-cloud-rt/rtdry.cpp +++ b/examples/2024-XZhang-cloud-rt/rtdry.cpp @@ -203,9 +203,9 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) { for (; i <= ie; ++i) { air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); AirParcelHelper::distribute_to_conserved(this, k, j, i, air); - // pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", - // grav); - pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav, 1.e-3); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", + grav); + //pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav, 1.e-4); } } } diff --git a/examples/2024-XZhang-cloud-rt/rtdry.yaml b/examples/2024-XZhang-cloud-rt/rtdry.yaml new file mode 100644 index 00000000..7f5f5bbc --- /dev/null +++ b/examples/2024-XZhang-cloud-rt/rtdry.yaml @@ -0,0 +1,95 @@ +opacity-sources: + - name: H2-vis + class: FreedmanMean + parameters: {scale: 0.5, met: 1.7} + + - name: H2-ir + class: FreedmanMean + parameters: {scale: 1., met: 1.7} + + - name: premix + class: AbsorberCK + data: ../data/ExoFMS/output.nc + + - name: heliosmix + class: HeliosCK + data: ../data/ck_data_01242024/ck/PM_ck_HELIOSK_cond_11_nOPT_wcia.txt + parameters: {band_id: 1} + +bands: [vistck] +#bands: [ird, irt] + #bands: [visd, vist] + +ird: + units: cm-1 + grid-type: regular + wavenumber-range: [100., 10000.] + num-bins: 1 + opacity: [H2-ir] + rt-solver: Disort + parameters: {btemp: 0., ttemp: 0} #assuming Tint = 350. K + flags: [broad_band, thermal_emission] + +irt: + units: cm-1 + grid-type: regular + wavenumber-range: [100., 10000.] + num-bins: 1 + opacity: [H2-ir] + rt-solver: Toon + flags: [broad_band, thermal_emission] + +irtck: + units: cm-1 + wavenumber-range: [100., 10000.] + grid-type: cktable + opacity: [premix] + rt-solver: Toon + flags: [thermal_emission] + +visd: + units: cm-1 + grid-type: regular + wavenumber-range: [10000., 50000.] + num-bins: 1 + opacity: [H2-vis] + rt-solver: Disort + parameters: {temp0: 1000., dist_au: 1., umu0: 1., phi0: 0.} + flags: [broad_band, stellar_beam] + +vist: + units: cm-1 + grid-type: regular + wavenumber-range: [10000., 50000.] + num-bins: 1 + opacity: [H2-vis] + rt-solver: Toon + parameters: {temp0: 1000., dist_au: 1., umu0: 1., phi0: 0.} + flags: [broad_band, stellar_beam] + +vistck: + units: cm-1 + wavenumber-range: [100., 10000.] + grid-type: cktable + opacity: [premix] + rt-solver: Toon + parameters: {temp0: 1000., dist_au: 1., umu0: 1., phi0: 0.} + flags: [broad_band, stellar_beam] + +Disort-flags: + ibcnd: false + usrtau: false + usrang: false + onlyfl: true + spher: false + intensity_correction: true + old_intensity_correction: false + general_source: false + output_uum: false + quiet: true + print-input: false + print-fluxes: false + print-intensity: false + print-transmissivity: false + print-phase-function: true + diff --git a/examples/2024-XZhang-cloud-rt/svp_functions.hpp b/examples/2024-XZhang-cloud-rt/svp_functions.hpp index 3652e7da..e90a6c3d 100644 --- a/examples/2024-XZhang-cloud-rt/svp_functions.hpp +++ b/examples/2024-XZhang-cloud-rt/svp_functions.hpp @@ -137,10 +137,63 @@ void enroll_vapor_Function_fe(Thermodynamics::SVPFunc1Container &svp_func1, } } + +// CH4 svp +double sat_vapor_p_CH4(double T) { + double betal = 10.15, gammal = 1.87, tr = 90.67, + pr = 11690.; + return SatVaporPresIdeal(T / tr, pr, betal, gammal); +} + +void enroll_vapor_Function_CH4(Thermodynamics::SVPFunc1Container &svp_func1, + std::vector const &cloud_index_set) { + Application::Logger app("snap"); + app->Log("Enrolling CH4 vapor pressures"); + + auto pindex = IndexMap::GetInstance(); + if (!pindex->HasVapor("CH4")) return; + int iCH4 = pindex->GetVaporId("CH4"); + + svp_func1[iCH4][0] = [](AirParcel const &qfrac, int, int) { + return sat_vapor_p_CH4(qfrac.w[IDN]); + }; + + for (int n = 1; n < cloud_index_set[iCH4].size(); ++n) { + svp_func1[iCH4][n] = NullSatVaporPres1; + } +} + +// H2S svp +double sat_vapor_p_H2S(double T) { + double betal = 23.67, gammal = 3.96, tr = 187.7, + pr = 23200.; + return SatVaporPresIdeal(T / tr, pr, betal, gammal); +} + +void enroll_vapor_Function_H2S(Thermodynamics::SVPFunc1Container &svp_func1, + std::vector const &cloud_index_set) { + Application::Logger app("snap"); + app->Log("Enrolling H2S vapor pressures"); + + auto pindex = IndexMap::GetInstance(); + if (!pindex->HasVapor("H2S")) return; + int iH2S = pindex->GetVaporId("H2S"); + + svp_func1[iH2S][0] = [](AirParcel const &qfrac, int, int) { + return sat_vapor_p_H2S(qfrac.w[IDN]); + }; + + for (int n = 1; n < cloud_index_set[iH2S].size(); ++n) { + svp_func1[iH2S][n] = NullSatVaporPres1; + } +} + void Thermodynamics::enrollVaporFunctions() { enroll_vapor_Function_H2O(svp_func1_, cloud_index_set_); enroll_vapor_Function_NH3(svp_func1_, cloud_index_set_); enroll_vapor_Function_NH4SH(svp_func1_, cloud_index_set_); enroll_vapor_Function_mgsio3(svp_func1_, cloud_index_set_); enroll_vapor_Function_fe(svp_func1_, cloud_index_set_); + enroll_vapor_Function_H2S(svp_func1_, cloud_index_set_); + enroll_vapor_Function_CH4(svp_func1_, cloud_index_set_); } diff --git a/examples/2024-XZhang-minichem/jup-dry.yaml b/examples/2024-XZhang-minichem/jup-dry.yaml new file mode 100644 index 00000000..04f3b210 --- /dev/null +++ b/examples/2024-XZhang-minichem/jup-dry.yaml @@ -0,0 +1,15 @@ +reference-state: + Tref: 0. + Pref: 1.e5 + +species: + - name: dry + composition: {H: 1.5, He: 0.15} + cv_R: 2.5 + +geometry: + type: cartesian + +dynamics: + equation-of-state: + type: ideal-moist diff --git a/examples/2024-XZhang-minichem/jup_onetra.cpp b/examples/2024-XZhang-minichem/jup_onetra.cpp new file mode 100644 index 00000000..9cb472e5 --- /dev/null +++ b/examples/2024-XZhang-minichem/jup_onetra.cpp @@ -0,0 +1,290 @@ +/* ------------------------------------------------------------------------------------- + * SNAP Example Program + * + * Contributer: + * Cheng Li, University of Michigan + * + * Year: 2023 + * Contact: chengcli@umich.edu + * Reference: Test Jupiter CRM + * ------------------------------------------------------------------------------------- + */ +// C++ headers +#include +#include +#include +#include +#include + +// athena +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// application +#include +#include + +// canoe +#include +#include + +// climath +#include +#include + +// snap +#include +#include + +// special includes +#include + +// minichem +//#include + +Real grav, P0, T0, Tmin, prad, hrate; +int iH2O; +//'OH','H2','H2O','H','CO','CO2','O','CH4','C2H2','NH3','N2','HCN', 'He' +std::vector vmass = {17.01, 2.02, 18.02, 1.01, 28.01, 44.01, 16., + 16.05, 26.04, 17.04, 28.02, 27.03, 4.}; +//MiniChem *mc; + +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { + AllocateUserOutputVariables(5 + NVAPOR); + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + SetUserOutputVariableName(2, "thetav"); + SetUserOutputVariableName(3, "mse"); + SetUserOutputVariableName(4, "pres"); + for (int n = 1; n <= NVAPOR; ++n) { + std::string name = "rh" + std::to_string(n); + SetUserOutputVariableName(4 + n, name.c_str()); + } +} + +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + 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(this, k, j, i); + user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); + // theta_v + user_out_var(2, k, j, i) = + user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); + // mse + user_out_var(3, k, j, i) = + pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i); + user_out_var(4, k, j, i) = phydro->w(IPR, k, j, i); + // relative humidity + for (int n = 1; n <= NVAPOR; ++n) + user_out_var(4 + n, k, j, i) = + pthermo->RelativeHumidity(this, n, k, j, i); + } +} + +void Forcing(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, const AthenaArray &prim_scalar, + AthenaArray const &bcc, AthenaArray &du, + AthenaArray &cons_scalar) { + int is = pmb->is, js = pmb->js, ks = pmb->ks; + int ie = pmb->ie, je = pmb->je, ke = pmb->ke; + auto pthermo = Thermodynamics::GetInstance(); + auto phydro = pmb->phydro; + std::vector vmr(NCHEMISTRY); + std::vector vmr_sp(NCHEMISTRY - 1); + std::vector mmr(NCHEMISTRY); + Real sumVMR; + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); + + air.ToMoleFraction(); + Real cv = + get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass()); + + if (w(IPR, k, j, i) < prad) { + du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv * + (1. + 1.E-4 * sin(2. * M_PI * rand() / RAND_MAX)); + } +/* + // minichem + Real temp = pthermo->GetTemp(w.at(k, j, i)); + Real pres = w(IPR, k, j, i); + + // change mmr to vmr and normalize + for (int n = 0; n < NCHEMISTRY; ++n) + vmr[n] = prim_scalar(NCLOUD + n, k, j, i) / vmass[n]; + sumVMR = std::accumulate(vmr.begin(), vmr.end(), static_cast(0)); + for (auto &value : vmr) value /= sumVMR; + + // call minichem for sp species + for (int n = 0; n < NCHEMISTRY - 1; ++n) vmr_sp[n] = vmr[n]; +// if (temp > 200.) mc->Run(temp, pres, dt, vmr_sp.data(), "NCHO"); + for (int n = 0; n < NCHEMISTRY - 1; ++n) vmr[n] = vmr_sp[n]; + + // change vmr to mmr and normalize + for (int n = 0; n < NCHEMISTRY; ++n) mmr[n] = vmr[n] * vmass[n]; + sumVMR = std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); + for (auto &value : mmr) value /= sumVMR; + + for (int n = 0; n < NCHEMISTRY; ++n) { + cons_scalar(NCLOUD + n, k, j, i) += + phydro->w(IDN, k, j, i) * + (mmr[n] - prim_scalar(NCLOUD + n, k, j, i)); + } +*/ + + } +} + +void Mesh::InitUserMeshData(ParameterInput *pin) { + grav = -pin->GetReal("hydro", "grav_acc1"); + + P0 = pin->GetReal("problem", "P0"); + T0 = pin->GetReal("problem", "T0"); + + Tmin = pin->GetReal("problem", "Tmin"); + prad = pin->GetReal("problem", "prad"); + hrate = pin->GetReal("problem", "hrate") / 86400.; + + // index + auto pindex = IndexMap::GetInstance(); + iH2O = pindex->GetVaporId("H2O"); +// EnrollUserExplicitSourceFunction(Forcing); + + // minichem +/* mc = new MiniChem(); + mc->SetDataFile("chem_data/mini_chem_data_NCHO.txt"); + mc->SetSpeciesFile("chem_data/mini_chem_sp_NCHO.txt"); + mc->SetNetworkDir("chem_data/1x/"); + mc->SetMetallicityStr("1x"); + mc->Initialize(); +*/ +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + srand(Globals::my_rank + time(0)); + + Application::Logger app("main"); + app->Log("ProblemGenerator: jupiter_crm"); + + auto pthermo = Thermodynamics::GetInstance(); + + // mesh limits + Real x1min = pmy_mesh->mesh_size.x1min; + Real x1max = pmy_mesh->mesh_size.x1max; + + // request temperature and pressure + app->Log("request T", T0); + app->Log("request P", P0); + + // thermodynamic constants + Real gamma = pin->GetReal("hydro", "gamma"); + Real Rd = pthermo->GetRd(); + Real cp = gamma / (gamma - 1.) * Rd; + + // set up an adiabatic atmosphere + int max_iter = 400, iter = 0; + Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); + + AirParcel air(AirParcel::Type::MoleFrac); + + // estimate surface temperature and pressure + Real Ts = T0 - grav / cp * x1min; + Real Ps = P0 * pow(Ts / T0, cp / Rd); + Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6; + Real xtra = pin->GetReal("problem", "tracer.ppb") / 1.E9; + + while (iter++ < max_iter) { + // read in vapors + air.w[iH2O] = xH2O; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // stop at just above P0 + for (int i = is; i <= ie; ++i) { + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); + if (air.w[IPR] < P0) break; + } + + // make up for the difference + Ts += T0 - air.w[IDN]; + if (std::abs(T0 - air.w[IDN]) < Ttol) break; + + app->Log("Iteration #", iter); + app->Log("T", air.w[IDN]); + } + + if (iter > max_iter) { + throw RuntimeError("ProblemGenerator", "maximum iteration reached"); + } + + // construct atmosphere from bottom up + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + air.SetZero(); + air.w[iH2O] = xH2O; + air.w[IPR] = Ps; + air.w[IDN] = Ts; + + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); + + int i = is; + for (; i <= ie; ++i) { + air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); + if (air.w[IDN] < Tmin) break; + AirParcelHelper::distribute_to_primitive(this, k, j, i, air); + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5); + } + + // Replace adiabatic atmosphere with isothermal atmosphere if temperature + // is too low + for (; i <= ie; ++i) { + air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); + AirParcelHelper::distribute_to_primitive(this, k, j, i, air); + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + } + } + + // minichem initialize conserved variables + // interpolate from ce table + std::vector vmr_ic(NCHEMISTRY); + std::vector mmr(NCHEMISTRY); + std::string ic_file = "chem_data/IC/mini_chem_IC_FastChem_1x.txt"; + + for (int k = ks; k <= ke; k++) + for (int j = js; j <= je; j++) + for (int i = is; i <= ie; i++) { + double mu; + double T_in = pthermo->GetTemp(this, k, j, i); + double P_in = phydro->w(IPR, k, j, i); +// interp_ce_table(NCHEMISTRY, T_in, P_in, vmr_ic.data(), &mu, ic_file); + + // change vmr to mmr and normalize +// for (int n = 0; n < NCHEMISTRY; ++n) mmr[n] = vmr_ic[n] * vmass[n]; +// Real sumVMR = +// std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); +// for (auto &value : mmr) value /= sumVMR; + if (P_in > 6.e7) { + for (int n = 0; n < NCHEMISTRY; ++n) { + mmr[n] = xtra; + pscalars->s(NCLOUD + n, k, j, i) += phydro->w(IDN, k, j, i) * mmr[n]; + } + } + } +} diff --git a/examples/2024-XZhang-minichem/mini.inp b/examples/2024-XZhang-minichem/mini.inp index 2e061720..3b4f92ed 100644 --- a/examples/2024-XZhang-minichem/mini.inp +++ b/examples/2024-XZhang-minichem/mini.inp @@ -68,7 +68,8 @@ nx2 = 16 nx3 = 1 -gamma = 1.42 +Rd = 3936.04 +gamma = 1.4 grav_acc1 = -24.79 implicit_flag = 1 @@ -79,18 +80,6 @@ cloud = H2O(c), H2O(p) microphysics_config = water_ammonia.yaml - -Rd = 3777. # mu = 2.3175 g/mol -eps1 = 8.18 8.18 8.18 -beta1 = 0. 24.845 24.845 -rcp1 = 0.15 0.33 0.33 -Ttriple1 = 273.16 -Ptriple1 = 611.7 - -sa.relax = 1.0 -sa.max_iter = 5 -sa.ftol = 1.e-2 - #packages = bot_heating, top_cooling #bot_heating.flux = 8. @@ -106,6 +95,8 @@ top_sponge_lyr.width = 20.E3 #relax_time = 1.E4 +config_file = jup-dry.yaml + qH2O.ppmv = 2500. #minichem diff --git a/examples/2024-XZhang-minichem/minichem.cpp b/examples/2024-XZhang-minichem/minichem.cpp index bf2683d7..444731fb 100644 --- a/examples/2024-XZhang-minichem/minichem.cpp +++ b/examples/2024-XZhang-minichem/minichem.cpp @@ -1,14 +1,3 @@ -/* ------------------------------------------------------------------------------------- - * SNAP Example Program - * - * Contributer: - * Cheng Li, University of Michigan - * - * Year: 2023 - * Contact: chengcli@umich.edu - * Reference: Test Jupiter CRM - * ------------------------------------------------------------------------------------- - */ // C++ headers #include #include @@ -16,6 +5,9 @@ #include #include +// yaml +#include + // athena #include #include @@ -33,30 +25,26 @@ #include #include +// snap +#include +#include + // canoe -#include -#include -#include +#include + #include +#include // climath #include #include -// snap -#include -#include - -// special includes -#include - // minichem #include Real grav, P0, T0, Tmin, prad, hrate, xH2O; Real Tbot, xtra, ptra, Ttop_tau, Tbot_tau, Trabot_tau; std::string met; -int iH2O; bool use_mini, use_mini_ic, use_tra_ic, fix_bot_tra; //'OH','H2','H2O','H','CO','CO2','O','CH4','C2H2','NH3','N2','HCN', 'He' @@ -65,14 +53,16 @@ std::vector vmass = {17.01, 2.02, 18.02, 1.01, 28.01, 44.01, 16., MiniChem *mc; void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { - if (NVAPOR > 0) { - AllocateUserOutputVariables(5 + NVAPOR); + int nvapor = 0; + + if (nvapor > 0) { + AllocateUserOutputVariables(5 + nvapor); SetUserOutputVariableName(0, "temp"); SetUserOutputVariableName(1, "theta"); SetUserOutputVariableName(2, "pres"); SetUserOutputVariableName(3, "thetav"); SetUserOutputVariableName(4, "mse"); - for (int n = 1; n <= NVAPOR; ++n) { + for (int n = 1; n <= nvapor; ++n) { std::string name = "rh" + std::to_string(n); SetUserOutputVariableName(4 + n, name.c_str()); } @@ -85,16 +75,24 @@ void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { } 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(P0 / pres, chi); + + auto temp_a = temp.accessor(); + auto theta_a = theta.accessor(); + auto pres_a = pres.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(this, k, j, i); - user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i); - user_out_var(2, k, j, i) = phydro->w(IPR, k, j, 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]; + user_out_var(2, k, j, i) = pres_a[k][j][i]; - if (NVAPOR > 0) { + /*if (NVAPOR > 0) { // theta_v user_out_var(3, k, j, i) = user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i); @@ -105,7 +103,7 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { for (int n = 1; n <= NVAPOR; ++n) user_out_var(4 + n, k, j, i) = pthermo->RelativeHumidity(this, n, k, j, i); - } + }*/ } } @@ -115,82 +113,86 @@ void Forcing(MeshBlock *pmb, Real const time, Real const dt, AthenaArray &cons_scalar) { int is = pmb->is, js = pmb->js, ks = pmb->ks; int ie = pmb->ie, je = pmb->je, ke = pmb->ke; - auto pthermo = Thermodynamics::GetInstance(); auto phydro = pmb->phydro; + auto temp = get_temp(pmb->pimpl->peos, w); + auto cv = get_cv(pmb->pimpl->peos, w); + + auto temp_a = temp.accessor(); + auto cv_a = cv.accessor(); + for (int k = ks; k <= ke; ++k) for (int j = js; j <= je; ++j) for (int i = is; i <= ie; ++i) { - auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); - - Real tem = pthermo->GetTemp(w.at(k,j,i)); - air.ToMoleFraction(); - Real cv = - get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass()); - - if (w(IPR, k, j, i) < prad) { -// body cooling + auto tem = temp_a[k][j][i]; + auto cv = cv_a[k][j][i]; + if (w(IPR, k, j, i) < prad) { + // body cooling du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv * (1. + 1.E-4 * sin(2. * M_PI * rand() / RAND_MAX)); -// newtonian at the top - if (tem < Tmin) { - du(IEN, k, j, i) += dt * (Tmin - tem)/Ttop_tau * w(IDN, k, j, i) * cv; - } + // newtonian at the top + if (tem < Tmin) { + du(IEN, k, j, i) += + dt * (Tmin - tem) / Ttop_tau * w(IDN, k, j, i) * cv; + } } -// relax T at the bottom - if(i == is) du(IEN, k, j, i) += dt * (Tbot - tem)/Tbot_tau * w(IDN, k, j, i) * cv; - } - - if (NCHEMISTRY > 0) { + // relax T at the bottom + if (i == is) + du(IEN, k, j, i) += + dt * (Tbot - tem) / Tbot_tau * w(IDN, k, j, i) * cv; + } + + if (NCHEM > 0) { if (use_mini) { - std::vector vmr(NCHEMISTRY); - std::vector vmr_sp(NCHEMISTRY - 1); - std::vector mmr(NCHEMISTRY); - Real sumVMR; - for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) - for (int i = is; i <= ie; ++i) { - // minichem - Real temp = pthermo->GetTemp(w.at(k, j, i)); - Real pres = w(IPR, k, j, i); - - // change mmr to vmr and normalize - for (int n = 0; n < NCHEMISTRY; ++n) - vmr[n] = prim_scalar(NCLOUD + n, k, j, i) / vmass[n]; - sumVMR = - std::accumulate(vmr.begin(), vmr.end(), static_cast(0)); - for (auto &value : vmr) value /= sumVMR; - - // call minichem for sp species - for (int n = 0; n < NCHEMISTRY - 1; ++n) vmr_sp[n] = vmr[n]; - if (temp > 200.) mc->Run(temp, pres, dt, vmr_sp.data(), "NCHO"); - for (int n = 0; n < NCHEMISTRY - 1; ++n) vmr[n] = vmr_sp[n]; - - // change vmr to mmr and normalize - for (int n = 0; n < NCHEMISTRY; ++n) mmr[n] = vmr[n] * vmass[n]; - sumVMR = - std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); - for (auto &value : mmr) value /= sumVMR; - - for (int n = 0; n < NCHEMISTRY; ++n) { - cons_scalar(NCLOUD + n, k, j, i) += - phydro->w(IDN, k, j, i) * - (mmr[n] - prim_scalar(NCLOUD + n, k, j, i)); + std::vector vmr(NCHEM); + std::vector vmr_sp(NCHEM - 1); + std::vector mmr(NCHEM); + Real sumVMR; + auto temp_a = temp.accessor(); + + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + // minichem + Real temp = temp_a[k][j][i]; + Real pres = w(IPR, k, j, i); + + // change mmr to vmr and normalize + for (int n = 0; n < NCHEM; ++n) + vmr[n] = prim_scalar(n, k, j, i) / vmass[n]; + sumVMR = + std::accumulate(vmr.begin(), vmr.end(), static_cast(0)); + for (auto &value : vmr) value /= sumVMR; + + // call minichem for sp species + for (int n = 0; n < NCHEM - 1; ++n) vmr_sp[n] = vmr[n]; + if (temp > 200.) mc->Run(temp, pres, dt, vmr_sp.data(), "NCHO"); + for (int n = 0; n < NCHEM - 1; ++n) vmr[n] = vmr_sp[n]; + + // change vmr to mmr and normalize + for (int n = 0; n < NCHEM; ++n) mmr[n] = vmr[n] * vmass[n]; + sumVMR = + std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); + for (auto &value : mmr) value /= sumVMR; + + for (int n = 0; n < NCHEM; ++n) { + cons_scalar(n, k, j, i) += + phydro->w(IDN, k, j, i) * (mmr[n] - prim_scalar(n, k, j, i)); + } } - } - } else if (fix_bot_tra) { - for (int k = ks; k <= ke; k++) - for (int j = js; j <= je; j++) - for (int i = is; i <= ie; i++) { - if (phydro->w(IPR, k, j, i) > ptra) { - for (int n = 0; n < NCHEMISTRY; ++n) { - cons_scalar(NCLOUD + n, k, j, i) += - dt/Trabot_tau * phydro->w(IDN, k, j, i) * - (xtra - prim_scalar(NCLOUD + n, k, j, i)); + } else if (fix_bot_tra) { + for (int k = ks; k <= ke; k++) + for (int j = js; j <= je; j++) + for (int i = is; i <= ie; i++) { + if (phydro->w(IPR, k, j, i) > ptra) { + for (int n = 0; n < NCHEM; ++n) { + cons_scalar(n, k, j, i) += dt / Trabot_tau * + phydro->w(IDN, k, j, i) * + (xtra - prim_scalar(n, k, j, i)); + } } } - } - } + } } } @@ -213,11 +215,6 @@ void Mesh::InitUserMeshData(ParameterInput *pin) { xtra = pin->GetReal("problem", "tracer.ppb") / 1.E9; ptra = pin->GetReal("problem", "tracer.pre"); - if (NVAPOR > 0) { - // index - auto pindex = IndexMap::GetInstance(); - iH2O = pindex->GetVaporId("H2O"); - } EnrollUserExplicitSourceFunction(Forcing); if (use_mini) { @@ -231,134 +228,174 @@ void Mesh::InitUserMeshData(ParameterInput *pin) { } } -void MeshBlock::ProblemGenerator(ParameterInput *pin) { - srand(Globals::my_rank + time(0)); +torch::Tensor setup_moist_adiabatic_profile(std::string infile) { + auto config = YAML::LoadFile(infile); + auto Ps = config["problem"]["Ps"].as(1.e5); + auto Ts = config["problem"]["Ts"].as(300.); + auto Tmin = config["problem"]["Tmin"].as(200.); + auto grav = -config["forcing"]["const-gravity"]["grav1"].as(); + + // initialize the block + auto block = snap::MeshBlock(snap::MeshBlockOptions::from_yaml(infile)); + + // useful modules + auto phydro = block->phydro; + auto pcoord = phydro->pcoord; + auto peos = phydro->peos; + auto m = block->named_modules()["hydro.eos.thermo"]; + auto thermo_y = std::dynamic_pointer_cast(m); + + // dimensions and indices + int nc3 = pcoord->x3v.size(0); + int nc2 = pcoord->x2v.size(0); + int nc1 = pcoord->x1v.size(0); + int ny = thermo_y->options.species().size() - 1; + int nvar = peos->nvar(); + + // construct an adiabatic atmosphere + kintera::ThermoX thermo_x(thermo_y->options); + + auto temp = Ts * torch::ones({nc3, nc2}, torch::kFloat64); + auto pres = Ps * torch::ones({nc3, nc2}, torch::kFloat64); + auto xfrac = torch::zeros({nc3, nc2, 1 + ny}, torch::kFloat64); + auto w = torch::zeros({nvar, nc3, nc2, nc1}, torch::kFloat64); + + // read in compositions + for (int i = 1; i <= ny; ++i) { + auto name = thermo_y->options.species()[i]; + auto xmixr = config["problem"]["x" + name].as(0.); + xfrac.select(2, i) = xmixr; + } + + // dry air mole fraction + xfrac.select(2, 0) = 1. - xfrac.narrow(-1, 1, ny).sum(-1); - Application::Logger app("main"); - app->Log("ProblemGenerator: jupiter_crm"); + // adiabatic extrapolate half a grid to cell center + int is = pcoord->is(); + int ie = pcoord->ie(); + auto dz = pcoord->dx1f[is].item(); + thermo_x->extrapolate_ad(temp, pres, xfrac, grav, dz / 2.); - auto pthermo = Thermodynamics::GetInstance(); + int i = is; + int nvapor = thermo_x->options.vapor_ids().size(); + int ncloud = thermo_x->options.cloud_ids().size(); + for (; i <= ie; ++i) { + auto conc = thermo_x->compute("TPX->V", {temp, pres, xfrac}); - // mesh limits - Real x1min = pmy_mesh->mesh_size.x1min; - Real x1max = pmy_mesh->mesh_size.x1max; + w[IPR].select(2, i) = pres; + w[IDN].select(2, i) = thermo_x->compute("V->D", {conc}); - // request temperature and pressure - app->Log("request T", T0); - app->Log("request P", P0); + auto result = thermo_x->compute("X->Y", {xfrac}); + w.narrow(0, snap::ICY, ny).select(3, i) = + thermo_x->compute("X->Y", {xfrac}); + + if ((temp < Tmin).any().item()) break; + dz = pcoord->dx1f[i].item(); + thermo_x->extrapolate_ad(temp, pres, xfrac, grav, dz); + } + + // isothermal extrapolation + for (; i <= ie; ++i) { + auto mu = (thermo_x->mu * xfrac).sum(-1); + dz = pcoord->dx1f[i].item(); + pres *= exp(-grav * mu * dz / (kintera::constants::Rgas * temp)); + auto conc = thermo_x->compute("TPX->V", {temp, pres, xfrac}); + w[IPR].select(2, i) = pres; + w[IDN].select(2, i) = thermo_x->compute("V->D", {conc}); + w.narrow(0, snap::ICY, ny).select(3, i) = + thermo_x->compute("X->Y", {xfrac}); + } + + // add noise + w[IVX] += 0.01 * torch::rand_like(w[IVX]); + w[IVY] += 0.01 * torch::rand_like(w[IVY]); + + return w; +} + +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + // auto infile = pin->GetString("problem", "config_file"); + // auto w = setup_moist_adiabatic_profile(infile) + + srand(Globals::my_rank + time(0)); // thermodynamic constants Real gamma = pin->GetReal("hydro", "gamma"); - Real Rd = pthermo->GetRd(); + Real Rd = pin->GetReal("hydro", "Rd"); Real cp = gamma / (gamma - 1.) * Rd; + std::cout << "Rd = " << Rd << ", cp = " << cp << std::endl; // set up an adiabatic atmosphere int max_iter = 400, iter = 0; Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01); - AirParcel air(AirParcel::Type::MoleFrac); - // estimate surface temperature and pressure + auto x1min = pcoord->x1f(is); + auto x1max = pcoord->x1f(ie + 1); + Real Ts = T0 - grav / cp * x1min; Real Ps = P0 * pow(Ts / T0, cp / Rd); - if (NVAPOR > 0) xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6; - - while (iter++ < max_iter) { - // read in vapors - if (NVAPOR > 0) air.w[iH2O] = xH2O; - air.w[IPR] = Ps; - air.w[IDN] = Ts; - - // stop at just above P0 - for (int i = is; i <= ie; ++i) { - pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav); - if (air.w[IPR] < P0) break; - } - - // make up for the difference - Ts += T0 - air.w[IDN]; - if (std::abs(T0 - air.w[IDN]) < Ttol) break; - app->Log("Iteration #", iter); - app->Log("T", air.w[IDN]); - } - - if (iter > max_iter) { - throw RuntimeError("ProblemGenerator", "maximum iteration reached"); - } - - // construct atmosphere from bottom up + // Loop over the grids and set initial condition for (int k = ks; k <= ke; ++k) - for (int j = js; j <= je; ++j) { - air.SetZero(); - if (NVAPOR > 0) air.w[iH2O] = xH2O; - air.w[IPR] = Ps; - air.w[IDN] = Ts; - - // half a grid to cell center - pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav); - - int i = is; - for (; i <= ie; ++i) { - air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); - if (air.w[IDN] < Tmin) break; - AirParcelHelper::distribute_to_primitive(this, k, j, i, air); - AirParcelHelper::distribute_to_conserved(this, k, j, i, air); - //pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5); - pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 0.); - } - - // Replace adiabatic atmosphere with isothermal atmosphere if temperature - // is too low - for (; i <= ie; ++i) { - air.w[IVX] = 0.001 * (1. * rand() / RAND_MAX - 0.5); - AirParcelHelper::distribute_to_primitive(this, k, j, i, air); - AirParcelHelper::distribute_to_conserved(this, k, j, i, air); - pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav); + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + Real x1 = pcoord->x1v(i); + Real x2 = pcoord->x2v(j); + Real x3 = pcoord->x3v(k); + Real temp = T0 - grav * x1 / cp; + phydro->w(IPR, k, j, i) = P0 * pow(temp / T0, cp / Rd); + phydro->w(IDN, k, j, i) = phydro->w(IPR, k, j, i) / (Rd * temp); + phydro->w(IVX, k, j, i) = 0.001 * (1. * rand() / RAND_MAX - 0.5); } - } - Tbot = pthermo->GetTemp(this, ks, js, is);; - - if (NCHEMISTRY > 0) { - - if (use_mini_ic) { - // minichem initialize conserved variables - // interpolate from ce table - std::vector vmr_ic(NCHEMISTRY); - std::vector mmr(NCHEMISTRY); - std::string ic_file = "chem_data/IC/mini_chem_IC_FastChem_" + met + ".txt"; - - for (int k = ks; k <= ke; k++) - for (int j = js; j <= je; j++) - for (int i = is; i <= ie; i++) { - double mu; - double T_in = pthermo->GetTemp(this, k, j, i); - double P_in = phydro->w(IPR, k, j, i); - interp_ce_table(NCHEMISTRY, T_in, P_in, vmr_ic.data(), &mu, ic_file); - - // change vmr to mmr and normalize - for (int n = 0; n < NCHEMISTRY; ++n) mmr[n] = vmr_ic[n] * vmass[n]; - Real sumVMR = - std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); - for (auto &value : mmr) value /= sumVMR; - - for (int n = 0; n < NCHEMISTRY; ++n) { - pscalars->s(NCLOUD + n, k, j, i) += - phydro->w(IDN, k, j, i) * mmr[n]; + auto temp = get_temp(pimpl->peos, phydro->w); + auto temp_a = temp.accessor(); + + Tbot = temp_a[ks][js][is]; + + std::cout << "ptra = " << ptra << " Pa" << std::endl; + + if (NCHEM > 0) { + if (use_mini_ic) { + // minichem initialize conserved variables + // interpolate from ce table + std::vector vmr_ic(NCHEM); + std::vector mmr(NCHEM); + std::string ic_file = + "chem_data/IC/mini_chem_IC_FastChem_" + met + ".txt"; + + for (int k = ks; k <= ke; k++) + for (int j = js; j <= je; j++) + for (int i = is; i <= ie; i++) { + double mu; + double T_in = temp_a[k][j][i]; + double P_in = phydro->w(IPR, k, j, i); + interp_ce_table(NCHEM, T_in, P_in, vmr_ic.data(), &mu, ic_file); + + // change vmr to mmr and normalize + for (int n = 0; n < NCHEM; ++n) mmr[n] = vmr_ic[n] * vmass[n]; + Real sumVMR = + std::accumulate(mmr.begin(), mmr.end(), static_cast(0)); + for (auto &value : mmr) value /= sumVMR; + + for (int n = 0; n < NCHEM; ++n) { + pscalars->s(n, k, j, i) += phydro->w(IDN, k, j, i) * mmr[n]; + } } - } - } else if (use_tra_ic) { - for (int k = ks; k <= ke; k++) - for (int j = js; j <= je; j++) - for (int i = is; i <= ie; i++) { - if (phydro->w(IPR, k, j, i) > ptra) { - for (int n = 0; n < NCHEMISTRY; ++n) { - pscalars->s(NCLOUD + n, k, j, i) += - phydro->w(IDN, k, j, i) * xtra; + } else if (use_tra_ic) { + for (int k = ks; k <= ke; k++) + for (int j = js; j <= je; j++) + for (int i = is; i <= ie; i++) { + if (phydro->w(IPR, k, j, i) > ptra) { + for (int n = 0; n < NCHEM; ++n) { + pscalars->s(n, k, j, i) += phydro->w(IDN, k, j, i) * xtra; + } } } - } - } + } } + + peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie, + js, je, ks, ke); } diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 2430bf57..8346a142 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -42,6 +42,7 @@ endif() if (${TASK} STREQUAL "rtdry" OR ${TASK} STREQUAL "rtcld" OR + ${TASK} STREQUAL "ice_h2s" OR ${TASK} STREQUAL "rtcldexo") add_subdirectory(2024-XZhang-cloud-rt) endif() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2b35c535..d736eb73 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -19,6 +19,7 @@ include_directories( SYSTEM ${HARP_INCLUDE_DIR} SYSTEM ${KINTERA_INCLUDE_DIR} SYSTEM ${SNAP_INCLUDE_DIR} + SYSTEM ${FMT_INCLUDE_DIR} ) set(namel canoe) @@ -27,7 +28,17 @@ string(TOUPPER ${namel} nameu) string(TOLOWER ${CMAKE_BUILD_TYPE} buildl) string(TOUPPER ${CMAKE_BUILD_TYPE} buildu) -file(GLOB src_files *.cpp) +file(GLOB src_files + *.cpp + climath/*.c + astro/*.cpp + astro/Jupiter/*.cpp + astro/Saturn/*.cpp + utils/*.cpp + tasklist/*.cpp + exo3/*.cpp + exchanger/*.cpp + ) add_library(${namel}_${buildl} OBJECT @@ -39,15 +50,27 @@ set_target_properties(${namel}_${buildl} COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}} ) -add_subdirectory(climath) -add_subdirectory(astro) -add_subdirectory(utils) +target_link_libraries(${namel}_${buildl} + PUBLIC + ${NETCDF_LIBRARIES} + ${PNETCDF_LIBRARY} + ${TORCH_LIBRARY} + ${TORCH_CPU_LIBRARY} + ${C10_LIBRARY} + ${DISORT_LIBRARY} + ${HARP_LIBRARY} + ${KINTERA_LIBRARY} + ${SNAP_LIBRARY} + ${MPI_CXX_LIBRARIES} + ${GLOG_LIBRARIES} + ${OpenMP_libomp_LIBRARY} + ${EXTRA_LIBS} + fmt::fmt + yaml-cpp::yaml-cpp + ) + add_subdirectory(snap) add_subdirectory(outputs) -add_subdirectory(tasklist) -add_subdirectory(exo3) -add_subdirectory(exchanger) -add_subdirectory(special) if(UNIX AND NOT APPLE) set(EXTRA_LIBS dl stdc++fs) @@ -57,33 +80,10 @@ set(CANOE_LIBRARY_${buildu} "canoe_${buildl}" "athenapp_${buildl}" "application_${buildl}" - "climath_${buildl}" - "astro_${buildl}" - "utils_${buildl}" "snap_${buildl}" "outputs_${buildl}" - "tasklist_${buildl}" - "exo3_${buildl}" - "exchanger_${buildl}" $<$:minichem::minichem++> $<$:exofmsrt::exofmsrt++> - ${YAML_CPP_LIBRARIES} - ${CPPDISORT_LIBRARY_${buildu}} - ${PYTHON_LIBRARY_RELEASE} - ${NETCDF_LIBRARIES} - ${PNETCDF_LIBRARY} - ${MPI_CXX_LIBRARIES} - ${GLOG_LIBRARIES} - ${OpenMP_libomp_LIBRARY} - ${CANTERA_LIBRARY} - ${TORCH_LIBRARY} - ${TORCH_CPU_LIBRARY} - ${C10_LIBRARY} - ${DISORT_LIBRARY} - ${HARP_LIBRARY} - ${KINTERA_LIBRARY} - ${SNAP_LIBRARY} - ${EXTRA_LIBS} CACHE INTERNAL "canoe library") set(CANOE_INCLUDE_DIR diff --git a/src/interface/eos.hpp b/src/interface/eos.hpp index ae48da76..67b0ae85 100644 --- a/src/interface/eos.hpp +++ b/src/interface/eos.hpp @@ -4,6 +4,8 @@ #include // snap +#include + #include // canoe @@ -24,6 +26,19 @@ inline torch::Tensor get_temp(snap::IdealMoist peos, return peos->compute("W->T", {w}); } +inline torch::Tensor get_cv(snap::IdealMoist peos, + AthenaArray const& hydro_w) { + auto pthermo = peos->pthermo; + auto mud = kintera::species_weights[0]; + auto Rd = kintera::constants::Rgas / mud; + auto cvd = kintera::species_cref_R[0] * Rd; + int ny = pthermo->options.vapor_ids().size() + + pthermo->options.cloud_ids().size() - 1; + auto w = get_all(hydro_w); + auto fsig = peos->f_sig(w.narrow(0, snap::ICY, ny)); + return cvd * fsig; +} + //! \brief Effective polytropic index //! //! Eq.63 in Li2019