(feat) Added mesh update capability to the solver interface

New surface mesh is imposed in the viscous solver in 2D and stagnation point is recomputed at the first coupling iteration

blast/coupler.py

@@ -85,6 +85,7 @@ class Coupler:
             # Write inviscid Cp distribution.
             if couplIter == 0:
                 self.isol.initializeViscousSolver()
+                self.isol.updateStagnation()
                 if write:
                     self.isol.writeCp(sfx='_inviscid'+self.filesfx)
 
@@ -142,11 +143,7 @@ class Coupler:
         # Inviscid solver
 
         self.isol.reset()
-        # Viscous solver
-        for s in self.isol.sec:
-            for reg in s:
-                reg.reset()
-        
+
         # Blowing velocity
         for b in self.isol.iBnd:
             for i in range(len(b.blowingVel)):
@@ -156,22 +153,28 @@ class Coupler:
                 for i in range(len(side.blowingVel)):
                     side.blowingVel[i] = 0.
         self.isol.setBlowingVelocity()
+
+        # Viscous solver
+        if self.isol.vinit:
+            for s in self.isol.sec:
+                for reg in s:
+                    reg.reset()
         
-        for iSec in range(self.isol.cfg['nSections']):
-            for i, reg in enumerate(self.isol.sec[iSec]):
-                if reg.name == 'wake':
-                    iReg = 1
-                elif reg.name == 'lower' or reg.name == 'upper':
-                    iReg = 0
-                else:
-                    raise RuntimeError('Invalid region', reg.name)
-                loc = np.zeros(reg.getnNodes())
-                for iPoint in range(reg.getnNodes()):
-                    loc[iPoint] = reg.loc[iPoint]
-                self.isol.xxExt[iSec][i] = loc
-                self.isol.deltaStarExt[iSec][i] = np.zeros(reg.getnNodes())
-                self.isol.vtExt[iSec][i] = np.zeros(reg.getnNodes())
-                reg.setxxExt(self.isol.xxExt[iSec][i])
-                reg.setVtExt(self.isol.vtExt[iSec][i])
-                reg.setDeltaStarExt(self.isol.deltaStarExt[iSec][i])
+            for iSec in range(self.isol.cfg['nSections']):
+                for i, reg in enumerate(self.isol.sec[iSec]):
+                    if reg.name == 'wake':
+                        iReg = 1
+                    elif reg.name == 'lower' or reg.name == 'upper':
+                        iReg = 0
+                    else:
+                        raise RuntimeError('Invalid region', reg.name)
+                    loc = np.zeros(reg.getnNodes())
+                    for iPoint in range(reg.getnNodes()):
+                        loc[iPoint] = reg.loc[iPoint]
+                    self.isol.xxExt[iSec][i] = loc
+                    self.isol.deltaStarExt[iSec][i] = np.zeros(reg.getnNodes())
+                    self.isol.vtExt[iSec][i] = np.zeros(reg.getnNodes())
+                    reg.setxxExt(self.isol.xxExt[iSec][i])
+                    reg.setVtExt(self.isol.vtExt[iSec][i])
+                    reg.setDeltaStarExt(self.isol.deltaStarExt[iSec][i])
 

blast/interfaces/DSolversInterface.py

@@ -2,7 +2,6 @@ import numpy as np
 import blast
 import tbox
 from blast.interfaces.DDataStructure import Group
-import time
 
 class SolversInterface:
     def __init__(self, _vSolver, _cfg):
@@ -17,25 +16,22 @@ class SolversInterface:
         self.vBnd = [[Group('vAirfoil'+str(k) if iReg == 0 else 'vWake'+str(k), self.getnDim()) for iReg in range(2)] for k in range(self.cfg['nSections'])]
         
         # Initialize interpolator and meshes
