From 8e2402bc391cd533c837428cbb6ad964661080b0 Mon Sep 17 00:00:00 2001 From: Linfeng Li Date: Thu, 7 Dec 2023 14:27:28 -0500 Subject: [PATCH 1/6] lavaplanet input and cpp files --- cmake/lava.cmake | 22 +++ examples/2023-Linf-lava/CMakeLists.txt | 11 ++ examples/2023-Linf-lava/balance_atm.txt | 81 +++++++++ examples/2023-Linf-lava/lava.cpp | 229 ++++++++++++++++++++++++ examples/2023-Linf-lava/lava.inp | 133 ++++++++++++++ examples/2023-Linf-lava/lava.yaml | 152 ++++++++++++++++ examples/CMakeLists.txt | 1 + 7 files changed, 629 insertions(+) create mode 100644 cmake/lava.cmake create mode 100644 examples/2023-Linf-lava/CMakeLists.txt create mode 100644 examples/2023-Linf-lava/balance_atm.txt create mode 100644 examples/2023-Linf-lava/lava.cpp create mode 100644 examples/2023-Linf-lava/lava.inp create mode 100644 examples/2023-Linf-lava/lava.yaml diff --git a/cmake/lava.cmake b/cmake/lava.cmake new file mode 100644 index 00000000..0c9485e2 --- /dev/null +++ b/cmake/lava.cmake @@ -0,0 +1,22 @@ +# configure file for hydrogen-water world + +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(NETCDF ON) +set(PNETCDF ON) +set(MPI ON) +set(DISORT ON) +set(PYTHON_BINDINGS ON) +set(TASKLIST ImplicitHydroTasks) +set(RSOLVER lmars) diff --git a/examples/2023-Linf-lava/CMakeLists.txt b/examples/2023-Linf-lava/CMakeLists.txt new file mode 100644 index 00000000..4a6cf64d --- /dev/null +++ b/examples/2023-Linf-lava/CMakeLists.txt @@ -0,0 +1,11 @@ + +# 1. Compile lava +if (${NVAPOR} EQUAL 1 AND ${NCLOUD} EQUAL 2 AND ${NPHASE_LEGACY} EQUAL 3) + setup_problem(lava) +endif() + +# 2. Copy input files to run directory +file(GLOB inputs *.inp *.yaml) +foreach(input ${inputs}) + file(COPY ${input} DESTINATION ${CMAKE_BINARY_DIR}/bin) +endforeach() diff --git a/examples/2023-Linf-lava/balance_atm.txt b/examples/2023-Linf-lava/balance_atm.txt new file mode 100644 index 00000000..1171dd98 --- /dev/null +++ b/examples/2023-Linf-lava/balance_atm.txt @@ -0,0 +1,81 @@ +IDX HGT PRE TEM SiO +1 0.4793 9.882340930 2147.3 1000000.0 +2 1.4374 9.651159599 2145.0 1000000.0 +3 2.3943 9.425386380 2141.7 1000000.0 +4 3.3506 9.204894760 2142.9 1000000.0 +5 4.3068 8.989561183 2140.8 1000000.0 +6 5.2628 8.779264985 2142.0 1000000.0 +7 6.2168 8.573888326 2132.1 1000000.0 +8 7.1704 8.373316121 2139.9 1000000.0 +9 8.1230 8.177435977 2127.6 1000000.0 +10 9.0741 7.986138131 2133.5 1000000.0 +11 10.0251 7.799315387 2127.0 1000000.0 +12 10.9756 7.616863059 2131.0 1000000.0 +13 11.9251 7.438678907 2122.8 1000000.0 +14 12.8740 7.264663084 2128.6 1000000.0 +15 13.8229 7.094718079 2122.4 1000000.0 +16 14.7704 6.928748661 2122.2 1000000.0 +17 15.7173 6.766661829 2120.2 1000000.0 +18 16.6636 6.608366755 2119.2 1000000.0 +19 17.6087 6.453774738 2114.8 1000000.0 +20 18.5528 6.302799149 2114.9 1000000.0 +21 19.4954 6.155355390 2107.6 1000000.0 +22 20.4378 6.011360839 2114.6 1000000.0 +23 21.3794 5.870734806 2103.8 1000000.0 +24 22.3192 5.733398491 2106.5 1000000.0 +25 23.2583 5.599274936 2100.4 1000000.0 +26 24.1958 5.468288984 2100.0 1000000.0 +27 25.1321 5.340367237 2094.5 1000000.0 +28 26.0669 5.215438010 2093.3 1000000.0 +29 27.0010 5.093431301 2091.2 1000000.0 +30 27.9331 4.974278740 2084.6 1000000.0 +31 28.8636 4.857913560 2084.0 1000000.0 +32 29.7923 4.744270555 2076.5 1000000.0 +33 30.7189 4.633286043 2074.5 1000000.0 +34 31.6439 4.524897834 2069.7 1000000.0 +35 32.5670 4.419045192 2065.7 1000000.0 +36 33.4884 4.315668800 2062.1 1000000.0 +37 34.4083 4.214710731 2059.0 1000000.0 +38 35.3266 4.116114412 2055.0 1000000.0 +39 36.2415 4.019824594 2044.0 1000000.0 +40 37.1541 3.925787320 2044.5 1000000.0 +41 38.0651 3.833949896 2036.6 1000000.0 +42 38.9736 3.744260859 2033.4 1000000.0 +43 39.8803 3.656669952 2028.7 1000000.0 +44 40.7841 3.571128091 2020.6 1000000.0 +45 41.6859 3.487587344 2019.5 1000000.0 +46 42.5860 3.406000897 2012.8 1000000.0 +47 43.4827 3.326323033 2004.6 1000000.0 +48 44.3772 3.248509103 2002.8 1000000.0 +49 45.2682 3.172515503 1989.0 1000000.0 +50 46.1562 3.098299651 1989.0 1000000.0 +51 47.0403 3.025819958 1972.0 1000000.0 +52 47.9226 2.955035810 1980.9 1000000.0 +53 48.8009 2.885907542 1953.8 1000000.0 +54 49.6745 2.818396418 1959.8 1000000.0 +55 50.5444 2.752464606 1937.7 1000000.0 +56 51.4107 2.688075163 1943.3 1000000.0 +57 52.2756 2.625192006 1931.3 1000000.0 +58 53.1364 2.563779899 1924.9 1000000.0 +59 53.9919 2.503804428 1907.8 1000000.0 +60 54.8411 2.445231986 1896.4 1000000.0 +61 55.6879 2.388029752 1897.2 1000000.0 +62 56.5299 2.332165670 1874.9 1000000.0 +63 57.3662 2.277608439 1871.9 1000000.0 +64 58.1971 2.224327485 1850.7 1000000.0 +65 59.0233 2.172292953 1850.4 1000000.0 +66 59.8453 2.121475684 1832.3 1000000.0 +67 60.6627 2.071847203 1829.7 1000000.0 +68 61.4756 2.023379699 1812.0 1000000.0 +69 62.2791 1.976046014 1787.7 1000000.0 +70 63.0805 1.929819624 1802.3 1000000.0 +71 63.8754 1.884674625 1759.1 1000000.0 +72 64.6611 1.840585720 1760.6 1000000.0 +73 65.4419 1.797528204 1737.3 1000000.0 +74 66.2142 1.755477948 1722.9 1000000.0 +75 66.9785 1.714411390 1701.0 1000000.0 +76 67.7369 1.674305517 1696.8 1000000.0 +77 68.4864 1.635137857 1660.6 1000000.0 +78 69.2213 1.596886460 1632.0 1000000.0 +79 69.9472 1.559529892 1619.8 1000000.0 +80 70.6546 1.523047221 1549.2 1000000.0 diff --git a/examples/2023-Linf-lava/lava.cpp b/examples/2023-Linf-lava/lava.cpp new file mode 100644 index 00000000..4b7b862e --- /dev/null +++ b/examples/2023-Linf-lava/lava.cpp @@ -0,0 +1,229 @@ +// athena +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// application +#include +#include + +// canoe +#include +#include + +// climath +#include +#include + +// snap +#include + +// harp +#include + +// special includes +#include + +// utils +#include + +//Real grav, P0, T0, Tmin; +//int iH2O; +// +// +// +// +//void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { +// AllocateUserOutputVariables(4 + NVAPOR); +// SetUserOutputVariableName(0, "temp"); +// SetUserOutputVariableName(1, "theta"); +// SetUserOutputVariableName(2, "thetav"); +// SetUserOutputVariableName(3, "mse"); +// for (int n = 1; n <= NVAPOR; ++n) { +// std::string name = "rh" + std::to_string(n); +// SetUserOutputVariableName(3 + 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); +// for (int n = 1; n <= NVAPOR; ++n) +// user_out_var(3 + n, k, j, i) = +// pthermo->RelativeHumidity(this, n, 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); +// +// Real cv = pthermo->GetCvMass(air, 0); +// +// // 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)); +// // } +// +// // if (air.w[IDN] < Tmin) { +// // u(IEN,k,j,i) += w(IDN,k,j,i)*cv*(Tmin - temp)/sponge_tau*dt; +// // } +// } +//} +// +//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"); +// +// // index +// auto pindex = IndexMap::GetInstance(); +// iH2O = pindex->GetVaporId("H2O"); +// EnrollUserExplicitSourceFunction(Forcing); +//} + +void MeshBlock::ProblemGenerator(ParameterInput *pin){ // sets initial conditions + std::string input_atm_path = "/home/linfel/canoe_1/examples/2023-Linf-lava/balance_atm.txt"; + DataVector atm_data = read_data_vector(input_atm_path); + + // Access the data using keys + for (const auto& entry : atm_data) { + const std::string& key = entry.first; + const std::vector& values = entry.second; + + // Print the key (column header) + std::cout << key << ": "; + + // Print the values in the vector + for (double value : values) { + std::cout << value << " "; + } + + std::cout << std::endl; + } +} + + +//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; +// +// 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), +// Thermodynamics::Method::PseudoAdiabat, 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., +// Thermodynamics::Method::ReversibleAdiabat, grav); +// +// int i = is; +// for (; i <= ie; ++i) { +// if (air.w[IDN] < Tmin) break; +// air.w[IVX] = 1. * sin(2. * M_PI * rand() / RAND_MAX); +// AirParcelHelper::distribute_to_conserved(this, k, j, i, air); +// pthermo->Extrapolate(&air, pcoord->dx1f(i), +// Thermodynamics::Method::PseudoAdiabat, 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), +// Thermodynamics::Method::Isothermal, grav); +// } +// +// peos->ConservedToPrimitive(phydro->u, phydro->w, pfield->b, phydro->w, +// pfield->bcc, pcoord, is, ie, j, j, k, k); +// +// pimpl->prad->CalFlux(this, k, j, is, ie + 1); +// } +//} diff --git a/examples/2023-Linf-lava/lava.inp b/examples/2023-Linf-lava/lava.inp new file mode 100644 index 00000000..afad01bc --- /dev/null +++ b/examples/2023-Linf-lava/lava.inp @@ -0,0 +1,133 @@ + +problem = Cloud radiative dynamics + Hydrogen Atm + Ocean below +configure = -Dpnetcdf_path=/home3/xzhang11/pnetcdf + + +problem_id = hywater2d # problem ID: basename of output filenames + + +file_type = rst +dt = 20.E8 + + +file_type = hst # History data dump +dt = 1.E6 # time increment between outputs + + +file_type = pnetcdf # Binary data dump +variable = prim # variables to be output +dt = 1.E5 # time increment between outputs + + +file_type = pnetcdf +variable = uov +dt = 1.E5 + + +file_type = pnetcdf +variable = rad +dt = 1.E5 + +# +#file_type = pnetcdf +#variable = diag +#dt = 1.E5 + +