Merge pull request #36 from WiltonLoch/coupldyn_yac_3d

Extended the YAC integration for 3d data and added a related example

examples/CMakeLists.txt

@@ -10,7 +10,6 @@ add_subdirectory(boxmodelcollisions/long/src EXCLUDE_FROM_ALL)
 add_subdirectory(boxmodelcollisions/lowlist/src EXCLUDE_FROM_ALL)
 add_subdirectory(constthermo2d/src EXCLUDE_FROM_ALL)
 add_subdirectory(divfreemotion/src EXCLUDE_FROM_ALL)
-add_subdirectory(divfreemotion_yac/src EXCLUDE_FROM_ALL)
 add_subdirectory(rainshaft1d/src EXCLUDE_FROM_ALL)
 add_subdirectory(speedtest/src EXCLUDE_FROM_ALL)
-add_subdirectory(yac_examples/ EXCLUDE_FROM_ALL)
+add_subdirectory(yac/ EXCLUDE_FROM_ALL)

examples/yac_examples/CMakeLists.txt → examples/yac/CMakeLists.txt

@@ -4,4 +4,6 @@ if(NOT DEFINED CMAKE_MINIMUM_REQUIRED_VERSION)
 endif()
 
 # subdirectories for specific yac examples of CLEO
-add_subdirectory(yac1_fromfile/src EXCLUDE_FROM_ALL)
+add_subdirectory(fromfile/src EXCLUDE_FROM_ALL)
+add_subdirectory(yac_3d/src EXCLUDE_FROM_ALL)
+add_subdirectory(divfreemotion/src EXCLUDE_FROM_ALL)

examples/divfreemotion_yac/yac_data_reader.py → examples/yac/divfreemotion/yac_data_reader.py

@@ -98,23 +98,25 @@ def read_metadata(file):
 
 # Read binary data from files
 press_values = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_press.dat")
-temp_values = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_temp.dat")
-qvap_values = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_qvap.dat")
+temp_values  = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_temp.dat")
+qvap_values  = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_qvap.dat")
 qcond_values = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_qcond.dat")
-uvel = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_uvel.dat")
-wvel = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_wvel.dat")
+uvel         = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_uvel.dat")
+wvel         = thermodynamicvar_from_binary("../build/share/df2d_dimlessthermo_wvel.dat")
 
 # Pack all wind velocities edge data into one united field for YAC exchange
 united_edge_data = []