-        initialTime = time.time()
-        self.setMesh()
+        self.getWing()
         if self.cfg['interpolator'] == 'Matching':
             from blast.interfaces.interpolators.DMatchingInterpolator import MatchingInterpolator as interp
             # Initialize viscous structures for matching computations
             for iSec in range(len(self.vBnd)):
                 for iReg, reg in enumerate(self.vBnd[iSec]):
                     reg.initStructures(self.iBnd[iReg].nPoints)
-                    self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord
-                    self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord
-                    self.vBnd[iSec][iReg].connectListNodes = self.iBnd[iReg].connectListNodes
-                    self.vBnd[iSec][iReg].connectListElems = self.iBnd[iReg].connectListElems
+                    self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord.copy()
+                    self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord.copy()
+                    self.vBnd[iSec][iReg].connectListNodes = self.iBnd[iReg].connectListNodes.copy()
+                    self.vBnd[iSec][iReg].connectListElems = self.iBnd[iReg].connectListElems.copy()
 
         elif self.cfg['interpolator'] == 'Rbf':
             from blast.interfaces.interpolators.DRbfInterpolator import RbfInterpolator as interp
-            print('RBF interpolator initialized.')
-            print('Loading viscous mesh... ', end='')
+            print('Initializing RBF interpolator.')
             self.getVWing()
-            print('done in {:.2f}s.'.format(time.time()-initialTime))
         else:
             raise RuntimeError('Incorrect interpolator specified in config.')
 
@@ -99,6 +95,102 @@ class SolversInterface:
         else:
             raise RuntimeError('Invalid direction', dir)
     
+    def updateMesh(self):
+        """Collect the mesh on blowing boundaries and update the mesh in the viscous solver.
+        This function does not update connectivity lists
+        """
+        if self.getnDim() != 2:
+            raise RuntimeError('SolverInterface::updateMesh - Only 2D meshes are supported.')
+        if self.cfg['interpolator'] != 'Matching':
+            raise RuntimeError('SolverInterface::updateMesh - Only Matching interpolator is supported.')
+
+        meshList = self.getMesh()
+
+        # Update the mesh in data structures
+        for ireg, reg in enumerate(self.iBnd):
+            for i, row in enumerate(reg.connectListNodes):
+                for idim in range(self.getnDim()):
+                    reg.nodesCoord[i, idim] = meshList[ireg][row, idim]
+
+        # Recompute elements data
+        for ireg, reg in enumerate(self.iBnd):
+            for i in range(len(reg.nodesCoord[:,0])-1):
+                for idim in range(self.getnDim()):
+                    reg.elemsCoord[i, idim] = (reg.nodesCoord[i, idim] + reg.nodesCoord[i+1, idim])/2
+        
+        # Update the mesh in the data structures
+        for iSec in range(len(self.vBnd)):
+            for iReg, reg in enumerate(self.vBnd[iSec]):
+                self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord.copy()
+                self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord.copy()
+
+        # Set mesh in the viscous solver
+        for iSec in range(self.cfg['nSections']):
+            # Compute section chord
+            xMin = np.min(self.vBnd[iSec][0].nodesCoord[:,0])
+            xMax = np.max(self.vBnd[iSec][0].nodesCoord[:,0])
+            chord = xMax - xMin
+            for reg in self.sec[iSec]:
+                if reg.name == 'wake':
+                    iReg = 1
+                elif reg.name == 'lower' or reg.name == 'upper':
+                    iReg = 0
+                else:
+                    raise RuntimeError('Invalid region', reg.name)
+                
+                xIdx = 0
+                if self.getnDim() == 2:
+                    yIdx = 1
+                    zIdx = 2
+                elif self.getnDim() == 3:
+                    yIdx = 2
+                    zIdx = 1
+                else:
+                    raise RuntimeError('Incorrect number of dimensions', self.getnDim())
+
+                nodeRows = tbox.std_vector_int()
+                if self.cfg['interpolator'] == 'Matching':
+                    for val in self.vBnd[iSec][iReg].getNodesRow(reg.name):
+                        nodeRows.push_back(int(val))
+                else:
+                    for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))):
+                        nodeRows.push_back(int(0))
+
+                connectElems = tbox.std_vector_int()
+                if self.cfg['interpolator'] == 'Matching':
+                    for val in self.vBnd[iSec][iReg].getConnectList(reg.name):
+                        connectElems.push_back(int(val))
+                else:
+                    for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))-1):
+                        connectElems.push_back(int(0))
+                
+                # x, y, z
+                xv = tbox.std_vector_double()
+                for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx):
+                    xv.push_back(val)
+                yv = tbox.std_vector_double()
+                for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx):
+                    yv.push_back(val)
+                zv = tbox.std_vector_double()
+                for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx):
+                    zv.push_back(val)
+
+                # Set mesh
+                reg.setMesh(xv,
+                            yv,
+                            zv,
+                            chord,
+                            xMin,
+                            nodeRows,
+                            connectElems)
+    
+    def updateStagnation(self):
+        """Update the stagnation points in the viscous solver.
+        """
+        for iSec in range(self.cfg['nSections']):
+            for reg in self.vBnd[iSec]:
+                reg.computeStagPoint()
+    
     def initializeViscousSolver(self):
         """Initialize the viscous sections in the viscous solver.
         """
