wip

Author: chengcli

examples/2024-XZhang-minichem/mini.inp

@@ -68,7 +68,8 @@ nx2         = 16
 nx3         = 1
 
 <hydro>
-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)
 <chemistry>
 microphysics_config = water_ammonia.yaml
 
-<thermodynamics>
-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
-
 <forcing>
 #packages = bot_heating, top_cooling
 #bot_heating.flux = 8.

examples/2024-XZhang-minichem/minichem.cpp

@@ -47,8 +47,6 @@ Real Tbot, xtra, ptra, Ttop_tau, Tbot_tau, Trabot_tau;
 std::string met;
 bool use_mini, use_mini_ic, use_tra_ic, fix_bot_tra;
 
-Real p0 = 1.e5;
-
 //'OH','H2','H2O','H','CO','CO2','O','CH4','C2H2','NH3','N2','HCN', 'He'
 std::vector<double> 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.};
@@ -81,7 +79,7 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
   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 theta = temp * pow(P0 / pres, chi);
 
   auto temp_a = temp.accessor<Real, 3>();
   auto theta_a = theta.accessor<Real, 3>();
@@ -316,20 +314,25 @@ torch::Tensor setup_moist_adiabatic_profile(std::string infile) {
 }
 
 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");
-  auto mud = kintera::species_weights[0];
-  auto Rd = kintera::constants::Rgas / mud;
+  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);
 
   // 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);
 
@@ -340,17 +343,19 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
         Real x1 = pcoord->x1v(i);
         Real x2 = pcoord->x2v(j);
         Real x3 = pcoord->x3v(k);
-        Real temp = Ts - grav * x1 / cp;
-        phydro->w(IPR, k, j, i) = p0 * pow(temp / Ts, cp / Rd);
+        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) = phydro->w(IVY, k, j, i) = 0.;
+        phydro->w(IVX, k, j, i) = 0.001 * (1. * rand() / RAND_MAX - 0.5);
       }
 
   auto temp = get_temp(pimpl->peos, phydro->w);
   auto temp_a = temp.accessor<Real, 3>();
 
   Tbot = temp_a[ks][js][is];
 
+  std::cout << "ptra = " << ptra << " Pa" << std::endl;
+
   if (NCHEM > 0) {
     if (use_mini_ic) {
       // minichem initialize conserved variables