More analysis breakup partii

docs/source/testcases/collisions.rst

@@ -84,6 +84,11 @@ To compare the results of the collision test cases including breakup, first buil
 breakup excutables ``colls_testikstraub``, ``colls_straub_schlottke``, ``colls_straub_fixednfrags`` and
 ``constcoalbu_fixednfrags`` as above, as well as the long executable ``colls_long``.
 
+To run and plot the breakup comparison with the marshall palmer distribution you first need to
+set ``use_marshall_parmer`` and ``plot_marshall_parmer``
+to ``True`` in the scripts for ``initconds_compare_breakup.py`` and ``plot_compare_breakup``,
+and in ``run_compare_breakup.sh`` you need to set ``label=bucomp_marshpam``
+
 The generate the initial conditons similarly but with the ``initconds_compare_breakup.py`` script
 and ``config_compare_breakup.yaml``.
 E.g.
@@ -113,6 +118,6 @@ E.g.
 
   $ python ./scripts/collisions/plot_compare_breakup.py`` \
       /home/m/m300950/CLEO \
-      /work/mh0731/m300950/validating_cleo/build/share/collisions/dimlessGBxboundaries_bucomp.dat \
+      /work/mh0731/m300950/validating_cleo/build/share/collisions/dimlessGBxboundaries_bucomp[_marshpam].dat \
       /work/mh0731/m300950/validating_cleo/build/bin/collisions \
       /work/mh0731/m300950/validating_cleo/build/bin/collisions

scripts/collisions/initconds_compare_breakup.py

@@ -36,7 +36,13 @@
 import collisions.initconds_compare_breakup as iccb
 
 
-def main(path2pySD, path2build, original_config, isfigures=[False, False]):
+def main(
+    path2pySD,
+    path2build,
+    original_config,
+    use_marshall_parmer=False,
+    isfigures=[False, False],
+):
     if not (
         path2build.is_dir()
         and all([(path2build / f"{d}").is_dir() for d in ["tmp", "share", "bin"]])
@@ -50,7 +56,10 @@ def main(path2pySD, path2build, original_config, isfigures=[False, False]):
             path2dir.mkdir(exist_ok=True)
 
     config_filenames = iccb.generate_configurations(
-        path2pySD, path2build, original_config
+        path2pySD,
+        path2build,
+        original_config,
+        use_marshall_parmer,
     )
 
     iccb.gridbox_boundaries(path2pySD, config_filenames[0], isfigures=isfigures)
@@ -93,4 +102,12 @@ def main(path2pySD, path2build, original_config, isfigures=[False, False]):
     if args.is_savefigs == "TRUE":
         isfigures[1] = True
 
-    main(path2pySD, path2build, original_config, isfigures=isfigures)
+    use_marshall_parmer = False
+
+    main(
+        path2pySD,
+        path2build,
+        original_config,
+        use_marshall_parmer=use_marshall_parmer,
+        isfigures=isfigures,
+    )

scripts/collisions/plot_compare_breakup.py

@@ -30,28 +30,43 @@
 import collisions.plot_compare_breakup as breakupplt
 
 
-def main(path2CLEO, grid_filename, path2bin, path4figs):
-    datasets = {}
-    setups = {}
-    for r in range(9):
-        datasets[r] = path2bin / f"sol_bucomp_{r}.zarr"
-        setups[r] = path2bin / f"setup_bucomp_{r}.txt"
-
-    fig1, fig2, fig3 = breakupplt.plot_compare_breakup(
-        path2CLEO, grid_filename, datasets, setups
-    )
+def main(path2CLEO, grid_filename, path2bin, path4figs, plot_marshall_parmer=False):
+    if not plot_marshall_parmer:
+        datasets = {}
+        setups = {}
+        for r in range(9):
+            datasets[r] = path2bin / f"sol_bucomp_{r}.zarr"
+            setups[r] = path2bin / f"setup_bucomp_{r}.txt"
+
+        fig1, fig2, fig3 = breakupplt.plot_compare_breakup(
+            path2CLEO, grid_filename, datasets, setups
+        )
+
+        savename = path4figs / "dejong2023_fig7a_as_timeseries.png"
+        fig1.savefig(savename, dpi=400, bbox_inches="tight", facecolor="w")
+        print("Figure .png saved as: " + str(savename))
 
-    savename = path4figs / "dejong2023_fig7a_as_timeseries.png"
-    fig1.savefig(savename, dpi=400, bbox_inches="tight", facecolor="w")
-    print("Figure .png saved as: " + str(savename))
+        savename = path4figs / "dejong2023_fig7a.pdf"
+        fig2.savefig(savename, dpi=400, bbox_inches="tight", facecolor="w")
+        print("Figure .pdf saved as: " + str(savename))
 
-    savename = path4figs / "dejong2023_fig7a.pdf"
-    fig2.savefig(savename, dpi=400, bbox_inches="tight", facecolor="w")
-    print("Figure .pdf saved as: " + str(savename))
+        savename = path4figs / "dejong2023_fig8.png"
+        fig3.savefig(savename, dpi=400, bbox_inches="tight", facecolor="w")
+        print("Figure .png saved as: " + str(savename))
+    else:
+        datasets = {}
+        setups = {}
+        for r in [2, 3, 4, 5]:
+            datasets[r] = path2bin / f"sol_bucomp_marshpam_{r}.zarr"
+            setups[r] = path2bin / f"setup_bucomp_marshpam_{r}.txt"
 
-    savename = path4figs / "dejong2023_fig8.pdf"
-    fig3.savefig(savename, dpi=400, bbox_inches="tight", facecolor="w")
-    print("Figure .pdf saved as: " + str(savename))
+        fig4 = breakupplt.plot_compare_marshall_parmer_breakup(
+            path2CLEO, grid_filename, datasets, setups
+        )
+
+        savename = path4figs / "dejong2023_fig9.pdf"
+        fig4.savefig(savename, dpi=400, bbox_inches="tight", facecolor="w")
+        print("Figure .pdf saved as: " + str(savename))
 
 
 if __name__ == "__main__":
@@ -68,4 +83,12 @@ def main(path2CLEO, grid_filename, path2bin, path4figs):
     )
     args = parser.parse_args()
 
