diff --git a/.gitignore b/.gitignore index de5fba0..edeff61 100644 --- a/.gitignore +++ b/.gitignore @@ -28,5 +28,6 @@ src_cs/f2py/get_f2py.py idwarp.egg-info input_files/*/ input_files/*.tar.gz +input_files/*.json testflo_report.out doc/_build diff --git a/idwarp/MultiUnstructuredMesh.py b/idwarp/MultiUnstructuredMesh.py index 1d2244d..e5978da 100644 --- a/idwarp/MultiUnstructuredMesh.py +++ b/idwarp/MultiUnstructuredMesh.py @@ -1,14 +1,12 @@ #!/usr/bin/python """ -Multiple USMesh instances contained in this MultiUSMesh object. -The MultiUnstructuredMesh module is used for interacting with multiple -unstructured (or structured!) meshes - typically used in a 3D CFD -program. +The MultiUnstructuredMesh module is used for creating multiple +USMesh instances, typically from a structured overset mesh. It contains the following classes: -MultiUSMesh: General class for working with multiple unstructured meshes +MultiUSMesh: General class for working with multi-component meshes Copyright (c) 2017 by John Jasa All rights reserved. Not to be used for commercial purposes. @@ -28,8 +26,10 @@ # ============================================================================= import os import numpy as np +from random import randint from mpi4py import MPI from .MExt import MExt +from . import USMesh, USMesh_C from baseclasses.utils import Error try: @@ -91,9 +91,15 @@ def __init__(self, CGNSFile, optionsDict, comm=None, dtype="d", debug=False): # Store original file name self.CGNSFile = CGNSFile + # Get rank of the current proc + self.myID = self.comm.Get_rank() + # Store scalar type self.dtype = dtype + # Set a random prefix to avoid I/O clashes between instances + prefix = f"tmp{randint(1, 1e5)}" + # Only the root processor will take the combined CGNS file # and explode it by instance. if self.myID == 0: @@ -124,7 +130,7 @@ def __init__(self, CGNSFile, optionsDict, comm=None, dtype="d", debug=False): # Save temporary grid files with the exploded zones for grid, zoneName in zip(grids, zoneNames): - grid.writeToCGNS("_" + zoneName + ".cgns") + grid.writeToCGNS(f"{prefix}_{zoneName}.cgns") # Store the number of blocks in each zone self.cgnsBlockIntervals.append([blockCounter, blockCounter + len(grid.blocks)]) @@ -182,20 +188,20 @@ def __init__(self, CGNSFile, optionsDict, comm=None, dtype="d", debug=False): # Assign the name of the temporary CGNS file to the options. # This is the file that contains the mesh o a single component. # Remember that we should use the temporary grid file. - optionsDict[zoneName]["gridFile"] = "_" + zoneName + ".cgns" + optionsDict[zoneName]["gridFile"] = f"{prefix}_{zoneName}.cgns" # Initialize an IDWarp instance with the current options if self.dtype == "d": - currMesh = self.warp.USMesh(options=optionsDict[zoneName], comm=self.comm) + currMesh = USMesh(options=optionsDict[zoneName], comm=self.comm) elif self.dtype == "D": - currMesh = self.warp.USMesh_C(options=optionsDict[zoneName], comm=self.comm) + currMesh = USMesh_C(options=optionsDict[zoneName], comm=self.comm) else: # We have a background mesh # Regenerate the temporary filename for the background grid - bgFile = "_" + zoneName + ".cgns" + bgFile = f"{prefix}_{zoneName}.cgns" # ------------------------------------------------------ # READING BACKGROUND MESHES @@ -217,26 +223,25 @@ def __init__(self, CGNSFile, optionsDict, comm=None, dtype="d", debug=False): if self.myID == 0: # Make a copy of the background mesh file - os.system("cp " + bgFile + " tmp_bg_file.cgns") + os.system(f"cp {bgFile} {prefix}_bg_file.cgns") # Create a temporary BC file - with open("tmp_bcdata.dat", "w") as fid: + with open(f"{prefix}_bcdata.dat", "w") as fid: fid.write("1 iLow BCwall wall\n") # Use CGNS utils to modify the BCs - os.system("cgns_utils overwritebc tmp_bg_file.cgns tmp_bcdata.dat") + os.system(f"cgns_utils overwriteBC {prefix}_bg_file.cgns {prefix}_bcdata.dat") # Create dummy set of options just to load the CGNS file dummyOptions = { - "gridFile": "tmp_bg_file.cgns", - "warpType": "unstructured", + "gridFile": f"{prefix}_bg_file.cgns", } # Initialize an IDWarp instance with the current options if self.dtype == "d": - currMesh = self.warp.USMesh(options=dummyOptions, comm=self.comm) + currMesh = USMesh(options=dummyOptions, comm=self.comm) elif self.dtype == "D": - currMesh = self.warp.USMesh_C(options=dummyOptions, comm=self.comm) + currMesh = USMesh_C(options=dummyOptions, comm=self.comm) # Initialize a dummy surface in the background mesh """ @@ -268,8 +273,8 @@ def __init__(self, CGNSFile, optionsDict, comm=None, dtype="d", debug=False): if self.myID == 0: # Make a copy of the background mesh file - os.system("rm tmp_bg_file.cgns") - os.system("rm tmp_bcdata.dat") + os.system(f"rm {prefix}_bg_file.cgns") + os.system(f"rm {prefix}_bcdata.dat") # Store the ID of this zone self.backgroundInstanceIDs.append(zoneNumber) @@ -282,7 +287,7 @@ def __init__(self, CGNSFile, optionsDict, comm=None, dtype="d", debug=False): # Now the root proc can remove the temporary grid files if self.myID == 0: for zoneName in zoneNames: - os.system("rm _" + zoneName + ".cgns") + os.system(f"rm {prefix}_{zoneName}.cgns") # ------------------------------------------------------ # Initialize other fields for completness @@ -617,19 +622,13 @@ def setSurfaceDefinition(self, pts, conn=None, faceSizes=None, cgnsBlockIDs=None print("Done") print("") - ''' - def getWarpGrid(self, warp_to_common=True): + def getWarpGrid(self): """ Return the current grids. This function is typically unused. See getSolverGrid for the more useful interface functionality. This only returns the nearfield meshes. - warp_to_common is a flag indicating if we should do an scatter - operation to update commonGridVec with the values in Xv. - If you want to get grid coordinates prior to initializing the full - warping procedure, then you can set this value to false. - Returns ------- volNodesList, list of 1D numpy arrays, real: These are the local volume nodes (in a flat 1D array) @@ -656,7 +655,7 @@ def getWarpGrid(self, warp_to_common=True): # Get volume nodes. # volNodes is a flattened vector that contains the background # mesh volume nodes that belong to the current proc. - volNodes = currMesh.getCommonGrid(warp_to_common) + volNodes = currMesh.getCommonGrid() # Store the nodes of the current instance in the list volNodesList.append(volNodes) @@ -673,7 +672,6 @@ def getWarpGrid(self, warp_to_common=True): # RETURNS return volNodesList, numCoorTotal - ''' def getdXs(self): """Return the current values in dXs. This is the result from a @@ -721,6 +719,7 @@ def warpMesh(self): if self.myID == 0: print("") print("Starting IDWarpMulti mesh warping") + print("") # Set mesh counter meshCounter = 1 @@ -730,27 +729,21 @@ def warpMesh(self): # Print log if self.myID == 0: - print("") - print(" warping mesh", meshCounter, "of", len(self.meshes)) + print(" Warping mesh", meshCounter, "of", len(self.meshes)) # Warp current instance mesh.warpMesh() - # Print log - if self.myID == 0: - print("") - print(" Done") - # Increment counter meshCounter = meshCounter + 1 # Print log if self.myID == 0: print("") - print("IDWarpMulti successfully warped all instances!") + print("IDWarpMulti finished warping all instances") print("") - def warpDeriv(self, dXv): + def warpDeriv(self, dXv, solverVec=True): """Compute the warping derivative (dXv/dXs^T)*Vec (where vec is the dXv argument to this function. @@ -759,7 +752,7 @@ def warpDeriv(self, dXv): Parameters ---------- - solverdXv : numpy array + dXv : numpy array Vector of size external solver_grid. This is typically obtained from the external solver's dRdx^T * psi calculation. @@ -785,6 +778,7 @@ def warpDeriv(self, dXv): if self.myID == 0: print("") print("Starting IDWarpMulti reverse AD") + print("") # --------------------------------------------------- # RUN REVERSE AD ON EACH INSTANCE @@ -794,7 +788,6 @@ def warpDeriv(self, dXv): # Print log if self.myID == 0: - print("") print(" Working on mesh", instanceID + 1, "of", len(self.meshes)) # Get current seeds @@ -802,17 +795,12 @@ def warpDeriv(self, dXv): # Run reverse AD. # This will update the surface seeds inside the mesh object - mesh.warpDeriv(curr_dXv) - - # Print log - if self.myID == 0: - print("") - print(" Done") + mesh.warpDeriv(curr_dXv, solverVec=solverVec) # Print log if self.myID == 0: print("") - print("IDWarpMulti successfully finished reverse AD on all instances!") + print("IDWarpMulti finished reverse AD on all instances") print("") # Get derivative seeds @@ -850,6 +838,7 @@ def warpDerivFwd(self, dXs): if self.myID == 0: print("") print("Starting IDWarpMulti forward AD") + print("") # --------------------------------------------------- # SPLIT SURFACE SEEDS ACROSS ALL INSTANCES AND GATHER @@ -863,7 +852,6 @@ def warpDerivFwd(self, dXs): # Print log if self.myID == 0: - print("") print(" Working on mesh", instanceID + 1, "of", len(self.meshes)) # Get current surface node forward AD seeds @@ -875,39 +863,16 @@ def warpDerivFwd(self, dXs): # Add seeds to the full volume vector dXv[self.cgnsVolNodeMasks[instanceID]] = curr_dXvWarp - # Print log - if self.myID == 0: - print("") - print(" Done") - # Print log if self.myID == 0: print("") - print("IDWarpMulti successfully finished forward AD on all instances!") + print("IDWarpMulti finished forward AD on all instances") print("") # --------------------------------------------------- # RETURNS return dXv - ''' - def verifyWarpDeriv(self, dXv=None, solverVec=True, dofStart=0, - dofEnd=10, h=1e-6, randomSeed=314): - """Run an internal verification of the solid warping - derivatives""" - - if dXv is None: - np.random.seed(randomSeed) # 'Random' seed to ensure runs are same - dXvWarp = np.random.random(self.warp.griddata.warpmeshdof) - else: - if solverVec: - dXvWarp = np.zeros(self.warp.griddata.warpmeshdof, self.dtype) - self.warp.solver_to_warp_grid(dXv, dXvWarp) - else: - dXvWarp = dXv - - self.warp.verifywarpderiv(dXvWarp, dofStart, dofEnd, h) - ''' # ========================================================================== # Output Functionality # ========================================================================== diff --git a/idwarp/MultiUnstructuredMesh_C.py b/idwarp/MultiUnstructuredMesh_C.py index 81d9bec..20d09d7 100644 --- a/idwarp/MultiUnstructuredMesh_C.py +++ b/idwarp/MultiUnstructuredMesh_C.py @@ -1,7 +1,7 @@ #!/usr/bin/python """ -This is the complex version of the Unstructured mesh warping. See -UnstructuredMesh.py for more info. +This is the complex version of the MultiUSMesh class. +See MultiUnstructuredMesh.py for more info. """ # ============================================================================= @@ -18,7 +18,7 @@ class MultiUSMesh_C(MultiUSMesh): """ - Represents a (Complex) Unstructured mesh + Represents a (complex) multi-component mesh """ def __init__(self, *args, **kwargs): diff --git a/idwarp/UnstructuredMesh.py b/idwarp/UnstructuredMesh.py index 81a7ed2..5c34525 100755 --- a/idwarp/UnstructuredMesh.py +++ b/idwarp/UnstructuredMesh.py @@ -399,7 +399,7 @@ def warpDeriv(self, dXv, solverVec=True): Parameters ---------- - solverdXv : numpy array + dXv : numpy array Vector of size external solver_grid. This is typically obtained from the external solver's dRdx^T * psi calculation. diff --git a/idwarp/__init__.py b/idwarp/__init__.py index 3638417..a5446cb 100644 --- a/idwarp/__init__.py +++ b/idwarp/__init__.py @@ -1,8 +1,8 @@ -__version__ = "2.5.0" +__version__ = "2.6.0" from .UnstructuredMesh import USMesh -from .MultiUnstructuredMesh import MultiUSMesh from .UnstructuredMesh_C import USMesh_C +from .MultiUnstructuredMesh import MultiUSMesh from .MultiUnstructuredMesh_C import MultiUSMesh_C __all__ = ["USMesh", "MultiUSMesh", "USMesh_C", "MultiUSMesh_C"] diff --git a/tests/ref/test_onera_m6.ref b/tests/ref/test_onera_m6.ref new file mode 100644 index 0000000..4693290 --- /dev/null +++ b/tests/ref/test_onera_m6.ref @@ -0,0 +1,23 @@ +{ + "Test_onera_m6 - Sum of dXs:": 300303.80115047796, + "Test_onera_m6 - Sums of Initial Volume Coords:": { + "__ndarray__": [ + 274789.83833522064, + 75732.81670770115 + ], + "dtype": "float64", + "shape": [ + 2 + ] + }, + "Test_onera_m6 - Sums of Warped Volume Coords:": { + "__ndarray__": [ + 309457.45451071725, + 84883.14062531918 + ], + "dtype": "float64", + "shape": [ + 2 + ] + } +} \ No newline at end of file diff --git a/tests/test_MultiUSMesh.py b/tests/test_MultiUSMesh.py new file mode 100644 index 0000000..679b568 --- /dev/null +++ b/tests/test_MultiUSMesh.py @@ -0,0 +1,192 @@ +# ============================================================================= +# Standard Python modules +# ============================================================================= +import os + +# ============================================================================= +# External Python modules +# ============================================================================= +import numpy as np +import unittest +from mpi4py import MPI +from parameterized import parameterized_class + +# ============================================================================= +# Extension modules +# ============================================================================= +from idwarp import MultiUSMesh, MultiUSMesh_C +from baseclasses import BaseRegTest +from baseclasses.utils import readJSON + +baseDir = os.path.dirname(os.path.abspath(__file__)) # Path to current directory + + +def eval_warp(handler, test_name, gridFile, optionsDict, isComplex, N_PROCS): + # Create warping object + if isComplex: + # Check if the complex verision of the code has been built + try: + from idwarp import libidwarp_cs # noqa: F401 + + h = 1e-200 + mesh = MultiUSMesh_C(gridFile, optionsDict) + except ImportError: + raise unittest.SkipTest("Skipping because you do not have complex idwarp compiled") + else: + try: + from idwarp import libidwarp # noqa: F401 + + mesh = MultiUSMesh(gridFile, optionsDict) + except ImportError: + raise unittest.SkipTest("Skipping because you do not have real idwarp compiled") + + # Get mesh data from JSON file to avoid needing ADflow for the test + jsonFile = os.path.join(baseDir, "../input_files/onera_m6.json") + meshData = readJSON(jsonFile) + + # Set the mesh surface from data + mesh.setExternalMeshIndices(meshData["meshInd"]) + mesh.setSurfaceDefinition(meshData["pts"], meshData["conn"], meshData["faceSizes"], meshData["cgnsBlockIDs"]) + + # Get the sum of the initial volume coordinates + volNodesList = mesh.getWarpGrid()[0] + sums = sumVolCoords(volNodesList) + + # Test the initial sums + handler.root_add_val(f"{test_name} - Sums of Initial Volume Coords:", sums, tol=1e-14) + + # Get initial surface mesh + coords0 = mesh.getSurfaceCoordinates() + + # Displace the surface coordinates with a shearing sweep + new_coords = coords0.copy() + for i in range(len(coords0)): + span = coords0[i, 1] + new_coords[i, 0] += 0.25 * span + mesh.setSurfaceCoordinates(new_coords) + + # Warp the volume mesh + mesh.warpMesh() + + # Get the sum of the warped volume coordinates + volNodesList = mesh.getWarpGrid()[0] + sums = sumVolCoords(volNodesList) + + # Test the warped sums + handler.root_add_val(f"{test_name} - Sums of Warped Volume Coords:", sums, tol=1e-14) + + # Create a volume seed vector to test the mesh warping derivatives + # We divide by the number of procs because we are not partitioning the volume mesh like ADflow would + dXv = np.linspace(0, 1.0, len(meshData["meshInd"])) / N_PROCS + + if not isComplex: + # Compute the warping derivatives (surface seeds) + dXs = mesh.warpDeriv(dXv) + + # Sum the derivatives + val = MPI.COMM_WORLD.reduce(np.sum(dXs.flatten()), op=MPI.SUM) + + else: + # Add a complex perturbation on all surface nodes simultaneously + for i in range(len(coords0)): + new_coords[i, :] += 1j * h + + # Set the complex surface coordinates + mesh.setSurfaceCoordinates(new_coords) + + # Warp the volume mesh + mesh.warpMesh() + + # Sum the derivatives + solverGrid = mesh.getSolverGrid() + deriv = np.imag(solverGrid) / h + val = MPI.COMM_WORLD.reduce(np.dot(dXv, deriv), op=MPI.SUM) + + handler.root_add_val(f"{test_name} - Sum of dXs:", val, tol=1e-14) + + +def sumVolCoords(volNodesList): + """ + Helper function to take the sum of the volume coordinates for each USMesh instance. + """ + + numMeshes = len(volNodesList) + dtype = volNodesList[0].dtype + sums = np.zeros(numMeshes, dtype=dtype) + + for i in range(numMeshes): + sums[i] = MPI.COMM_WORLD.reduce(np.sum(volNodesList[i]), op=MPI.SUM) + + return sums + + +# Set up parameterized variables +test_params = [ + {"N_PROCS": 1, "isComplex": False}, + {"N_PROCS": 2, "name": "parallel", "isComplex": False}, + {"N_PROCS": 1, "name": "complex", "isComplex": True}, + {"N_PROCS": 2, "name": "parallel_complex", "isComplex": True}, +] + + +@parameterized_class(test_params) +class Test_MultiUSMesh(unittest.TestCase): + """ + testflo does not run the ``train_`` functions by default. + Use ``testflo -m "train_*" -n 1`` to regenerate the .ref files. + """ + + N_PROCS = 1 + + def train_onera_m6(self): + # Skip training complex tests, parallel tests, and the extra test that runs because of a bug in parameterized + if not hasattr(self, "isComplex") or self.isComplex or self.N_PROCS != 1: + raise unittest.SkipTest() + else: + self.test_onera_m6(train=True) + + def test_onera_m6(self, train=False): + ref_file = os.path.join(baseDir, "ref/test_onera_m6.ref") + with BaseRegTest(ref_file, train=train) as handler: + + # Test the ONERA M6 overset mesh + test_name = "Test_onera_m6" + gridFile = os.path.join(baseDir, "../input_files/onera_m6.cgns") + + # This instance tests default options + nearTipOptions = { + "gridFile": gridFile, + } + + # This instance tests non-default options + nearWingOptions = { + "gridFile": gridFile, + "aExp": 4.0, + "bExp": 6.0, + "LdefFact": 100.0, + "alpha": 0.3, + "errTol": 1e-5, + "useRotations": False, + "zeroCornerRotations": False, + "cornerAngle": 45.0, + "bucketSize": 4, + } + + optionsDict = { + "near_tip_vol_L3": nearTipOptions, + "near_wing_vol_L3": nearWingOptions, + } + + # Call warping test function + eval_warp( + handler=handler, + test_name=test_name, + gridFile=gridFile, + optionsDict=optionsDict, + isComplex=self.isComplex, + N_PROCS=self.N_PROCS, + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_usmesh.py b/tests/test_USMesh.py similarity index 100% rename from tests/test_usmesh.py rename to tests/test_USMesh.py diff --git a/tests/writeMeshData.py b/tests/writeMeshData.py new file mode 100644 index 0000000..fc0bbc4 --- /dev/null +++ b/tests/writeMeshData.py @@ -0,0 +1,40 @@ +""" +This script was used to generate the onera_m6.json input file needed for the MultiUSMesh test. +It could also be useful for creating other tests in the future. +We use ADflow to get mesh information and then write it as a dictionary in a JSON file. +Unlike USMesh, MultiUSMesh needs surface information before generating internal surfaces. +By saving the mesh data, we avoid adding ADflow as a testing dependency. +""" + +import os +from adflow import ADFLOW +from baseclasses.utils import writeJSON + +baseDir = os.path.dirname(os.path.abspath(__file__)) +gridFile = os.path.join(baseDir, "../input_files/onera_m6.cgns") + +# Set up ADflow +aeroOptions = { + "gridFile": gridFile, + "outputDirectory": baseDir, + "MGCycle": "sg", # This is needed to avoid problems with coarse grid initialization +} +CFDSolver = ADFLOW(options=aeroOptions) + +# Get the mesh data from ADflow +meshInd = CFDSolver.getSolverMeshIndices() +conn, faceSizes, cgnsBlockIDs = CFDSolver.getSurfaceConnectivity( + CFDSolver.meshFamilyGroup, includeZipper=False, includeCGNS=True +) +pts = CFDSolver.getSurfaceCoordinates(CFDSolver.meshFamilyGroup, includeZipper=False) + +# Save the ADflow output +meshData = { + "meshInd": meshInd, + "conn": conn, + "faceSizes": faceSizes, + "cgnsBlockIDs": cgnsBlockIDs, + "pts": pts, +} +jsonFile = os.path.join(baseDir, "onera_m6.json") +writeJSON(jsonFile, meshData)