diff --git a/.gitignore b/.gitignore index 3c5e416b..03530987 100644 --- a/.gitignore +++ b/.gitignore @@ -28,3 +28,9 @@ tests/Testing # external external/*.tar.gz + +# swap files +*.swp + +# log files +*.log diff --git a/cmake/lava.cmake b/cmake/lava.cmake new file mode 100644 index 00000000..995b9b69 --- /dev/null +++ b/cmake/lava.cmake @@ -0,0 +1,22 @@ +# configure file for lava 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/cmake/lstraka.cmake b/cmake/lstraka.cmake new file mode 100644 index 00000000..848b186c --- /dev/null +++ b/cmake/lstraka.cmake @@ -0,0 +1,22 @@ +# configuration for lstraka hydrodynamcis + +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/2019-Li-snap/CMakeLists.txt b/examples/2019-Li-snap/CMakeLists.txt index de44b373..9856e18c 100644 --- a/examples/2019-Li-snap/CMakeLists.txt +++ b/examples/2019-Li-snap/CMakeLists.txt @@ -4,6 +4,7 @@ # 1. Compile straka, robert and bryan setup_problem(straka) +setup_problem(lstraka) setup_problem(robert) setup_problem(bryan) diff --git a/examples/2019-Li-snap/balance_atm.txt b/examples/2019-Li-snap/balance_atm.txt new file mode 100644 index 00000000..1171dd98 --- /dev/null +++ b/examples/2019-Li-snap/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/2019-Li-snap/lstraka.cpp b/examples/2019-Li-snap/lstraka.cpp new file mode 100644 index 00000000..049e4cc1 --- /dev/null +++ b/examples/2019-Li-snap/lstraka.cpp @@ -0,0 +1,325 @@ +/* ------------------------------------------------------------------------------------- + * SNAP Example Program + * + * Contributer: + * Cheng Li, University of Michigan + * + * Year: 2021 + * Contact: chengcli@umich.edu + * Reference: Straka et al., 1993 + * ------------------------------------------------------------------------------------- + */ + +// @sect3{Include files} + +// First we include some Athena++ headers files. They contain declarations +// of various components of the solver system. +// Athena++ is able to run both single precision (float) and double precision +// (double) applications. A macro Real is defined to indicate the +// precision and it is define in the following header file. +#include +// Then the basic multi-dimension array that stores fluid dynamic data is +// declared in the AthenaArray class. +#include +#include +// The model executes according to the parameters specified in an input file. +// The input file and the parameters within are managed by the ParameterInput +// class. +#include +// Coordinate related information are stored in the Coordinates class. +#include +// Everying regarding the equation of state is treated by the EquationOfState +// class and declared in the following file. +#include +// Evolution of hydrodynamic fields like pressure, density and velocities are +// managed by the Hydro class. +#include +// Athena++ can alo simulate magnetohydrodynamics (MHD). +// However, since this example is a hydro-only application, the MHD file is only +// included as a place holder. +#include +// Dynamic fields are evolving on a mesh and this is the file working with +// meshes and the partition of meshes among processors +#include + +// The following header file contains the definition of the MeshBlock::Impl +// class +#include +#include + +// Finally, the Thermodynamics class works with thermodynamic aspects of the +// problem such as the temperature, potential temperature, condensation of +// vapor, etc. +#include + +// harp +//#include + +// Functions in the math library are protected by a specific namespace because +// math functions are usually concise in the names, such as min, +// sqr. It is better to protect them in a namespace to avoid +// conflicts. +#include +#include + +// special includes +//#include + +// utils +#include // read_data_vector + +// Now, we get to the real program. +// First, we define several global variables specific to this appliciation +// For example, here K is the kinematic visocisty, p0 +// is the surface pressure, cp is the specific heat capacity, and +// Rd is the ideal gas constant of dry air. The Prandtl number is 1. +Real K, p0, cp, Rd; +int iSiO; + +// @sect3{User-defined output variables} + +// The hydrodynamic solver evolves density, pressure and velocities with time, +// meaning that the (potential) temperature is a diagnostic quantity for output +// only. This block of code allocates the memory for the outputs of temperature +// and potential temperature. +void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) { + AllocateUserOutputVariables(2 + NVAPOR); + SetUserOutputVariableName(0, "temp"); + SetUserOutputVariableName(1, "theta"); + for (int n = 1; n <= NVAPOR; ++n) { + std::string name = "rh" + std::to_string(n); + SetUserOutputVariableName(1 + n, name.c_str()); + } +} + +// Since temperature and potential temperature are only used for output purpose, +// they do not need to be calculated every dynamic time step. The subroutine +// below loops over all grids and calculates the value of temperature and +// potential temperature before the output time step. Particularly, the pointer +// to the Thermodynamics class pthermo calculates the temperature +// and potential temperature using its own member function +// Thermodynamics::GetTemp and PotentialTemp. pthermo, is a member +// of a higher level management class MeshBlock. So, you can use it directly +// inside a member function of class MeshBlock. In fact, all physics modules are +// managed by MeshBlock. As long as you have a pointer to a MeshBlock, you can +// access all physics in the simulation. +void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) { + auto pthermo = Thermodynamics::GetInstance(); + + // Loop over the grids. js,je,is,ie are members of the MeshBlock + // class. They are integer values representing the start index and the end + // index of the grids in each dimension. + int k = ks; + for (int j = js; j <= je; ++j) + for (int i = is; i <= ie; ++i) { + // user_out_var stores the actual data. + // phydro is a pointer to the Hydro class, which has a member + // w that stores density, pressure, and velocities at each + // grid. + user_out_var(0, j, i) = pthermo->GetTemp(this, k, j, i); + user_out_var(1, j, i) = pthermo->PotentialTemp(this, p0, k, j, i); + for (int n = 1; n <= NVAPOR; ++n){ + user_out_var(1 + n, k, j, i) = pthermo->RelativeHumidity(this, n, k, j, i); + } + } +} + +// @sect3{Forcing function} + +// The sinking bubble is forced by the gravity and the viscous and thermal +// dissipation. The gravitational forcing is taken care of by other parts of the +// program and we only need to write a few lines of code to facilitate +// dissipation. The first step is to define a function that takes a pointer to +// the MeshBlock as an argument such that we can access all physics via this +// pointer. The name of this function is not important. It can be anything. But +// the order and types of the arguments must be (MeshBlock *, Real const, +// Real const, AthenaArray const&, AthenaArray const&, +// AthenaArray &). They are called the signature of the function. +void Diffusion(MeshBlock *pmb, Real const time, Real const dt, + AthenaArray const &w, AthenaArray const &r, + AthenaArray const &bcc, AthenaArray &u, + AthenaArray &s) { + // pcoord is a pointer to the Coordinates class and it is a + // member of the MeshBlock class. We use the pointer to the MeshBlock class, + // pmb to access pcoord and use its member function + // to get the spacing of the grid. + Real dx = pmb->pcoord->dx1f(pmb->is); + Real dy = pmb->pcoord->dx2f(pmb->js); + + // Similarly, we use pmb to find the pointer to the + // Thermodynamics class, pthermo. + auto pthermo = Thermodynamics::GetInstance(); + + // Loop over the grids. + int k = pmb->ks; + for (int j = pmb->js; j <= pmb->je; ++j) + for (int i = pmb->is; i <= pmb->ie; ++i) { + // Similar to what we have done in MeshBlock::UserWorkBeforeOutput, we use + // the Thermodynamics class to calculate temperature and potential + // temperature. + Real temp = pthermo->GetTemp(pmb, pmb->ks, j, i); + Real theta = pthermo->PotentialTemp(pmb, p0, pmb->ks, j, i); + + // The thermal diffusion is applied to the potential temperature field, + // which is not exactly correct. But this is the setting of the test + // program. + Real theta_ip1_j = pthermo->PotentialTemp(pmb, p0, k, j + 1, i); + Real theta_im1_j = pthermo->PotentialTemp(pmb, p0, k, j - 1, i); + Real theta_i_jp1 = pthermo->PotentialTemp(pmb, p0, k, j, i + 1); + Real theta_i_jm1 = pthermo->PotentialTemp(pmb, p0, k, j, i - 1); + + // Add viscous dissipation to the velocities. Now you encounter another + // variable called u. u stands for `conserved + // variables`, which are density, momentums and total energy (kinetic + + // internal). In contrast, w stands for `primitive + // variables`, which are density, velocities, and pressure. Their relation + // is handled by the EquationOfState class. The `conserved variables` are + // meant for internal calculation. Solving for `conserved variables` + // guarantees the conservation properties. However, `conserved variables` + // are not easy to use for diagnostics. Therefore, another group of + // variables called the `primitive variables` are introduced for external + // calculations, like calculating transport fluxes, radiative fluxes and + // interacting with other physical components of the system. In this case, + // the diffusion is calculated by using the `primitive variabes` and the + // result is added to the `conserved variables` to ensure conservation + // properties. + u(IM1, j, i) += dt * K * w(IDN, j, i) / (dx * dy) * + (w(IVX, j, i - 1) + w(IVX, j, i + 1) + w(IVX, j - 1, i) + + w(IVX, j + 1, i) - 4. * w(IVX, j, i)); + u(IM2, j, i) += dt * K * w(IDN, j, i) / (dx * dy) * + (w(IVY, j, i - 1) + w(IVY, j, i + 1) + w(IVY, j - 1, i) + + w(IVY, j + 1, i) - 4. * w(IVY, j, i)); + + // Adding thermal dissipation is similar to adding viscous dissipation. + u(IEN, j, i) += + dt * K * w(IDN, j, i) / (dx * dy) * cp * temp / theta * + (theta_ip1_j + theta_im1_j + theta_i_jp1 + theta_i_jm1 - 4. * theta); + } +} + +// @sect3{Program specific variables} + +// This is the place where program specific variables are initialized. +// Not that the function is a member function of the Mesh class rather than the +// MeshBlock class we have been working with. The difference between class Mesh +// and class MeshBlock is that class Mesh is an all-encompassing class that +// manages multiple MeshBlocks while class MeshBlock manages all physics +// modules. During the instantiation of the classes. class Mesh is instantiated +// first and then it instantiates all MeshBlocks inside it. Therefore, this +// subroutine runs before any MeshBlock. +void Mesh::InitUserMeshData(ParameterInput *pin) { + // The program specific forcing parameters are set here. + // They come from the input file, which is parsed by the ParameterInput class + Real gamma = pin->GetReal("hydro", "gamma"); + K = pin->GetReal("problem", "K"); + p0 = pin->GetReal("problem", "p0"); + Rd = pin->GetReal("thermodynamics", "Rd"); + cp = gamma / (gamma - 1.) * Rd; + + // index + auto pindex = IndexMap::GetInstance(); + iSiO = pindex->GetVaporId("SiO"); + // This line code enrolls the forcing function we wrote in + // section Forcing function + EnrollUserExplicitSourceFunction(Diffusion); +} + +// @sect3{Initial condition} + +// This is the final part of the program, in which we set the initial condition. +// It is called the `problem generator` in a general Athena++ application. +void MeshBlock::ProblemGenerator(ParameterInput *pin) { + + auto pthermo = Thermodynamics::GetInstance(); + + // Get the value of gravity. Positive is upward and negative is downward. + // We take the negative to get the absolute value. + //Real grav = -phydro->hsrc.GetG1(); + Real grav = -pin->GetReal("hydro", "grav_acc1"); + // These lines read the parameter values in the input file, which + // are organized into sections. The problem specific parameters are usually + // placed under section `problem`. + //Real Ts = pin->GetReal("problem", "Ts"); + Real xc = pin->GetReal("problem", "xc"); + Real xr = pin->GetReal("problem", "xr"); + Real zc = pin->GetReal("problem", "zc"); + Real zr = pin->GetReal("problem", "zr"); + Real dT = pin->GetReal("problem", "dT"); + + // read initial atmosphere conditions + std::string input_atm_path = "/home/linfel/canoe_1/examples/2019-Li-snap/balance_atm.txt"; + DataVector atm = read_data_vector(input_atm_path); + size_t atm_size = atm["HGT"].size(); + + // convert km -> m, hPa -> Pa, ppm -> fraction + double *hgtptr = atm["HGT"].data(); + double *preptr = atm["PRE"].data(); + double *sioptr = atm["SiO"].data(); + for (int i = 0; i < atm_size; ++i) { + hgtptr[i] *= 1000.0; + preptr[i] *= 100.0; + sioptr[i] /= 1.E8; + } + + AirParcel air(AirParcel::Type::MoleFrac); + + // Loop over the grids. The purpose is to set the temperature, pressure and + // density fields at each cell-centered grid. + for (int k = ks; k <= ke; ++k) + for (int j = js; j <= je; ++j) { + // half a grid to cell center + pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., + Thermodynamics::Method::ReversibleAdiabat, grav); + + for (int i = is; i <= ie; ++i) { + // Get the Cartesian coordinates in the vertical and horizontal + // directions. `x1` is usually the vertical direction and `x2` is the + // horizontal direction. The meaning of `x1` and `x2` may change with the + // coordinate system. x1v and x2v retrieve the + // coordinate at cell centers. + Real x1 = pcoord->x1v(i); + Real x2 = pcoord->x2v(j); + // Distance to the center of a cold air bubble at (xc,zc). + Real L = sqrt(sqr((x2 - xc) / xr) + sqr((x1 - zc) / zr)); + // Adiabatic temperature gradient. + air.SetZero(); + //Real temp = Ts - grav * x1 / cp; + Real temp = interp1(x1, atm["TEM"].data(), + atm["HGT"].data(), atm_size); + // Once we know the temperature, we can calculate the adiabatic and + // hydrostatic pressure as $p=p_0(\frac{T}{T_s})^{c_p/R_d}$, where $p_0$ + // is the surface pressure and $T_s$ is the surface temperature. + //phydro->w(IPR, j, i) = p0 * pow(temp / Ts, cp / Rd); + air.w[IPR] = interp1(x1, atm["PRE"].data(), + atm["HGT"].data(), atm_size); + + if (L <= 1.) temp += dT * (cos(M_PI * L) + 1.) / 2.; + air.w[IDN] = temp; + // Set density using ideal gas law. + //air.w[IDN] = air.w[IPR] / (Rd * temp); + air.w[iSiO] = interp1(x1, atm["SiO"].data(), + atm["HGT"].data(), atm_size); + + // Initialize velocities to zero. + air.w[IVX] = air.w[IVY] = 0.; + //phydro->w(IVX, j, i) = phydro->w(IVY, j, i) = 0.; + + AirParcelHelper::distribute_to_conserved(this, k, j, i, air); + pthermo->Extrapolate(&air, pcoord->dx1f(i), + Thermodynamics::Method::PseudoAdiabat, grav, + 1.e-5); + } + 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); + } + + // We have set all `primitive variables`. The last step is to calculate the + // `conserved variables` based on the `primitive variables`. This is done by + // calling the member function EquationOfState::PrimitiveToConserved. Note + // that pfield->bcc is a placeholder because there is not MHD + // involved in this example. + //peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie, + // js, je, ks, ke); +} diff --git a/examples/2019-Li-snap/lstraka.inp b/examples/2019-Li-snap/lstraka.inp new file mode 100644 index 00000000..12fe0cef --- /dev/null +++ b/examples/2019-Li-snap/lstraka.inp @@ -0,0 +1,85 @@ + +problem = Dense sinking bubble test +reference = Straka et al., 1993 +configure = --prob=straka --nghost=3 -netcdf + + +problem_id = straka # problem ID: basename of output filenames + + +file_type = hst # History data dump +dt = 300. # time increment between outputs + + +file_type = netcdf # Binary data dump +variable = prim # variables to be output +dt = 300. # time increment between outputs + + +file_type = netcdf +variable = uov +dt = 300. + +