From 44904a461bb8b41c5951fb99f115003bb687778d Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Sun, 29 Sep 2024 10:36:26 +0200
Subject: [PATCH] (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 | 47 +++---
blast/interfaces/DSolversInterface.py | 192 +++++++++++++++++++++++--
blast/interfaces/dart/DartInterface.py | 26 +++-
blast/utils.py | 4 +-
4 files changed, 227 insertions(+), 42 deletions(-)
diff --git a/blast/coupler.py b/blast/coupler.py
index 9c3d903..2878d67 100644
--- a/blast/coupler.py
+++ b/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])
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index 50eae4b..b7542bc 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/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')
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index 0334541..7913ab5 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/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.
"""
diff --git a/blast/utils.py b/blast/utils.py
index 4da3ce6..db3d523 100644
--- a/blast/utils.py
+++ b/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':
--
GitLab