updating the README.md and addiding comments to the code, addressing …

…issue NOAA-EMC#55

unst_msh_gen/README.md

@@ -1,32 +1,32 @@
 # Project Name
 Unstructured mesh generation for WW3 using JIGSSAW.
-## Description
+# Description
 Mesh generation script capable of creating unstructured meshes for WW3 global modeling. The tool leverages JIGSAWPY (https://github.com/dengwirda/jigsaw-python) for efficient triangulation.
 Main changes include:
 Implementation of the ocn_ww3.py to create uniform unstructured mesh for global WW3 model.
 This tool is under active development, with future work focused on variable unstructured mesh generation.
 
-## Installation
+# Installation
 
-1- Install jigsawpy (https://github.com/dengwirda/jigsaw-python)
-2- you need following packages:
-	numpy
-        scipy
-	packaging
-	netcdf4
+## 1- Install jigsawpy (https://github.com/dengwirda/jigsaw-python)
 
-3- clone the repo
-        $git clone https://github.com/NOAA-EMC/WW3-tools
-	$cd WW3-tools/unst_msh_gen
+## 2- you need following packages:
+- numpy
+- scipy
+- packaging
+- netcdf4
 
-3- get the DEM and make sure it is in the WW3-tools/unst_msh_gen directory
-	$wget https://github.com/dengwirda/dem/releases/download/v0.1.1/RTopo_2_0_4_GEBCO_v2023_60sec_pixel.zip
-4- unzip the DEMi 
-	unzip *.zip
+## 3- clone the repo
+- $git clone https://github.com/NOAA-EMC/WW3-tools
+- $cd WW3-tools/unst_msh_gen
+
+## 4- get the DEM and make sure it is in the WW3-tools/unst_msh_gen directory
+- $wget https://github.com/dengwirda/dem/releases/download/v0.1.1/RTopo_2_0_4_GEBCO_v2023_60sec_pixel.zip
+- $unzip *.zip
 
-## Usage
-5- run the script inside of WW3-tools/unst_msh_gen:
-	$python3 ocn_ww3.py --black_sea [option]
+# Usage
+## 5- run the script inside of WW3-tools/unst_msh_gen:
+- $python3 ocn_ww3.py --black_sea [option]
 
 		option=  3: default which will have the Black Sea and the connections.
 			 2: will have the Balck sea as a seperate basin.
@@ -35,37 +35,43 @@ This tool is under active development, with future work focused on variable unst
 NOTE: the output will be gmsh format which will be used by WW3.
 
 NOTE: for different resolution (uniform) in km, you should change the following:
-	opts.hfun_hmax
-	hmax
-	hshr
-	hmin
+- opts.hfun_hmax
+- hmax
+- hshr
+- hmin
 
 NOTE: The output mesh will have -180:180 longitude, you can convert this by unisg ShiftMesh.py script, to 0:360 longitude.
 	input_file_path: your jigsaw mesh in gmsh format with -18:180 long
 	output_file_path: shifted mesh in gmsh format with 0:360 long
 
 
-6- Using variable mode:
-	$python3 window_mask.py --windows '[{"min_lon": -98, "max_lon": -64, "min_lat": 24, "max_lat": 44.5, "hshr": 5}'
+## 6- Using variable mode:
+- To create a mesh with finer resolution near the US coastlines you can define different region in json format (east coast, west coast and golf od Mexico, Purto Rico, and Hawaii):
+
+- $python3 window_mask.py --windows '[{"min_lon": -98, "max_lon": -64, "min_lat": 24, "max_lat": 44.5, "hshr": 5}, {"min_lon": -158, "max_lon": -155, "min_lat": 19, "max_lat": 22, "hshr": 5}, {"min_lon": -128, "max_lon": -64, "min_lat": 34.5, "max_lat": 48.5, "hshr": 5}, {"min_lon": -67.4, "max_lon": -64.1, "min_lat": 17.2, "max_lat": 18.3, "hshr": 5}]'
+
+NOTE: hshr is the shoreline resolution which can be smaller than the hmin which is defined globally.
 
-NOTE: You can define multiple windows with different shoreline resolution (hshr) in json format.
 NOTE: You can define different background mesh based on lat location in the window_mask.py:
 
-	 # Define the boundaries and scaling values
-   	 upper_bound = 50
-    	 middle_bound = -20
-    	 lower_bound = -90
-    	 scale_north = 9
-    	 scale_middle = 20
-    	 scale_south_upper = 30
-    	 scale_south_lower = 9
+NOTE: for the background mesh you can define different resolution (for eaxmple, the globe is divided to three regions based on latitude and corresponding resolutions: 
+
+- upper_bound = 50          # for lat > 50N
+- middle_bound = -20        # for -20S < lat < 50N
+- lower_bound = -90         # for lat < -20S
+
+- scale_north = 9           # mesh resolution for lat < 50N
+- scale_middle = 20         # mesh resolution for -20S< lat <50N
+- scale_south_upper = 30    # mesh resolution linearly decreasing from 30km at -20S to 
+- scale_south_lower = 9     # 9km mesh resolution at -90S
 
 NOTE: The output will be wmask.nc file which will have the mesh spacing info.
 
 
-        $python3 ocn_ww3.py --black_sea [option] --mask_file="wmask.nc"
+NOTE:To create the variable mesh based on specified mesh spacing file in "wmask.nc" you can use the following comand:
+
+- $python3 ocn_ww3.py --black_sea [option] --mask_file="wmask.nc"
 
-NOTE: This will create the variable mesh based on specified mesh spacing file.
 
 ## Contributing
 This is ongoing effort with the great help of Darren Engwirda, JIGSAW developer.

unst_msh_gen/ocn_ww3.py

@@ -51,9 +51,9 @@ def create_msh():
 
     print("*create-msh...")
 
-    opts.geom_file = "geom.msh"
-    opts.hfun_file = "spac.msh"
-    opts.jcfg_file = "opts.jig"
+    opts.geom_file = "geom.msh"  #saves the geometry info for jigsaw
+    opts.hfun_file = "spac.msh"  #saves the final mesh spacing info 
+    opts.jcfg_file = "opts.jig"  #jigsaw ctlr file
 
     geom.mshID = "ellipsoid-mesh"
     geom.radii = np.full(
@@ -68,26 +68,26 @@ def create_msh():
     # solve |dh/dx| constraints in spacing
     jigsawpy.cmd.marche(opts, spac)
 
-    opts.mesh_file = "uglo_100km_BlkS.msh"
+    opts.mesh_file = "uglo_100km_BlkS.msh"  #jigsaw format mesh file
 
     opts.hfun_scal = "absolute"
-    opts.hfun_hmax = 100.           # uniform at 30.km
+    opts.hfun_hmax = 100.           # global maximum mesh resolution (similar to hmax)
     opts.mesh_dims = +2             # 2-dim. simplexes
-    opts.optm_iter = +64
+    opts.optm_iter = +64            # number of itereation for the optimization
     opts.optm_cost = "skew-cos"
 
     jigsawpy.cmd.jigsaw(opts, mesh)
 
 
 def create_siz():
 
-#-- create mesh spacing function for the globe
+    #-- create mesh spacing function for the globe: for uniform mesh hmax = hshr = hmin
 
-    hmax = 100.0 # maximum spacing [km]
+    hmax = 100.0 # maximum spacing [km] 
     hshr = 100  # shoreline spacing
     nwav = 400.  # number of cells per sqrt(g*H)
     hmin = 100.0  # minimum spacing
-    dhdx = 0.05  # allowable spacing gradient
+    dhdx = 0.05  # allowable spacing gradient: for more gradual transition use lower value
 
     data = nc.Dataset(
         "RTopo_2_0_4_GEBCO_v2023_60sec_pixel.nc", "r")
@@ -527,7 +527,7 @@ def write_gmsh_mesh(filename, node_data, tri):
     point = jigsawpy.R3toS2(geom.radii, point)  # to [lon,lat] in deg
     point*= 180. / np.pi
     depth = np.reshape(-1*mesh.value, (mesh.value.size, 1))
-    depth[depth <= 0] = 2
+    depth[depth <= 0] = 50
     point = np.hstack((point, depth))  # append elev. as 3rd coord.
     cells = [("triangle", mesh.tria3["index"])]
     tri_data=cells[0][1]+1

unst_msh_gen/window_mask.py

@@ -29,14 +29,14 @@ def create_mask_file(windows):
     # Initialize scal with a default condition
     #scal = np.where(ymat > 50, 9, np.where(ymat < -20, 30, 20))
 
-    # Define the boundaries and scaling values
-    upper_bound = 50
-    middle_bound = -20
-    lower_bound = -90
-    scale_north = 9
-    scale_middle = 20
-    scale_south_upper = 30
-    scale_south_lower = 9
+    # Define the boundaries and scaling values: The globe is divided to three regions based on latitude and corresponding resolutions
+    upper_bound = 50       # for lat > 50N
+    middle_bound = -20     # for -20S < lat < 50N
+    lower_bound = -90      # for lat < -20S
+    scale_north = 9        # mesh resolution for lat < 50N 
+    scale_middle = 20      # mesh resolution for -20S< lat <50N
+    scale_south_upper = 30 # mesh resolution linearly decreasing from 30km at -20S to
+    scale_south_lower = 9  # 9km mesh resolution at -90S
 
     # Calculate the scaling using conditions
     scal = np.where(ymat > upper_bound,