Merge pull request #4 from akturner/seaice/two_dimensional_hex_square…

…_testcase

Improvements to square test case

testing_and_setup/seaice/testcases/square/create_grids.py

@@ -0,0 +1,212 @@
+from netCDF4 import Dataset
+import math
+import os
+import matplotlib.pyplot as plt
+import numpy as np
+
+mpas_tools_dir = "/home/akt/Work/MPAS-Seaice/MPAS-Tools/feature_branches/master/MPAS-Tools/"
+
+#-------------------------------------------------------------------------------
+
+def create_grid_hex(nx, ny, dc):
+
+    Lx = float(nx) * dc
+    Ly = float(ny) * dc
+
+    gridname = "grid_hex_%4.4ix%4.4i.nc" %(nx,ny)
+
+    fileout = open("namelist.input","w")
+
+    line = "&periodic_grid\n"
+    fileout.write(line)
+
+    line = "   nx = %i,\n" %(nx+2)
+    fileout.write(line)
+
+    line = "   ny = %i,\n" %(ny+2)
+    fileout.write(line)
+
+    line = "   dc = %f,\n" %(dc)
+    fileout.write(line)
+
+    line = "   nVertLevels = 1,\n"
+    fileout.write(line)
+
+    line = "   nTracers = 1,\n"
+    fileout.write(line)
+
+    line = "   nproc = 1,\n"
+    fileout.write(line)
+
+    line = "/\n"
+    fileout.write(line)
+
+    fileout.close()
+
+    os.system("%s/mesh_tools/periodic_hex/periodic_grid" %(mpas_tools_dir))
+    os.system("python %s/mesh_tools/periodic_hex/mark_periodic_boundaries_for_culling.py -f grid.nc" %(mpas_tools_dir))
+    os.system("%s/mesh_tools/mesh_conversion_tools/MpasCellCuller.x grid.nc grid_culled.nc" %(mpas_tools_dir))
+
+    filein = Dataset("grid_culled.nc","a")
+
+    filein.Lx = Lx
+    filein.Ly = Ly
+
+    filein.nx = nx
+    filein.ny = ny
+
+    filein.dc = dc
+
+    filein.close()
+
+    cmd = "mv grid_culled.nc %s" %(gridname)
+    os.system(cmd)
+
+#-------------------------------------------------------------------------------
+
+def create_grid_quad(nx, ny, dc):
+
+    Lx = float(nx) * dc
+    Ly = float(ny) * dc
+
+    gridname = "grid_quad_%4.4ix%4.4i.nc" %(nx,ny)
+
+    vertexDegree = 4
+
+    fileGrid = Dataset("grid_in.nc","w",format="NETCDF3_CLASSIC")
+
+    fileGrid.on_a_sphere = "NO"
+
+    nCells = nx * ny
+    nVertices = (nx+1) * (ny+1)
+
+    fileGrid.createDimension("nCells", nCells)
+    fileGrid.createDimension("nVertices", nVertices)
+    fileGrid.createDimension("vertexDegree", vertexDegree)
+
+    xCell = np.zeros(nCells)
+    yCell = np.zeros(nCells)
+    zCell = np.zeros(nCells)
+
+    for j in range(0,ny):
+        for i in range(0,nx):
+            iCell = i + j * nx
+            xCell[iCell] = (float(i) + 0.5) * dc
+            yCell[iCell] = (float(j) + 0.5) * dc
+
+    xVertex = np.zeros(nVertices)
+    yVertex = np.zeros(nVertices)
+    zVertex = np.zeros(nVertices)
+    cellsOnVertex = np.zeros((nVertices,vertexDegree),dtype="i")
+
+    for j in range(0,ny+1):
+        for i in range(0,nx+1):
+            iVertex = i + j * (nx+1)
+            xVertex[iVertex] = float(i) * dc
+            yVertex[iVertex] = float(j) * dc
+
+            ic = i - 1 ; jc = j - 1
+            iCell = ic + jc * nx
+            if (ic >= 0 and jc >= 0):
+                cellsOnVertex[iVertex,0] = iCell
+            else:
+                cellsOnVertex[iVertex,0] = -1
+
+            ic = i ; jc = j - 1
+            iCell = ic + jc * nx
+            if (ic < nx and jc >= 0):
+                cellsOnVertex[iVertex,1] = iCell
+            else:
+                cellsOnVertex[iVertex,1] = -1
+
+            ic = i ; jc = j
+            iCell = ic + jc * nx
+            if (ic < nx and jc < ny):
+                cellsOnVertex[iVertex,2] = iCell
+            else:
+                cellsOnVertex[iVertex,2] = -1
+
+            ic = i - 1 ; jc = j
+            iCell = ic + jc * nx
+            if (ic >= 0 and jc < ny):
+                cellsOnVertex[iVertex,3] = iCell
+            else:
+                cellsOnVertex[iVertex,3] = -1
+
+            #print(iVertex+1,cellsOnVertex[iVertex,:]+1)
+
+    var = fileGrid.createVariable("xCell","d",dimensions=["nCells"])
+    var[:] = xCell[:]
+    var = fileGrid.createVariable("yCell","d",dimensions=["nCells"])
+    var[:] = yCell[:]
+    var = fileGrid.createVariable("zCell","d",dimensions=["nCells"])
+    var[:] = zCell[:]
+
+    var = fileGrid.createVariable("xVertex","d",dimensions=["nVertices"])
+    var[:] = xVertex[:]
+    var = fileGrid.createVariable("yVertex","d",dimensions=["nVertices"])
+    var[:] = yVertex[:]
+    var = fileGrid.createVariable("zVertex","d",dimensions=["nVertices"])
+    var[:] = zVertex[:]
+
+    cellsOnVertex[:] = cellsOnVertex[:] + 1
+    var = fileGrid.createVariable("cellsOnVertex","i",dimensions=["nVertices","vertexDegree"])
+    var[:] = cellsOnVertex[:]
+
+    fileGrid.close()
+
+    os.system("%s/mesh_tools/mesh_conversion_tools/MpasMeshConverter.x grid_in.nc %s" %(mpas_tools_dir, gridname))
+
+    filein = Dataset(gridname,"a")
+
+    filein.Lx = Lx
+    filein.Ly = Ly
+
+    filein.nx = nx
+    filein.ny = ny
+
+    filein.dc = dc
+
+    filein.close()
+
+#-------------------------------------------------------------------------------
+
+nGrid = 4
+
+# hex
+dc = 16000.0
+nx = 82
+ny = 94
+
+nxs = []
+nys = []
+dcs = []
+for i in range(0,nGrid):
+    nxs.append(nx)
+    nys.append(ny)
+    dcs.append(dc)
+    nx = nx*2
+    ny = ny*2
+    dc = dc/2
+
+for nx, ny, dc in zip(nxs, nys, dcs):
+   create_grid_hex(nx, ny, dc)
+
+# quad
+dc = 16000.0
+nx = 80
+ny = 80
+
+nxs = []
+nys = []
+dcs = []
+for i in range(0,nGrid):
+    nxs.append(nx)
+    nys.append(ny)
+    dcs.append(dc)
+    nx = nx*2
+    ny = ny*2
+    dc = dc/2
+
+for nx, ny, dc in zip(nxs, nys, dcs):
+   create_grid_quad(nx, ny, dc)

