Merge pull request #28 from yoctoyotta1024/add_pysdm

Add pysdm as git submodule

Author: yoctoyotta1024

cleo_1dkid/cleo_deps/CMakeLists.txt

@@ -23,12 +23,12 @@ include(FetchContent)
 FetchContent_Declare(
   cleo
   DOWNLOAD_EXTRACT_TIMESTAMP TRUE
-  URL https://github.com/yoctoyotta1024/CLEO/archive/refs/tags/v0.63.0.zip
-  GIT_TAG 5af4642cc6305d93577ad3e4c027890fbd52fded
+  URL https://github.com/yoctoyotta1024/CLEO/archive/refs/tags/v0.64.0.zip
+  GIT_TAG 4cf4985b741c098c12baa0c0261ae602a221c31a
 )
 FetchContent_MakeAvailable(cleo)
 
-message(STATUS "CLEO version v0.63.0, installation in: ${CLEO_BINARY_DIR}")
+message(STATUS "CLEO version v0.64.0, installation in: ${CLEO_BINARY_DIR}")
 
 # print default compiler flags
 message(STATUS "CLEO primary CMAKE_CXX_FLAGS: ${CMAKE_CXX_FLAGS}")

cleo_1dkid/libs/cleo_sdm/cleo_sdm.py

@@ -16,7 +16,7 @@
 -----
 File Description:
 class and functions for handing setup and running of CLEO via python bindings,
-functions adapted from CLEO v0.63.0 python_bindings example.
+functions adapted from CLEO v0.64.0 python_bindings example.
 NOTE: To use the wrapper, you must first export "CLEO_PYTHON_BINDINGS".
 E.g. if python bindings are built in $HOME/superdrops-in-action/build/, do:
 export CLEO_PYTHON_BINDINGS=$HOME/superdrops-in-action/build/_deps/cleo-build/cleo_python_bindings/

cleo_1dkid/libs/cleo_sdm/initconds/create_config.py

@@ -62,6 +62,11 @@
     type=float,
     help="alpha value for superdroplet initial conditions",
 )
+parser.add_argument(
+    "--numconc",
+    type=float,
+    help="number concentration /cm^3 for superdroplet initial conditions",
+)
 args = parser.parse_args()
 
 from cleopy import editconfigfile
@@ -74,6 +79,7 @@
 zarrbasedir = args.zarrbasedir
 nsupers_pergbx = args.nsupers_pergbx
 alpha = args.alpha
+numconc_perm3 = float(args.numconc) * 1e6  # [m^-3]
 
 cnfg = yaml.safe_load(open(args.src_config_filename))
 initnsupers = int(cnfg["domain"]["ngbxs"]) * nsupers_pergbx
@@ -101,6 +107,8 @@
     "maxnsupers": int(maxnsupers),
     "initnsupers": int(initnsupers),
     "newnsupers": int(newnsupers),
+    "numconc": int(numconc_perm3),
+    "NUMCONC_a": int(numconc_perm3),
 }
 
 print("--- create_config configuration arguments ---")

cleo_1dkid/libs/cleo_sdm/initconds/initial_pressure_profile.py

@@ -52,7 +52,7 @@ def get_initial_pressure_profile(
         is_approx_drhod_dz=False,
     )
     zfull = np.arange(z_min + z_delta / 2, z_max + z_delta / 2, z_delta)
-    press_prof = settings.press(zfull)
+    press_prof = settings.press0(zfull)
 
     gbxbounds = rgrid.read_dimless_gbxboundaries_binary(grid_filename, isprint=False)
 

cleo_1dkid/libs/test_case_1dkid/kid_dynamics.py

