From bb8fc364fe38aeb2881e1ea91f47ee1850aa0ab6 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Thu, 24 Mar 2022 19:40:34 -0400 Subject: [PATCH 01/31] Added copy/rotate utility --- utilities/copyRotateSectorGrid.py | 190 ++++++++++++++++++++++++++++++ utilities/interpolate.py | 2 +- 2 files changed, 191 insertions(+), 1 deletion(-) create mode 100755 utilities/copyRotateSectorGrid.py diff --git a/utilities/copyRotateSectorGrid.py b/utilities/copyRotateSectorGrid.py new file mode 100755 index 00000000..938f4722 --- /dev/null +++ b/utilities/copyRotateSectorGrid.py @@ -0,0 +1,190 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +"""This utility executes a copy/rotate operation on a sector grid. + +Requires a path to a "from" folder that contains g.* and conn.yaml files. + +The utility will output the copy/rotated grid in the "to" folder. + +""" + +import argparse +import peregrinepy as pg +import numpy as np +import os +from match_periodic_faces import * + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Interpolate from one grid to another") + parser.add_argument( + "-from", + "--fromDir", + action="store", + metavar="", + dest="fromDir", + default="./from", + help="Directory containing the g.*.h5 and q.*.h5 files to interpolate from. Default is ./from", + type=str, + ) + parser.add_argument( + "-to", + "--toDir", + action="store", + metavar="", + dest="toDir", + default="./to", + help="Directory containing the g.*.h5 files to interpolate to. Default is ./to", + type=str, + ) + parser.add_argument( + "-nseg", + "--numSegments", + action="store", + metavar="", + dest="nseg", + default="2", + help="Number of output segments (must be > 1)", + type=int, + ) + parser.add_argument( + "-sectorAngle", + "--sectorAngle", + action="store", + metavar="", + dest="sectorAngle", + help="Angle of sector in degrees.", + type=float, + ) + + args = parser.parse_args() + fromDir = args.fromDir + toDir = args.toDir + nseg = args.nseg + sectorAngle = args.sectorAngle + + if nseg <= 1: + raise ValueError( + "nseg must be >= 1 (it corresponds to the total number of output segments)" + ) + + nblks = len( + [f for f in os.listdir(fromDir) if f.startswith("g.") and f.endswith(".h5")] + ) + + fromGrid = pg.multiBlock.grid(nblks) + pg.readers.readGrid(fromGrid, fromDir) + pg.readers.readConnectivity(fromGrid, fromDir) + + toGrid = pg.multiBlock.grid(nblks * nseg) + + # Copy the original grid to new grid's first sector + # Also collect "front" and "back" faces. + backside = [] + frontside = [] + for i, fromBlk in enumerate(fromGrid): + toGrid[i] = fromBlk + + for face in fromBlk.faces: + if face.bcType == "b1": + if np.mean(fromBlk.array["z"]) > 0.0: + backside.append(fromBlk.nblki) + elif np.mean(fromBlk.array["z"] < 0.0): + frontside.append(fromBlk.nblki) + else: + print("HELP") + + # Now copy/rotate sector by sector + # We will update the front and back side + # faces and the new sector faces on this + # pass + for i in range(nseg - 1): + angle = sectorAngle * (i + 1) * np.pi / 180.0 + for j in range(nblks): + fromBlk = fromGrid[j] + + rotNblki = (i + 1) * nblks + j + rotBlk = toGrid[rotNblki] + + # copy/rotate block coordinates + rotBlk.ni = fromBlk.ni + rotBlk.nj = fromBlk.nj + rotBlk.nk = fromBlk.nk + + rotBlk.array["x"] = fromBlk.array["z"] * np.sin(angle) + fromBlk.array[ + "x" + ] * np.cos(angle) + rotBlk.array["y"] = fromBlk.array["y"] + rotBlk.array["z"] = fromBlk.array["z"] * np.cos(angle) - fromBlk.array[ + "x" + ] * np.sin(angle) + + # copy connectivity + for toFace, fromFace in zip(rotBlk.faces, fromBlk.faces): + + # treat boundary faces + if fromFace.neighbor is None: + toFace.bcFam = fromFace.bcFam + toFace.bcType = fromFace.bcType + toFace.orientation = None + toFace.neighbor = None + continue + + # treat the rotated faces + if fromFace.bcType == "b1": + toFace.bcFam = None + toFace.orientation = fromFace.orientation + # back side faces + if fromBlk.nblki in backside: + toFace.bcType = "b0" + toFace.neighbor = fromFace.neighbor + nblks * i + # front side faces + elif fromBlk.nblki in frontside: + if i == nseg - 1: + toFace.bcType = "b1" + toFace.neighbor = fromFace.neighbor + else: + toFace.bcType = "b0" + toFace.neighbor = fromFace.neighbor + nblks * (i + 2) + # treat the internal faces + elif fromFace.bcType == "b0": + toFace.bcType = "b0" + toFace.bcFam = None + toFace.orientation = fromFace.orientation + toFace.neighbor = fromFace.neighbor + nblks * (i + 1) + + # At this point, the original sector front and back side faces are out of + # date, as is the last sector's front side + if abs(nseg * sectorAngle - 360.0) < 1e-10: + is360 = True + else: + is360 = False + + for i, fromBlk in enumerate(fromGrid): + # backside + if fromBlk.nblki in backside: + rotBlk = toGrid[i] + for toFace, fromFace in zip(rotBlk.faces, fromBlk.faces): + if fromFace.bcType == "b1": + toFace.neighbor = fromFace.neighbor + nblks * (nseg - 1) + if is360: + toFace.bcType = "b0" + + # original and far frontside + elif fromBlk.nblki in frontside: + # original + rotBlk = toGrid[i] + for toFace, fromFace in zip(rotBlk.faces, fromBlk.faces): + if fromFace.bcType == "b1": + toFace.bcType = "b0" + toFace.neighbor = fromFace.neighbor + nblks + # far frontside + rotBlk = toGrid[i + nblks * (nseg - 1)] + for toFace, fromFace in zip(rotBlk.faces, fromBlk.faces): + if is360: + toFace.bcType = "b0" + + pg.writers.writeGrid(toGrid, toDir) + pg.writers.writeConnectivity(toGrid, toDir) + # match_faces(toGrid) diff --git a/utilities/interpolate.py b/utilities/interpolate.py index 09b5b00c..b137fa70 100755 --- a/utilities/interpolate.py +++ b/utilities/interpolate.py @@ -3,7 +3,7 @@ """This utility interpolates a PEREGRINE restart from one grid to another. -Requires a path to folder that contains g.* and dtms.* grid and restart files respectively. +Requires a path to folder that contains g.*.h5 and q.*.h5 grid and output files respectively. Also must know the number of species in the case. The utility will output the interpolated restart file in the "to" folder. From 5a2edaae884f132bc0fc0c7676a711f4f437a587 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Thu, 24 Mar 2022 19:47:34 -0400 Subject: [PATCH 02/31] Seems to be working. --- utilities/copyRotateSectorGrid.py | 4 ++-- utilities/verifyGrid.py | 28 ++++++++-------------------- 2 files changed, 10 insertions(+), 22 deletions(-) diff --git a/utilities/copyRotateSectorGrid.py b/utilities/copyRotateSectorGrid.py index 938f4722..a36ef647 100755 --- a/utilities/copyRotateSectorGrid.py +++ b/utilities/copyRotateSectorGrid.py @@ -120,7 +120,7 @@ "x" ] * np.sin(angle) - # copy connectivity + # transfer connectivity for toFace, fromFace in zip(rotBlk.faces, fromBlk.faces): # treat boundary faces @@ -141,7 +141,7 @@ toFace.neighbor = fromFace.neighbor + nblks * i # front side faces elif fromBlk.nblki in frontside: - if i == nseg - 1: + if i == nseg - 2: toFace.bcType = "b1" toFace.neighbor = fromFace.neighbor else: diff --git a/utilities/verifyGrid.py b/utilities/verifyGrid.py index 4624f6c4..4a53ccb0 100755 --- a/utilities/verifyGrid.py +++ b/utilities/verifyGrid.py @@ -1,21 +1,25 @@ #!/usr/bin/env python3 -"""This utility goes through a grid face by face, verifying that all the block's connectivities agree, +""" +This utility goes through a grid face by face, +verifying that all the block's connectivities agree, and that the coordinates of matching faces are identical. Inputs are the path to the grid files, and path to the conn.inp file. -It can handle b 1 periodicity in the conn.inp, however it will not compare the x,y,z coordinate locations +It can handle b1 periodicity in the conn.yaml, +however it will not compare the x,y,z coordinate locations of the faces. Output will print any discrepencies to the screen """ -import os import argparse -import peregrinepy as pg +import os + import numpy as np +import peregrinepy as pg def verify(mb): @@ -37,22 +41,6 @@ def verify(mb): 5: 2, 6: 2, } - orientToSmallFaceMapping = { - 1: 2, - 2: 4, - 3: 6, - 4: 1, - 5: 3, - 6: 5, - } - orientToLargeFaceMapping = { - 1: 1, - 2: 3, - 3: 5, - 4: 2, - 5: 4, - 6: 6, - } largeIndexMapping = {0: "k", 1: "k", 2: "j"} needToTranspose = { From 8552a9786a5ebde93eacabca774d39d3b8f3c8a0 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Thu, 24 Mar 2022 20:01:46 -0400 Subject: [PATCH 03/31] Removed un needed import --- utilities/copyRotateSectorGrid.py | 1 - 1 file changed, 1 deletion(-) diff --git a/utilities/copyRotateSectorGrid.py b/utilities/copyRotateSectorGrid.py index a36ef647..f6284103 100755 --- a/utilities/copyRotateSectorGrid.py +++ b/utilities/copyRotateSectorGrid.py @@ -13,7 +13,6 @@ import peregrinepy as pg import numpy as np import os -from match_periodic_faces import * if __name__ == "__main__": From 1c2e0a99062339f10f53b9af8ac5e2559d30f03f Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Fri, 25 Mar 2022 07:19:06 -0400 Subject: [PATCH 04/31] Got annulus creator working. --- src/peregrinepy/grid/create.py | 77 ++++++++++++++++++---------------- 1 file changed, 41 insertions(+), 36 deletions(-) diff --git a/src/peregrinepy/grid/create.py b/src/peregrinepy/grid/create.py index ecec9672..cd966f6b 100644 --- a/src/peregrinepy/grid/create.py +++ b/src/peregrinepy/grid/create.py @@ -247,8 +247,6 @@ def multiBlockCube( def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): - raise ValueError("The annulus is jacked, needs to be updated.") - """Function to populate the coordinate arrays of a provided peregrinepy.grid.grid_block in the shape of an annulus with prescribed location, extents, and discretization. If the input multiBlock object is a restart block the shape and size of the flow data arrays are also updated. @@ -306,10 +304,20 @@ def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): n12 = (p2 - p1) / np.linalg.norm(p2 - p1) n13 = (p3 - p1) / np.linalg.norm(p3 - p1) - shape = (blk.ni + 2, blk.nj + 2, blk.nk + 2) - blk.array["x"] = np.zeros(shape) - blk.array["y"] = np.zeros(shape) - blk.array["z"] = np.zeros(shape) + blk.ni = dimensions[0] + blk.nj = dimensions[1] + blk.nk = dimensions[2] + if blk.blockType == "solver": + ng = blk.ng + else: + ng = 0 + + blk.initGridArrays() + + if blk.blockType == "solver": + s_i = np.s_[ng:-ng, ng:-ng, ng:-ng] + else: + s_i = np.s_[:, :, :] dx = np.linalg.norm(p2 - p1) / (blk.ni - 1) dr = thickness / (blk.nj - 1) @@ -319,17 +327,17 @@ def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): for i in range(blk.ni): p_ij = np.append(p3 + dx * i * n12 + dr * j * n13, 1) - blk.array["x"][i, j, 0] = p_ij[0] - blk.array["y"][i, j, 0] = p_ij[1] - blk.array["z"][i, j, 0] = p_ij[2] + blk.array["x"][s_i][i, j, 0] = p_ij[0] + blk.array["y"][s_i][i, j, 0] = p_ij[1] + blk.array["z"][s_i][i, j, 0] = p_ij[2] - xflat = np.reshape(blk.array["x"][:, :, 0], (blk.ni * blk.nj, 1)) - yflat = np.reshape(blk.array["y"][:, :, 0], (blk.ni * blk.nj, 1)) - zflat = np.reshape(blk.array["z"][:, :, 0], (blk.ni * blk.nj, 1)) + xflat = np.reshape(blk.array["x"][s_i][:, :, 0], (blk.ni * blk.nj, 1)) + yflat = np.reshape(blk.array["y"][s_i][:, :, 0], (blk.ni * blk.nj, 1)) + zflat = np.reshape(blk.array["z"][s_i][:, :, 0], (blk.ni * blk.nj, 1)) pts = np.hstack((xflat, yflat, zflat)) p = pts - p1 - shape = blk.array["x"][:, :, 0].shape + shape = blk.array["x"][s_i][:, :, 0].shape for k in range(1, blk.nk): # See http://paulbourke.net/geometry/rotate/ @@ -354,16 +362,17 @@ def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): q[:, 1] += p1[1] q[:, 2] += p1[2] - blk.array["x"][:, :, k] = np.reshape(q[:, 0], shape) - blk.array["y"][:, :, k] = np.reshape(q[:, 1], shape) - blk.array["z"][:, :, k] = np.reshape(q[:, 2], shape) + blk.array["x"][s_i][:, :, k] = np.reshape(q[:, 0], shape) + blk.array["y"][s_i][:, :, k] = np.reshape(q[:, 1], shape) + blk.array["z"][s_i][:, :, k] = np.reshape(q[:, 2], shape) def multiBlockAnnulus( mb, p1, p2, p3, sweep, thickness, mbDims, dimsPerBlock, periodic=False ): - """Function to populate the coordinate arrays of a peregrinepy.multiBlock.grid (or one of its descendants) in the shape + """ + Function to populate the coordinate arrays of a peregrinepy.multiBlock.grid (or one of its descendants) in the shape of an annulus with prescribed location, extents, and discretization split into as manj blocks as mb.nblks. Will also update connectivity of interblock faces. If the input multiBlock object is a restart block the shape and size of the flow data arrays are also updated. @@ -379,7 +388,7 @@ def multiBlockAnnulus( p2 : list, tuple List/tuple of length 3 containing the location of the end of the annulus to be created, i.e. - the center of the end of the the whole. + the axial center of the end of the the annulus. p3 : list, tuple List/tuple of length 3 containing the location of a point orthogonal to the line (p1,p2) marking @@ -401,11 +410,14 @@ def multiBlockAnnulus( List/tuple of length 3 containing number of blocks in axial direction, radial direction, and theta direction. NOTE: product of mbDims must equal mb.nblks! - dimensions : list, tuple + dimsPerBlock : list, tuple List/tuple of length 3 containing discretization (ni,nj,nk) in each dimension of every block (all will be uniform). Where the "x,i,xi" direction is along the annulus axis, the "y,j,eta" direction is along the radial direction, and the "z,k,zeta" direction is along the theta direction. + periodic : bool + Whether domain is periodic about the rotational axis. + Returns ------- None @@ -470,10 +482,6 @@ def multiBlockAnnulus( blkNum = k * mbDims[1] * mbDims[0] + j * mbDims[0] + i blk = mb[blkNum] - blk.nblki = blkNum - blk.ni = dimsPerBlock[0] - blk.nj = dimsPerBlock[1] - blk.nk = dimsPerBlock[2] newp1 = p1 + dx * i * n12 newp2 = p1 + dx * (i + 1) * n12 @@ -482,30 +490,27 @@ def multiBlockAnnulus( annulus(blk, newp1, newp2, newp3, dtheta, dr, dimsPerBlock) # Update connectivity - conn = blk.connectivity cubicConnectivity(blk, mbDims, blkNum, i, j, k) # k faces if connect: if k == 0: + face = blk.getFace(5) if float(sweep) == 360.0: - conn["5"]["bc"] = "b0" + face.bcType = "b0" else: - conn["5"]["bc"] = "b1" - conn["5"]["neighbor"] = blkNum + mbDims[0] * mbDims[1] * ( - mbDims[2] - 1 - ) - conn["5"]["orientation"] = "123" + face.bcType = "b1" + face.neighbor = blkNum + mbDims[0] * mbDims[1] * (mbDims[2] - 1) + face.orientation = "123" elif k == mbDims[2] - 1: + face = blk.getFace(6) if float(sweep) == 360.0: - conn["6"]["bc"] = "b0" + face.bcType = "b0" else: - conn["6"]["bc"] = "b1" - conn["6"]["neighbor"] = blkNum - mbDims[0] * mbDims[1] * ( - mbDims[2] - 1 - ) - conn["6"]["orientation"] = "123" + face.bcType = "b1" + face.neighbor = blkNum - mbDims[0] * mbDims[1] * (mbDims[2] - 1) + face.orientation = "123" for blk in mb: if blk.blockType == "solver" and blk._isInitialized: From ff1096cb9d1aa88f45c870796ce70eea4aaf5038 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Fri, 25 Mar 2022 17:22:12 -0400 Subject: [PATCH 05/31] Forgot fstring --- utilities/verifyGrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utilities/verifyGrid.py b/utilities/verifyGrid.py index 4a53ccb0..71a03b84 100755 --- a/utilities/verifyGrid.py +++ b/utilities/verifyGrid.py @@ -186,7 +186,7 @@ def extractFace(blk, nface): cross = np.cross(vI, vJ) if np.dot(vK, cross) < 0.0: - print("Warning, block {blk.nbkli} is left handed. This must be fixed.") + print(f"Warning, block {blk.nbkli} is left handed. This must be fixed.") warn = True if warn: From 0b6d19f63ab0d26ff810bacaef84f17cf6d9020e Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Fri, 25 Mar 2022 17:22:36 -0400 Subject: [PATCH 06/31] type --- utilities/verifyGrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utilities/verifyGrid.py b/utilities/verifyGrid.py index 71a03b84..1694fd1c 100755 --- a/utilities/verifyGrid.py +++ b/utilities/verifyGrid.py @@ -186,7 +186,7 @@ def extractFace(blk, nface): cross = np.cross(vI, vJ) if np.dot(vK, cross) < 0.0: - print(f"Warning, block {blk.nbkli} is left handed. This must be fixed.") + print(f"Warning, block {blk.nblki} is left handed. This must be fixed.") warn = True if warn: From 4162e3caf234a37e2a84de2fb47a36b5a11e44e6 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Sun, 3 Apr 2022 16:17:57 -0400 Subject: [PATCH 07/31] Added periodic info to grid creation, topo face, conn writer, and bcFam reader. --- src/peregrinepy/grid/create.py | 193 ++++++++++--------- src/peregrinepy/multiBlock/topologyFace.py | 7 +- src/peregrinepy/readers/readBcs.py | 8 +- src/peregrinepy/readers/readConnectivity.py | 6 + src/peregrinepy/writers/writeConnectivity.py | 16 +- 5 files changed, 135 insertions(+), 95 deletions(-) diff --git a/src/peregrinepy/grid/create.py b/src/peregrinepy/grid/create.py index cd966f6b..58002b58 100644 --- a/src/peregrinepy/grid/create.py +++ b/src/peregrinepy/grid/create.py @@ -11,106 +11,111 @@ import numpy as np -def cubicConnectivity(blk, mbDims, blkNum, i, j, k): +def cubicConnectivity( + blk, mbDims, blkNum, i, j, k, periodicI=False, periodicJ=False, periodicK=False +): # i faces - if i == 0 and mbDims[0] != 1: - face = blk.getFace(1) - face.bcType = "adiabaticNoSlipWall" - face.neighbor = None - face.orientation = None - - face = blk.getFace(2) - face.bcType = "b0" - face.neighbor = blkNum + 1 - face.orientation = "123" - - elif i == mbDims[0] - 1 and mbDims[0] != 1: - face = blk.getFace(1) - face.bcType = "b0" - face.neighbor = blkNum - 1 - face.orientation = "123" - - face = blk.getFace(2) - face.bcType = "adiabaticNoSlipWall" - face.neighbor = None - face.orientation = None - - elif mbDims[0] != 1: - face = blk.getFace(1) + # face 1 + face = blk.getFace(1) + if i == 0: + if periodicI: + face.bcType = "periodicTrans" + face.neighbor = blkNum + (mbDims[0] - 1) + face.orientation = "123" + face.isPeriodicLow = True + else: + face.bcType = "adiabaticNoSlipWall" + face.neighbor = None + face.orientation = None + else: face.bcType = "b0" face.neighbor = blkNum - 1 face.orientation = "123" - face = blk.getFace(2) + # face 2 + face = blk.getFace(2) + if i == mbDims[0] - 1: + if periodicI: + face.bcType = "periodicTrans" + face.neighbor = blkNum - (mbDims[0] - 1) + face.orientation = "123" + face.isPeriodicLow = False + else: + face.bcType = "adiabaticNoSlipWall" + face.neighbor = None + face.orientation = None + else: face.bcType = "b0" face.neighbor = blkNum + 1 face.orientation = "123" # j faces - if j == 0 and mbDims[1] != 1: - face = blk.getFace(3) - face.bcType = "adiabaticNoSlipWall" - face.neighbor = None - face.orientation = None - - face = blk.getFace(4) - face.bcType = "b0" - face.neighbor = blkNum + mbDims[0] - face.orientation = "123" - - elif j == mbDims[1] - 1 and mbDims[1] != 1: - face = blk.getFace(3) - face.bcType = "b0" - face.neighbor = blkNum - mbDims[0] - face.orientation = "123" - - face = blk.getFace(4) - face.bcType = "adiabaticNoSlipWall" - face.neighbor = None - face.orientation = None - - elif mbDims[1] != 1: - face = blk.getFace(3) + # face 3 + face = blk.getFace(3) + if j == 0: + if periodicJ: + face.bcType = "periodicTrans" + face.neighbor = blkNum + mbDims[0] * (mbDims[1] - 1) + face.orientation = "123" + face.isPeriodicLow = True + else: + face.bcType = "adiabaticNoSlipWall" + face.neighbor = None + face.orientation = None + else: face.bcType = "b0" face.neighbor = blkNum - mbDims[0] face.orientation = "123" - face = blk.getFace(4) + # face 4 + face = blk.getFace(4) + if j == mbDims[1] - 1: + if periodicJ: + face.bcType = "periodicTrans" + face.neighbor = blkNum - mbDims[0] * (mbDims[1] - 1) + face.orientation = "123" + face.isPeriodicLow = False + else: + face.bcType = "adiabaticNoSlipWall" + face.neighbor = None + face.orientation = None + else: face.bcType = "b0" face.neighbor = blkNum + mbDims[0] face.orientation = "123" # k faces - if k == 0 and mbDims[2] != 1: - face = blk.getFace(5) - face.bcType = "adiabaticNoSlipWall" - face.neighbor = None - face.orientation = None - - face = blk.getFace(6) - face.bcType = "b0" - face.neighbor = blkNum + mbDims[0] * mbDims[1] - face.orientation = "123" - - elif k == mbDims[2] - 1 and mbDims[2] != 1: - face = blk.getFace(5) - face.bcType = "b0" - face.neighbor = blkNum - mbDims[0] * mbDims[1] - face.orientation = "123" - - face = blk.getFace(6) - face.bcType = "adiabaticNoSlipWall" - face.neighbor = None - face.orientation = None - - elif mbDims[2] != 1: - face = blk.getFace(5) + # face 5 + face = blk.getFace(5) + if k == 0: + if periodicK: + face.bcType = "periodicTrans" + face.neighbor = blkNum + mbDims[0] * mbDims[1] * (mbDims[2] - 1) + face.orientation = "123" + face.isPeriodicLow = True + else: + face.bcType = "adiabaticNoSlipWall" + face.neighbor = None + face.orientation = None + else: face.bcType = "b0" face.neighbor = blkNum - mbDims[0] * mbDims[1] face.orientation = "123" - face = blk.getFace(6) + # face 6 + face = blk.getFace(6) + if k == mbDims[2] - 1: + if periodicK: + face.bcType = "periodicTrans" + face.neighbor = blkNum - mbDims[0] * mbDims[1] * (mbDims[2] - 1) + face.orientation = "123" + face.isPeriodicLow = False + else: + face.bcType = "b0" + face.neighbor = blkNum + mbDims[0] * mbDims[1] + face.orientation = "123" + else: face.bcType = "b0" face.neighbor = blkNum + mbDims[0] * mbDims[1] face.orientation = "123" @@ -174,6 +179,7 @@ def multiBlockCube( lengths=[1, 1, 1], mbDims=[1, 1, 1], dimsPerBlock=[10, 10, 10], + periodic=[False, False, False], ): """Function to populate the coordinate arrays of a peregrinepy.multiBlock.grid (or one of its descendants) in the shape of a cube @@ -198,6 +204,9 @@ def multiBlockCube( dimsPerBlock : list, tuple List/tuple of length 3 containing discretization (nx,nj,nk) in each dimension of each block to be created. + periodic: list[bool] + Whether any of the axes are periodic [I,J,K] + Returns ------- None @@ -237,7 +246,21 @@ def multiBlockCube( cube(blk, origin, lengths, dimensions) # Update connectivity - cubicConnectivity(blk, mbDims, blkNum, i, j, k) + cubicConnectivity( + blk, mbDims, blkNum, i, j, k, periodic[0], periodic[1], periodic[2] + ) + + # Set the peiodic info + for face in blk.faces: + if face.bcType == "periodicTrans" and face.nface in [1, 2]: + face.periodicAxis = np.array([1.0, 0.0, 0.0]) + face.periodicSpan = lengths[0] + elif face.bcType == "periodicTrans" and face.nface in [3, 4]: + face.periodicAxis = np.array([0.0, 1.0, 0.0]) + face.periodicSpan = lengths[1] + elif face.bcType == "periodicTrans" and face.nface in [5, 6]: + face.periodicAxis = np.array([0.0, 0.0, 1.0]) + face.periodicSpan = lengths[2] for blk in mb: if blk.blockType == "solver" and blk._isInitialized: @@ -490,27 +513,25 @@ def multiBlockAnnulus( annulus(blk, newp1, newp2, newp3, dtheta, dr, dimsPerBlock) # Update connectivity - cubicConnectivity(blk, mbDims, blkNum, i, j, k) + cubicConnectivity(blk, mbDims, blkNum, i, j, k, periodicK=connect) - # k faces + # Update the k faces if necessary if connect: + if k == 0: face = blk.getFace(5) if float(sweep) == 360.0: face.bcType = "b0" else: - face.bcType = "b1" - face.neighbor = blkNum + mbDims[0] * mbDims[1] * (mbDims[2] - 1) - face.orientation = "123" - + face.bcType = "periodicRot" elif k == mbDims[2] - 1: face = blk.getFace(6) if float(sweep) == 360.0: face.bcType = "b0" else: - face.bcType = "b1" - face.neighbor = blkNum - mbDims[0] * mbDims[1] * (mbDims[2] - 1) - face.orientation = "123" + face.bcType = "periodicRot" + face.periodicAxis = n12 + face.periodicSpan = sweep for blk in mb: if blk.blockType == "solver" and blk._isInitialized: diff --git a/src/peregrinepy/multiBlock/topologyFace.py b/src/peregrinepy/multiBlock/topologyFace.py index a325d7f5..307e826e 100644 --- a/src/peregrinepy/multiBlock/topologyFace.py +++ b/src/peregrinepy/multiBlock/topologyFace.py @@ -11,6 +11,10 @@ def __init__(self, nface): self._neighbor = None self._orientation = None + self.periodicAxis = None + self.periodicSpan = None + self.isPeriodicLow = None + @property def nface(self): return self._nface @@ -44,7 +48,8 @@ def bcType(self, value): validBcTypes = ( # Interior, periodic "b0", - "b1", + "periodicTrans", + "periodicRot", # Inlets "constantVelocitySubsonicInlet", "supersonicInlet", diff --git a/src/peregrinepy/readers/readBcs.py b/src/peregrinepy/readers/readBcs.py index 42611d83..61c27334 100644 --- a/src/peregrinepy/readers/readBcs.py +++ b/src/peregrinepy/readers/readBcs.py @@ -41,10 +41,10 @@ def readBcs(mb, pathToFile): for face in blk.faces: bcFam = face.bcFam bcType = face.bcType + # Not all bc types need inputs from bcFam, so we check here if bcFam is None: if bcType not in ( "b0", - "b1", "supersonicExit", "adiabaticNoSlipWall", "adiabaticSlipWall", @@ -65,6 +65,12 @@ def readBcs(mb, pathToFile): print(f"Warning, no values found for {bcsIn[bcFam]}") continue + # Add the information from any periodic faces + if bcType in ["periodicTrans", "periodicRot"]: + face.periodicAxis = bcsIn[bcFam]["periodicAxis"] + face.periodicSpan = bcsIn[bcFam]["periodicSpan"] + continue + # If we are a solver face, we need to create the kokkos arrays if face.faceType != "solver": continue diff --git a/src/peregrinepy/readers/readConnectivity.py b/src/peregrinepy/readers/readConnectivity.py index f7010b53..88e2cf28 100644 --- a/src/peregrinepy/readers/readConnectivity.py +++ b/src/peregrinepy/readers/readConnectivity.py @@ -49,4 +49,10 @@ def readConnectivity(mb, pathToFile, parallel=False): face.neighbor = fdict["neighbor"] face.orientation = fdict["orientation"] + try: + face.isPeriodicLow = fdict["isPeriodicLow"] + assert type(face.isPeriodicLow) is bool + except KeyError: + pass + mb.totalBlocks = conn["Total_Blocks"] diff --git a/src/peregrinepy/writers/writeConnectivity.py b/src/peregrinepy/writers/writeConnectivity.py index 8eac9d95..71a34ad8 100644 --- a/src/peregrinepy/writers/writeConnectivity.py +++ b/src/peregrinepy/writers/writeConnectivity.py @@ -26,13 +26,15 @@ def writeConnectivity(mb, path="./"): for blk in mb: conn[f"Block{blk.nblki}"] = {} bdict = conn[f"Block{blk.nblki}"] - for fc in blk.faces: - bdict[f"Face{fc.nface}"] = {} - fdict = bdict[f"Face{fc.nface}"] - fdict["bcType"] = fc.bcType - fdict["bcFam"] = fc.bcFam - fdict["neighbor"] = fc.neighbor - fdict["orientation"] = fc.orientation + for face in blk.faces: + bdict[f"Face{face.nface}"] = {} + fdict = bdict[f"Face{face.nface}"] + fdict["bcType"] = face.bcType + fdict["bcFam"] = face.bcFam + fdict["neighbor"] = face.neighbor + fdict["orientation"] = face.orientation + if face.isPeriodicLow is not None: + fdict["isPeriodicLow"] = face.isPeriodicLow with open(f"{path}/conn.yaml", "w") as f: yaml.dump(conn, f, sort_keys=False) From e91ff8376cec6183bd1b395f2188a1016e70e2be Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Mon, 4 Apr 2022 07:24:45 -0400 Subject: [PATCH 08/31] Add defaults to annulus --- src/peregrinepy/grid/create.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/peregrinepy/grid/create.py b/src/peregrinepy/grid/create.py index 58002b58..aa0b9893 100644 --- a/src/peregrinepy/grid/create.py +++ b/src/peregrinepy/grid/create.py @@ -391,7 +391,15 @@ def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): def multiBlockAnnulus( - mb, p1, p2, p3, sweep, thickness, mbDims, dimsPerBlock, periodic=False + mb, + p1=[0, 0, 0], + p2=[0, 1, 0], + p3=[0, 0, 1], + sweep=45, + thickness=0.1, + mbDims=[1, 1, 1], + dimsPerBlock=[10, 10, 10], + periodic=False, ): """ From 1ac713aa4b542036abac52d77364c91782666564 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Tue, 5 Apr 2022 10:00:31 -0400 Subject: [PATCH 09/31] bcType will now tell us whether face is periodic Low or High Added rotation matrix when setting periodic axis and span. --- src/peregrinepy/grid/create.py | 48 ++++++++++---------- src/peregrinepy/multiBlock/topologyFace.py | 47 +++++++++++++++++-- src/peregrinepy/readers/readConnectivity.py | 6 --- src/peregrinepy/writers/writeConnectivity.py | 2 - 4 files changed, 68 insertions(+), 35 deletions(-) diff --git a/src/peregrinepy/grid/create.py b/src/peregrinepy/grid/create.py index aa0b9893..4ed718b0 100644 --- a/src/peregrinepy/grid/create.py +++ b/src/peregrinepy/grid/create.py @@ -20,10 +20,9 @@ def cubicConnectivity( face = blk.getFace(1) if i == 0: if periodicI: - face.bcType = "periodicTrans" + face.bcType = "periodicTransLow" face.neighbor = blkNum + (mbDims[0] - 1) face.orientation = "123" - face.isPeriodicLow = True else: face.bcType = "adiabaticNoSlipWall" face.neighbor = None @@ -37,10 +36,9 @@ def cubicConnectivity( face = blk.getFace(2) if i == mbDims[0] - 1: if periodicI: - face.bcType = "periodicTrans" + face.bcType = "periodicTransHigh" face.neighbor = blkNum - (mbDims[0] - 1) face.orientation = "123" - face.isPeriodicLow = False else: face.bcType = "adiabaticNoSlipWall" face.neighbor = None @@ -55,10 +53,9 @@ def cubicConnectivity( face = blk.getFace(3) if j == 0: if periodicJ: - face.bcType = "periodicTrans" + face.bcType = "periodicTransLow" face.neighbor = blkNum + mbDims[0] * (mbDims[1] - 1) face.orientation = "123" - face.isPeriodicLow = True else: face.bcType = "adiabaticNoSlipWall" face.neighbor = None @@ -72,10 +69,9 @@ def cubicConnectivity( face = blk.getFace(4) if j == mbDims[1] - 1: if periodicJ: - face.bcType = "periodicTrans" + face.bcType = "periodicTransHigh" face.neighbor = blkNum - mbDims[0] * (mbDims[1] - 1) face.orientation = "123" - face.isPeriodicLow = False else: face.bcType = "adiabaticNoSlipWall" face.neighbor = None @@ -90,10 +86,9 @@ def cubicConnectivity( face = blk.getFace(5) if k == 0: if periodicK: - face.bcType = "periodicTrans" + face.bcType = "periodicTransLow" face.neighbor = blkNum + mbDims[0] * mbDims[1] * (mbDims[2] - 1) face.orientation = "123" - face.isPeriodicLow = True else: face.bcType = "adiabaticNoSlipWall" face.neighbor = None @@ -107,10 +102,9 @@ def cubicConnectivity( face = blk.getFace(6) if k == mbDims[2] - 1: if periodicK: - face.bcType = "periodicTrans" + face.bcType = "periodicTransHigh" face.neighbor = blkNum - mbDims[0] * mbDims[1] * (mbDims[2] - 1) face.orientation = "123" - face.isPeriodicLow = False else: face.bcType = "b0" face.neighbor = blkNum + mbDims[0] * mbDims[1] @@ -252,15 +246,21 @@ def multiBlockCube( # Set the peiodic info for face in blk.faces: - if face.bcType == "periodicTrans" and face.nface in [1, 2]: - face.periodicAxis = np.array([1.0, 0.0, 0.0]) + if face.bcType.startswith("periodicTrans") and face.nface in [1, 2]: face.periodicSpan = lengths[0] - elif face.bcType == "periodicTrans" and face.nface in [3, 4]: - face.periodicAxis = np.array([0.0, 1.0, 0.0]) + face.periodicAxis = np.array([1.0, 0.0, 0.0]) + elif face.bcType.startswith("periodicTrans") and face.nface in [ + 3, + 4, + ]: face.periodicSpan = lengths[1] - elif face.bcType == "periodicTrans" and face.nface in [5, 6]: - face.periodicAxis = np.array([0.0, 0.0, 1.0]) + face.periodicAxis = np.array([0.0, 1.0, 0.0]) + elif face.bcType.startswith("periodicTrans") and face.nface in [ + 5, + 6, + ]: face.periodicSpan = lengths[2] + face.periodicAxis = np.array([0.0, 0.0, 1.0]) for blk in mb: if blk.blockType == "solver" and blk._isInitialized: @@ -393,8 +393,8 @@ def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): def multiBlockAnnulus( mb, p1=[0, 0, 0], - p2=[0, 1, 0], - p3=[0, 0, 1], + p2=[1, 0, 0], + p3=[0, 1, 0], sweep=45, thickness=0.1, mbDims=[1, 1, 1], @@ -531,15 +531,17 @@ def multiBlockAnnulus( if float(sweep) == 360.0: face.bcType = "b0" else: - face.bcType = "periodicRot" + face.bcType = "periodicRotLow" + face.periodicSpan = sweep + face.periodicAxis = n12 elif k == mbDims[2] - 1: face = blk.getFace(6) if float(sweep) == 360.0: face.bcType = "b0" else: - face.bcType = "periodicRot" - face.periodicAxis = n12 + face.bcType = "periodicRotHigh" face.periodicSpan = sweep + face.periodicAxis = n12 for blk in mb: if blk.blockType == "solver" and blk._isInitialized: diff --git a/src/peregrinepy/multiBlock/topologyFace.py b/src/peregrinepy/multiBlock/topologyFace.py index 307e826e..91c6c2f4 100644 --- a/src/peregrinepy/multiBlock/topologyFace.py +++ b/src/peregrinepy/multiBlock/topologyFace.py @@ -11,9 +11,9 @@ def __init__(self, nface): self._neighbor = None self._orientation = None - self.periodicAxis = None + self.periodicOrigin = [0, 0, 0] self.periodicSpan = None - self.isPeriodicLow = None + self._periodicAxis = None @property def nface(self): @@ -48,8 +48,10 @@ def bcType(self, value): validBcTypes = ( # Interior, periodic "b0", - "periodicTrans", - "periodicRot", + "periodicTransLow", + "periodicTransHigh", + "periodicRotLow", + "periodicRotHigh", # Inlets "constantVelocitySubsonicInlet", "supersonicInlet", @@ -93,6 +95,43 @@ def orientation(self, value): raise TypeError(f"orientation must be a str not {type(value)}.") self._orientation = value + # Periodic stuff + @property + def periodicAxis(self): + return self._periodicAxis + + @periodicAxis.setter + def periodicAxis(self, axis): + import numpy as np + + axis = axis / np.linalg.norm(np.array(axis)) + self._periodicAxis = axis + + if not self.bcType.startswith("periodicRot"): + return + elif self.periodicSpan is None: + raise AttributeError("Please set periodicSpan before setting periodicAxis") + + # Compute rotation matrix + rot = np.zeros((3, 3)) + th = self.periodicSpan * np.pi / 180.0 + ct = np.cos(th) + st = np.sin(th) + ux, uy, uz = tuple(axis) + rot[0, 0] = ct + ux ** 2 * (1 - ct) + rot[0, 1] = ux * uy * (1 - ct) * uz * st + rot[0, 2] = ux * uz * (1 - ct) + uy * st + + rot[1, 0] = uy * ux * (1 - ct) + uz * st + rot[1, 1] = ct + uy ** 2 * (1 - ct) + rot[1, 2] = uy * uz * (1 - ct) - ux * st + + rot[2, 0] = uz * ux * (1 - ct) - uy * st + rot[2, 1] = uz * uy * (1 - ct) + ux * st + rot[2, 2] = ct + uz ** 2 * (1 - ct) + + self.periodicRotMatrix = rot + @property def neighborNface(self): diff --git a/src/peregrinepy/readers/readConnectivity.py b/src/peregrinepy/readers/readConnectivity.py index 88e2cf28..f7010b53 100644 --- a/src/peregrinepy/readers/readConnectivity.py +++ b/src/peregrinepy/readers/readConnectivity.py @@ -49,10 +49,4 @@ def readConnectivity(mb, pathToFile, parallel=False): face.neighbor = fdict["neighbor"] face.orientation = fdict["orientation"] - try: - face.isPeriodicLow = fdict["isPeriodicLow"] - assert type(face.isPeriodicLow) is bool - except KeyError: - pass - mb.totalBlocks = conn["Total_Blocks"] diff --git a/src/peregrinepy/writers/writeConnectivity.py b/src/peregrinepy/writers/writeConnectivity.py index 71a34ad8..c8104b1b 100644 --- a/src/peregrinepy/writers/writeConnectivity.py +++ b/src/peregrinepy/writers/writeConnectivity.py @@ -33,8 +33,6 @@ def writeConnectivity(mb, path="./"): fdict["bcFam"] = face.bcFam fdict["neighbor"] = face.neighbor fdict["orientation"] = face.orientation - if face.isPeriodicLow is not None: - fdict["isPeriodicLow"] = face.isPeriodicLow with open(f"{path}/conn.yaml", "w") as f: yaml.dump(conn, f, sort_keys=False) From b9bf8895d63ec6c30c7ae3d7a9d4c076187acf11 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Tue, 5 Apr 2022 10:11:56 -0400 Subject: [PATCH 10/31] Added periodic templates --- .../bcs/bcfam_templates/walls/periodicRot.yaml | 7 +++++++ .../bcfam_templates/walls/periodicTrans.yaml | 7 +++++++ src/peregrinepy/readers/readBcs.py | 17 ++++++++++++----- 3 files changed, 26 insertions(+), 5 deletions(-) create mode 100644 src/peregrinepy/bcs/bcfam_templates/walls/periodicRot.yaml create mode 100644 src/peregrinepy/bcs/bcfam_templates/walls/periodicTrans.yaml diff --git a/src/peregrinepy/bcs/bcfam_templates/walls/periodicRot.yaml b/src/peregrinepy/bcs/bcfam_templates/walls/periodicRot.yaml new file mode 100644 index 00000000..0205f46f --- /dev/null +++ b/src/peregrinepy/bcs/bcfam_templates/walls/periodicRot.yaml @@ -0,0 +1,7 @@ +--- + +periodic: + bcType: periodicRot + bcVals: + periodicSpan: 45 # degrees + periodicAxis: [1, 0, 0] diff --git a/src/peregrinepy/bcs/bcfam_templates/walls/periodicTrans.yaml b/src/peregrinepy/bcs/bcfam_templates/walls/periodicTrans.yaml new file mode 100644 index 00000000..a71eff69 --- /dev/null +++ b/src/peregrinepy/bcs/bcfam_templates/walls/periodicTrans.yaml @@ -0,0 +1,7 @@ +--- + +periodic: + bcType: periodicTrans + bcVals: + periodicSpan: 0.1 + periodicAxis: [1, 0, 0] diff --git a/src/peregrinepy/readers/readBcs.py b/src/peregrinepy/readers/readBcs.py index 61c27334..3a19132f 100644 --- a/src/peregrinepy/readers/readBcs.py +++ b/src/peregrinepy/readers/readBcs.py @@ -55,10 +55,17 @@ def readBcs(mb, pathToFile): continue # Make sure the type in the input file matches the type in the connectivity + # Because periodics have a high/low designation they wont match exactly if bcsIn[bcFam]["bcType"] != face.bcType: - raise KeyError( - f'ERROR, block {blk.nblki} face {face.nface} does not match the bcType between input *{bcsIn[bcFam]["bcType"]}* and connectivity *{face.bcType}*.' - ) + if bcsIn[bcFam]["bcType"].startswith("periodic"): + if not face.bcType.startswith(bcsIn[bcFam]["bcType"]): + raise KeyError( + f'ERROR, block {blk.nblki} face {face.nface} does not match the bcType between input *{bcsIn[bcFam]["bcType"]}* and connectivity *{face.bcType}*.' + ) + else: + raise KeyError( + f'ERROR, block {blk.nblki} face {face.nface} does not match the bcType between input *{bcsIn[bcFam]["bcType"]}* and connectivity *{face.bcType}*.' + ) # If there are no values to set, continue if "bcVals" not in bcsIn[bcFam]: @@ -66,9 +73,9 @@ def readBcs(mb, pathToFile): continue # Add the information from any periodic faces - if bcType in ["periodicTrans", "periodicRot"]: - face.periodicAxis = bcsIn[bcFam]["periodicAxis"] + if bcType.startswith("periodic"): face.periodicSpan = bcsIn[bcFam]["periodicSpan"] + face.periodicAxis = bcsIn[bcFam]["periodicAxis"] continue # If we are a solver face, we need to create the kokkos arrays From 3f953e26d5f67fb33d33df0c100059c7f9b73024 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Tue, 5 Apr 2022 14:20:12 -0400 Subject: [PATCH 11/31] GridPro translator and erify grid now using periodic info --- src/peregrinepy/multiBlock/topologyFace.py | 22 ++++- src/peregrinepy/readers/readBcs.py | 4 +- utilities/gridPro2pg.py | 48 +++++++++- utilities/verifyGrid.py | 102 ++++++++++++++------- 4 files changed, 140 insertions(+), 36 deletions(-) diff --git a/src/peregrinepy/multiBlock/topologyFace.py b/src/peregrinepy/multiBlock/topologyFace.py index 91c6c2f4..a5490f8b 100644 --- a/src/peregrinepy/multiBlock/topologyFace.py +++ b/src/peregrinepy/multiBlock/topologyFace.py @@ -112,7 +112,7 @@ def periodicAxis(self, axis): elif self.periodicSpan is None: raise AttributeError("Please set periodicSpan before setting periodicAxis") - # Compute rotation matrix + # Compute rotation matrix for positive and negative rotatoin rot = np.zeros((3, 3)) th = self.periodicSpan * np.pi / 180.0 ct = np.cos(th) @@ -130,7 +130,25 @@ def periodicAxis(self, axis): rot[2, 1] = uz * uy * (1 - ct) + ux * st rot[2, 2] = ct + uz ** 2 * (1 - ct) - self.periodicRotMatrix = rot + self.periodicRotMatrixUp = rot + + rot = np.zeros((3, 3)) + ct = np.cos(-th) + st = np.sin(-th) + ux, uy, uz = tuple(axis) + rot[0, 0] = ct + ux ** 2 * (1 - ct) + rot[0, 1] = ux * uy * (1 - ct) * uz * st + rot[0, 2] = ux * uz * (1 - ct) + uy * st + + rot[1, 0] = uy * ux * (1 - ct) + uz * st + rot[1, 1] = ct + uy ** 2 * (1 - ct) + rot[1, 2] = uy * uz * (1 - ct) - ux * st + + rot[2, 0] = uz * ux * (1 - ct) - uy * st + rot[2, 1] = uz * uy * (1 - ct) + ux * st + rot[2, 2] = ct + uz ** 2 * (1 - ct) + + self.periodicRotMatrixDown = rot @property def neighborNface(self): diff --git a/src/peregrinepy/readers/readBcs.py b/src/peregrinepy/readers/readBcs.py index 3a19132f..b7725241 100644 --- a/src/peregrinepy/readers/readBcs.py +++ b/src/peregrinepy/readers/readBcs.py @@ -74,8 +74,8 @@ def readBcs(mb, pathToFile): # Add the information from any periodic faces if bcType.startswith("periodic"): - face.periodicSpan = bcsIn[bcFam]["periodicSpan"] - face.periodicAxis = bcsIn[bcFam]["periodicAxis"] + face.periodicSpan = bcsIn[bcFam]["bcVals"]["periodicSpan"] + face.periodicAxis = bcsIn[bcFam]["bcVals"]["periodicAxis"] continue # If we are a solver face, we need to create the kokkos arrays diff --git a/utilities/gridPro2pg.py b/utilities/gridPro2pg.py index 1f0d9235..2049791e 100755 --- a/utilities/gridPro2pg.py +++ b/utilities/gridPro2pg.py @@ -102,9 +102,9 @@ # Read in bcFam.yaml file so we know what the bcType is for each label. with open(args.bcFam, "r") as f: bcFamDict = yaml.load(f, Loader=yaml.FullLoader) + gpSurfaceToPgBcType = { "pdc:INTERBLK": "b0", - "pdc:PERIODIC": "b1", "pdc:WALL": "adiabaticNoSlipWall", "pdc:user8": "adiabaticSlipWall", } @@ -182,12 +182,24 @@ pgConns.append(line[2:-1]) mb = mbg(nblks) +gridIsPeriodic = False for temp, blk in zip(pgConns, mb): for face in blk.faces: faceData = temp[(int(face.nface) - 1) * 4 : (int(face.nface) - 1) * 4 + 4] + # Set the named BCs in the bcFam.yaml if faceData[0] in bcFamDict: face.bcFam = f"{faceData[0]}" face.bcType = bcFamDict[faceData[0]]["bcType"] + # Need to treat periodics + elif faceData[0] == "PERIODIC": + gridIsPeriodic = True + # get the periodic from bcFam + key = [key for key in bcFamDict.keys() if key.startswith("periodic")][0] + face.bcFam = key + # We will set High/Low bcType after we have all the point data + face.bcType = bcFamDict[key]["bcType"] + "Low" + face.periodicSpan = bcFamDict[key]["bcVals"]["periodicSpan"] + face.periodicAxis = bcFamDict[key]["bcVals"]["periodicAxis"] else: face.bcType = f"{faceData[0]}" face.bcFam = None @@ -244,6 +256,40 @@ gpPtyFile.close() gpBlkFile.close() +# Now we go through the grid to set periodic High/Low +if gridIsPeriodic: + from verifyGrid import extractFace + + for blk in mb: + for face in blk.faces: + if face.bcType.startswith("periodic"): + myX, myY, myZ = extractFace(blk, face.nface) + myFaceCenter = np.array([np.mean(i) for i in (myX, myY, myZ)]) + nBlk = mb.getBlock(face.neighbor) + nX, nY, nZ = extractFace(nBlk, face.neighborNface) + nFaceCenter = np.array([np.mean(i) for i in (nX, nY, nZ)]) + else: + continue + + if face.bcType == "periodicTransLow": + faceUp = myFaceCenter + face.periodicAxis * face.periodicSpan + faceDown = myFaceCenter - face.periodicAxis * face.periodicSpan + elif face.bcType == "periodicRotLow": + faceUp = np.matmul(face.periodicRotMatrixUp, myFaceCenter) + faceDown = np.matmul(face.periodicRotMatrixDown, myFaceCenter) + else: + raise ValueError("Didnt find either periodic Trans or Rot") + + if np.allclose(faceUp, nFaceCenter): + # Keep the face low + pass + elif np.allclose(faceDown, nFaceCenter): + face.bcType = face.bcType.replace("Low", "High") + else: + raise ValueError( + f"Not a periodic match between block {blk.nblki} face {face.nface} and block {nBlk.nblki} face {face.neighborNface}." + ) + if verify(mb): pass diff --git a/utilities/verifyGrid.py b/utilities/verifyGrid.py index 1694fd1c..91bb1711 100755 --- a/utilities/verifyGrid.py +++ b/utilities/verifyGrid.py @@ -7,9 +7,8 @@ Inputs are the path to the grid files, and path to the conn.inp file. -It can handle b1 periodicity in the conn.yaml, -however it will not compare the x,y,z coordinate locations -of the faces. +If you have periodicity in the grid, you must have the periodic +information populated in the grid, i.e. periodicSpan and periodicAxis Output will print any discrepencies to the screen @@ -22,7 +21,23 @@ import peregrinepy as pg -def verify(mb): +faceToOrientIndexMapping = { + 1: 0, + 2: 0, + 3: 1, + 4: 1, + 5: 2, + 6: 2, +} + +largeIndexMapping = {0: "k", 1: "k", 2: "j"} +needToTranspose = { + "k": {"k": [1, 2, 4, 5], "j": [1, 4]}, + "j": {"k": [1, 2, 4, 5], "j": [1, 4]}, +} + + +def extractFace(blk, nface): faceSliceMapping = { 1: {"i": 0, "j": slice(None), "k": slice(None)}, @@ -33,30 +48,16 @@ def verify(mb): 6: {"i": slice(None), "j": slice(None), "k": -1}, } - faceToOrientIndexMapping = { - 1: 0, - 2: 0, - 3: 1, - 4: 1, - 5: 2, - 6: 2, - } - - largeIndexMapping = {0: "k", 1: "k", 2: "j"} - needToTranspose = { - "k": {"k": [1, 2, 4, 5], "j": [1, 4]}, - "j": {"k": [1, 2, 4, 5], "j": [1, 4]}, - } + face_i = faceSliceMapping[nface] - def extractFace(blk, nface): + x = blk.array["x"][face_i["i"], face_i["j"], face_i["k"]] + y = blk.array["y"][face_i["i"], face_i["j"], face_i["k"]] + z = blk.array["z"][face_i["i"], face_i["j"], face_i["k"]] - face_i = faceSliceMapping[nface] + return x, y, z - x = blk.array["x"][face_i["i"], face_i["j"], face_i["k"]] - y = blk.array["y"][face_i["i"], face_i["j"], face_i["k"]] - z = blk.array["z"][face_i["i"], face_i["j"], face_i["k"]] - return x, y, z +def verify(mb): warn = False for blk in mb: @@ -68,10 +69,6 @@ def extractFace(blk, nface): if neighbor is None: continue - if bc == "b1": - periodic = True - else: - periodic = False (face_x, face_y, face_z) = extractFace(blk, face.nface) @@ -116,6 +113,35 @@ def extractFace(blk, nface): face_y = np.flip(face_y, 1) face_z = np.flip(face_z, 1) + # Here we translate the coordinates to periodic face + if bc == "periodicTransLow": + face_x += face.periodicAxis[0] * face.periodicSpan + face_y += face.periodicAxis[1] * face.periodicSpan + face_z += face.periodicAxis[2] * face.periodicSpan + elif bc == "periodicTransHigh": + face_x -= face.periodicAxis[0] * face.periodicSpan + face_y -= face.periodicAxis[1] * face.periodicSpan + face_z -= face.periodicAxis[2] * face.periodicSpan + elif bc == "periodicRotLow": + shape = face_x.shape + points = np.column_stack( + (face_x.ravel(), face_y.ravel(), face_z.ravel()) + ) + points = np.matmul(face.periodicRotMatrixUp, points.T).T + face_x = points[:, 0].reshape(shape) + face_y = points[:, 1].reshape(shape) + face_z = points[:, 2].reshape(shape) + elif bc == "periodicRotHigh": + shape = face_x.shape + points = np.column_stack( + (face_x.ravel(), face_y.ravel(), face_z.ravel()) + ) + points = np.matmul(face.periodicRotMatrixDown, points.T).T + face_x = points[:, 0].reshape(shape) + face_y = points[:, 1].reshape(shape) + face_z = points[:, 2].reshape(shape) + + # Compare the face points to the neighbor face try: diff_x = ( np.mean( @@ -146,21 +172,21 @@ def extractFace(blk, nface): f"Error when comparing block {blk.nblki} and block {blk2.nblki} connection" ) - if diff_x > 1e-06 and not periodic: + if diff_x > 1e-06: print( f"Warning, the x coordinates of face {nface} on block {blk.nblki} are not matching the x coordinates of face {nface2} of block {blk2.nblki}" ) print(f"Off by average of {diff_x}") warn = True - if diff_y > 1e-06 and not periodic: + if diff_y > 1e-06: print( f"Warning, the y coordinates of face {nface} on block {blk.nblki} are not matching the y coordinates of face {nface2} of block {blk2.nblki}" ) print(f"Off by average of {diff_y}") warn = True - if diff_z > 1e-06 and not periodic: + if diff_z > 1e-06: print( f"Warning, the z coordinates of face {nface} on block {blk.nblki} are not matching the z coordinates of face {nface2} of block {blk2.nblki}" ) @@ -217,17 +243,31 @@ def extractFace(blk, nface): help="Path to conn.yaml", type=str, ) + parser.add_argument( + "-bcFamPath", + action="store", + metavar="", + dest="bcFamPath", + default="./", + help="""If your grid has periodics, we need to the periodic data from bcFams.""", + type=str, + ) args = parser.parse_args() gp = args.gridPath cp = args.connPath + bcFamPath = args.bcFamPath nblks = len([i for i in os.listdir(gp) if i.startswith("g.") and i.endswith(".h5")]) assert nblks > 0 mb = pg.multiBlock.grid(nblks) pg.readers.readGrid(mb, gp) pg.readers.readConnectivity(mb, cp) + try: + pg.readers.readBcs(mb, bcFamPath) + except FileNotFoundError: + print("No bcFam.yaml file provided, assuming no periodics.") if verify(mb): print("Grid is valid!") From c0056cd0940413803fb6e8caf7541ab8c6890972 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Tue, 5 Apr 2022 14:40:32 -0400 Subject: [PATCH 12/31] cutGrid not handles periodics. --- utilities/cutGrid.py | 121 +++++++++++++++---------------------------- 1 file changed, 41 insertions(+), 80 deletions(-) diff --git a/utilities/cutGrid.py b/utilities/cutGrid.py index 5224860a..d23daa2f 100755 --- a/utilities/cutGrid.py +++ b/utilities/cutGrid.py @@ -82,13 +82,35 @@ def findInteriorNeighbor(mb, incompleteBlocks, foundFaces): # Now we go block by block, face by face and look for matching faces for face in blk.faces: - if face.bcType != "b0": + if face.bcType != "b0" and not face.bcType.startswith("periodic"): continue if foundFaces[index][face.nface - 1]: continue meanX = np.mean([i for i in corners["x"][index][face.nface - 1]]) meanY = np.mean([i for i in corners["y"][index][face.nface - 1]]) meanZ = np.mean([i for i in corners["z"][index][face.nface - 1]]) + # Translate periodics + if face.bcType == "periodicTransLow": + meanX += face.periodicAxis[0] * face.periodicSpan + meanY += face.periodicAxis[1] * face.periodicSpan + meanZ += face.periodicAxis[2] * face.periodicSpan + elif face.bcType == "periodicTransHigh": + meanX -= face.periodicAxis[0] * face.periodicSpan + meanY -= face.periodicAxis[1] * face.periodicSpan + meanZ -= face.periodicAxis[2] * face.periodicSpan + elif face.bcType == "periodicRotLow": + mean = np.array([meanX, meanY, meanZ]) + points = np.matmul(face.periodicRotMatrixUp, mean) + meanX = points[0] + meanY = points[1] + meanZ = points[2] + elif face.bcType == "periodicRotHigh": + mean = np.array([meanX, meanY, meanZ]) + points = np.matmul(face.periodicRotMatrixDown, mean) + meanX = points[0] + meanY = points[1] + meanZ = points[2] + for testIndex, nblki in enumerate(incompleteBlocks): testBlk = mb.getBlock(nblki) for testFace in testBlk.faces: @@ -128,84 +150,6 @@ def findInteriorNeighbor(mb, incompleteBlocks, foundFaces): corners[var].pop(index) -def findPeriodicNeighbor(mb, incompleteBlocks, foundFaces): - if len(incompleteBlocks) == []: - return - - corners = extractCorners(mb, incompleteBlocks, foundFaces) - assert all( - [ - len(corners[var]) == len(incompleteBlocks) == len(foundFaces) - for var in corners - ] - ) - - for index, nblki in enumerate(incompleteBlocks): - blk = mb.getBlock(nblki) - - # Is this block complete? - if all(foundFaces[index]): - continue - - # Now we go block by block, face by face and look for matching faces - for face in blk.faces: - if foundFaces[index][face.nface - 1]: - continue - assert face.bcType.startswith("b") and face.bcType != "b0" - Xs = corners["x"][index][face.nface - 1] - Ys = corners["y"][index][face.nface - 1] - Zs = corners["z"][index][face.nface - 1] - faceCenterX = np.mean(Xs) - faceCenterY = np.mean(Ys) - faceCenterZ = np.mean(Zs) - - p0 = np.array([Xs[0], Ys[0], Zs[0]]) - p1 = np.array([Xs[1], Ys[1], Zs[1]]) - p2 = np.array([Xs[3], Ys[3], Zs[3]]) - pFace = np.array([faceCenterX, faceCenterY, faceCenterZ]) - cross = np.cross(p1 - p0, p2 - p0) - cross = cross / np.linalg.norm(cross) - - for testIndex, nblki in enumerate(incompleteBlocks): - testBlk = mb.getBlock(nblki) - for testFace in testBlk.faces: - if foundFaces[testIndex][testFace.nface - 1]: - continue - assert testFace.bcType.startswith("b") and testFace.bcType != "b0" - if blk.nblki == testBlk.nblki and face.nface == testFace.nface: - continue - testXs = corners["x"][testIndex][testFace.nface - 1] - testYs = corners["y"][testIndex][testFace.nface - 1] - testZs = corners["z"][testIndex][testFace.nface - 1] - testFaceCenterX = np.mean(testXs) - testFaceCenterY = np.mean(testYs) - testFaceCenterZ = np.mean(testZs) - - # See if test face is in place with face so we dont accidentially pick it up - pTest = np.array( - [testFaceCenterX, testFaceCenterY, testFaceCenterZ] - ) - testVector = pTest - pFace - testVector = testVector / np.linalg.norm(testVector) - dot = np.abs(np.dot(cross, testVector)) - if dot == 1.0: - face.neighbor = testBlk.nblki - testFace.neighbor = blk.nblki - foundFaces[index][face.nface - 1] = True - foundFaces[testIndex][testFace.nface - 1] = True - - # Filter out all complete blocks - completeIndex = [] - for index, found in enumerate(foundFaces): - if all(found): - completeIndex.append(index) - for index in sorted(completeIndex, reverse=True): - incompleteBlocks.pop(index) - foundFaces.pop(index) - for var in corners: - corners[var].pop(index) - - def cutBlock(mb, nblki, cutAxis, cutIndex, incompleteBlocks, foundFaces): oldBlk = mb.getBlock(nblki) @@ -273,6 +217,10 @@ def cutBlock(mb, nblki, cutAxis, cutIndex, incompleteBlocks, foundFaces): newSplitFace.orientation = oldSplitFace.orientation newSplitFace.bcFam = oldSplitFace.bcFam newSplitFace.bcType = oldSplitFace.bcType + # if the split face is a periodic, they need the perodic info + if oldSplitFace.bcType.startswith("periodic"): + newSplitFace.periodicSpan = oldSplitFace.periodicSpan + newSplitFace.periodicAxis = oldSplitFace.periodicAxis # If this face is a boundary, then we can add it to the found faces if oldSplitFace.neighbor is None: newSplitFace.neighbor = None @@ -360,7 +308,6 @@ def performCutOperations(mb, cutOps): cutBlock(mb, cutNblki, cutAxis, index, incompleteBlocks, foundFaces) findInteriorNeighbor(mb, incompleteBlocks, foundFaces) - findPeriodicNeighbor(mb, incompleteBlocks, foundFaces) assert incompleteBlocks == [] assert foundFaces == [] @@ -417,6 +364,15 @@ def performCutOperations(mb, cutOps): help="If this option is specified, the utility will automatically cut down the grid such that the max block size (ni*nj*nk) < targetBlockSize", type=int, ) + parser.add_argument( + "-bcFamPath", + action="store", + metavar="", + dest="bcFamPath", + default="./", + help="""If your grid has periodics, we need to the periodic data from bcFams.""", + type=str, + ) args = parser.parse_args() @@ -425,6 +381,7 @@ def performCutOperations(mb, cutOps): cutOps = args.cutOps maxBlockSize = args.maxBlockSize + bcFamPath = args.bcFamPath if cutOps and maxBlockSize > 0: raise ValueError("Cannot specify both cutOps and maxBlockSize") @@ -433,6 +390,10 @@ def performCutOperations(mb, cutOps): mb = pg.multiBlock.grid(nblks) pg.readers.readGrid(mb, fromDir) pg.readers.readConnectivity(mb, fromDir) + try: + pg.readers.readBcs(mb, bcFamPath) + except FileNotFoundError: + print("No bcFam.yaml file provided, assuming no periodics.") if cutOps: cutOperations = [] From e3bd99acd0a068d276a44aae78bba70df1d56c94 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Tue, 5 Apr 2022 15:15:42 -0400 Subject: [PATCH 13/31] Updated copyRotateSectorGrid to use periodic scheme --- utilities/copyRotateSectorGrid.py | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/utilities/copyRotateSectorGrid.py b/utilities/copyRotateSectorGrid.py index f6284103..2742265b 100755 --- a/utilities/copyRotateSectorGrid.py +++ b/utilities/copyRotateSectorGrid.py @@ -86,13 +86,10 @@ toGrid[i] = fromBlk for face in fromBlk.faces: - if face.bcType == "b1": - if np.mean(fromBlk.array["z"]) > 0.0: - backside.append(fromBlk.nblki) - elif np.mean(fromBlk.array["z"] < 0.0): - frontside.append(fromBlk.nblki) - else: - print("HELP") + if face.bcType == "periodicRotLow": + backside.append(fromBlk.nblki) + elif face.bcType == "periodicRotHigh": + frontside.append(fromBlk.nblki) # Now copy/rotate sector by sector # We will update the front and back side @@ -131,21 +128,24 @@ continue # treat the rotated faces - if fromFace.bcType == "b1": - toFace.bcFam = None + if fromFace.bcType.startswith("periodic"): + toFace.bcFam = fromFace.bcFam toFace.orientation = fromFace.orientation # back side faces if fromBlk.nblki in backside: toFace.bcType = "b0" toFace.neighbor = fromFace.neighbor + nblks * i + toFace.bcFam = None # front side faces elif fromBlk.nblki in frontside: if i == nseg - 2: - toFace.bcType = "b1" + toFace.bcType = "periodicRotHigh" toFace.neighbor = fromFace.neighbor + toFace.bcFam = fromFace.bcFam else: toFace.bcType = "b0" toFace.neighbor = fromFace.neighbor + nblks * (i + 2) + toFace.bcFam = None # treat the internal faces elif fromFace.bcType == "b0": toFace.bcType = "b0" @@ -165,25 +165,27 @@ if fromBlk.nblki in backside: rotBlk = toGrid[i] for toFace, fromFace in zip(rotBlk.faces, fromBlk.faces): - if fromFace.bcType == "b1": + if fromFace.bcType == "periodicRotLow": toFace.neighbor = fromFace.neighbor + nblks * (nseg - 1) if is360: toFace.bcType = "b0" + toFace.bcFam = None # original and far frontside elif fromBlk.nblki in frontside: # original rotBlk = toGrid[i] for toFace, fromFace in zip(rotBlk.faces, fromBlk.faces): - if fromFace.bcType == "b1": + if fromFace.bcType == "periodicRotHigh": toFace.bcType = "b0" toFace.neighbor = fromFace.neighbor + nblks + toFace.bcFam = None # far frontside rotBlk = toGrid[i + nblks * (nseg - 1)] for toFace, fromFace in zip(rotBlk.faces, fromBlk.faces): if is360: toFace.bcType = "b0" + toFace.bcFam = None pg.writers.writeGrid(toGrid, toDir) pg.writers.writeConnectivity(toGrid, toDir) - # match_faces(toGrid) From a08631e8586883f132ece2e25030033dd485bb34 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Tue, 5 Apr 2022 15:18:11 -0400 Subject: [PATCH 14/31] Reducing warning output from readBcs --- src/peregrinepy/readers/readBcs.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/peregrinepy/readers/readBcs.py b/src/peregrinepy/readers/readBcs.py index b7725241..affccf8c 100644 --- a/src/peregrinepy/readers/readBcs.py +++ b/src/peregrinepy/readers/readBcs.py @@ -37,6 +37,7 @@ def readBcs(mb, pathToFile): print("No bcFams.yaml found, assuming all defaults.") return + warn = [] for blk in mb: for face in blk.faces: bcFam = face.bcFam @@ -69,7 +70,9 @@ def readBcs(mb, pathToFile): # If there are no values to set, continue if "bcVals" not in bcsIn[bcFam]: - print(f"Warning, no values found for {bcsIn[bcFam]}") + if bcsIn[bcFam] not in warn: + print(f"Warning, no values found for {bcsIn[bcFam]}") + warn.append(bcsIn[bcFam]) continue # Add the information from any periodic faces From 1ef629eb8592e10271b5b79e47d11d683518dbca Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Tue, 5 Apr 2022 18:50:58 -0400 Subject: [PATCH 15/31] Updated the grid unification scheme to handle periodics. --- src/peregrinepy/grid/create.py | 2 +- src/peregrinepy/grid/unifySolverGrid.py | 115 ++++++++++++----------- src/peregrinepy/multiBlock/solverFace.py | 2 +- 3 files changed, 60 insertions(+), 59 deletions(-) diff --git a/src/peregrinepy/grid/create.py b/src/peregrinepy/grid/create.py index 4ed718b0..d982831e 100644 --- a/src/peregrinepy/grid/create.py +++ b/src/peregrinepy/grid/create.py @@ -534,7 +534,7 @@ def multiBlockAnnulus( face.bcType = "periodicRotLow" face.periodicSpan = sweep face.periodicAxis = n12 - elif k == mbDims[2] - 1: + if k == mbDims[2] - 1: face = blk.getFace(6) if float(sweep) == 360.0: face.bcType = "b0" diff --git a/src/peregrinepy/grid/unifySolverGrid.py b/src/peregrinepy/grid/unifySolverGrid.py index b0f069b5..14da8533 100644 --- a/src/peregrinepy/grid/unifySolverGrid.py +++ b/src/peregrinepy/grid/unifySolverGrid.py @@ -1,5 +1,3 @@ -from mpi4py.MPI import DOUBLE as MPIDOUBLE -from mpi4py.MPI import Request from .. import mpiComm @@ -13,62 +11,65 @@ def unifySolverGrid(mb): for _ in range(3): mpiComm.communicate(mb, ["x", "y", "z"]) - comm, rank, size = mpiComm.mpiUtils.getCommRankSize() - - for var in ["x", "y", "z"]: - for blk in mb: - # Need to update host data - blk.updateHostView(var) - for _ in range(3): - reqs = [] - # Post non-blocking recieves - for blk in mb: - for face in blk.faces: - bc = face.bcType - if bc != "b1": - continue + for blk in mb: + for face in blk.faces: + bc = face.bcType + if not bc.startswith("periodic"): + continue + for i, sR in enumerate(face.sliceR3): + x = blk.array["x"][sR] + y = blk.array["y"][sR] + z = blk.array["z"][sR] - ssize = face.array["recvBuffer3"].size - reqs.append( - comm.Irecv( - [face.array["recvBuffer3"][:], ssize, MPIDOUBLE], - source=face.commRank, - tag=face.tagR, - ) + # Translate periodics + if face.bcType == "periodicTransLow": + x[:] -= face.periodicAxis[0] * face.periodicSpan + y[:] -= face.periodicAxis[1] * face.periodicSpan + z[:] -= face.periodicAxis[2] * face.periodicSpan + elif face.bcType == "periodicTransHigh": + x[:] += face.periodicAxis[0] * face.periodicSpan + y[:] += face.periodicAxis[1] * face.periodicSpan + z[:] += face.periodicAxis[2] * face.periodicSpan + elif face.bcType == "periodicRotLow": + print(face.nface, face.periodicRotMatrixDown) + tempx = ( + face.periodicRotMatrixDown[0, 0] * x[:] + + face.periodicRotMatrixDown[0, 1] * y[:] + + face.periodicRotMatrixDown[0, 2] * z[:] ) - - # Post non-blocking sends - for blk in mb: - for face in blk.faces: - bc = face.bcType - if bc != "b1": - continue - for i, sS in enumerate(face.sliceS3): - face.array["sendBuffer3"][i] = face.orient( - blk.array[var][sS] - blk.array[var][face.s1_] - ) - ssize = face.array["sendBuffer3"].size - comm.Send( - [face.array["sendBuffer3"], ssize, MPIDOUBLE], - dest=face.commRank, - tag=face.tagS, + tempy = ( + face.periodicRotMatrixDown[1, 0] * x[:] + + face.periodicRotMatrixDown[1, 1] * y[:] + + face.periodicRotMatrixDown[1, 2] * z[:] ) + tempz = ( + face.periodicRotMatrixDown[2, 0] * x[:] + + face.periodicRotMatrixDown[2, 1] * y[:] + + face.periodicRotMatrixDown[2, 2] * z[:] + ) + x[:] = tempx[:] + y[:] = tempy[:] + z[:] = tempz[:] + elif face.bcType == "periodicRotHigh": + tempx = ( + face.periodicRotMatrixUp[0, 0] * x[:] + + face.periodicRotMatrixUp[0, 1] * y[:] + + face.periodicRotMatrixUp[0, 2] * z[:] + ) + tempy = ( + face.periodicRotMatrixUp[1, 0] * x[:] + + face.periodicRotMatrixUp[1, 1] * y[:] + + face.periodicRotMatrixUp[1, 2] * z[:] + ) + tempz = ( + face.periodicRotMatrixUp[2, 0] * x[:] + + face.periodicRotMatrixUp[2, 1] * y[:] + + face.periodicRotMatrixUp[2, 2] * z[:] + ) + x[:] = tempx[:] + y[:] = tempy[:] + z[:] = tempz[:] - # wait and assign - reqs = iter(reqs) - for blk in mb: - for face in blk.faces: - bc = face.bcType - if bc != "b1": - continue - Request.Wait(reqs.__next__()) - for i, sR in enumerate(face.sliceR3): - blk.array[var][sR] = ( - blk.array[var][face.s1_] + face.array["recvBuffer3"][i] - ) - - comm.Barrier() - - for blk in mb: - # Push back up the device - blk.updateDeviceView(var) + for blk in mb: + # Push back up the device + blk.updateDeviceView(["x", "y", "z"]) diff --git a/src/peregrinepy/multiBlock/solverFace.py b/src/peregrinepy/multiBlock/solverFace.py index 277bca77..d95b21bd 100644 --- a/src/peregrinepy/multiBlock/solverFace.py +++ b/src/peregrinepy/multiBlock/solverFace.py @@ -96,7 +96,7 @@ def bcType(self, value): def setBcFunc(self): bcType = self.bcType - if bcType in ["b0", "b1"]: + if bcType == "b0" or bcType.startswith("periodic"): self.bcFunc = null else: for bcmodule in [bcs.inlets, bcs.exits, bcs.walls]: From c21a3b2d54d35ea3b5f80d7716e3d903fdf8bb15 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 09:24:57 -0400 Subject: [PATCH 16/31] Added periodic bc to kokkos compute module --- src/bindings.cpp | 2 + .../bindingsBoundaryConditions.cpp | 17 ++ src/compute/boundaryConditions/periodics.cpp | 151 ++++++++++++++++++ src/compute/compute.hpp | 9 ++ src/compute/face_.hpp | 3 + src/peregrinepy/multiBlock/solverFace.py | 8 +- src/peregrinepy/multiBlock/topologyFace.py | 58 +++---- 7 files changed, 219 insertions(+), 29 deletions(-) create mode 100644 src/compute/boundaryConditions/periodics.cpp diff --git a/src/bindings.cpp b/src/bindings.cpp index 7aa0d784..e5370c6d 100644 --- a/src/bindings.cpp +++ b/src/bindings.cpp @@ -169,6 +169,8 @@ PYBIND11_MODULE(compute, m) { .def_readwrite("recvBuffer4", &face_::recvBuffer4) .def_readwrite("tempRecvBuffer3", &face_::tempRecvBuffer3) .def_readwrite("tempRecvBuffer4", &face_::tempRecvBuffer4) + .def_readwrite("periodicRotMatrixUp", &face_::periodicRotMatrixUp) + .def_readwrite("periodicRotMatrixDown", &face_::periodicRotMatrixDown) .def_readwrite("intervalDt", &face_::intervalDt) .def_readwrite("cubicSplineAlphas", &face_::cubicSplineAlphas) diff --git a/src/compute/boundaryConditions/bindingsBoundaryConditions.cpp b/src/compute/boundaryConditions/bindingsBoundaryConditions.cpp index 2b8bb529..f506cd39 100644 --- a/src/compute/boundaryConditions/bindingsBoundaryConditions.cpp +++ b/src/compute/boundaryConditions/bindingsBoundaryConditions.cpp @@ -101,4 +101,21 @@ void bindBoundaryConditions(py::module_ &m) { py::arg("thtrdat_ object"), py::arg("terms"), py::arg("tme")); + + // |----> periodics.cpp + py::module periodics = bcs.def_submodule("periodics", "periodics boundary conditions module"); + periodics.def("periodicRotHigh", &periodicRotHigh, "High rotational periodic face", + py::arg("block_ object"), + py::arg("face_ object"), + py::arg("eos pointer"), + py::arg("thtrdat_ object"), + py::arg("terms"), + py::arg("tme")); + periodics.def("periodicRotLow", &periodicRotLow, "Low rotational periodic face", + py::arg("block_ object"), + py::arg("face_ object"), + py::arg("eos pointer"), + py::arg("thtrdat_ object"), + py::arg("terms"), + py::arg("tme")); } diff --git a/src/compute/boundaryConditions/periodics.cpp b/src/compute/boundaryConditions/periodics.cpp new file mode 100644 index 00000000..56af2ec5 --- /dev/null +++ b/src/compute/boundaryConditions/periodics.cpp @@ -0,0 +1,151 @@ +#include "Kokkos_Core.hpp" +#include "block_.hpp" +#include "compute.hpp" +#include "face_.hpp" +#include "kokkos_types.hpp" +#include "thtrdat_.hpp" + +void periodicRotHigh(block_ b, + face_ face, + const std::function &eos, + const thtrdat_ th, + const std::string terms, + const double /*tme*/) { + //-------------------------------------------------------------------------------------------| + // Apply BC to face, slice by slice. + //-------------------------------------------------------------------------------------------| + const int ng = b.ng; + int s0, s1, s2, plus; + setHaloSlices(s0, s1, s2, plus, b.ni, b.nj, b.nk, ng, face._nface); + + if (terms.compare("euler") == 0) { + + threeDsubview q1 = getHaloSlice(b.q, face._nface, s1); + MDRange2 range_face = MDRange2({0, 0}, {q1.extent(0), q1.extent(1)}); + for (int g = 0; g < b.ng; g++) { + s0 -= plus * g; + + threeDsubview q0 = getHaloSlice(b.q, face._nface, s0); + + Kokkos::parallel_for( + "Rotate periodic euler terms", range_face, + KOKKOS_LAMBDA(const int i, const int j) { + // rotate velocity vector up to high face + double tempU, tempV, tempW; + double u = q0(i, j, 1); + double v = q0(i, j, 2); + double w = q0(i, j, 3); + tempU = face.periodicRotMatrixUp(0,0)*u + face.periodicRotMatrixUp(0,1)*v + face.periodicRotMatrixUp(0,2)*w; + tempV = face.periodicRotMatrixUp(1,0)*u + face.periodicRotMatrixUp(1,1)*v + face.periodicRotMatrixUp(1,2)*w; + tempW = face.periodicRotMatrixUp(2,0)*u + face.periodicRotMatrixUp(2,1)*v + face.periodicRotMatrixUp(2,2)*w; + + q0(i, j, 1) = tempU; + q0(i, j, 2) = tempV; + q0(i, j, 3) = tempW; + }); + } + + } else if (terms.compare("viscous") == 0) { + + threeDsubview dqdx1 = getHaloSlice(b.dqdx, face._nface, s1); + MDRange3 range_face = + MDRange3({0, 0, 0}, {static_cast(dqdx1.extent(0)), + static_cast(dqdx1.extent(1)), b.ne}); + for (int g = 0; g < b.ng; g++) { + s0 -= plus * g; + threeDsubview dqdx0 = getHaloSlice(b.dqdx, face._nface, s0); + threeDsubview dqdy0 = getHaloSlice(b.dqdy, face._nface, s0); + threeDsubview dqdz0 = getHaloSlice(b.dqdz, face._nface, s0); + + Kokkos::parallel_for( + "Constant pressure subsonic exit euler terms", range_face, + KOKKOS_LAMBDA(const int i, const int j, const int l) { + // rotate gradients vector up to high face + double tempdx, tempdy, tempdz; + double dx = dqdz0(i, j, 1); + double dy = dqdy0(i, j, 2); + double dz = dqdz0(i, j, 3); + tempdx = face.periodicRotMatrixUp(0,0)*dx + face.periodicRotMatrixUp(0,1)*dx + face.periodicRotMatrixUp(0,2)*dx; + tempdy = face.periodicRotMatrixUp(1,0)*dy + face.periodicRotMatrixUp(1,1)*dy + face.periodicRotMatrixUp(1,2)*dy; + tempdz = face.periodicRotMatrixUp(2,0)*dz + face.periodicRotMatrixUp(2,1)*dz + face.periodicRotMatrixUp(2,2)*dz; + + dqdx0(i, j, l) = tempdx; + dqdy0(i, j, l) = tempdy; + dqdz0(i, j, l) = tempdz; + }); + } + } +} + + +void periodicRotLow(block_ b, + face_ face, + const std::function &eos, + const thtrdat_ th, + const std::string terms, + const double /*tme*/) { + //-------------------------------------------------------------------------------------------| + // Apply BC to face, slice by slice. + //-------------------------------------------------------------------------------------------| + const int ng = b.ng; + int s0, s1, s2, plus; + setHaloSlices(s0, s1, s2, plus, b.ni, b.nj, b.nk, ng, face._nface); + + if (terms.compare("euler") == 0) { + + threeDsubview q1 = getHaloSlice(b.q, face._nface, s1); + MDRange2 range_face = MDRange2({0, 0}, {q1.extent(0), q1.extent(1)}); + for (int g = 0; g < b.ng; g++) { + s0 -= plus * g; + + threeDsubview q0 = getHaloSlice(b.q, face._nface, s0); + + Kokkos::parallel_for( + "Rotate periodic euler terms", range_face, + KOKKOS_LAMBDA(const int i, const int j) { + // rotate velocity vector up to high face + double tempU, tempV, tempW; + double u = q0(i, j, 1); + double v = q0(i, j, 2); + double w = q0(i, j, 3); + tempU = face.periodicRotMatrixDown(0,0)*u + face.periodicRotMatrixDown(0,1)*v + face.periodicRotMatrixDown(0,2)*w; + tempV = face.periodicRotMatrixDown(1,0)*u + face.periodicRotMatrixDown(1,1)*v + face.periodicRotMatrixDown(1,2)*w; + tempW = face.periodicRotMatrixDown(2,0)*u + face.periodicRotMatrixDown(2,1)*v + face.periodicRotMatrixDown(2,2)*w; + + q0(i, j, 1) = tempU; + q0(i, j, 2) = tempV; + q0(i, j, 3) = tempW; + }); + } + + } else if (terms.compare("viscous") == 0) { + + threeDsubview dqdx1 = getHaloSlice(b.dqdx, face._nface, s1); + MDRange3 range_face = + MDRange3({0, 0, 0}, {static_cast(dqdx1.extent(0)), + static_cast(dqdx1.extent(1)), b.ne}); + for (int g = 0; g < b.ng; g++) { + s0 -= plus * g; + threeDsubview dqdx0 = getHaloSlice(b.dqdx, face._nface, s0); + threeDsubview dqdy0 = getHaloSlice(b.dqdy, face._nface, s0); + threeDsubview dqdz0 = getHaloSlice(b.dqdz, face._nface, s0); + + Kokkos::parallel_for( + "Constant pressure subsonic exit euler terms", range_face, + KOKKOS_LAMBDA(const int i, const int j, const int l) { + // rotate gradients vector up to high face + double tempdx, tempdy, tempdz; + double dx = dqdz0(i, j, 1); + double dy = dqdy0(i, j, 2); + double dz = dqdz0(i, j, 3); + tempdx = face.periodicRotMatrixDown(0,0)*dx + face.periodicRotMatrixDown(0,1)*dx + face.periodicRotMatrixDown(0,2)*dx; + tempdy = face.periodicRotMatrixDown(1,0)*dy + face.periodicRotMatrixDown(1,1)*dy + face.periodicRotMatrixDown(1,2)*dy; + tempdz = face.periodicRotMatrixDown(2,0)*dz + face.periodicRotMatrixDown(2,1)*dz + face.periodicRotMatrixDown(2,2)*dz; + + dqdx0(i, j, l) = tempdx; + dqdy0(i, j, l) = tempdy; + dqdz0(i, j, l) = tempdz; + }); + } + } +} diff --git a/src/compute/compute.hpp b/src/compute/compute.hpp index 2fef5023..9e5ff28e 100644 --- a/src/compute/compute.hpp +++ b/src/compute/compute.hpp @@ -89,6 +89,15 @@ void supersonicExit( block_ b, face_ face, const std::function &eos, const thtrdat_ th, const std::string terms, const double tme); +// |------> periodics +void periodicRotHigh( + block_ b, face_ face, + const std::function &eos, + const thtrdat_ th, const std::string terms, const double tme); +void periodicRotLow( + block_ b, face_ face, + const std::function &eos, + const thtrdat_ th, const std::string terms, const double tme); // ./chemistry // |------> CH4_O2_Stanford_Skeletal diff --git a/src/compute/face_.hpp b/src/compute/face_.hpp index b42d43b4..2cf1ca61 100644 --- a/src/compute/face_.hpp +++ b/src/compute/face_.hpp @@ -25,6 +25,9 @@ struct face_ { double intervalDt; int currentInterval = -1; + // Periodic rotation matricies + twoDview periodicRotMatrixUp, periodicRotMatrixDown; + }; #endif diff --git a/src/peregrinepy/multiBlock/solverFace.py b/src/peregrinepy/multiBlock/solverFace.py index d95b21bd..8327a6e1 100644 --- a/src/peregrinepy/multiBlock/solverFace.py +++ b/src/peregrinepy/multiBlock/solverFace.py @@ -59,6 +59,8 @@ def __init__(self, nface, ng): "tempRecvBuffer4": None, "qBcVals": None, "QBcVals": None, + "periodicRotMatrixUp": None, + "periodicRotMatrixDown": None, } ) self.mirror = frozenDict( @@ -71,6 +73,8 @@ def __init__(self, nface, ng): "tempRecvBuffer4": None, "qBcVals": None, "QBcVals": None, + "periodicRotMatrixUp": None, + "periodicRotMatrixDown": None, } ) # Boundary function @@ -96,10 +100,10 @@ def bcType(self, value): def setBcFunc(self): bcType = self.bcType - if bcType == "b0" or bcType.startswith("periodic"): + if bcType == "b0" or bcType.startswith("peiodicTrans"): self.bcFunc = null else: - for bcmodule in [bcs.inlets, bcs.exits, bcs.walls]: + for bcmodule in [bcs.inlets, bcs.exits, bcs.walls, bcs.periodics]: try: self.bcFunc = getattr(bcmodule, bcType) break diff --git a/src/peregrinepy/multiBlock/topologyFace.py b/src/peregrinepy/multiBlock/topologyFace.py index a5490f8b..72cecd9d 100644 --- a/src/peregrinepy/multiBlock/topologyFace.py +++ b/src/peregrinepy/multiBlock/topologyFace.py @@ -11,7 +11,6 @@ def __init__(self, nface): self._neighbor = None self._orientation = None - self.periodicOrigin = [0, 0, 0] self.periodicSpan = None self._periodicAxis = None @@ -113,42 +112,47 @@ def periodicAxis(self, axis): raise AttributeError("Please set periodicSpan before setting periodicAxis") # Compute rotation matrix for positive and negative rotatoin - rot = np.zeros((3, 3)) + rotUp = np.zeros((3, 3)) th = self.periodicSpan * np.pi / 180.0 ct = np.cos(th) st = np.sin(th) ux, uy, uz = tuple(axis) - rot[0, 0] = ct + ux ** 2 * (1 - ct) - rot[0, 1] = ux * uy * (1 - ct) * uz * st - rot[0, 2] = ux * uz * (1 - ct) + uy * st + rotUp[0, 0] = ct + ux ** 2 * (1 - ct) + rotUp[0, 1] = ux * uy * (1 - ct) * uz * st + rotUp[0, 2] = ux * uz * (1 - ct) + uy * st - rot[1, 0] = uy * ux * (1 - ct) + uz * st - rot[1, 1] = ct + uy ** 2 * (1 - ct) - rot[1, 2] = uy * uz * (1 - ct) - ux * st + rotUp[1, 0] = uy * ux * (1 - ct) + uz * st + rotUp[1, 1] = ct + uy ** 2 * (1 - ct) + rotUp[1, 2] = uy * uz * (1 - ct) - ux * st - rot[2, 0] = uz * ux * (1 - ct) - uy * st - rot[2, 1] = uz * uy * (1 - ct) + ux * st - rot[2, 2] = ct + uz ** 2 * (1 - ct) + rotUp[2, 0] = uz * ux * (1 - ct) - uy * st + rotUp[2, 1] = uz * uy * (1 - ct) + ux * st + rotUp[2, 2] = ct + uz ** 2 * (1 - ct) - self.periodicRotMatrixUp = rot - - rot = np.zeros((3, 3)) + rotDown = np.zeros((3, 3)) ct = np.cos(-th) st = np.sin(-th) ux, uy, uz = tuple(axis) - rot[0, 0] = ct + ux ** 2 * (1 - ct) - rot[0, 1] = ux * uy * (1 - ct) * uz * st - rot[0, 2] = ux * uz * (1 - ct) + uy * st - - rot[1, 0] = uy * ux * (1 - ct) + uz * st - rot[1, 1] = ct + uy ** 2 * (1 - ct) - rot[1, 2] = uy * uz * (1 - ct) - ux * st - - rot[2, 0] = uz * ux * (1 - ct) - uy * st - rot[2, 1] = uz * uy * (1 - ct) + ux * st - rot[2, 2] = ct + uz ** 2 * (1 - ct) - - self.periodicRotMatrixDown = rot + rotDown[0, 0] = ct + ux ** 2 * (1 - ct) + rotDown[0, 1] = ux * uy * (1 - ct) * uz * st + rotDown[0, 2] = ux * uz * (1 - ct) + uy * st + + rotDown[1, 0] = uy * ux * (1 - ct) + uz * st + rotDown[1, 1] = ct + uy ** 2 * (1 - ct) + rotDown[1, 2] = uy * uz * (1 - ct) - ux * st + + rotDown[2, 0] = uz * ux * (1 - ct) - uy * st + rotDown[2, 1] = uz * uy * (1 - ct) + ux * st + rotDown[2, 2] = ct + uz ** 2 * (1 - ct) + + if self.faceType == "topology": + self.periodicRotMatrixUp = rotUp + self.periodicRotMatrixDown = rotDown + elif self.faceType == "solver": + from ..misc import createViewMirrorArray + + createViewMirrorArray(self, "periodicRotMatrixUp", (3, 3)) + createViewMirrorArray(self, "periodicRotMatrixDown", (3, 3)) @property def neighborNface(self): From b5080fe23468b7c058284bd7bee718669d737b45 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 10:11:49 -0400 Subject: [PATCH 17/31] Fixed bug in periodic bc --- src/compute/boundaryConditions/periodics.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/compute/boundaryConditions/periodics.cpp b/src/compute/boundaryConditions/periodics.cpp index 56af2ec5..3156a9ce 100644 --- a/src/compute/boundaryConditions/periodics.cpp +++ b/src/compute/boundaryConditions/periodics.cpp @@ -62,9 +62,9 @@ void periodicRotHigh(block_ b, KOKKOS_LAMBDA(const int i, const int j, const int l) { // rotate gradients vector up to high face double tempdx, tempdy, tempdz; - double dx = dqdz0(i, j, 1); - double dy = dqdy0(i, j, 2); - double dz = dqdz0(i, j, 3); + double dx = dqdx0(i, j, l); + double dy = dqdy0(i, j, l); + double dz = dqdz0(i, j, l); tempdx = face.periodicRotMatrixUp(0,0)*dx + face.periodicRotMatrixUp(0,1)*dx + face.periodicRotMatrixUp(0,2)*dx; tempdy = face.periodicRotMatrixUp(1,0)*dy + face.periodicRotMatrixUp(1,1)*dy + face.periodicRotMatrixUp(1,2)*dy; tempdz = face.periodicRotMatrixUp(2,0)*dz + face.periodicRotMatrixUp(2,1)*dz + face.periodicRotMatrixUp(2,2)*dz; @@ -135,9 +135,9 @@ void periodicRotLow(block_ b, KOKKOS_LAMBDA(const int i, const int j, const int l) { // rotate gradients vector up to high face double tempdx, tempdy, tempdz; - double dx = dqdz0(i, j, 1); - double dy = dqdy0(i, j, 2); - double dz = dqdz0(i, j, 3); + double dx = dqdx0(i, j, l); + double dy = dqdy0(i, j, l); + double dz = dqdz0(i, j, l); tempdx = face.periodicRotMatrixDown(0,0)*dx + face.periodicRotMatrixDown(0,1)*dx + face.periodicRotMatrixDown(0,2)*dx; tempdy = face.periodicRotMatrixDown(1,0)*dy + face.periodicRotMatrixDown(1,1)*dy + face.periodicRotMatrixDown(1,2)*dy; tempdz = face.periodicRotMatrixDown(2,0)*dz + face.periodicRotMatrixDown(2,1)*dz + face.periodicRotMatrixDown(2,2)*dz; From 60c1681fa9dcb05ba07cbd5f76fe8ca2659162f2 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 10:12:54 -0400 Subject: [PATCH 18/31] Fixed bug in periodic --- src/compute/boundaryConditions/periodics.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/compute/boundaryConditions/periodics.cpp b/src/compute/boundaryConditions/periodics.cpp index 3156a9ce..2980113c 100644 --- a/src/compute/boundaryConditions/periodics.cpp +++ b/src/compute/boundaryConditions/periodics.cpp @@ -65,9 +65,9 @@ void periodicRotHigh(block_ b, double dx = dqdx0(i, j, l); double dy = dqdy0(i, j, l); double dz = dqdz0(i, j, l); - tempdx = face.periodicRotMatrixUp(0,0)*dx + face.periodicRotMatrixUp(0,1)*dx + face.periodicRotMatrixUp(0,2)*dx; - tempdy = face.periodicRotMatrixUp(1,0)*dy + face.periodicRotMatrixUp(1,1)*dy + face.periodicRotMatrixUp(1,2)*dy; - tempdz = face.periodicRotMatrixUp(2,0)*dz + face.periodicRotMatrixUp(2,1)*dz + face.periodicRotMatrixUp(2,2)*dz; + tempdx = face.periodicRotMatrixUp(0,0)*dx + face.periodicRotMatrixUp(0,1)*dy + face.periodicRotMatrixUp(0,2)*dz; + tempdy = face.periodicRotMatrixUp(1,0)*dx + face.periodicRotMatrixUp(1,1)*dy + face.periodicRotMatrixUp(1,2)*dz; + tempdz = face.periodicRotMatrixUp(2,0)*dx + face.periodicRotMatrixUp(2,1)*dy + face.periodicRotMatrixUp(2,2)*dz; dqdx0(i, j, l) = tempdx; dqdy0(i, j, l) = tempdy; @@ -138,9 +138,9 @@ void periodicRotLow(block_ b, double dx = dqdx0(i, j, l); double dy = dqdy0(i, j, l); double dz = dqdz0(i, j, l); - tempdx = face.periodicRotMatrixDown(0,0)*dx + face.periodicRotMatrixDown(0,1)*dx + face.periodicRotMatrixDown(0,2)*dx; - tempdy = face.periodicRotMatrixDown(1,0)*dy + face.periodicRotMatrixDown(1,1)*dy + face.periodicRotMatrixDown(1,2)*dy; - tempdz = face.periodicRotMatrixDown(2,0)*dz + face.periodicRotMatrixDown(2,1)*dz + face.periodicRotMatrixDown(2,2)*dz; + tempdx = face.periodicRotMatrixDown(0,0)*dx + face.periodicRotMatrixDown(0,1)*dy + face.periodicRotMatrixDown(0,2)*dz; + tempdy = face.periodicRotMatrixDown(1,0)*dx + face.periodicRotMatrixDown(1,1)*dy + face.periodicRotMatrixDown(1,2)*dz; + tempdz = face.periodicRotMatrixDown(2,0)*dx + face.periodicRotMatrixDown(2,1)*dy + face.periodicRotMatrixDown(2,2)*dz; dqdx0(i, j, l) = tempdx; dqdy0(i, j, l) = tempdy; From 0cc652dde48f6a7b00858b550f5c9d2075e67ee7 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 10:14:31 -0400 Subject: [PATCH 19/31] Fixed labels --- src/compute/boundaryConditions/exits.cpp | 2 +- src/compute/boundaryConditions/periodics.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/compute/boundaryConditions/exits.cpp b/src/compute/boundaryConditions/exits.cpp index f3172f48..98004e58 100644 --- a/src/compute/boundaryConditions/exits.cpp +++ b/src/compute/boundaryConditions/exits.cpp @@ -91,7 +91,7 @@ void constantPressureSubsonicExit(block_ b, threeDsubview dqdz0 = getHaloSlice(b.dqdz, face._nface, s0); Kokkos::parallel_for( - "Constant pressure subsonic exit euler terms", range_face, + "Constant pressure subsonic exit viscous terms", range_face, KOKKOS_LAMBDA(const int i, const int j, const int l) { // neumann all gradients dqdx0(i, j, l) = dqdx1(i, j, l); diff --git a/src/compute/boundaryConditions/periodics.cpp b/src/compute/boundaryConditions/periodics.cpp index 2980113c..83f5a64d 100644 --- a/src/compute/boundaryConditions/periodics.cpp +++ b/src/compute/boundaryConditions/periodics.cpp @@ -58,7 +58,7 @@ void periodicRotHigh(block_ b, threeDsubview dqdz0 = getHaloSlice(b.dqdz, face._nface, s0); Kokkos::parallel_for( - "Constant pressure subsonic exit euler terms", range_face, + "Periodic viscous terms", range_face, KOKKOS_LAMBDA(const int i, const int j, const int l) { // rotate gradients vector up to high face double tempdx, tempdy, tempdz; @@ -131,7 +131,7 @@ void periodicRotLow(block_ b, threeDsubview dqdz0 = getHaloSlice(b.dqdz, face._nface, s0); Kokkos::parallel_for( - "Constant pressure subsonic exit euler terms", range_face, + "Periodic viscous terms", range_face, KOKKOS_LAMBDA(const int i, const int j, const int l) { // rotate gradients vector up to high face double tempdx, tempdy, tempdz; From f51e00ea3c12453598c6488c374c16a15af418a2 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 16:56:44 -0400 Subject: [PATCH 20/31] Fixed spec of roationMatricies to face.array[""] --- src/peregrinepy/grid/unifySolverGrid.py | 44 ++++------------------ src/peregrinepy/multiBlock/topologyFace.py | 3 ++ 2 files changed, 11 insertions(+), 36 deletions(-) diff --git a/src/peregrinepy/grid/unifySolverGrid.py b/src/peregrinepy/grid/unifySolverGrid.py index 14da8533..989ac7b3 100644 --- a/src/peregrinepy/grid/unifySolverGrid.py +++ b/src/peregrinepy/grid/unifySolverGrid.py @@ -30,42 +30,14 @@ def unifySolverGrid(mb): x[:] += face.periodicAxis[0] * face.periodicSpan y[:] += face.periodicAxis[1] * face.periodicSpan z[:] += face.periodicAxis[2] * face.periodicSpan - elif face.bcType == "periodicRotLow": - print(face.nface, face.periodicRotMatrixDown) - tempx = ( - face.periodicRotMatrixDown[0, 0] * x[:] - + face.periodicRotMatrixDown[0, 1] * y[:] - + face.periodicRotMatrixDown[0, 2] * z[:] - ) - tempy = ( - face.periodicRotMatrixDown[1, 0] * x[:] - + face.periodicRotMatrixDown[1, 1] * y[:] - + face.periodicRotMatrixDown[1, 2] * z[:] - ) - tempz = ( - face.periodicRotMatrixDown[2, 0] * x[:] - + face.periodicRotMatrixDown[2, 1] * y[:] - + face.periodicRotMatrixDown[2, 2] * z[:] - ) - x[:] = tempx[:] - y[:] = tempy[:] - z[:] = tempz[:] - elif face.bcType == "periodicRotHigh": - tempx = ( - face.periodicRotMatrixUp[0, 0] * x[:] - + face.periodicRotMatrixUp[0, 1] * y[:] - + face.periodicRotMatrixUp[0, 2] * z[:] - ) - tempy = ( - face.periodicRotMatrixUp[1, 0] * x[:] - + face.periodicRotMatrixUp[1, 1] * y[:] - + face.periodicRotMatrixUp[1, 2] * z[:] - ) - tempz = ( - face.periodicRotMatrixUp[2, 0] * x[:] - + face.periodicRotMatrixUp[2, 1] * y[:] - + face.periodicRotMatrixUp[2, 2] * z[:] - ) + elif face.bcType.startswith("periodicRot"): + if face.bcType == "periodicRotLow": + rotM = face.array["periodicRotMatrixDown"] + elif face.bcType == "periodicRotHigh": + rotM = face.array["periodicRotMatrixUp"] + tempx = rotM[0, 0] * x[:] + rotM[0, 1] * y[:] + rotM[0, 2] * z[:] + tempy = rotM[1, 0] * x[:] + rotM[1, 1] * y[:] + rotM[1, 2] * z[:] + tempz = rotM[2, 0] * x[:] + rotM[2, 1] * y[:] + rotM[2, 2] * z[:] x[:] = tempx[:] y[:] = tempy[:] z[:] = tempz[:] diff --git a/src/peregrinepy/multiBlock/topologyFace.py b/src/peregrinepy/multiBlock/topologyFace.py index 72cecd9d..1980ca56 100644 --- a/src/peregrinepy/multiBlock/topologyFace.py +++ b/src/peregrinepy/multiBlock/topologyFace.py @@ -151,6 +151,9 @@ def periodicAxis(self, axis): elif self.faceType == "solver": from ..misc import createViewMirrorArray + self.array["periodicRotMatrixUp"] = rotUp + self.array["periodicRotMatrixDown"] = rotDown + createViewMirrorArray(self, "periodicRotMatrixUp", (3, 3)) createViewMirrorArray(self, "periodicRotMatrixDown", (3, 3)) From 201439ec42b65263bedd2f2f535796b6155f254b Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 18:25:52 -0400 Subject: [PATCH 21/31] Fixed typo --- src/peregrinepy/multiBlock/solverFace.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/peregrinepy/multiBlock/solverFace.py b/src/peregrinepy/multiBlock/solverFace.py index 8327a6e1..85653e12 100644 --- a/src/peregrinepy/multiBlock/solverFace.py +++ b/src/peregrinepy/multiBlock/solverFace.py @@ -100,7 +100,7 @@ def bcType(self, value): def setBcFunc(self): bcType = self.bcType - if bcType == "b0" or bcType.startswith("peiodicTrans"): + if bcType == "b0" or bcType.startswith("periodicTrans"): self.bcFunc = null else: for bcmodule in [bcs.inlets, bcs.exits, bcs.walls, bcs.periodics]: From 827353f27a9dd121da220706f1371e34a888f792 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 18:26:00 -0400 Subject: [PATCH 22/31] Periodic rotation axis must pass through origin --- src/peregrinepy/grid/create.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/peregrinepy/grid/create.py b/src/peregrinepy/grid/create.py index d982831e..0214dca8 100644 --- a/src/peregrinepy/grid/create.py +++ b/src/peregrinepy/grid/create.py @@ -314,7 +314,7 @@ def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): """ - p1 = np.array(p1) + p1 = np.array([0.0, 0.0, 0.0]) # All periodic axes go through origin!!! p2 = np.array(p2) p3 = np.array(p3) @@ -392,7 +392,6 @@ def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): def multiBlockAnnulus( mb, - p1=[0, 0, 0], p2=[1, 0, 0], p3=[0, 1, 0], sweep=45, @@ -413,10 +412,6 @@ def multiBlockAnnulus( mb : peregrinepy.multiBlock.grid (or one of its descendants) - p1 : list, tuple - List/tuple of length 3 containing the location of the origin of the annulus to be created, i.e. - the center of the beginning of the whole annulus. - p2 : list, tuple List/tuple of length 3 containing the location of the end of the annulus to be created, i.e. the axial center of the end of the the annulus. @@ -455,7 +450,7 @@ def multiBlockAnnulus( Updates elements in mb """ - p1 = np.array(p1) + p1 = np.array([0.0, 0.0, 0.0]) # All periodic axes go through origin!!! p2 = np.array(p2) p3 = np.array(p3) From 35e2d56dbb0a4fa763c344366de2c43fb81bd9d9 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 18:37:14 -0400 Subject: [PATCH 23/31] Fixed bug in cubic connectivity --- src/peregrinepy/grid/create.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/peregrinepy/grid/create.py b/src/peregrinepy/grid/create.py index 0214dca8..db00a214 100644 --- a/src/peregrinepy/grid/create.py +++ b/src/peregrinepy/grid/create.py @@ -106,9 +106,9 @@ def cubicConnectivity( face.neighbor = blkNum - mbDims[0] * mbDims[1] * (mbDims[2] - 1) face.orientation = "123" else: - face.bcType = "b0" - face.neighbor = blkNum + mbDims[0] * mbDims[1] - face.orientation = "123" + face.bcType = "adiabaticNoSlipWall" + face.neighbor = None + face.orientation = None else: face.bcType = "b0" face.neighbor = blkNum + mbDims[0] * mbDims[1] From 961f66af683edc7f9e40be7c4950c5fa56141742 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 19:35:47 -0400 Subject: [PATCH 24/31] Also update momentum --- src/compute/boundaryConditions/periodics.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/compute/boundaryConditions/periodics.cpp b/src/compute/boundaryConditions/periodics.cpp index 83f5a64d..b9478185 100644 --- a/src/compute/boundaryConditions/periodics.cpp +++ b/src/compute/boundaryConditions/periodics.cpp @@ -26,6 +26,7 @@ void periodicRotHigh(block_ b, s0 -= plus * g; threeDsubview q0 = getHaloSlice(b.q, face._nface, s0); + threeDsubview Q0 = getHaloSlice(b.Q, face._nface, s0); Kokkos::parallel_for( "Rotate periodic euler terms", range_face, @@ -39,9 +40,15 @@ void periodicRotHigh(block_ b, tempV = face.periodicRotMatrixUp(1,0)*u + face.periodicRotMatrixUp(1,1)*v + face.periodicRotMatrixUp(1,2)*w; tempW = face.periodicRotMatrixUp(2,0)*u + face.periodicRotMatrixUp(2,1)*v + face.periodicRotMatrixUp(2,2)*w; + // Update velocity q0(i, j, 1) = tempU; q0(i, j, 2) = tempV; q0(i, j, 3) = tempW; + + // Update momentum + Q0(i, j, 1) = tempU*Q0(i,j,0); + Q0(i, j, 2) = tempV*Q0(i,j,0); + Q0(i, j, 3) = tempW*Q0(i,j,0); }); } @@ -99,6 +106,7 @@ void periodicRotLow(block_ b, s0 -= plus * g; threeDsubview q0 = getHaloSlice(b.q, face._nface, s0); + threeDsubview Q0 = getHaloSlice(b.Q, face._nface, s0); Kokkos::parallel_for( "Rotate periodic euler terms", range_face, @@ -112,9 +120,15 @@ void periodicRotLow(block_ b, tempV = face.periodicRotMatrixDown(1,0)*u + face.periodicRotMatrixDown(1,1)*v + face.periodicRotMatrixDown(1,2)*w; tempW = face.periodicRotMatrixDown(2,0)*u + face.periodicRotMatrixDown(2,1)*v + face.periodicRotMatrixDown(2,2)*w; + // Update velocity q0(i, j, 1) = tempU; q0(i, j, 2) = tempV; q0(i, j, 3) = tempW; + + // Update momentum + Q0(i, j, 1) = tempU*Q0(i,j,0); + Q0(i, j, 2) = tempV*Q0(i,j,0); + Q0(i, j, 3) = tempW*Q0(i,j,0); }); } From c090e3df9516ce9526496dbfcc98e1b405f8acba Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 21:21:14 -0400 Subject: [PATCH 25/31] Fixed bug in rot matrix definition --- src/peregrinepy/multiBlock/topologyFace.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/peregrinepy/multiBlock/topologyFace.py b/src/peregrinepy/multiBlock/topologyFace.py index 1980ca56..6f24d9ca 100644 --- a/src/peregrinepy/multiBlock/topologyFace.py +++ b/src/peregrinepy/multiBlock/topologyFace.py @@ -103,13 +103,14 @@ def periodicAxis(self): def periodicAxis(self, axis): import numpy as np + axis = np.array(axis) axis = axis / np.linalg.norm(np.array(axis)) self._periodicAxis = axis if not self.bcType.startswith("periodicRot"): return elif self.periodicSpan is None: - raise AttributeError("Please set periodicSpan before setting periodicAxis") + raise AttributeError("Must set periodicSpan before setting periodicAxis") # Compute rotation matrix for positive and negative rotatoin rotUp = np.zeros((3, 3)) @@ -118,7 +119,7 @@ def periodicAxis(self, axis): st = np.sin(th) ux, uy, uz = tuple(axis) rotUp[0, 0] = ct + ux ** 2 * (1 - ct) - rotUp[0, 1] = ux * uy * (1 - ct) * uz * st + rotUp[0, 1] = ux * uy * (1 - ct) - uz * st rotUp[0, 2] = ux * uz * (1 - ct) + uy * st rotUp[1, 0] = uy * ux * (1 - ct) + uz * st @@ -132,9 +133,8 @@ def periodicAxis(self, axis): rotDown = np.zeros((3, 3)) ct = np.cos(-th) st = np.sin(-th) - ux, uy, uz = tuple(axis) rotDown[0, 0] = ct + ux ** 2 * (1 - ct) - rotDown[0, 1] = ux * uy * (1 - ct) * uz * st + rotDown[0, 1] = ux * uy * (1 - ct) - uz * st rotDown[0, 2] = ux * uz * (1 - ct) + uy * st rotDown[1, 0] = uy * ux * (1 - ct) + uz * st From 1ef6300e6339769ec9a29e666717b61d895752bc Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 21:21:32 -0400 Subject: [PATCH 26/31] Add tolerance to annulus axis creation --- src/peregrinepy/grid/create.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/peregrinepy/grid/create.py b/src/peregrinepy/grid/create.py index db00a214..36459618 100644 --- a/src/peregrinepy/grid/create.py +++ b/src/peregrinepy/grid/create.py @@ -318,7 +318,7 @@ def annulus(blk, p1, p2, p3, sweep, thickness, dimensions): p2 = np.array(p2) p3 = np.array(p3) - if np.dot(p2 - p1, p3 - p1) != 0.0: + if np.abs(np.dot(p2 - p1, p3 - p1)) > 1e-7: raise ValueError("Error: The line (p1,p2) is not orthogonal to (p1,p3)") if abs(sweep) < -360 or abs(sweep) > 360.0: @@ -459,7 +459,7 @@ def multiBlockAnnulus( "Error: multiBlock dimensions does not equal number of blocks!" ) - if np.dot(p2 - p1, p3 - p1) != 0.0: + if np.abs(np.dot(p2 - p1, p3 - p1)) > 1e-7: raise ValueError("Error: The line (p1,p2) is not orthogonal to (p1,p3)") if sweep < 0.0 or sweep > 360.0: From 77b3e0af06495d2c5644fa741d7823d856dbe559 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 21:22:11 -0400 Subject: [PATCH 27/31] Should iterate over halo slices not recieve slices --- src/peregrinepy/grid/unifySolverGrid.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/peregrinepy/grid/unifySolverGrid.py b/src/peregrinepy/grid/unifySolverGrid.py index 989ac7b3..be199903 100644 --- a/src/peregrinepy/grid/unifySolverGrid.py +++ b/src/peregrinepy/grid/unifySolverGrid.py @@ -16,10 +16,10 @@ def unifySolverGrid(mb): bc = face.bcType if not bc.startswith("periodic"): continue - for i, sR in enumerate(face.sliceR3): - x = blk.array["x"][sR] - y = blk.array["y"][sR] - z = blk.array["z"][sR] + for i, s0 in enumerate(face.s0_): + x = blk.array["x"][s0] + y = blk.array["y"][s0] + z = blk.array["z"][s0] # Translate periodics if face.bcType == "periodicTransLow": From b92b41421a6c48c9a6119ff06278a59507d322fe Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 21:30:37 -0400 Subject: [PATCH 28/31] Using underscore for functions that arent called by user. --- src/peregrinepy/multiBlock/solverBlock.py | 4 ++-- src/peregrinepy/multiBlock/solverFace.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/peregrinepy/multiBlock/solverBlock.py b/src/peregrinepy/multiBlock/solverBlock.py index e0f3494c..b61551d0 100644 --- a/src/peregrinepy/multiBlock/solverBlock.py +++ b/src/peregrinepy/multiBlock/solverBlock.py @@ -240,8 +240,8 @@ def setBlockCommunication(self): for face in self.faces: if face.neighbor is None: continue - face.setOrientFunc(self.ni, self.nj, self.nk, self.ne) - face.setCommBuffers(self.ni, self.nj, self.nk, self.ne, self.nblki) + face._setOrientFunc(self.ni, self.nj, self.nk, self.ne) + face._setCommBuffers(self.ni, self.nj, self.nk, self.ne, self.nblki) def updateDeviceView(self, vars): if type(vars) == str: diff --git a/src/peregrinepy/multiBlock/solverFace.py b/src/peregrinepy/multiBlock/solverFace.py index 85653e12..49605c9f 100644 --- a/src/peregrinepy/multiBlock/solverFace.py +++ b/src/peregrinepy/multiBlock/solverFace.py @@ -95,9 +95,9 @@ def __init__(self, nface, ng): @topologyFace.bcType.setter def bcType(self, value): super(solverFace, type(self)).bcType.fset(self, value) - self.setBcFunc() + self._setBcFunc() - def setBcFunc(self): + def _setBcFunc(self): bcType = self.bcType if bcType == "b0" or bcType.startswith("periodicTrans"): @@ -112,7 +112,7 @@ def setBcFunc(self): else: raise KeyError(f"{bcType} is not a valid bcType") - def setCommBuffers(self, ni, nj, nk, ne, nblki): + def _setCommBuffers(self, ni, nj, nk, ne, nblki): assert ( self.orient is not None ), "Must set orientFunc before commBuffers, or perhaps you are trying to set buffers for a non communication face." @@ -240,7 +240,7 @@ def setCommBuffers(self, ni, nj, nk, ne, nblki): self.tagR = int(nblki * 6 + self.nface) self.tagS = int(self.neighbor * 6 + self.neighborNface) - def setOrientFunc(self, ni, nj, nk, ne): + def _setOrientFunc(self, ni, nj, nk, ne): assert 0 not in [ ni, nj, From c41d12ca8fd7d15ad6d2caace5bb0d088d0807e8 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 21:31:14 -0400 Subject: [PATCH 29/31] Added mpi init to the conftest, allows for multiple tests to use it. Also added a periodic test case. --- tests/boundaryConditions/test_periodics.py | 147 +++++++++++++++++++++ tests/communication/test_orientation.py | 4 +- tests/conftest.py | 9 +- tests/shock/test_rusanov_sod.py | 2 - 4 files changed, 157 insertions(+), 5 deletions(-) create mode 100644 tests/boundaryConditions/test_periodics.py diff --git a/tests/boundaryConditions/test_periodics.py b/tests/boundaryConditions/test_periodics.py new file mode 100644 index 00000000..75eb66ce --- /dev/null +++ b/tests/boundaryConditions/test_periodics.py @@ -0,0 +1,147 @@ +import itertools + +import numpy as np +import peregrinepy as pg +import pytest + + +############################################## +# Test all inlet boundary conditions +############################################## + +pytestmark = pytest.mark.parametrize( + "adv,spdata", + list( + itertools.product( + ("secondOrderKEEP", "fourthOrderKEEP"), + (["Air"], "thtr_CH4_O2_Stanford_Skeletal.yaml"), + ) + ), +) + + +class TestPeriodics: + @classmethod + def setup_class(self): + # MPI.Init() + pass + + @classmethod + def teardown_class(self): + # MPI.Finalize() + pass + + def test_rotationalPeriodics(self, my_setup, adv, spdata): + + config = pg.files.configFile() + config["RHS"]["primaryAdvFlux"] = adv + config["thermochem"]["spdata"] = spdata + + mb = pg.multiBlock.generateMultiBlockSolver(1, config) + + axis = np.random.random(3) + axis /= np.linalg.norm(axis) + sweep = np.random.randint(1, 90) + p3 = np.random.random(3) + p3 /= np.linalg.norm(axis) + p3[0] = (-axis[1] * p3[1] - axis[2] * p3[2]) / axis[0] + + pg.grid.create.multiBlockAnnulus( + mb, + sweep=sweep, + p2=axis, + p3=p3, + periodic=True, + ) + blk = mb[0] + blk.getFace(5).commRank = 0 + blk.getFace(6).commRank = 0 + + mb.setBlockCommunication() + mb.initSolverArrays(config) + + mb.unifyGrid() + mb.computeMetrics(config["RHS"]["diffOrder"]) + + qshape = blk.array["q"][:, :, :, 0].shape + # NOTE: Nov, 2021 KAS: The currently un protected extrapolation of + # boundary conditions were making the constantMassFluxInlet test case + # behave poorly (negative species, etc.). So instead of random physical + # values everywhere we narrow the scope a bit. Maybe down the line + # we see how necessary it is to protect those BC extraplations. + p = np.random.uniform(low=101325 * 0.9, high=101325 * 1.1) + u = np.random.uniform(low=1, high=1000, size=qshape) + v = np.random.uniform(low=1, high=1000, size=qshape) + w = np.random.uniform(low=1, high=1000, size=qshape) + T = np.random.uniform(low=300 * 0.9, high=300 * 1.1) + + if blk.ns > 1: + Y = np.random.uniform(low=0.0, high=1.0, size=(blk.ns - 1)) + Y = Y / np.sum(Y) + + blk.array["q"][:, :, :, 0] = p + blk.array["q"][:, :, :, 1] = u + blk.array["q"][:, :, :, 2] = v + blk.array["q"][:, :, :, 3] = w + blk.array["q"][:, :, :, 4] = T + if blk.ns > 1: + blk.array["q"][:, :, :, 5::] = Y + blk.updateDeviceView("q") + + mb.eos(blk, mb.thtrdat, 0, "prims") + pg.consistify(mb) + + u = blk.array["q"][:, :, :, 1] + v = blk.array["q"][:, :, :, 2] + w = blk.array["q"][:, :, :, 3] + + dqdx = blk.array["dqdx"] + dqdy = blk.array["dqdy"] + dqdz = blk.array["dqdz"] + + nx = blk.array["knx"] + ny = blk.array["kny"] + nz = blk.array["knz"] + + ng = blk.ng + for g in range(ng): + + # face5 halo compares to interior on side 6 + s5 = np.s_[:, :, g] + s6c = np.s_[:, :, -2 * ng + g] + s6f = np.s_[:, :, -2 * ng - 1 + g] + + normals5 = np.column_stack((nx[s5].ravel(), ny[s5].ravel(), nz[s5].ravel())) + velo5 = np.column_stack((u[s5].ravel(), v[s5].ravel(), w[s5].ravel())) + + normals6 = np.column_stack( + (nx[s6f].ravel(), ny[s6f].ravel(), nz[s6f].ravel()) + ) + velo6 = np.column_stack((u[s6c].ravel(), v[s6c].ravel(), w[s6c].ravel())) + + # dot the velocities with the face normals kills two birds with one stone + # check that velo and gradients are rotated correctly and also halo + assert np.allclose( + np.sum(normals5 * velo5, axis=1), np.sum(normals6 * velo6, axis=1) + ) + + # check the gradients + for i in range(blk.ne): + dqdx5 = np.column_stack( + ( + dqdx[s5][:, :, i].ravel(), + dqdy[s5][:, :, i].ravel(), + dqdz[s5][:, :, i].ravel(), + ) + ) + + dqdx6 = np.column_stack( + ( + dqdx[s6c][:, :, i].ravel(), + dqdy[s6c][:, :, i].ravel(), + dqdz[s6c][:, :, i].ravel(), + ) + ) + assert np.allclose( + np.sum(normals5 * dqdx5, axis=1), np.sum(normals6 * dqdx6, axis=1) + ) diff --git a/tests/communication/test_orientation.py b/tests/communication/test_orientation.py index ff3ed259..a6c401bc 100644 --- a/tests/communication/test_orientation.py +++ b/tests/communication/test_orientation.py @@ -57,11 +57,11 @@ def __init__(self, adv, spdata): class TestOrientation: @classmethod def setup_class(self): - MPI.Init() + pass @classmethod def teardown_class(self): - MPI.Finalize() + pass ############################################## # Test for all positive i aligned orientations diff --git a/tests/conftest.py b/tests/conftest.py index 506603e8..5ee2772a 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,12 +1,19 @@ -import pytest +import mpi4py.rc + +mpi4py.rc.finalize = False +mpi4py.rc.initialize = False import kokkos +import pytest +from mpi4py import MPI @pytest.fixture(scope="session") def my_setup(request): kokkos.initialize() + MPI.Init() def fin(): kokkos.finalize() + MPI.Finalize() request.addfinalizer(fin) diff --git a/tests/shock/test_rusanov_sod.py b/tests/shock/test_rusanov_sod.py index 3e8edd6f..ded4b97e 100755 --- a/tests/shock/test_rusanov_sod.py +++ b/tests/shock/test_rusanov_sod.py @@ -1,10 +1,8 @@ import peregrinepy as pg import numpy as np from pathlib import Path -from pytest_easyMPI import mpi_parallel -@mpi_parallel(1) def test_rusanov_sod(my_setup): # Left State From 8c641a875e7781c44f493f15a9f655e9a3797960 Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 21:46:02 -0400 Subject: [PATCH 30/31] Cleaned up examles to use new create scheme (after periodic addition) --- examples/couetteFlow.py | 19 +------------------ examples/oneDAdvection.py | 10 ++-------- examples/threeDTaylorGreen.py | 7 +------ examples/threeDTaylorGreenSkewed.py | 5 +---- examples/twoDEulerVortex.py | 13 +++---------- 5 files changed, 8 insertions(+), 46 deletions(-) diff --git a/examples/couetteFlow.py b/examples/couetteFlow.py index 1827ead1..5865a009 100755 --- a/examples/couetteFlow.py +++ b/examples/couetteFlow.py @@ -41,6 +41,7 @@ def simulate(): mbDims=[1, 1, 1], dimsPerBlock=[2, ny, 2], lengths=[0.01, h, 0.01], + periodic=[True, False, False], ) mb.initSolverArrays(config) @@ -48,34 +49,16 @@ def simulate(): blk = mb[0] # face 1 - blk.getFace(1).bcFam = None - blk.getFace(1).bcType = "b1" - blk.getFace(1).neighbor = 0 - blk.getFace(1).orientation = "123" blk.getFace(1).commRank = 0 # face 2 - blk.getFace(2).bcFam = None - blk.getFace(2).bcType = "b1" - blk.getFace(2).neighbor = 0 - blk.getFace(2).orientation = "123" blk.getFace(2).commRank = 0 # face 3 - blk.getFace(3).bcFam = None blk.getFace(3).bcType = "adiabaticNoSlipWall" - blk.getFace(3).neighbor = None - blk.getFace(3).orientation = None - blk.getFace(3).commRank = 0 # face 4 isoT moving wall - blk.getFace(4).bcFam = "whoosh" blk.getFace(4).bcType = "adiabaticMovingWall" - blk.getFace(4).neighbor = None - blk.getFace(4).orientation = None - blk.getFace(4).commRank = 0 for face in [5, 6]: - blk.getFace(face).bcFam = None blk.getFace(face).bcType = "adiabaticSlipWall" - blk.getFace(face).commRank = 0 valueDict = {"u": wallSpeed, "v": 0.0, "w": 0.0} face4 = blk.getFace(4) diff --git a/examples/oneDAdvection.py b/examples/oneDAdvection.py index 5a4de696..ebe5553b 100755 --- a/examples/oneDAdvection.py +++ b/examples/oneDAdvection.py @@ -28,21 +28,15 @@ def simulate(): mbDims=[1, 1, 1], dimsPerBlock=[41, 2, 2], lengths=[1, 0.01, 0.01], + periodic=[True, False, False], ) mb.initSolverArrays(config) blk = mb[0] - for face in blk.faces: + for face in blk.faces[2::]: face.bcType = "adiabaticSlipWall" - blk.getFace(1).bcType = "b1" - blk.getFace(1).neighbor = 0 - blk.getFace(1).orientation = "123" blk.getFace(1).commRank = 0 - - blk.getFace(2).bcType = "b1" - blk.getFace(2).neighbor = 0 - blk.getFace(2).orientation = "123" blk.getFace(2).commRank = 0 mb.setBlockCommunication() diff --git a/examples/threeDTaylorGreen.py b/examples/threeDTaylorGreen.py index ac144ffe..9870dbab 100755 --- a/examples/threeDTaylorGreen.py +++ b/examples/threeDTaylorGreen.py @@ -33,6 +33,7 @@ def simulate(): mbDims=[1, 1, 1], dimsPerBlock=[NE, NN, NX], lengths=[2 * np.pi for _ in range(3)], + periodic=[True, True, True], ) mb.initSolverArrays(config) @@ -41,13 +42,7 @@ def simulate(): ng = blk.ng for face in blk.faces: - face.bcType = "b1" - face.neighbor = 0 - face.orientation = "123" face.commRank = 0 - pg.writers.writeGrid(mb) - pg.writers.writeConnectivity(mb) - mb.generateHalo() mb.setBlockCommunication() diff --git a/examples/threeDTaylorGreenSkewed.py b/examples/threeDTaylorGreenSkewed.py index da9cb31f..4d239bf5 100755 --- a/examples/threeDTaylorGreenSkewed.py +++ b/examples/threeDTaylorGreenSkewed.py @@ -39,6 +39,7 @@ def simulate(): mbDims=[1, 1, 1], dimsPerBlock=[NE, NN, NX], lengths=[2 * np.pi for _ in range(3)], + periodic=[True, True, True], ) mb.initSolverArrays(config) @@ -85,10 +86,6 @@ def simulate(): ) for face in blk.faces: - face.bcType = "b1" - face.bcFam = None - face.neighbor = 0 - face.orientation = "123" face.commRank = 0 mb.setBlockCommunication() diff --git a/examples/twoDEulerVortex.py b/examples/twoDEulerVortex.py index 46be8536..a8f25926 100755 --- a/examples/twoDEulerVortex.py +++ b/examples/twoDEulerVortex.py @@ -31,7 +31,8 @@ def simulate(): mb, mbDims=[1, 1, 1], dimsPerBlock=[NE, NN, 2], - lengths=[1, 1, 0.01], + lengths=[12, 12, 0.01], + periodic=[True, True, False], ) blk = mb[0] @@ -59,19 +60,11 @@ def simulate(): + Ay * np.sin(2 * np.pi * kappa) * np.sin(lamY * np.pi * E * delX / Lx) ) - for face in blk.faces: - face.bcType = "b1" - face.bcFam = None - face.neighbor = 0 - face.orientation = "123" + for face in blk.faces[0:4]: face.commRank = 0 for f in [5, 6]: face = blk.getFace(f) face.bcType = "adiabaticSlipWall" - face.bcFam = None - face.neighbor = None - face.orientation = "000" - face.commRank = None mb.initSolverArrays(config) mb.setBlockCommunication() From 3770aed08c6d42f16cc590bc69bd0259a02cc5dc Mon Sep 17 00:00:00 2001 From: Kyle Schau Date: Wed, 6 Apr 2022 23:23:47 -0400 Subject: [PATCH 31/31] Passing all tests on cpu and gpu --- src/compute/boundaryConditions/periodics.cpp | 8 ++-- src/peregrinepy/consistify.py | 1 - src/peregrinepy/grid/unifySolverGrid.py | 6 ++- tests/boundaryConditions/test_periodics.py | 45 ++++++++++++++++++-- tests/conftest.py | 2 +- 5 files changed, 51 insertions(+), 11 deletions(-) diff --git a/src/compute/boundaryConditions/periodics.cpp b/src/compute/boundaryConditions/periodics.cpp index b9478185..5955d855 100644 --- a/src/compute/boundaryConditions/periodics.cpp +++ b/src/compute/boundaryConditions/periodics.cpp @@ -7,8 +7,8 @@ void periodicRotHigh(block_ b, face_ face, - const std::function &eos, - const thtrdat_ th, + const std::function &/*eos*/, + const thtrdat_ /*th*/, const std::string terms, const double /*tme*/) { //-------------------------------------------------------------------------------------------| @@ -87,8 +87,8 @@ void periodicRotHigh(block_ b, void periodicRotLow(block_ b, face_ face, - const std::function &eos, - const thtrdat_ th, + const std::function &/*eos*/, + const thtrdat_ /*th*/, const std::string terms, const double /*tme*/) { //-------------------------------------------------------------------------------------------| diff --git a/src/peregrinepy/consistify.py b/src/peregrinepy/consistify.py index a92f10f8..381ccb0d 100644 --- a/src/peregrinepy/consistify.py +++ b/src/peregrinepy/consistify.py @@ -27,7 +27,6 @@ def consistify(mb): # Update spatial derivatives mb.dqdxyz(blk) - # TODO: can we get rid of this if check? if mb.config["RHS"]["diffusion"]: # Apply viscous boundary conditions for face in blk.faces: diff --git a/src/peregrinepy/grid/unifySolverGrid.py b/src/peregrinepy/grid/unifySolverGrid.py index be199903..3c586989 100644 --- a/src/peregrinepy/grid/unifySolverGrid.py +++ b/src/peregrinepy/grid/unifySolverGrid.py @@ -11,6 +11,10 @@ def unifySolverGrid(mb): for _ in range(3): mpiComm.communicate(mb, ["x", "y", "z"]) + # Device is up to date after communicate, so pull back down + for blk in mb: + blk.updateHostView(["x", "y", "z"]) + for blk in mb: for face in blk.faces: bc = face.bcType @@ -42,6 +46,6 @@ def unifySolverGrid(mb): y[:] = tempy[:] z[:] = tempz[:] + # Push back up the device for blk in mb: - # Push back up the device blk.updateDeviceView(["x", "y", "z"]) diff --git a/tests/boundaryConditions/test_periodics.py b/tests/boundaryConditions/test_periodics.py index 75eb66ce..190510c6 100644 --- a/tests/boundaryConditions/test_periodics.py +++ b/tests/boundaryConditions/test_periodics.py @@ -23,12 +23,10 @@ class TestPeriodics: @classmethod def setup_class(self): - # MPI.Init() pass @classmethod def teardown_class(self): - # MPI.Finalize() pass def test_rotationalPeriodics(self, my_setup, adv, spdata): @@ -91,6 +89,8 @@ def test_rotationalPeriodics(self, my_setup, adv, spdata): mb.eos(blk, mb.thtrdat, 0, "prims") pg.consistify(mb) + blk.updateHostView(["q", "dqdx", "dqdy", "dqdz"]) + u = blk.array["q"][:, :, :, 1] v = blk.array["q"][:, :, :, 2] w = blk.array["q"][:, :, :, 3] @@ -119,8 +119,6 @@ def test_rotationalPeriodics(self, my_setup, adv, spdata): ) velo6 = np.column_stack((u[s6c].ravel(), v[s6c].ravel(), w[s6c].ravel())) - # dot the velocities with the face normals kills two birds with one stone - # check that velo and gradients are rotated correctly and also halo assert np.allclose( np.sum(normals5 * velo5, axis=1), np.sum(normals6 * velo6, axis=1) ) @@ -145,3 +143,42 @@ def test_rotationalPeriodics(self, my_setup, adv, spdata): assert np.allclose( np.sum(normals5 * dqdx5, axis=1), np.sum(normals6 * dqdx6, axis=1) ) + + # Now check the other way + # face6 halo compares to interior on side 5 + s6 = np.s_[:, :, -(g + 1)] + s5c = np.s_[:, :, 2 * ng - g - 1] + s5f = np.s_[:, :, 2 * ng - g] + + normals6 = np.column_stack((nx[s6].ravel(), ny[s6].ravel(), nz[s6].ravel())) + velo6 = np.column_stack((u[s6].ravel(), v[s6].ravel(), w[s6].ravel())) + + normals5 = np.column_stack( + (nx[s5f].ravel(), ny[s5f].ravel(), nz[s5f].ravel()) + ) + velo5 = np.column_stack((u[s5c].ravel(), v[s5c].ravel(), w[s5c].ravel())) + + assert np.allclose( + np.sum(normals6 * velo6, axis=1), np.sum(normals5 * velo5, axis=1) + ) + + # check the gradients + for i in range(blk.ne): + dqdx6 = np.column_stack( + ( + dqdx[s6][:, :, i].ravel(), + dqdy[s6][:, :, i].ravel(), + dqdz[s6][:, :, i].ravel(), + ) + ) + + dqdx5 = np.column_stack( + ( + dqdx[s5c][:, :, i].ravel(), + dqdy[s5c][:, :, i].ravel(), + dqdz[s5c][:, :, i].ravel(), + ) + ) + assert np.allclose( + np.sum(normals6 * dqdx6, axis=1), np.sum(normals5 * dqdx5, axis=1) + ) diff --git a/tests/conftest.py b/tests/conftest.py index 5ee2772a..9cf31f33 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -9,8 +9,8 @@ @pytest.fixture(scope="session") def my_setup(request): - kokkos.initialize() MPI.Init() + kokkos.initialize() def fin(): kokkos.finalize()