@@ -289,4 +381,82 @@ class SolversInterface:
 
                 self.vBnd[iy][iReg].initStructures(nodesSec.shape[0])
                 self.vBnd[iy][iReg].nodesCoord = nodesSec
-                self.vBnd[iy][iReg].elemsCoord = cgSec
\ No newline at end of file
+                self.vBnd[iy][iReg].elemsCoord = cgSec
+
+    # Inviscid solver methods
+    def run()->int:
+        """Run inviscid solver
+        """
+        raise NotImplementedError('SolverInterface - run not implemented')
+    
+    def reset()->None:
+        """Reset solution
+        """
+        raise NotImplementedError('SolverInterface - reset not implemented')
+    
+    def save(self)->None:
+        """Save solution
+        """
+        raise NotImplementedError('SolverInterface - save not implemented')
+    
+    def writeCp(self)->None:
+        """Write pressure coefficient
+        """
+        raise NotImplementedError('SolverInterface - writeCp not implemented')
+    
+    # Inviscid solver properties
+    def getVerbose(self)->int:
+        """Get verbose
+        """
+        raise NotImplementedError('SolverInterface - getVerbose not implemented')
+
+    # Flow properties and aerodynamic coefficients
+    def getAoA(self)->float:
+        """Get angle of attack
+        """
+        raise NotImplementedError('SolverInterface - getAoA not implemented')
+
+    def getCl(self)->float:
+        """Get lift coefficient
+        """
+        raise NotImplementedError('SolverInterface - getCl not implemented')
+    
+    def getCd(self)->float:
+        """Get inviscid drag coefficient
+        """
+        raise NotImplementedError('SolverInterface - getCd not implemented')
+    
+    def getCm(self)->float:
+        """Get moment coefficient
+        """
+        raise NotImplementedError('SolverInterface - getCm not implemented')
+    
+    def getnDim(self)->int:
+        """Get number of dimensions
+        """
+        raise NotImplementedError('SolverInterface - getnDim not implemented')
+    
+    def getWing(self)->None:
+        """Get wing mesh
+        """
+        raise NotImplementedError('SolverInterface - getWing not implemented')
+    
+    def getMesh(self)->list:
+        """Get mesh
+        """
+        raise NotImplementedError('SolverInterface - getMesh not implemented')
+    
+    def getInviscidBC(self)->None:
+        """Get inviscid boundary condition
+        """
+        raise NotImplementedError('SolverInterface - getInviscidBC not implemented')
+    
+    def setBlowingVelocity(self)->None:
+        """Set blowing velocity
+        """
+        raise NotImplementedError('SolverInterface - setBlowingVelocity not implemented')
+    
+    def getnNodesDomain(self)->int:
+        """Get number of nodes in domain
+        """
+        raise NotImplementedError('SolverInterface - getnNodesDomain not implemented')

blast/interfaces/dart/DartInterface.py

