Merge pull request #30 from yoctoyotta1024/newexample

New Example for Time Varying 3-D Thermodynamics

examples/exampleplotting/plotssrc/pltsds.py

@@ -6,7 +6,7 @@
 Author: Clara Bayley (CB)
 Additional Contributors:
 -----
-Last Modified: Tuesday 21st November 2023
+Last Modified: Monday 25th March 2024
 Modified By: CB
 -----
 License: BSD 3-Clause "New" or "Revised" License
@@ -126,14 +126,17 @@ def plot_randomsample_superdrops(time, sddata, totnsupers,
 
 def plot_randomsample_superdrops_2dmotion(sddata, totnsupers,
                                           nsample, savename="",
-                                          arrows=False):
+                                          arrows=False, israndom=True):
   ''' plot timeseries of the attributes of a
   random sample of superdroplets '''
 
   fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6,6))
 
   minid, maxid = 0, int(totnsupers) # largest value of ids to sample
-  ids2plot = random.sample(list(range(minid, maxid, 1)), nsample)
+  if israndom:
+    ids2plot = random.sample(list(range(minid, maxid, 1)), nsample)
+  else:
+    ids2plot = np.linspace(0, maxid-1, nsample, dtype=int)
 
   mks = MarkerStyle('o', fillstyle='full')
   coordz = sdtracing.attribute_for_superdroplets_sample(sddata, "coord3",

examples/yac_examples/yac1_fromfile/src/config/yac1_fromfile_config.txt

@@ -6,7 +6,7 @@
  * Author: Clara Bayley (CB)
  * Additional Contributors:
  * -----
- * Last Modified: Friday 22nd March 2024
+ * Last Modified: Monday 25th March 2024
  * Modified By: CB
  * -----
  * License: BSD 3-Clause "New" or "Revised" License
@@ -33,15 +33,15 @@ maxchunk = 1250000                             # maximum no. of elements in chun
 ### SDM Runtime parameters ###
 # domain setup #
 nspacedims = 3                           # no. of spatial dimensions to model
-ngbxs = 1250                             # total number of Gbxs
-totnsupers = 1600                        # (initial) total no. of SDs
+ngbxs = 1875                             # total number of Gbxs
+totnsupers = 2400                        # (initial) total no. of SDs
 
 # timestepping #
 CONDTSTEP = 2                             # time between SD condensation events [s]
 COLLTSTEP = 2                             # time between SD collision events [s]
 MOTIONTSTEP = 3                           # time between SDM motion [s]
-COUPLTSTEP = 3600                         # time between dynamic couplings [s]
-OBSTSTEP = 180                            # time between SDM observations [s]
+COUPLTSTEP = 1800                         # time between dynamic couplings [s]
+OBSTSTEP = 1800                           # time between SDM observations [s]
 T_END = 7200                              # time span of integration from 0s to T_END [s]
 
 # microphysics #

examples/yac_examples/yac1_fromfile/src/gen_input_thermo.py

@@ -0,0 +1,122 @@
+'''
+Copyright (c) 2024 MPI-M, Clara Bayley
+
+
+----- CLEO -----
+File: gen_input_thermo.py
+Project: src
+Created Date: Monday 25th March 2024
+Author: Clara Bayley (CB)
+Additional Contributors:
+-----
+Last Modified: Monday 25th March 2024
+Modified By: CB
+-----
+License: BSD 3-Clause "New" or "Revised" License
+https://opensource.org/licenses/BSD-3-Clause
+-----
+File Description:
+Python functions used by yac1_fromfile example to make thermo and wind fields
+for CLEO to run example with 3-D time-varying thermodynamics.
+'''
+
+
+import sys
+import numpy as np
+
+sys.path.append("../../../..") # for imports from pySD package
+from pySD import cxx2py
+from pySD.thermobinary_src.create_thermodynamics import thermoinputsdict
+from pySD.gbxboundariesbinary_src import read_gbxboundaries as rgrid
+
+class TimeVarying3DThermo:
+  ''' create some sinusoidal thermodynamics that varies in time and is
+  hetergenous throughout 3D domain '''
+
+  def __init__(self, PRESSz0, TEMPz0, qvapz0, qcondz0,
+               WMAX, Zlength, Xlength, VMAX, Ylength):
+
+    ### parameters of profile ###
+    self.PRESSz0 = PRESSz0 # pressure at z=0m [Pa]
+    self.TEMPz0 = TEMPz0 # temperature at z=0m [K]
+    self.qvapz0 = qvapz0 # vapour mass mixing ratio at z=0m [Kg/Kg]
+    self.qcondz0 = qcondz0 # liquid mass mixing ratio at z=0m [Kg/Kg]
+    self.dimless_omega = np.pi/4.0 # ~ frequency of time modulation []
+
+    self.WMAX = WMAX  # max velocities constant [m/s]
+    self.Zlength = Zlength # wavelength of velocity modulation in z direction [m]
+    self.Xlength = Xlength # wavelength of velocity modulation in x direction [m]
+    self.VMAX = VMAX  # max horizontal (y) velocity
+    self.Ylength = Ylength # wavelength of velocity modulation in y direction [m]
+
+  def idealised_flowfield2D(self, gbxbounds, ndims):
+
+    zfaces, xcens_z, ycens_z = rgrid.coords_forgridboxfaces(gbxbounds, ndims, 'z')
+    zcens_x, xfaces, ycens_x = rgrid.coords_forgridboxfaces(gbxbounds, ndims, 'x')
+
+    ztilda = self.Zlength / np.pi
+    xtilda = self.Xlength / (2*np.pi)
+    wamp = 2 * self.WMAX
+
+    WVEL = wamp * np.sin(zfaces / ztilda) * np.sin(xcens_z / xtilda)
+    UVEL = wamp * xtilda / ztilda * np.cos(zcens_x/ztilda) * np.cos(xfaces/xtilda)
+
+    # modulation in y direction
+    WVEL *= (1.0 + 0.5 * np.cos(self.Ylength / np.pi * ycens_z))
+    UVEL *= (1.0 + 0.5 * np.cos(self.Ylength / np.pi * ycens_x))
+
+    return WVEL, UVEL
+
+  def gen_3dvvelocity(self, gbxbounds, ndims):
+
+    zcens_y, xcens_y, yfaces = rgrid.coords_forgridboxfaces(gbxbounds, ndims, 'y')
+    zxmod = zcens_y / self.Zlength + xcens_y / self.Xlength
+    VVEL = self.VMAX * (zxmod + np.cos(self.Ylength / np.pi * yfaces))
+
+    return VVEL
+
+  def generate_timevarying_3dwinds(self, gbxbounds, ndims, ntime, THERMODATA):
+
+    # time modulation factor for variables at each timestep
+    tmod = np.full(ntime, -0.5)
+    tmod = np.power(tmod, np.array(range(0, ntime, 1)))
+
+    WVEL, UVEL = self.idealised_flowfield2D(gbxbounds, ndims)
+    VVEL = self.gen_3dvvelocity(gbxbounds, ndims)
+
+    THERMODATA["WVEL"] = np.outer(tmod, WVEL).flatten()
+    THERMODATA["UVEL"] = np.outer(tmod, UVEL).flatten()
+    THERMODATA["VVEL"] = np.outer(tmod, VVEL).flatten()
+
+    return THERMODATA
+
+  def generate_3dsinusoidal_variable(self, gbxbounds, ndims, amp):
+
+    zfulls, xfulls, yfulls = rgrid.fullcoords_forallgridboxes(gbxbounds, ndims)
+
+    ztilda = self.Zlength / np.pi / 3.0
+    xtilda = self.Xlength / np.pi / 4.0
+    ytilda = self.Ylength / np.pi / 1.33
+
+    return amp + 0.25 * amp * (np.sin(zfulls / ztilda) * np.sin(xfulls / xtilda) + np.sin(yfulls/ ytilda))
+
+  def generate_thermo(self, gbxbounds, ndims, ntime):
+
+    PRESS = self.generate_3dsinusoidal_variable(gbxbounds, ndims, self.PRESSz0)
+    TEMP = self.generate_3dsinusoidal_variable(gbxbounds, ndims, self.TEMPz0)
+    qvap = self.generate_3dsinusoidal_variable(gbxbounds, ndims, self.qvapz0)
+    qcond = self.generate_3dsinusoidal_variable(gbxbounds, ndims, self.qcondz0)
+
+    tmod = np.cos(self.dimless_omega * np.arange(0.0, ntime, 1.0))
+    tmod = 1 + 0.5 * tmod
+
+    THERMODATA = {
+      "PRESS": np.outer(tmod, PRESS).flatten(),
+      "TEMP": np.outer(tmod, TEMP).flatten(),
+      "qvap": np.outer(tmod, qvap).flatten(),
+      "qcond": np.outer(tmod, qcond).flatten(),
+    }
+
+    THERMODATA = self.generate_timevarying_3dwinds(gbxbounds, ndims, ntime, THERMODATA)
+
+    return THERMODATA

examples/yac_examples/yac1_fromfile/src/plot_output_thermo.py

@@ -0,0 +1,115 @@
+'''
+Copyright (c) 2024 MPI-M, Clara Bayley
+
+
+----- CLEO -----
+File: plot_output_thermo.py
+Project: src
+Created Date: Monday 25th March 2024
+Author: Clara Bayley (CB)
+Additional Contributors:
+-----
+Last Modified: Monday 25th March 2024
+Modified By: CB
+-----
+License: BSD 3-Clause "New" or "Revised" License
+https://opensource.org/licenses/BSD-3-Clause
+-----
+File Description:
+Python functions used to make plots of CLEO's thermodynamics output for
+yac1_fromfile example.
+'''
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.gridspec import GridSpec
+from matplotlib import colors
+
+def plot_domain_thermodynamics_timeseries(time, gbxs, thermo, winds, savedir):
+  ''' plot 2-D cross-sections of domain along y axis for thermodynamics and
+  wind fields at a series of time slices. '''
+
+  labels = {
+    "press": "Pressure /hPa",
+    "temp": "Temperature /K",
+    "qvap": "q$_{v}$ Kg/Kg",
+    "qcond": "q$_{c}$ Kg/Kg",
+    "wvel": "w m/s",
+    "uvel": "u m/s",
+    "vvel": "v m/s",
+  }
+
+  cmaps = {
+    "press": ["PRGn", colors.CenteredNorm(vcenter=1015)],
+    "temp": ["RdBu_r", colors.CenteredNorm(vcenter=300)],
+    "qvap": ["BrBG", colors.CenteredNorm(vcenter=0.01)],
+    "qcond": ["BrBG", colors.CenteredNorm(vcenter=0.001)],
+    "wvel": ["coolwarm", colors.CenteredNorm(vcenter=0.0, halfrange=3.0)],
+    "uvel": ["coolwarm", colors.CenteredNorm(vcenter=0.0, halfrange=3.0)],
+    "vvel": ["coolwarm", colors.CenteredNorm(vcenter=0.0, halfrange=3.0)],
+  }
+
+  xxh, zzh = np.meshgrid(gbxs["xhalf"], gbxs["zhalf"], indexing="ij") # dims [xdims, zdims]
+  t2plts = np.linspace(time.mins[0], time.mins[-1], 5)
+
+  for key in labels.keys():
+    try:
+      data4d = thermo[key]
+    except:
+      data4d = winds[key]
+
+    plot_timeseries_domain_slices(key, labels[key], cmaps[key][0], cmaps[key][1],
+                                  time, t2plts, xxh, zzh, gbxs["yfull"],
+                                  data4d, savedir)
+
+def plot_timeseries_domain_slices(key, label, cmap, norm, time, t2plts,
+                                  xxh, zzh, yfull, data4d, savedir):
+  ''' plot 2-D cross-sections along y axis of 3-D data 'data3d' for several timeslices
+  given times and meshgrid for centres in x and z, 'xxh' and 'zzh'. '''
+
+  fig = plt.figure(figsize=(16, 9))
+  ncols = len(t2plts)
+  nrows = data4d.shape[1]
+  gs = GridSpec(nrows=nrows+1, ncols=ncols, figure=fig, height_ratios=[1]+[8]*nrows)
+  cax = fig.add_subplot(gs[0,:])
+
+  cax.set_title(label)
+  pcm = plot_domain_slices(fig, gs, nrows, ncols, time, t2plts,
+                           xxh, zzh, yfull, data4d, cmap, norm)
+  cbar = plt.colorbar(pcm, cax=cax, orientation='horizontal')
+
+  fig.tight_layout()
+
+  if savedir != "":
+    savename = savedir+"/timeseries_"+key+".png"
+    fig.savefig(savename, dpi=400,
+            bbox_inches="tight", facecolor='w', format="png")
+    print("Figure .png saved as: "+savename)
+    plt.show()
+
+def plot_domain_slices(fig, gs, nrows, ncols, time, t2plts,
+                       xxh, zzh, yfull, data4d, cmap, norm):
+  ''' plot 2-D cross-sections along y axis of 4-D data 'data4d' at a series of time slices
+  given meshgrid for centres in x and z, 'xxh' and 'zzh'. '''
+
+  for m in range(ncols):
+    t = np.argmin(abs(time.mins-t2plts[m]))
+    data3d = data4d[t, :, :, :]
+    for n in range(nrows):
+      ax = fig.add_subplot(gs[n+1, m])
+      ax.set_title("t={:.0f}min,\ny={:.2f}km".format(time.mins[t], yfull[n]/1000))
+      pcm = ax_colormap(ax, xxh, zzh, data3d[n, :, :], cmap, norm)
+
+  return pcm
+
+def ax_colormap(ax, xxh, zzh, data2d, cmap, norm):
+  ''' plot pcolormesh for 2-D data 'data2d' given meshgrid for
+  centres in x and z, 'xxh' and 'zzh'. '''
+
+  ax.set_aspect("equal")
+  pcm = ax.pcolormesh(xxh[:,:]/1000, zzh[:,:]/1000, data2d, cmap=cmap, norm=norm)
+  ax.set_xlabel("x /km")
+  ax.set_ylabel("z /km")
+
+  return pcm