testing_and_setup/seaice/testcases/square/create_ics.py

@@ -0,0 +1,140 @@
+from netCDF4 import Dataset
+import numpy as np
+import math
+
+#-------------------------------------------------------------
+
+def wind_velocity(x, y, Lx, Ly):
+
+    t = 0.0#21600.0
+    tau = 4.0 * 24.0 * 3600.0
+
+    u = 5.0 + (math.sin((2.0 * math.pi * t) / tau) - 3.0) * math.sin((2.0 * math.pi * x) / Lx) * math.sin((math.pi * y) / Ly)
+    v = 5.0 + (math.sin((2.0 * math.pi * t) / tau) - 3.0) * math.sin((2.0 * math.pi * y) / Ly) * math.sin((math.pi * x) / Lx)
+
+    return u, v
+
+#-------------------------------------------------------------
+
+def ocean_currents(x, y, Lx, Ly):
+
+    u =  0.1 * ((2.0 * y - Ly) / Ly)
+    v = -0.1 * ((2.0 * x - Lx) / Lx)
+
+    return u, v
+
+#-------------------------------------------------------------
+
+def ice_concentration(x, y, Lx, Ly):
+
+    conc = max(min(x / Lx, 1.0),0.0)
+
+    return conc
+
+#-------------------------------------------------------------
+
+def create_ic(gridfile, icfile):
+
+    # load grid file
+    grid = Dataset(gridfile, "r")
+
+    xCell = grid.variables["xCell"][:]
+    yCell = grid.variables["yCell"][:]
+
+    xVertex = grid.variables["xVertex"][:]
+    yVertex = grid.variables["yVertex"][:]
+
+    nCells = len(grid.dimensions["nCells"])
+    nVertices = len(grid.dimensions["nVertices"])
+
+    # calculate output variables
+    uAirVelocity = np.empty(nCells)
+    vAirVelocity = np.empty(nCells)
+
+    uOceanVelocity = np.empty(nCells)
+    vOceanVelocity = np.empty(nCells)
+
+    iceConcentration = np.empty((nCells,1,1))
+    iceVolume        = np.empty((nCells,1,1))
+
+    fVertex = np.empty((nVertices,1,1))
+
+    xmin = np.amin(xVertex)
+    xmax = np.amax(xVertex)
+    Lx = xmax - xmin
+
+    print(xmin, xmax, Lx)
+
+    ymin = np.amin(yVertex)
+    ymax = np.amax(yVertex)
+    Ly = ymax - ymin
+
+    print(ymin, ymax, Ly)
+
+    for iCell in range(0,nCells):
+
+        x = xCell[iCell] - xmin
+        y = yCell[iCell] - ymin
+
+        uAirVelocity[iCell],   vAirVelocity[iCell]   = wind_velocity (x, y, Lx, Ly)
+        uOceanVelocity[iCell], vOceanVelocity[iCell] = ocean_currents(x, y, Lx, Ly)
+
+        iceConcentration[iCell,0,0] = ice_concentration(x, y, Lx, Ly)
+
+        iceVolume[iCell,0,0] = 2.0 * iceConcentration[iCell,0,0]
+
+    for iVertex in range(0,nVertices):
+        fVertex[iVertex] = 1.46e-4
+
+    # create output file
+    ic = Dataset(icfile, "w", format="NETCDF3_64BIT")
+    # format: 'NETCDF4', 'NETCDF4_CLASSIC', 'NETCDF3_CLASSIC', 'NETCDF3_64BIT_OFFSET' or 'NETCDF3_64BIT_DATA
+
+    ic.createDimension("nCells", nCells)
+    ic.createDimension("nVertices", nVertices)
+    ic.createDimension("nCategories", 1)
+    ic.createDimension("ONE", 1)
+    ic.createDimension("Time", None)
+
+    uAirVelocityVar = ic.createVariable("uAirVelocity", 'd', ('nCells'))
+    vAirVelocityVar = ic.createVariable("vAirVelocity", 'd', ('nCells'))
+    uAirVelocityVar[:] = uAirVelocity[:]
+    vAirVelocityVar[:] = vAirVelocity[:]
+
+    uOceanVelocityVar = ic.createVariable("uOceanVelocity", 'd', ('nCells'))
+    vOceanVelocityVar = ic.createVariable("vOceanVelocity", 'd', ('nCells'))
+    uOceanVelocityVar[:] = uOceanVelocity[:]
+    vOceanVelocityVar[:] = vOceanVelocity[:]
+
+    iceConcentrationVar = ic.createVariable("iceAreaCategory", 'd', ('nCells', 'nCategories', 'ONE'))
+    iceConcentrationVar[:,:,:] = iceConcentration[:,:,:]
+
+    iceVolumeVar = ic.createVariable("iceVolumeCategory", 'd', ('nCells', 'nCategories', 'ONE'))
+    iceVolumeVar[:,:,:] = iceVolume[:,:,:]
+
+    fVertexVar = ic.createVariable("fVertex", 'd', ('nVertices'))
+    fVertexVar[:] = fVertex[:]
+
+    ic.close()
+    grid.close()
+
+#-------------------------------------------------------------
+
+gridTypes = ["hex","quad"]
+
+grids = {"hex": ["0082x0094",
+                 "0164x0188",
+                 "0328x0376",
+                 "0656x0752"],
+         "quad":["0080x0080",
+                 "0160x0160",
+                 "0320x0320",
+                 "0640x0640"]}
+
+for gridType in gridTypes:
+    for grid in grids[gridType]:
+
+        gridfile = "grid_%s_%s.nc" %(gridType,grid)
+        icfile = "ic_%s_%s.nc" %(gridType,grid)
+
+        create_ic(gridfile, icfile)