@@ -72,7 +72,7 @@ def __init__(
         self.mpdata = MPDATA(
             nz=self.settings.nz,
             dt=self.settings.dt,
-            qv_of_zZ_at_t0=lambda zZ: self.settings.qv(
+            qv_of_zZ_at_t0=lambda zZ: self.settings.qv0(
                 zZ * self.settings.dz + self.settings.z_min
             ),
             g_factor_of_zZ=lambda zZ: self.settings.rhod(
@@ -87,13 +87,9 @@ def __init__(
         assert self.settings.z_min == z_min
         assert self.settings.z_max == z_max
         assert self.settings.dz == z_delta
-        zfull = np.arange(z_min + z_delta / 2, z_max + z_delta / 2, z_delta)
+        self.zfull = np.arange(z_min + z_delta / 2, z_max + z_delta / 2, z_delta)
         self.zhalf = np.arange(z_min, z_max + z_delta, z_delta)
-        assert np.all((self.zhalf[1:] + self.zhalf[:-1]) / 2 == zfull)
-
-        self.rhod_prof = self.settings.rhod(zfull)
-        self.temp_prof = self.settings.temp(zfull)
-        self.press_prof = self.settings.press(zfull)
+        assert np.all((self.zhalf[1:] + self.zhalf[:-1]) / 2 == self.zfull)
 
         key = f"nr={self.mpdata.nr}, dz={self.settings.dz}, dt={self.settings.dt}, options={options}"
         print(f"Simulating {self.settings.nt} timesteps using {key}")
@@ -107,15 +103,15 @@ def set_thermo(self, time, thermo):
         Returns:
             Thermodynamics: Updated thermodynamic state.
         """
-        thermo.temp[:] = self.temp_prof
-        thermo.rho[:] = self.rhod_prof
-        thermo.press[:] = self.press_prof
         thermo.massmix_ratios["qvap"][:] = self.mpdata["qvap"].advectee.get()
         thermo.massmix_ratios["qcond"][:] = self.mpdata["qcond"].advectee.get()
         thermo.massmix_ratios["qice"][:] = self.mpdata["qice"].advectee.get()
         thermo.massmix_ratios["qrain"][:] = self.mpdata["qrain"].advectee.get()
         thermo.massmix_ratios["qsnow"][:] = self.mpdata["qsnow"].advectee.get()
         thermo.massmix_ratios["qgrau"][:] = self.mpdata["qgrau"].advectee.get()
+        thermo.rhod[:] = self.settings.rhod(self.zfull)
+        thermo.temp[:] = self.settings.temp(self.zfull)
+        thermo.press[:] = self.settings.press(self.zfull, thermo.massmix_ratios["qvap"])
 
         if thermo.wvel.size != 0:
             z_half_reps = np.repeat(

cleo_1dkid/libs/test_case_1dkid/settings.py

@@ -56,14 +56,14 @@ def __init__(
         self.z_max = z_max
         self.t_max = t_max
 
-        self.qv = interp1d(
+        self.qv0 = interp1d(
             (0, 740, 3260), (0.015, 0.0138, 0.0024), fill_value="extrapolate"
         )
         self._th = interp1d(
             (0, 740, 3260), (297.9, 297.9, 312.66), fill_value="extrapolate"
         )
 
-        self.thd = lambda z: formulae.th_dry(self._th(z), self.qv(z))
+        self.thd = lambda z, qv: formulae.th_dry(self._th(z), qv)
 
         z_points = np.arange(0, self.z_max + self.dz / 2, self.dz / 2)
 
@@ -104,7 +104,7 @@ def z2exner(z, return_rho=False, return_press=False, return_temp=False):
                     return exner
 
             self.temp = lambda z: z2exner(zpos(z), return_temp=True)
-            self.press = lambda z: z2exner(zpos(z), return_press=True)
+            self.press = lambda z, qv: z2exner(zpos(z), return_press=True)
             self.rhod = lambda z: z2exner(zpos(z), return_rho=True)
             if is_exner_novapour_uniformrho:
                 rhod0 = 1.0 * np.ones_like(z_points)  # [kg / m^3]
@@ -114,30 +114,31 @@ def z2exner(z, return_rho=False, return_press=False, return_temp=False):
             # note: not in the paper,
             # https://github.com/BShipway/KiD/tree/master/src/physconst.f90#L43
             def drhod_dz(z, rhod):
-                T = formulae.temperature(rhod[0], self.thd(z))
-                p = formulae.pressure(rhod[0], T, self.qv(z))
-                drhod_dz = formulae.drho_dz(const.g, p, T, self.qv(z), const.lv)
+                T = formulae.temperature(rhod[0], self.thd(z, self.qv0(z)))
+                p = formulae.pressure(rhod[0], T, self.qv0(z))
+                drhod_dz = formulae.drho_dz(const.g, p, T, self.qv0(z), const.lv)
                 if not is_approx_drhod_dz:  # to resolve issue #335
-                    qv = self.qv(z)
-                    dqv_dz = Derivative(self.qv)(z)
+                    qv = self.qv0(z)
+                    dqv_dz = Derivative(self.qv0)(z)
                     drhod_dz = drhod_dz / (1 + qv) - rhod * dqv_dz / (1 + qv) ** 2
                 return drhod_dz
 
-            rhod0 = formulae.rho_d(p_surf, self.qv(0), self._th(0))
+            rhod0 = formulae.rho_d(p_surf, self.qv0(0), self._th(0))
             rhod_solution = solve_ivp(
                 fun=drhod_dz,
                 t_span=(0, self.z_max),
                 y0=np.asarray((rhod0,)),
                 t_eval=z_points,
+                max_step=self.dz / 2,
             )
             assert rhod_solution.success
 
             self.rhod = lambda z: interp1d(z_points, rhod_solution.y[0])(zpos(z))
             self.temp = lambda z: formulae.temperature(
-                self.rhod(zpos(z)), self.thd(zpos(z))
+                self.rhod(zpos(z)), self.thd(zpos(z), self.qv0(zpos(z)))
             )
-            self.press = lambda z: formulae.pressure(
-                self.rhod(zpos(z)), self.temp(zpos(z)), self.qv(zpos(z))
+            self.press = lambda z, qv: formulae.pressure(
+                self.rhod(zpos(z)), self.temp(zpos(z)), qv
             )
 
         rhod_w_const = wmax_const * si.m / si.s * si.kg / si.m**3
@@ -147,6 +148,7 @@ def drhod_dz(z, rhod):
         )
 
         self.nz_vec = arakawa_c.z_vector_coord((self.nz,))
+        self.press0 = lambda z: self.press(z, self.qv0(zpos(z)))
 
     @property
     def nz(self):

cleo_1dkid/libs/thermo/output_thermodynamics.py

@@ -91,8 +91,8 @@ class OutputThermodynamics:
           time (s).
         temp (OutputVariable):
           Temperature (K).
-        rho (OutputVariable):
-          Density of moist air (kg/m3).
+        rhod (OutputVariable):
+          Density of dry air (kg/m3).
         press (OutputVariable):
           Pressure (Pa).
         qvap (OutputVariable):
@@ -121,7 +121,7 @@ def __init__(self, shape, zhalf=None, xhalf=None, yhalf=None):
 
         self.time = OutputVariable("time", "s", [shape[0]])
         self.temp = OutputVariable("temp", "K", shape)
-        self.rho = OutputVariable("rho", "kg m-3", shape)
+        self.rhod = OutputVariable("rhod", "kg m-3", shape)
         self.press = OutputVariable("press", "Pa", shape)
         self.qvap = OutputVariable("qvap", "kg/kg", shape)
         self.qcond = OutputVariable("qcond", "kg/kg", shape)
@@ -165,7 +165,7 @@ def output_thermodynamics(self, time: float, thermo: Thermodynamics):
         """
         self.time.write(time)
         self.temp.write(thermo.temp)
-        self.rho.write(thermo.rho)
+        self.rhod.write(thermo.rhod)
         self.press.write(thermo.press)
         self.qvap.write(thermo.massmix_ratios["qvap"])
         self.qcond.write(thermo.massmix_ratios["qcond"])
@@ -212,7 +212,7 @@ def finalize(self):
         """
         self.time.finalize()
         self.temp.finalize()
-        self.rho.finalize()
+        self.rhod.finalize()
         self.press.finalize()
         self.qvap.finalize()
         self.qcond.finalize()

cleo_1dkid/libs/thermo/thermodynamics.py

@@ -31,8 +31,8 @@ class Thermodynamics:
     Parameters:
       temp (np.ndarray):
         Temperature (K).
-      rho (np.ndarray):
-        Density of moist air (kg/m3).
+      rhod (np.ndarray):
+        Density of dry air (kg/m3).
       press (np.ndarray):
         Pressure (Pa).
       qvap (np.ndarray):
@@ -51,8 +51,8 @@ class Thermodynamics:
     Attributes:
       temp (np.ndarray):
         Temperature (K).
-      rho (np.ndarray):
-        Density of moist air (kg/m3).
+      rhod (np.ndarray):
+        Density of dry air (kg/m3).
       press (np.ndarray):
         Pressure (Pa).
       massmix_ratios (dict):
@@ -71,7 +71,7 @@ class Thermodynamics:
     def __init__(
         self,
         temp: np.ndarray,
-        rho: np.ndarray,
+        rhod: np.ndarray,
         press: np.ndarray,
         qvap: np.ndarray,
         qcond: np.ndarray,
@@ -88,7 +88,7 @@ def __init__(
         Parameters:
             press (np.ndarray): Pressure (Pa).
             temp (np.ndarray): Temperature (K).
-            rho (np.ndarray): Density of moist air (kg/m3)
+            rhod (np.ndarray): Density of moist air (kg/m3)
             qvap (np.ndarray): Specific water vapor content (kg/kg)
             qcond (np.ndarray): Specific cloud water content (kg/kg)
             qice (np.ndarray): Specific cloud ice content (kg/kg)
@@ -100,7 +100,7 @@ def __init__(
             vvel (np.ndarray): meridional wind velocity, 'y', (m/s)
         """
         self.temp = deepcopy(temp)
-        self.rho = deepcopy(rho)
+        self.rhod = deepcopy(rhod)
         self.press = deepcopy(press)
         self.massmix_ratios = {
             "qvap": deepcopy(qvap),
@@ -116,7 +116,7 @@ def __init__(
 
     def print_state(self):
         print(self.temp)
-        print(self.rho)
+        print(self.rhod)
         print(self.press)
         print(self.massmix_ratios)
 

cleo_1dkid/scripts/bash/inputfiles_cleo_1dkid.sh

@@ -24,13 +24,14 @@ isfigures=FALSE
 python=/work/bm1183/m300950/bin/envs/superdrops-in-action/bin/python
 path2initcondsscripts=${path2cleo1dkid}/libs/cleo_sdm/initconds
 
+runtype=condevap_only
 nsupers_pergbxs=(256) # for superdroplet initial conditions
-alphas=(0 0.5 1.0) # for superdroplet initial conditions alpha sampling
+alphas=(0.5) # for superdroplet initial conditions alpha sampling
+numconc=150
 
 ### src_configs is list of absolute paths to source config files seperated by spaces
 ### e.g. src_configs=("$HOME/config1" "$HOME/config2"), following lists are
-src_configs=("${path2cleo1dkid}/share/cleo_initial_conditions/1dkid/condevap_only/config.yaml"
-      "${path2cleo1dkid}/share/cleo_initial_conditions/1dkid/fullscheme/config.yaml")
+src_configs=("${path2cleo1dkid}/share/cleo_initial_conditions/1dkid/${runtype}/config.yaml")
 
 ### IDs of ensemble members (diff superdrop initial conditions) for each src_configs to create
 if [[ "${start_id}" == "" || "${end_id}" == "" ]]
@@ -45,13 +46,11 @@ cleoconstants_filepath="${path2build}/_deps/cleo-src/libs/"
 grid_filename="${path2build}/share/dimlessGBxboundaries.dat"
 
 ### same files for all src_configs, different for run_ids
-initsupers_directory="${path2build}/share"
+initsupers_directory="${path2build}/share/share_${numconc}cm3"
 
 ### different for all src_configs and different for run_ids
-configs_directory=("${path2build}/tmp/condevap_only"
-                 "${path2build}/tmp/fullscheme")
-bin_directory=("${path2build}/bin/condevap_only"
-                "${path2build}/bin/fullscheme")
+configs_directory=("${path2build}/tmp_${numconc}cm3/${runtype}")
+bin_directory=("${path2build}/bin_${numconc}cm3/${runtype}")
 ### ---------------------------------------------------- ###
 ### ---------------------------------------------------- ###
 ### ---------------------------------------------------- ###
@@ -104,7 +103,7 @@ do
         initsupers_filename="${initsupers_directory}/dimlessSDsinit_${label}.dat"
         setup_filename="${bin_directory[i]}/setup_${label}.txt"
         zarrbasedir="${bin_directory[i]}/sol_${label}.zarr"
-        echo "---- src ${i}, run number: ${m} ----"
+        echo "---- src ${i}, run number: ${m}, numconc: ${numconc}cm^-3  ----"
         echo "---- nsupers ${nsupers_pergbxs[k]}, alpha ${alphas[l]} ----"
         echo "path to build directory: ${path2build}"
         echo "src config file: ${src_configfile}"
@@ -115,6 +114,7 @@ do
         echo "-- ${zarrbasedir}"
         echo "-- ${nsupers_pergbxs[k]}"
         echo "-- ${alphas[l]}"
+        echo "-- ${numconc}"
 
         echo "python create_config.py ${src_configfile} ${dest_configfile} [...]"
         ${python} ${path2initcondsscripts}/create_config.py ${src_configfile} ${dest_configfile} \
@@ -124,7 +124,8 @@ do
           --setup_filename="${setup_filename}" \
           --zarrbasedir="${zarrbasedir}" \
           --nsupers_pergbx="${nsupers_pergbxs[k]}" \
-          --alpha="${alphas[l]}"
+          --alpha="${alphas[l]}" \
+          --numconc="${numconc}"
       done
     done
   done

cleo_1dkid/scripts/bash/run_cleo_1dkid.sh

@@ -4,8 +4,7 @@
 #SBATCH --nodes=1
 #SBATCH --ntasks-per-node=1
 #SBATCH --cpus-per-task=128
-#SBATCH --mem=20G
-#SBATCH --time=00:60:00
+#SBATCH --time=02:00:00
 #SBATCH --mail-user=clara.bayley@mpimet.mpg.de
 #SBATCH --mail-type=FAIL
 #SBATCH --account=bm1183
@@ -26,17 +25,16 @@ path2cleopythonbindings="${path2build}/_deps/cleo-build/cleo_python_bindings"
 python="/work/bm1183/m300950/bin/envs/superdrops-in-action/bin/python"
 pythonlibs="/work/bm1183/m300950/bin/envs/superdrops-in-action/lib/python3.13/site-packages"
 
+runtype=condevap_only
 nsupers_pergbxs=(256) # for superdroplet initial conditions
-alphas=(0 0.5 1.0) # for superdroplet initial conditions alpha sampling
+alphas=(0.5) # for superdroplet initial conditions alpha sampling
+numconc=150
 
 ### loop over configs_directory for all different for run_ids
-configs_directory=("${path2build}/tmp/condevap_only"
-                 "${path2build}/tmp/fullscheme")
-run_labels=("condevap_only" "fullscheme")
-bin_directory=("${path2build}/bin/condevap_only"
-                "${path2build}/bin/fullscheme") # Note! must match paths in config's 'outputdata'
-fig_directory=("${path2build}/bin/condevap_only"
-                "${path2build}/bin/fullscheme")
+configs_directory=("${path2build}/tmp_${numconc}cm3/${runtype}")
+run_labels=("${runtype}")
+bin_directory=("${path2build}/bin_${numconc}cm3/${runtype}") # Note! must match paths in config's 'outputdata'
+fig_directory=("${path2build}/bin_${numconc}cm3/${runtype}")
 
 ### IDs of ensemble members (diff superdrop initial conditions) for each src_configs to create
 if [[ "${start_id}" == "" || "${end_id}" == "" ]]
@@ -93,7 +91,7 @@ do
         run_name="${run_labels[i]}_${label}"
         binpath="${bin_directory[i]}"
         figpath="${fig_directory[i]}"
-        echo "---- src ${i}, run number: ${m} ----"
+        echo "---- src ${i}, run number: ${m}, numconc: ${numconc}cm^-3 ----"
         echo "---- nsupers ${nsupers_pergbxs[k]}, alpha ${alphas[l]} ----"
         echo "--run_name=${run_name}"
         echo "--config_filename=${config_filename}"