Coarsen and overwritebcs modifications (#44)

Co-authored-by: josh <[email protected]> Co-authored-by: Anil Yildirim <[email protected]> Co-authored-by: Sabet Seraj <[email protected]>
mdolab · Feb 4, 2022 · d3948bb · d3948bb
1 parent 29d69b6
commit d3948bb
Show file tree

Hide file tree

Showing 6 changed files with 323 additions and 48 deletions.
diff --git a/cgnsutilities/cgns_utils.py b/cgnsutilities/cgns_utils.py
@@ -232,9 +232,8 @@ def get_parser():
         "overwriteBC", help="Overwrite boundary condition information", formatter_class=argparse.RawTextHelpFormatter
     )
     p_sub.add_argument("gridFile", help="Name of input CGNS file")
-    p_sub.add_argument(
-        "bcFile",
-        help="""File containing additional bc info. The file must consist of one or more lines contaning the following data:
+    bcFile_txt = """
+The file must consist of one or more lines contaning the following data:
 <blockID> <faceID> <BCType> <family> [dataset]
 
 where:
@@ -257,9 +256,7 @@ def get_parser():
     BCSetType     - bc dataset type. This is in most cases the same type as the BCType specified
     DirNeuArr     - can have one of two options: Dirichlet or Neumann
     DataArrNameN  - name of first property specified. This can be a range of things. Refer to ICEM or ADflow for supported BC properties
-    dataArrN      - the actual data for the property
-
-Note that only scalar values are supported at the moment for dataArrN.
+    dataArrN      - the actual data for the property. Either a scalar or a flattened nodal array. If an array is passed, the solver will convert the 1D array to the (possibly) 2D BC face.
 
 To find blockID of any mesh using Tecplot,
 1. Load the mesh with Advanced options > One Tecplot zone per non-poly CGNS zone/solution
@@ -270,10 +267,26 @@ def get_parser():
     7 kLow bcwallviscous wing
     4 jHigh bcsymmetryplane sym
     5 khigh bcoutflowsubsonic turb_inlet BCDataSet_1 BCInFlowSubsonic Dirichlet PressureStagnation 1234.0 TemperatureStagnation 4556.0
-""",
+"""
+    p_sub.add_argument(
+        "bcFile",
+        help="File containing additional bc info." + bcFile_txt,
     )
     p_sub.add_argument("outFile", nargs="?", default=None, help="Optional output file")
 
+    # ------------ Options for 'writebcinfo' mode --------------------
+    p_sub = subparsers.add_parser(
+        "writebcinfo",
+        help="Writes boundary condition information to a file.",
+        formatter_class=argparse.RawTextHelpFormatter,
+    )
+    p_sub.add_argument("gridFile", help="Name of input CGNS file")
+    p_sub.add_argument(
+        "bcOutFile",
+        default=None,
+        help="A file containing bc info." + bcFile_txt,
+    )
+
     # ------------ Options for 'rebunch' mode --------------------
     p_bunch = subparsers.add_parser("rebunch", help="Rebunch offwall spacing (experimental)")
     p_bunch.add_argument("gridFile", help="Name of input CGNS file")
@@ -873,7 +886,11 @@ def main():
     elif args.mode == "overwriteBC":
         curGrid.overwriteBCs(args.bcFile)
 
-    elif args.mode == "removeBC":
+    elif args.mode == "writebcinfo":
+        curGrid.writeBCs(args.bcOutFile)
+        sys.exit(0)
+
+    elif args.mode == "removebc":
         curGrid.removeBCs()
 
     elif args.mode == "rebunch":

diff --git a/cgnsutilities/cgnsutilities.py b/cgnsutilities/cgnsutilities.py
@@ -382,6 +382,13 @@ def removeBCs(self):
     def overwriteBCs(self, bcFile):
         """Overwrite BCs with information given in the file"""
 
+        def isFloat(string):
+            try:
+                float(string)
+                return True
+            except ValueError:
+                return False
+
         with open(bcFile, "r") as f:
             for line in f:
                 if line.strip():
@@ -399,13 +406,27 @@ def overwriteBCs(self, bcFile):
                         DirNeu = aux[6]
                         bocoDataSet = BocoDataSet(bocoSetName, BC[bocoDataSetType.lower()])
 
-                        for i in range(7, len(aux), 2):
+                        i = 7
+                        while i < len(aux):
                             arrayName = aux[i]
+                            i += 1
                             dType = CGNSDATATYPES["RealDouble"]
+                            nDims = 1
+
+                            dataArr = []
+
+                            for j in range(i, len(aux)):
+                                if isFloat(aux[j]):
+                                    dataArr.append(aux[j])
+                                    i += 1
+                                else:
+                                    break
+
+                            dataArr = numpy.array(dataArr, dtype=numpy.float64)
+
                             nDims = 1
                             dataDims = numpy.ones(3, dtype=numpy.int32, order="F")
-                            dataArr = numpy.zeros(1, dtype=numpy.float64, order="F")
-                            dataArr[0] = float(aux[i + 1])
+                            dataDims[0] = dataArr.size
 
                             bcDataArr = BocoDataSetArray(arrayName, dType, nDims, dataDims, dataArr)
                             if DirNeu == "Dirichlet":
@@ -419,6 +440,17 @@ def overwriteBCs(self, bcFile):
 
                     self.blocks[blockID].overwriteBCs(face, bocoType, family, dataSet)
 
+    def writeBCs(self, bcFile):
+        """write BCs to file"""
+        if bcFile is None:
+            bcFile = self.name + ".bc"
+
+        with open(bcFile, "w") as fid:
+
+            for idx_block, block in enumerate(self.blocks):
+                blockID = idx_block + 1
+                block.writeBCs(blockID, fid)
+
     def autoOversetBC(self, sym, connectSelf, tol):
         """This is essentially a simplified version of autoBC that flags all
         kMin faces as walls and all kMax faces as BCOverset"""
@@ -1631,34 +1663,56 @@ def coarsen(self):
         conditions and B2B if necessary"""
         # We will coarsen one direction at a time. We do this to check if the block
         # is already 1-cell wide, which can't be coarsened any further
-        if self.dims[0] != 2:
-            self.coords = self.coords[0::2, :, :, :]
-            # Update BCs and B2B
-            for boco in self.bocos:
-                boco.coarsen(0)
-            for b2b in self.B2Bs:
-                b2b.coarsen(0)
-
-        if self.dims[1] != 2:
-            self.coords = self.coords[:, 0::2, :, :]
-            # Update BCs and B2B
-            for boco in self.bocos:
-                boco.coarsen(1)
-            for b2b in self.B2Bs:
-                b2b.coarsen(1)
-
-        if self.dims[2] != 2:
-            self.coords = self.coords[:, :, 0::2, :]
-            # Update BCs and B2B
-            for boco in self.bocos:
-                boco.coarsen(2)
-            for b2b in self.B2Bs:
-                b2b.coarsen(2)
-
-        # Update dimensions
-        self.dims[0] = self.coords.shape[0]
-        self.dims[1] = self.coords.shape[1]
-        self.dims[2] = self.coords.shape[2]
+
+        # the new dimensions are half rounded up of the old dimensions
+        new_dims = copy.deepcopy(self.dims)
+        for i in range(3):
+            if self.dims[i] > 2:
+
+                if self.dims[i] % 2 == 0:
+                    print(
+                        f"INFO: unevenly coarsing block {self.name.decode('utf-8', 'ignore')} along dimension {i} (size {self.dims[i]}) "
+                    )
+
+                # The RHS takes odd numbers to *1/2 rounded up and even numbers to *1/2
+                # for example
+                # old dim: 0 1 2 3 4 5 6 7 8 9
+                # new dim: 1 1 1 2 2 3 3 4 4 5
+                new_dims[i] = int(numpy.ceil((self.dims[i]) / 2))
+
+        new_coords = numpy.zeros((new_dims[0], new_dims[1], new_dims[2], 3))
+
+        # Loop over all directions
+        s = slice(None)
+        fine_slicer = [s] * 3
+        for idx_dim in range(3):
+            # can this direction be coarsened?
+            if self.dims[idx_dim] > 2:
+                # set the slice in that direction to take every other point
+                fine_slicer[idx_dim] = slice(None, None, 2)
+
+        new_coords = self.coords[tuple(fine_slicer)]
+
+        # set the last point to be the same so we don't create any gaps if
+        # the number of points in that direction isn't odd
+        for idx_dim in range(3):
+            end_fine_slicer = copy.deepcopy(fine_slicer)
+            end_coarse_slicer = [s] * 3
+
+            end_fine_slicer[idx_dim] = -1
+            end_coarse_slicer[idx_dim] = -1
+
+            new_coords[tuple(end_coarse_slicer)] = self.coords[tuple(end_fine_slicer)]
+
+        # coarsen the bc and connectivity for the blk as well
+        for boco in self.bocos:
+            boco.coarsen()
+        for b2b in self.B2Bs:
+            b2b.coarsen()
+
+        self.coords = new_coords
+
+        self.dims = new_dims
 
     def refine(self, axes):
         """Refine the block uniformly. We will also update the
@@ -2135,6 +2189,66 @@ def overwriteBCs(self, face, bocoType, family, dataSet):
         self.addBoco(Boco("boco_%d" % self.bocoCounter, BC[bocoType.lower()], ptRange, family, dataSet))
         self.bocoCounter += 1
 
+    def writeBCs(self, blk_num, file_handle):
+        """write the bc data to a file"""
+
+        d = self.dims
+
+        for boco in self.bocos:
+            # check the point range to see what face it is on
+            ptRange = boco.ptRange
+            if (ptRange[0] == [1, 1]).all():
+                face = "ilow"
+            elif (ptRange[0] == [d[0], d[0]]).all():
+                face = "ihigh"
+
+            elif (ptRange[1] == [1, 1]).all():
+                face = "jlow"
+            elif (ptRange[1] == [d[1], d[1]]).all():
+                face = "jhigh"
+
+            elif (ptRange[2] == [1, 1]).all():
+                face = "klow"
+            elif (ptRange[2] == [d[2], d[2]]).all():
+                face = "khigh"
+            else:
+                raise ValueError("Face could not be determined to be one of (iLow, iHigh, jLow, jHigh, kLow, or kHigh)")
+
+            data_arr_str = ""
+            for data_arr in boco.dataSets:
+
+                # use the BC dictionary in reverse to find the bc type string
+                for bctype in BC:
+                    if data_arr.type == BC[bctype]:
+                        bctype_str = bctype
+                        break
+
+                data_arr_str += " " + data_arr.name.decode("utf-8", "ignore")
+                data_arr_str += " " + bctype_str
+
+                if data_arr.dirichletArrays:
+                    data_arr_str += " Dirichlet"
+
+                    for d_arr in data_arr.dirichletArrays:
+                        data_arr_str += " " + d_arr.name.decode("utf-8", "ignore")
+                        data_arr_str += " " + " ".join([f"{data}" for data in d_arr.dataArr])
+
+                if data_arr.neumannArrays:
+                    data_arr_str += " Neumann"
+
+                    for n_arr in data_arr.neumannArrays:
+                        data_arr_str += " " + n_arr.name
+                        data_arr_str += " " + " ".join([f"{data}" for data in n_arr.dataArr])
+
+            # use the BC dictionary in reverse to find the bc type string
+            for bctype in BC:
+                if boco.type == BC[bctype]:
+                    bctype_str = bctype
+                    break
+
+            fam_name = boco.family.strip().decode("utf-8", "ignore")
+            file_handle.write(f"{blk_num} {face} {bctype_str} {fam_name} {data_arr_str}\n")
+
     def isFaceInPtRange(self, face, ptRange):
         """
         Identifies if the provided face matches a given point range.
@@ -2323,10 +2437,39 @@ def addBocoDataSet(self, bocoDataSet):
         """Add a boundary condition dataset to this bc"""
         self.dataSets.append(bocoDataSet)
 
-    def coarsen(self, direction):
-        """Coarsen the range of the BC along the specified direction"""
-        for j in range(2):
-            self.ptRange[direction, j] = (self.ptRange[direction, j] - 1) // 2 + 1
+    def coarsen(self):
+        """Coarsen the range of the BC"""
+
+        for idim in range(3):
+
+            self.ptRange[idim, 0] = int(numpy.floor((self.ptRange[idim, 0]) / 2)) + 1
+            if self.ptRange[idim, 1] > 2:
+
+                # coarsen the data set if it is an array
+                if self.dataSets:
+                    for data_set in self.dataSets:
+                        for dir_arr in data_set.dirichletArrays:
+                            if numpy.prod(dir_arr.dataDimensions) == 1:
+                                # one value is being used for all points
+                                # thus, there is no need to coarsen the data
+                                continue
+
+                            slicer = [slice(None)] * 3
+                            slicer[idim] = slice(None, None, 2)
+
+                            data_mat = dir_arr.dataArr.reshape(self.ptRange[:, 1])
+                            new_data_mat = data_mat[tuple(slicer)]
+
+                            # make sure the last data point is always copied
+                            slicer[idim] = -1
+
+                            new_data_mat[tuple(slicer)] = data_mat[tuple(slicer)]
+
+                            dir_arr.dataDimensions[0] = numpy.prod(new_data_mat.shape)
+
+                            dir_arr.dataArr = new_data_mat.flatten()
+
+                self.ptRange[idim, 1] = int(numpy.ceil((self.ptRange[idim, 1]) / 2))
 
     def refine(self, axes):
         """refine the range of the BC"""
@@ -2400,12 +2543,19 @@ def __init__(self, connectName, donorName, ptRange, donorRange, transform):
         self.donorRange = donorRange
         self.transform = transform
 
-    def coarsen(self, direction):
+    def coarsen(self):
         """Coarsen the range of the B2B along the specified direction"""
-        donorDir = abs(self.transform[direction]) - 1
-        for j in range(2):
-            self.ptRange[direction, j] = (self.ptRange[direction, j] - 1) // 2 + 1
-            self.donorRange[donorDir, j] = (self.donorRange[donorDir, j] - 1) // 2 + 1
+        for idim in range(3):
+
+            donorDir = abs(self.transform[idim]) - 1
+
+            for j in range(2):
+
+                if self.ptRange[idim, j] > 2:
+                    self.ptRange[idim, j] = (self.ptRange[idim, j] - 1) // 2 + 1
+
+                if self.donorRange[donorDir, j] > 2:
+                    self.donorRange[donorDir, j] = (self.donorRange[donorDir, j] + 1) // 2
 
     def refine(self, axes):
         """refine the range of the B2B"""

diff --git a/examples/block_4x2x3.cgns b/examples/block_4x2x3.cgns
diff --git a/examples/hotwall_boco.info b/examples/hotwall_boco.info
@@ -0,0 +1,14 @@
+1 klow bcwallviscous wall 
+1 khigh bcfarfield Far 
+1 jhigh bcsymmetryplane Sym 
+2 klow bcwallviscous wall 
+2 khigh bcfarfield Far 
+2 jlow bcsymmetryplane Sym 
+3 klow bcwallviscous wall 
+3 khigh bcfarfield Far 
+3 jhigh bcsymmetryplane Sym 
+4 klow bcwallviscous wall 
+4 khigh bcfarfield Far 
+4 jhigh bcsymmetryplane Sym 
+5 klow bcwallviscousisothermal hotwall BCDataSet_t bcwallviscousisothermal Dirichlet Temperature 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
+5 khigh bcfarfield Far 
diff --git a/examples/overwriteBCs_array.sh b/examples/overwriteBCs_array.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# This script demonstrates how to run and use the overwritebc option to add array bc data
+CGNS_INFILE="717_wl_L2.cgns"
+CGNS_OUTFILE="717_wl_L2_overwriteBCs_array.cgns"
+BC_FILE="hotwall_boco.info"
+
+
+# Run the command
+cgns_utils overwriteBC $CGNS_INFILE $BC_FILE $CGNS_OUTFILE