@@ -76,7 +76,8 @@ class DartInterface(SolversInterface):
                     for j in range(size):
                         data[i,3+j] = vData[vElems[i%vElems.size()].nodes[0].row][j]
                 i += 1
-            print('writing data file ' + 'Cp' +sfx + '.dat')
+            if self.getVerbose() > 1:
+                print('writing data file ' + 'Cp' +sfx + '.dat')
             np.savetxt('Cp'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
         
         elif self.getnDim() == 3:
@@ -90,7 +91,8 @@ class DartInterface(SolversInterface):
                 np.savetxt(self.iobj['msh'].name+'_Cp_surface'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
                 import modules.dartflo.dart.utils as invUtils
                 if self.cfg['Format'] == 'vtk':
-                    print('Saving Cp files in vtk format. Msh {:s}, Tag {:.0f}'.format(self.iobj['msh'].name, self.cfg['saveTag']))
+                    if self.getVerbose() > 1:
+                        print('Saving Cp files in vtk format. Msh {:s}, Tag {:.0f}'.format(self.iobj['msh'].name, self.cfg['saveTag']))
                     invUtils.write_slices(self.iobj['msh'].name, self.cfg['writeSections'], self.cfg['saveTag'], sfx=sfx)
     
     def save(self, sfx='inviscid'):
@@ -154,17 +156,27 @@ class DartInterface(SolversInterface):
             for i, blw in enumerate(self.iBnd[n].blowingVel):
                 self.blw[n].setU(i, blw)
 
-    def setMesh(self):
+    def getWing(self):
         if self.getnDim() == 2:
-            self.getAirfoil()
+            self.__getWing2D()
         elif self.getnDim() == 3:
-            self.getWing()
+            self.__getWing3D()
         else:
             raise RuntimeError('dartInterface::Could not resolve how to set\
             the mesh. Either ''nDim'' is incorrect or ''msh'' was not set for\
             3D computations')
+    
+    def getMesh(self):
+        """Return the mesh on blowing boundaries
+        """
+        meshList = [np.zeros((b.nodes.size(), self.getnDim())) for b in self.blw]
+        for ib, b in enumerate(self.blw):
+            for inod, n in enumerate(b.nodes):
+                for idim in range(self.getnDim()):
+                    meshList[ib][inod, idim] = n.pos[idim]
+        return meshList
 
-    def getAirfoil(self):
+    def __getWing2D(self):
         """Create data structure for information transfer.
         """
         self.iBnd[0].initStructures(self.blw[0].nodes.size())
@@ -332,7 +344,7 @@ class DartInterface(SolversInterface):
             elemCoordW[i,3] = i
         self.iBnd[1].elemsCoord = elemCoordW
 
-    def getWing(self):
+    def __getWing3D(self):
         """Obtain the nodes' location and row and cg of all elements.
         """
 

blast/utils.py

@@ -44,7 +44,7 @@ def initBL(Re, Minf, CFL0, nSections, xtrF = [None, None], span=0, verb=None):
         if xtrF[i] is None:
             xtrF[i] = -1
         if xtrF[i] != -1 and not(0<= xtrF[i] <= 1):
-            raise RuntimeError('Incorrect forced transition location') 
+            raise RuntimeError('Incorrect forced transition location')
 
     solver = blast.Driver(Re, Minf, CFL0, nSections, xtrF_top=xtrF[0], xtrF_bot=xtrF[1], _span=span, _verbose=verbose)
     return solver
@@ -70,7 +70,7 @@ def initBlast(iconfig, vconfig, iSolver='DART'):
     if 'sfx' not in vconfig:
         vconfig['sfx'] = ''
     # Viscous solver
-    vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], xtrF=vconfig['xtrF'])
+    vSolver = initBL(vconfig['Re'], vconfig['Minf'], vconfig['CFL0'], vconfig['nSections'], xtrF=vconfig['xtrF'], verb=vconfig['Verb'])
 
     # Inviscid solver
     if iSolver == 'DART':