-    main(args.path2CLEO, args.grid_filename, args.path2bin, args.path4figs)
+    plot_marshall_parmer = False
+
+    main(
+        args.path2CLEO,
+        args.grid_filename,
+        args.path2bin,
+        args.path4figs,
+        plot_marshall_parmer=plot_marshall_parmer,
+    )

scripts/collisions/run_compare_breakup.sh

@@ -21,13 +21,16 @@ buildtype=${2:-openmp}                                          # "serial", "thr
 
 SCRIPT_DIR=/home/m/m300950/validating_cleo/scripts
 runscript=${SCRIPT_DIR}/run.sh
+
+label="bucomp"
+# label="bucomp_marshpam"  # only need runs 2,3,4 and 5 (straub efficiency)
 ### ---------------------------------------------------- ###
 
 ### ------------ run long executable ---------- ###
 for run in 0; do
   runcmd="${runscript} ${path2build} ${buildtype} \
             /work/mh0731/m300950/validating_cleo/build/collisions/colls_long \
-            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_bucomp_${run}.yaml"
+            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_${label}_${run}.yaml"
   echo ${runcmd}
   eval ${runcmd}
 done;
@@ -37,7 +40,7 @@ done;
 for run in 1; do
   runcmd="${runscript} ${path2build} ${buildtype} \
             /work/mh0731/m300950/validating_cleo/build/collisions/colls_testikstraub \
-            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_bucomp_${run}.yaml"
+            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_${label}_${run}.yaml"
   echo ${runcmd}
   eval ${runcmd}
 done;
@@ -47,7 +50,7 @@ done;
 for run in 2; do
   runcmd="${runscript} ${path2build} ${buildtype} \
             /work/mh0731/m300950/validating_cleo/build/collisions/colls_straub_schlottke \
-            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_bucomp_${run}.yaml"
+            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_${label}_${run}.yaml"
   echo ${runcmd}
   eval ${runcmd}
 done;
@@ -57,7 +60,7 @@ done;
 for run in 3 4 5; do
   runcmd="${runscript} ${path2build} ${buildtype} \
             /work/mh0731/m300950/validating_cleo/build/collisions/colls_straub_fixednfrags \
-            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_bucomp_${run}.yaml"
+            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_${label}_${run}.yaml"
   echo ${runcmd}
   eval ${runcmd}
 done;
@@ -67,7 +70,7 @@ done;
 for run in 6 7 8; do
   runcmd="${runscript} ${path2build} ${buildtype} \
             /work/mh0731/m300950/validating_cleo/build/collisions/colls_constcoalbu_fixednfrags \
-            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_bucomp_${run}.yaml"
+            /work/mh0731/m300950/validating_cleo/build/tmp/collisions/config_${label}_${run}.yaml"
   echo ${runcmd}
   eval ${runcmd}
 done;

src/collisions/config_compare_breakup.yaml

@@ -28,6 +28,7 @@ python_initconds:
     dryradius: 1.0e-16                                        # all SDs have negligible solute [m]
     volexpr0: XXX                                             # peak of volume exponential distribution [m]
     numconc: XXX                                              # = 2^23 = total no. conc of real droplets [m^-3]
+    use_marshall_parmer: XXX                                  # use Marshall-Palmer distribution instead
 
 ### SDM Runtime Parameters ###
 domain:

src/collisions/initconds_compare_breakup.py

