test_shipway_hill.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PySDM_examples.Shipway_and_Hill_2012 import Simulation, Settings\n",
    "from PySDM.physics import si\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "settings = Settings(n_sd_per_gridbox=2048, precip=False, p0=1007 * si.hPa, particles_per_volume_STP=1e10, z_max=2000, dz=2000/256, t_max=1000, dt=0.1, rho_times_w_1=1 * si.m / si.s * si.kg / si.m**3, kappa=0.1)\n",
    "simulation = Simulation(settings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = simulation.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "products = results.products"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "products.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ylim((100,2000))\n",
    "plt.pcolormesh(products[\"t\"], products[\"z\"], products[\"nc\"])\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(products[\"z\"], products[\"nc\"][:,-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PySDM_ds = xr.Dataset(\n",
    "    data_vars={\n",
    "        \"N_aer\" : ([\"t\", \"zc\"], products[\"na\"].T),\n",
    "        \"N_liq\" : ([\"t\", \"zc\"], products[\"nc\"].T),\n",
    "        \"q_v\" : ([\"t\", \"zc\"], products[\"qv\"].T / (products[\"qv\"].T + 1)),\n",
    "        \"q_liq\" : ([\"t\", \"zc\"], (products[\"qc\"].T / 1000) / (products[\"qc\"].T / 1000 + 1)),\n",
    "    },\n",
    "    coords={\n",
    "        \"t\" : products[\"t\"],\n",
    "        \"zc\" : products[\"z\"],\n",
    "    }\n",
    ")\n",
    "PySDM_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PySDM_ds.to_netcdf(\"kinematic_simulations/PySDM_output_2.nc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ylim(0, 2000)\n",
    "plt.xlabel(\"t (s)\")\n",
    "plt.ylabel(\"z (m)\")\n",
    "plt.title(\"PySDM N_liq\")\n",
    "plt.pcolormesh(PySDM_ds[\"t\"], PySDM_ds[\"zc\"], (PySDM_ds[\"N_liq\"]).T)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "KiD_ARG_ds = xr.open_dataset(\"kinematic_simulations/KiD_output_ARG_2.nc\", group=\"profiles\")\n",
    "KiD_ARG_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ARG_diff_ds = (KiD_ARG_ds.interp(PySDM_ds.coords) - PySDM_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "KiD_emulated_ds = xr.open_dataset(\"kinematic_simulations/KiD_output_emulated_2.nc\", group=\"profiles\")\n",
    "KiD_emulated_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emulated_diff_ds = (KiD_emulated_ds.interp(PySDM_ds.coords) - PySDM_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ylim(0, 2000)\n",
    "plt.xlabel(\"t (s)\")\n",
    "plt.ylabel(\"z (m)\")\n",
    "plt.title(\"ARG Scheme N_liq\")\n",
    "plt.pcolormesh(KiD_ARG_ds[\"t\"], KiD_ARG_ds[\"zc\"], KiD_ARG_ds[\"N_liq\"].T)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ylim(0, 2000)\n",
    "plt.xlabel(\"t (s)\")\n",
    "plt.ylabel(\"z (m)\")\n",
    "plt.title(\"ARG Scheme N_liq Deviation from PySDM\")\n",
    "plt.pcolormesh(ARG_diff_ds[\"t\"], ARG_diff_ds[\"zc\"], ARG_diff_ds[\"N_liq\"].T, cmap=\"coolwarm\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ylim(0, 2000)\n",
    "plt.xlabel(\"t (s)\")\n",
    "plt.ylabel(\"z (m)\")\n",
    "plt.title(\"Emulated N_liq\")\n",
    "plt.pcolormesh(KiD_emulated_ds[\"t\"], KiD_emulated_ds[\"zc\"], KiD_emulated_ds[\"N_liq\"].T)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ylim(0, 2000)\n",
    "plt.xlabel(\"t (s)\")\n",
    "plt.ylabel(\"z (m)\")\n",
    "plt.title(\"Emulated N_liq Deviation from PySDM\")\n",
    "plt.pcolormesh(emulated_diff_ds[\"t\"], emulated_diff_ds[\"zc\"], emulated_diff_ds[\"N_liq\"].T, cmap=\"coolwarm\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_diff_ds = np.abs(emulated_diff_ds) - np.abs(ARG_diff_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.ylim(0, 2000)\n",
    "plt.xlabel(\"t (s)\")\n",
    "plt.ylabel(\"z (m)\")\n",
    "plt.title(\"Difference of Absolute Deviations (Emulated - ARG)\")\n",
    "plt.pcolormesh(diff_diff_ds[\"t\"], diff_diff_ds[\"zc\"], diff_diff_ds[\"N_liq\"].T, cmap=\"RdYlGn_r\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"N_liq (m^-3)\")\n",
    "plt.ylabel(\"z (m)\")\n",
    "plt.ylim(0, 2000)\n",
    "plt.plot(PySDM_ds[\"q_liq\"][-1,:], PySDM_ds[\"zc\"], label=\"PySDM\")\n",
    "plt.plot(KiD_ARG_ds[\"q_liq\"][-1,:], KiD_ARG_ds[\"zc\"], label=\"ARG\")\n",
    "plt.plot(KiD_emulated_ds[\"q_liq\"][-1,:], KiD_emulated_ds[\"zc\"], ls=\"--\", label=\"Emulated\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel(\"N_aer (m^-3)\")\n",
    "plt.ylabel(\"z (m)\")\n",
    "plt.ylim(0, 2000)\n",
    "plt.plot(PySDM_ds[\"N_aer\"][-1,:], PySDM_ds[\"zc\"], label=\"PySDM\")\n",
    "plt.plot(KiD_ARG_ds[\"N_aer\"][-1,:], KiD_ARG_ds[\"zc\"], label=\"ARG\")\n",
    "plt.plot(KiD_emulated_ds[\"N_aer\"][-1,:], KiD_emulated_ds[\"zc\"], ls=\"--\", label=\"Emulated\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "CliMA",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}