-for lat_index in range(0, len(lat) * 2 - 1):
-    if (lat_index % 2 == 0):
-        lower_index = (lat_index // 2) * (len(lon) - 1)
-        upper_index = (lat_index // 2 + 1) * (len(lon) - 1)
-        united_edge_data.extend(uvel[lower_index:upper_index])
-    else:
-        lower_index = ((lat_index - 1) // 2) * len(lon)
-        upper_index = ((lat_index + 1) // 2) * len(lon)
-        united_edge_data.extend(wvel[lower_index:upper_index])
+for timestep in range(0, 3):
+    timestep_offset = timestep * (len(edge_centers_lon) // 2 )
+    for lat_index in range(0, len(lat) * 2 - 1):
+        if (lat_index % 2 == 0):
+            lower_index = timestep_offset + (lat_index // 2) * (len(lon) - 1)
+            upper_index = timestep_offset + (lat_index // 2 + 1) * (len(lon) - 1)
+            united_edge_data.extend(uvel[lower_index:upper_index])
+        else:
+            lower_index = timestep_offset + ((lat_index - 1) // 2) * len(lon)
+            upper_index = timestep_offset + ((lat_index + 1) // 2) * len(lon)
+            united_edge_data.extend(wvel[lower_index:upper_index])
 
 # Use YAC to exchange the values for first timestep
 press.put(press_values[0:900])

examples/yac/yac_3d/src/CMakeLists.txt

@@ -0,0 +1,29 @@
+# set cmake version
+if(NOT DEFINED CMAKE_MINIMUM_REQUIRED_VERSION)
+  cmake_minimum_required(VERSION 3.18.0)
+  # cmake_minimum_required(VERSION 3.21.1) # if using Kokkos c++ with NVC++ compiler
+endif()
+
+# set project name and print directory of this CMakeLists.txt (source directory of project)
+project("yac_3d")
+message("PROJECT_SOURCE_DIR: ${PROJECT_SOURCE_DIR}")
+
+# Set libraries from CLEO to link with executable
+set(CLEOLIBS gridboxes initialise observers runcleo superdrops zarr)
+
+# create primary executable for CLEO in 2-D coupled to thermodynamics from file setup
+add_executable(yac_3d EXCLUDE_FROM_ALL "main_yac1.cpp")
+
+# Add directories and link libraries to target
+target_link_libraries(yac_3d PRIVATE coupldyn_yac cartesiandomain "${CLEOLIBS}")
+target_link_libraries(yac_3d PUBLIC Kokkos::kokkos)
+target_include_directories(yac_3d PRIVATE "${CMAKE_SOURCE_DIR}/libs")
+
+# set specific C++ compiler options for target (optional)
+#target_compile_options(yac1 PRIVATE)
+
+# set compiler properties for target(s)
+set_target_properties(yac_3d PROPERTIES
+  CMAKE_CXX_STANDARD_REQUIRED ON
+  CMAKE_CXX_EXTENSIONS ON
+  CXX_STANDARD 20)

examples/yac/yac_3d/src/config/yac1_fromfile_config.txt

@@ -0,0 +1,67 @@
+/*
+ * ----- CLEO -----
+ * File: yac1_fromfile_config.txt
+ * Project: config
+ * Created Date: Tuesday 21st November 2023
+ * Author: Clara Bayley (CB)
+ * Additional Contributors:
+ * -----
+ * Last Modified: Tuesday 26th March 2024
+ * Modified By: CB
+ * -----
+ * License: BSD 3-Clause "New" or "Revised" License
+ * https://opensource.org/licenses/BSD-3-Clause
+ * -----
+ * Copyright (c) 2023 MPI-M, Clara Bayley
+ * -----
+ * File Description:
+ * configuration input parameters for CLEO yac test 1, with 3D time vayring
+ * thermodyanamics read from a file.
+ */
+
+### Initialisation parameters ###
+constants_filename = ../libs/cleoconstants.hpp              # name of file for values of physical constants
+initsupers_filename = ./share/yac1_dimlessSDsinit.dat       # binary filename for initialisation of SDs
+grid_filename = ./share/yac1_dimlessGBxboundaries.dat       # binary filename for initialisation of GBxs / GbxMaps
+
+### Output Data parameters ###
+setuptxt = ./bin/yac1_setup.txt                # .txt filename to copy configuration to
+stats_filename = ./bin/yac1_stats.txt          # .txt file to output runtime statistics to
+zarrbasedir = ./bin/yac1_sol.zarr              # zarr store base directory
+maxchunk = 1250000                             # maximum no. of elements in chunks of zarr store array
+
+### SDM Runtime parameters ###
+# domain setup #
+nspacedims = 3                           # no. of spatial dimensions to model
+ngbxs = 2250                             # total number of Gbxs
+totnsupers = 2880                        # (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 = 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 #
+cond_iters = 2                            # no. iterations of Newton Raphson Method before testing for convergence
+cond_SUBTSTEP = 0.1                       # smallest timestep in cases where substepping occurs [s]
+cond_rtol = 0.0                           # relative tolerance for implicit euler integration
+cond_atol = 0.01                          # abolute tolerance for implicit euler integration
+
+# superdroplets #
+doAlterThermo = false                   # enable condensation to alter the thermodynamic state
+
+### Coupled Dynamics Solver Parameters ###
+# type of coupling #
+thermosolver = fromfile                                   # dynamics solver to configure
+
+### read in dynamics from file ###
+press_filename = ./share/yac1_dimlessthermo_press.dat          # binary filename for pressure
+temp_filename = ./share/yac1_dimlessthermo_temp.dat            # binary filename for temperature
+qvap_filename = ./share/yac1_dimlessthermo_qvap.dat            # binary filename for vapour mixing ratio
+qcond_filename = ./share/yac1_dimlessthermo_qcond.dat          # binary filename for liquid mixing ratio
+wvel_filename = ./share/yac1_dimlessthermo_wvel.dat            # binary filename for vertical (z) velocity
+uvel_filename = ./share/yac1_dimlessthermo_uvel.dat            # binary filename for horizontal x velocity
+vvel_filename = ./share/yac1_dimlessthermo_vvel.dat            # binary filename for horizontal y velocity

examples/yac/yac_3d/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