@@ -22,7 +22,9 @@
 """
 
 
-def generate_configurations(path2pySD, path2build, original_config):
+def generate_configurations(
+    path2pySD, path2build, original_config, use_marshall_parmer
+):
     import shutil
     import sys
 
@@ -37,14 +39,22 @@ def generate_configurations(path2pySD, path2build, original_config):
     constants_filename = (
         path2build / "_deps" / "cleo-src" / "libs" / "cleoconstants.hpp"
     )
-    grid_filename = sharepath / "dimlessGBxboundaries_bucomp.dat"
+
+    if use_marshall_parmer:
+        grid_filename = sharepath / "dimlessGBxboundaries_bucomp_marshpam.dat"
+    else:
+        grid_filename = sharepath / "dimlessGBxboundaries_bucomp.dat"
 
     # parameters specific for each run
     nruns = (
         9  # [long, testikstraub+schlottke, straub+schlottke, 3x straub, 3x constcoalbu]
     )
-    volexpr0 = [30.531e-06] * nruns
-    numconc = [100e6] * nruns
+    if use_marshall_parmer:
+        volexpr0 = ["n/a"] * nruns
+        numconc = ["XXX"] * nruns
+    if not use_marshall_parmer:
+        volexpr0 = [30.531e-06] * nruns
+        numconc = [100e6] * nruns
     maxnsupers = [8192] * nruns
     COLLTSTEP = [1] * nruns
     OBSTSTEP = [10] + [60] * 5 + [10] * 3
@@ -54,11 +64,14 @@ def generate_configurations(path2pySD, path2build, original_config):
 
     config_filenames = []
     for r in range(nruns):
-        cf = tmppath / f"config_bucomp_{r}.yaml"
+        label = "bucomp"
+        if use_marshall_parmer:
+            label += "_marshpam"
+        cf = tmppath / f"config_{label}_{r}.yaml"
         shutil.copy(original_config, cf)
         config_filenames.append(cf)
 
-        initsupers_filename = str(sharepath / f"dimlessSDsinit_bucomp_{r}.dat")
+        initsupers_filename = str(sharepath / f"dimlessSDsinit_{label}_{r}.dat")
 
         params = {
             "savefigpath": str(savefigpath),
@@ -72,10 +85,11 @@ def generate_configurations(path2pySD, path2build, original_config):
             "constants_filename": str(constants_filename),
             "grid_filename": str(grid_filename),
             "initsupers_filename": initsupers_filename,
-            "setup_filename": str(binpath / f"setup_bucomp_{r}.txt"),
-            "zarrbasedir": str(binpath / f"sol_bucomp_{r}.zarr"),
+            "setup_filename": str(binpath / f"setup_{label}_{r}.txt"),
+            "zarrbasedir": str(binpath / f"sol_{label}_{r}.zarr"),
             "coaleff": coaleff[r],
             "nfrags": nfrags[r],
+            "use_marshall_parmer": use_marshall_parmer,
         }
         editconfigfile.edit_config_params(cf, params)
 
@@ -106,9 +120,11 @@ def gridbox_boundaries(path2pySD, config_filename, isfigures=[False, False]):
         savefigpath = Path(pyconfig["paths"]["savefigpath"])
 
         ### --- settings for 0-D Model gridbox boundaries --- ###
-        zgrid = [0, 10, 10]
-        xgrid = [0, 10, 10]
-        ygrid = [0, 10, 10]
+        if pyconfig["supers"]["use_marshall_parmer"]:
+            grid_spacing = [0, 100, 100]  # (!) must be consistent with ``volume'' below
+        else:
+            grid_spacing = [0, 10, 10]
+        zgrid = xgrid = ygrid = grid_spacing
         ### ---------------------------------------------------------------- ###
 
         ### -------------------- INPUT FILES GENERATION -------------------- ###
@@ -139,6 +155,8 @@ def initial_superdroplet_conditions(
     )
     from pySD import geninitconds as gic
 
+    from .probdists_marshallpalmer import get_marshall_palmer_generators
+
     config = yaml.safe_load(open(config_filename))
     pyconfig = config["python_initconds"]
 
@@ -162,14 +180,20 @@ def initial_superdroplet_conditions(
         ### --- settings for initial superdroplets --- ###
         # settings for superdroplet attributes
         dryradius = pyconfig["supers"]["dryradius"]
-        volexpr0 = pyconfig["supers"]["volexpr0"]
-        numconc = pyconfig["supers"]["numconc"]
-        rspan = [5e-6, 7e-5]
+        if pyconfig["supers"]["use_marshall_parmer"]:
+            volume = 100**3  # (!) must be consistent with grid above [m^3]
+            rspan = [5e-7, 4e-3]  # [m]
+            radiigen, xiprobdist = get_marshall_palmer_generators(rspan, volume)
+            numconc = xiprobdist.calc_numconc()
+        else:
+            rspan = [5e-6, 7e-5]
+            volexpr0 = pyconfig["supers"]["volexpr0"]
+            numconc = pyconfig["supers"]["numconc"]
+            xiprobdist = probdists.VolExponential(volexpr0, rspan)
+            radiigen = rgens.SampleLog10RadiiGen(rspan)
 
         # attribute generators
-        radiigen = rgens.SampleLog10RadiiGen(rspan)
         dryradiigen = rgens.MonoAttrGen(dryradius)
-        xiprobdist = probdists.VolExponential(volexpr0, rspan)
         coord3gen = None
         coord1gen = None
         coord2gen = None