From 1f95c6a3ca0c95a7d9187b480a30a5f5d525d7fe Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Thu, 27 Feb 2025 18:38:45 +0100
Subject: [PATCH] (feat) Modify interpolators for new data structures
---
.../interpolators/blMatchingInterpolator.py | 16 +-
.../interpolators/blRbfInterpolator.py | 275 +++++++++---------
2 files changed, 151 insertions(+), 140 deletions(-)
diff --git a/blast/interfaces/interpolators/blMatchingInterpolator.py b/blast/interfaces/interpolators/blMatchingInterpolator.py
index 5f90c89..4fd0ada 100644
--- a/blast/interfaces/interpolators/blMatchingInterpolator.py
+++ b/blast/interfaces/interpolators/blMatchingInterpolator.py
@@ -61,15 +61,11 @@ class MatchingInterpolator(Interpolator):
"""
if self.ndim == 2:
for ibody in range(len(iDict)):
- for ireg in range(len(iDict[ibody])):
- vDict[ibody][0][ireg].updateVars(iDict[ibody][ireg].V, iDict[ibody][ireg].M, iDict[ibody][ireg].Rho)
- elif self.ndim == 3:
- for ibody in range(len(iDict)):
- for isec, ysec in enumerate(self._sections):
- for ireg in range(2):
- print(iDict[ibody][ireg].nodesCoord[iDict[ireg].nodesCoord[:,1] == ysec])
- print(iDict[ibody][ireg].V[iDict[ireg].nodesCoord[:,1] == ysec])
- vDict[ibody][isec][ireg].updateVars(iDict[ibody][ireg].V[iDict[ibody][ireg].nodesCoord[:,1] == ysec], iDict[ibody][ireg].Rho[iDict[ibody][ireg].nodesCoord[:,1] == ysec])
+ for i, ireg in enumerate(iDict[ibody]):
+ vreg = vDict[ibody][0][i]
+ vreg.updateVariables(ireg.getMach('all'), ireg.getVelocity('all'), ireg.getDensity('all'))
+ else:
+ raise RuntimeError('Incorrect number of dimensions', self.ndim)
def viscousToInviscid(self, iDict, vDict):
"""
@@ -85,6 +81,6 @@ class MatchingInterpolator(Interpolator):
if self.ndim == 2:
for ibody in range(len(iDict)):
for iReg in range(len(iDict[ibody])):
- iDict[ibody][iReg].blowingVel = vDict[ibody][0][iReg].blowingVel
+ iDict[ibody][iReg].setBlowingVelocity(vDict[ibody][0][iReg].getBlowingVelocity())
else:
raise RuntimeError('Incorrect number of dimensions', self.ndim)
\ No newline at end of file
diff --git a/blast/interfaces/interpolators/blRbfInterpolator.py b/blast/interfaces/interpolators/blRbfInterpolator.py
index 35b8535..5f2d0a6 100644
--- a/blast/interfaces/interpolators/blRbfInterpolator.py
+++ b/blast/interfaces/interpolators/blRbfInterpolator.py
@@ -34,157 +34,172 @@ class RbfInterpolator(Interpolator):
def inviscidToViscous(self, iDict, vDict):
self.tms['tot'].start()
- for ibody in range(len(iDict)):
- if iDict[ibody][0].sidesIdentified and iDict[ibody][0].type == 'wing':
+ for ibody, inviscid_regions in enumerate(iDict):
+ inviscid_airfoil = inviscid_regions[0]
+ inviscid_wake = inviscid_regions[1]
+
+ if inviscid_airfoil.sidesIdentified and inviscid_airfoil.getType() == 'wing':
self.tms['getIdict'].start()
- xi_up = iDict[ibody][0].getUpperNodes()
- Mi_up = iDict[ibody][0].getUpperMach()
- rhoi_up = iDict[ibody][0].getUpperDensity()
- vi_up = iDict[ibody][0].getUpperVelocity()
-
- # lower side
- xi_lw = iDict[ibody][0].getLowerNodes()
- Mi_lw = iDict[ibody][0].getLowerMach()
- rhoi_lw = iDict[ibody][0].getLowerDensity()
- vi_lw = iDict[ibody][0].getLowerVelocity()
-
- # wake
- xi_wk = iDict[ibody][1].nodesCoord[:, :self.ndim]
- Mi_wk = iDict[ibody][1].M
- rhoi_wk = iDict[ibody][1].Rho
- vi_wk = iDict[ibody][1].V[:,:self.ndim]
- self.tms['getIdict'].stop()
+ # Extract inviscid data for interpolation
+ xi_up, Mi_up, rhoi_up, vi_up = (
+ inviscid_airfoil.getNodes('upper_fromLE')[:, :self.ndim],
+ inviscid_airfoil.getMach('upper_fromLE'),
+ inviscid_airfoil.getDensity('upper_fromLE'),
+ inviscid_airfoil.getVelocity('upper_fromLE')[:, :self.ndim],
+ )
+
+ xi_lw, Mi_lw, rhoi_lw, vi_lw = (
+ inviscid_airfoil.getNodes('lower_fromLE')[:, :self.ndim],
+ inviscid_airfoil.getMach('lower_fromLE'),
+ inviscid_airfoil.getDensity('lower_fromLE'),
+ inviscid_airfoil.getVelocity('lower_fromLE')[:, :self.ndim],
+ )
+
+ xi_wk, Mi_wk, rhoi_wk, vi_wk = (
+ inviscid_wake.getNodes('wake')[:, :self.ndim],
+ inviscid_wake.getMach('wake'),
+ inviscid_wake.getDensity('wake'),
+ inviscid_wake.getVelocity('wake')[:, :self.ndim],
+ )
+
+ self.tms['getIdict'].stop()
self.tms['getVdict'].start()
- # Viscous side
- xv_up = np.empty((0, self.ndim), dtype=float)
- xv_lw = np.empty((0, self.ndim), dtype=float)
- xv_wk = np.empty((0, self.ndim), dtype=float)
- for iSec in range(len(vDict[ibody])):
- xv_up = np.row_stack((xv_up, vDict[ibody][iSec][0].getUpperNodes()))
- xv_lw = np.row_stack((xv_lw, vDict[ibody][iSec][0].getLowerNodes()))
- xv_wk = np.row_stack((xv_wk, vDict[ibody][iSec][1].nodesCoord[:,:self.ndim]))
+
+ # Get viscous node
+ xv_up = np.vstack([sec[0].getNodes('upper_fromLE')[:, :self.ndim] for sec in vDict[ibody]])
+ xv_lw = np.vstack([sec[0].getNodes('lower_fromLE')[:, :self.ndim] for sec in vDict[ibody]])
+ xv_wk = np.vstack([sec[1].getNodes('wake')[:, :self.ndim] for sec in vDict[ibody]])
+
self.tms['getVdict'].stop()
- # Interpolation
self.tms['rbf'].start()
- Mv_up = self.__rbfToSection(xi_up, Mi_up, xv_up)
- rhov_up = self.__rbfToSection(xi_up, rhoi_up, xv_up)
- vv_up = np.zeros((xv_up.shape[0], self.ndim))
- for idim in range(self.ndim):
- vv_up[:,idim] = self.__rbfToSection(xi_up, vi_up[:,idim], xv_up)
-
- Mv_lw = self.__rbfToSection(xi_lw, Mi_lw, xv_lw)
- rhov_lw = self.__rbfToSection(xi_lw, rhoi_lw, xv_lw)
- vv_lw = np.zeros((xv_lw.shape[0], self.ndim))
- for idim in range(self.ndim):
- vv_lw[:,idim] = self.__rbfToSection(xi_lw, vi_lw[:,idim], xv_lw)
-
- Mv_wk = self.__rbfToSection(xi_wk, Mi_wk, xv_wk)
- rhov_wk = self.__rbfToSection(xi_wk, rhoi_wk, xv_wk)
- vv_wk = np.zeros((xv_wk.shape[0], self.ndim))
- for idim in range(self.ndim):
- vv_wk[:,idim] = self.__rbfToSection(xi_wk, vi_wk[:,idim], xv_wk)
- self.tms['rbf'].stop()
- # The viscous structure is sorted in selig format.
- # Here, upper and lower side are 'opposite' (upper is from TE to LE
- # while lower is from LE to TE). It is important that it always stays
- # this way for this concatenation operation.
+ # Perform vectorized RBF interpolation
+ Mv_up, rhov_up = self.__rbfToSection(xi_up, Mi_up, xv_up), self.__rbfToSection(xi_up, rhoi_up, xv_up)
+ Mv_lw, rhov_lw = self.__rbfToSection(xi_lw, Mi_lw, xv_lw), self.__rbfToSection(xi_lw, rhoi_lw, xv_lw)
+ Mv_wk, rhov_wk = self.__rbfToSection(xi_wk, Mi_wk, xv_wk), self.__rbfToSection(xi_wk, rhoi_wk, xv_wk)
+
+ vv_up = np.column_stack([self.__rbfToSection(xi_up, vi_up[:, dim], xv_up) for dim in range(self.ndim)])
+ vv_lw = np.column_stack([self.__rbfToSection(xi_lw, vi_lw[:, dim], xv_lw) for dim in range(self.ndim)])
+ vv_wk = np.column_stack([self.__rbfToSection(xi_wk, vi_wk[:, dim], xv_wk) for dim in range(self.ndim)])
+
+ self.tms['rbf'].stop()
self.tms['setting'].start()
- cpt_up = 0
- cpt_lw = 0
- cpt_wk = 0
- from matplotlib import pyplot as plt
+
+ # Assign values to each viscous section
+ cpt_up, cpt_lw, cpt_wk = 0, 0, 0
for sec in vDict[ibody]:
- M_airfoil = np.empty(sec[0].nPoints, dtype=float)
- M_wake = np.empty(sec[1].nPoints, dtype=float)
- rho_airfoil = np.empty(sec[0].nPoints, dtype=float)
- rho_wake = np.empty(sec[1].nPoints, dtype=float)
- v_airfoil = np.zeros((sec[0].nPoints, self.ndim), dtype=float)
- v_wake = np.zeros((sec[1].nPoints, self.ndim), dtype=float)
-
- # Mach number
- M_airfoil = np.concatenate((Mv_up[cpt_up:cpt_up+sec[0].getUpperNodes().shape[0]],
- Mv_lw[cpt_lw:cpt_lw+sec[0].getLowerNodes().shape[0]]))
- M_wake = Mv_wk[cpt_wk:cpt_wk+sec[1].nPoints]
- # Density
- rho_airfoil = np.concatenate((rhov_up[cpt_up:cpt_up+sec[0].getUpperNodes().shape[0]],
- rhov_lw[cpt_lw:cpt_lw+sec[0].getLowerNodes().shape[0]]))
- rho_wake = rhov_wk[cpt_wk:cpt_wk+sec[1].nPoints]
-
- # Velocity
- for idim in range(self.ndim):
- v_airfoil[:,idim] = np.concatenate((vv_up[cpt_up:cpt_up+sec[0].getUpperNodes().shape[0],idim],
- vv_lw[cpt_lw:cpt_lw+sec[0].getLowerNodes().shape[0],idim]))
- v_wake[:,idim] = vv_wk[cpt_wk:cpt_wk+sec[1].nPoints,idim]
- cpt_up += sec[0].getUpperNodes().shape[0]
- cpt_lw += sec[0].getLowerNodes().shape[0]
- cpt_wk += sec[1].nPoints
-
- sec[0].updateVars(v_airfoil, M_airfoil, rho_airfoil)
- sec[1].updateVars(v_wake, M_wake, rho_wake)
+ n_up = sec[0].getNodes('upper_fromLE').shape[0]
+ n_lw = sec[0].getNodes('lower_fromLE').shape[0]
+ n_wk = sec[1].getnNodes()
+
+ # Concatenate Mach number and density
+ M_airfoil = np.concatenate((Mv_up[cpt_up:cpt_up + n_up], Mv_lw[cpt_lw:cpt_lw + n_lw]))
+ rho_airfoil = np.concatenate((rhov_up[cpt_up:cpt_up + n_up], rhov_lw[cpt_lw:cpt_lw + n_lw]))
+ M_wake, rho_wake = Mv_wk[cpt_wk:cpt_wk + n_wk], rhov_wk[cpt_wk:cpt_wk + n_wk]
+
+ # Concatenate velocity
+ v_airfoil = np.column_stack([
+ np.concatenate((vv_up[cpt_up:cpt_up + n_up, dim], vv_lw[cpt_lw:cpt_lw + n_lw, dim]))
+ for dim in range(self.ndim)
+ ])
+ v_wake = vv_wk[cpt_wk:cpt_wk + n_wk]
+
+ # Update the variables
+ sec[0].updateVariables(M_airfoil, v_airfoil, rho_airfoil)
+ sec[1].updateVariables(M_wake, v_wake, rho_wake)
+
+ # Update counters
+ cpt_up += n_up
+ cpt_lw += n_lw
+ cpt_wk += n_wk
+
self.tms['setting'].stop()
+
else:
- for iSec in range(len(vDict[ibody])):
- for iReg, reg in enumerate(vDict[ibody][iSec]):
- v = np.zeros((reg.nodesCoord.shape[0], 3))
- M = np.zeros(reg.nodesCoord.shape[0])
- rho = np.zeros(reg.nodesCoord.shape[0])
- for iDim in range(self.ndim):
- v[:,iDim] = self.__rbfToSection(iDict[ibody][iReg].nodesCoord[:,:(self.ndim if iDict[ibody][iReg].name == 'iWing' else 1)], iDict[ibody][iReg].V[:,iDim], reg.nodesCoord[:,:(self.ndim if 'vAirfoil' in reg.name else 1)])
- M = self.__rbfToSection(iDict[ibody][iReg].nodesCoord[:,:(self.ndim if iDict[ibody][iReg].name == 'iWing' else 1)], iDict[ibody][iReg].M, reg.nodesCoord[:,:(self.ndim if 'vAirfoil' in reg.name else 1)])
- rho = self.__rbfToSection(iDict[ibody][iReg].nodesCoord[:,:(self.ndim if iDict[ibody][iReg].name == 'iWing' else 1)], iDict[ibody][iReg].Rho, reg.nodesCoord[:,:(self.ndim if 'vAirfoil' in reg.name else 1)])
- vDict[ibody][iSec][iReg].updateVars(v, M, rho)
+ # General case
+ for iSec, sections in enumerate(vDict[ibody]):
+ for iReg, reg in enumerate(sections):
+ n_nodes = reg.getnNodes()
+
+ slice_dim = self.ndim if 'vAirfoil' in reg._name else 1
+ xi = iDict[ibody][iReg].getNodes('all')[:, :slice_dim]
+ xv = reg.getNodes('all')[:, :slice_dim]
+
+ # Perform RBF interpolation
+ v = np.column_stack([
+ self.__rbfToSection(xi, iDict[ibody][iReg].getVelocity(reg='all')[:,dim], xv)
+ for dim in range(self.ndim)
+ ])
+ M = self.__rbfToSection(xi, iDict[ibody][iReg].getMach('all'), xv)
+ rho = self.__rbfToSection(xi, iDict[ibody][iReg].getDensity('all'), xv)
+
+ # Update the variables
+ reg.updateVariables(M, v, rho)
+ self.tms['tot'].stop()
def viscousToInviscid(self, iDict, vDict):
for ibody in range(len(iDict)):
if iDict[ibody][0].sidesIdentified:
- # Viscous side
- xv_up = np.zeros((0, self.ndim))
- xv_lw = np.zeros((0, self.ndim))
- xv_wk = np.zeros((0, self.ndim))
- blwv_up = np.zeros(0)
- blwv_lw = np.zeros(0)
- blwv_wk = np.zeros(0)
- for iSec, sec in enumerate(vDict[ibody]):
- xv_up = np.row_stack((xv_up, sec[0].getUpperElems()))
- xv_lw = np.row_stack((xv_lw, sec[0].getLowerElems()))
- xv_wk = np.row_stack((xv_wk, sec[1].elemsCoord[:,:self.ndim]))
- blwv_up = np.concatenate((blwv_up, sec[0].getUpperBlowingVel()))
- blwv_lw = np.concatenate((blwv_lw, sec[0].getLowerBlowingVel()))
- blwv_wk = np.concatenate((blwv_wk, sec[1].blowingVel))
-
- # Inviscid side
- xi_up = iDict[ibody][0].getUpperElems()
- xi_lw = iDict[ibody][0].getLowerElems()
- xi_wk = iDict[ibody][1].elemsCoord[:,:self.ndim]
-
- blwi_up = self.__rbfToSection(xv_up, blwv_up, xi_up[:,:self.ndim])
- blwi_lw = self.__rbfToSection(xv_lw, blwv_lw, xi_lw[:,:self.ndim])
- blwi_wk = self.__rbfToSection(xv_wk, blwv_wk, xi_wk[:,:self.ndim])
-
- iDict[ibody][0].blowingVel = np.concatenate((blwi_up, blwi_lw))
- iDict[ibody][1].blowingVel = blwi_wk
+ # Extract elements and blowing velocity for the viscous side
+ xv_up, xv_lw, xv_wk = [], [], []
+ blwv_up, blwv_lw, blwv_wk = [], [], []
+
+ for sec in vDict[ibody]:
+ xv_up.append(sec[0].getElms('upper_fromLE')[:, :self.ndim])
+ xv_lw.append(sec[0].getElms('lower_fromLE')[:, :self.ndim])
+ xv_wk.append(sec[1].getElms('wake')[:, :self.ndim])
+
+ blwv_up.append(sec[0].getBlowingVelocity('upper_fromLE'))
+ blwv_lw.append(sec[0].getBlowingVelocity('lower_fromLE'))
+ blwv_wk.append(sec[1].getBlowingVelocity('wake'))
+
+ xv_up, xv_lw, xv_wk = map(np.vstack, [xv_up, xv_lw, xv_wk])
+ blwv_up, blwv_lw, blwv_wk = map(np.concatenate, [blwv_up, blwv_lw, blwv_wk])
+
+ # Extract elements for the inviscid side
+ xi_up = iDict[ibody][0].getElms('upper_fromLE')[:, :self.ndim]
+ xi_lw = iDict[ibody][0].getElms('lower_fromLE')[:, :self.ndim]
+ xi_wk = iDict[ibody][1].getElms('wake')[:, :self.ndim]
+
+ # Perform RBF interpolation
+ blwi_up = self.__rbfToSection(xv_up, blwv_up, xi_up)
+ blwi_lw = self.__rbfToSection(xv_lw, blwv_lw, xi_lw)
+ blwi_wk = self.__rbfToSection(xv_wk, blwv_wk, xi_wk)
+
+ # Set interpolated blowing velocity for inviscid side
+ iDict[ibody][0].setBlowingVelocity(np.concatenate((blwi_up, blwi_lw)))
+ iDict[ibody][1].setBlowingVelocity(blwi_wk)
+
else:
if self.ndim == 2:
raise RuntimeError('2D not implemented')
- for iReg, reg in enumerate(iDict[ibody]):
- iDict[iReg].blowingVel = self.__rbfToSection(vDict[ibody][0][iReg].elemsCoord[:,:(self.ndim-1 if 'vAirfoil' in reg.name else 1)], vDict[0][iReg].blowingVel, reg.elemsCoord[:,:(self.ndim-1 if reg.name == 'iWing' else 1)])
+
elif self.ndim == 3:
for iReg, reg in enumerate(iDict[ibody]):
- viscElemsCoord = np.zeros((0,3))
- viscBlowing = np.zeros(0)
- for iSec, sec in enumerate(vDict[ibody]):
- viscElemsCoord = np.row_stack((viscElemsCoord, sec[iReg].elemsCoord[:,:3]))
- viscBlowing = np.concatenate((viscBlowing, sec[iReg].blowingVel))
- if len(self._sym) > 0:
+ # Collect elements and blowing velocity
+ viscElemsCoord, viscBlowing = [], []
+ for sec in vDict[ibody]:
+ viscElemsCoord.append(sec[iReg].getElms('all'))
+ viscBlowing.append(sec[iReg].getBlowingVelocity('all'))
+
+ # Apply symmetry transformations if needed
+ if self._sym:
for s in self._sym:
- dummy = viscElemsCoord.copy()
- dummy[:,1] = (dummy[:,1] - s)*(-1) + s
- viscElemsCoord = np.row_stack((viscElemsCoord, dummy))
- viscBlowing = np.concatenate((viscBlowing, viscBlowing))
- reg.blowingVel = self.__rbfToSection(viscElemsCoord, viscBlowing, reg.elemsCoord[:,:3])
+ mirrored = np.copy(np.vstack(viscElemsCoord))
+ mirrored[:, 1] = (mirrored[:, 1] - s) * (-1) + s
+ viscElemsCoord.append(mirrored)
+ viscBlowing.append(np.concatenate(viscBlowing))
+
+ # Convert to numpy arrays efficiently
+ viscElemsCoord = np.vstack(viscElemsCoord)
+ viscBlowing = np.concatenate(viscBlowing)
+
+ # Apply interpolation and update
+ blw = self.__rbfToSection(viscElemsCoord, viscBlowing, reg.getElms('all'))
+ reg.setBlowingVelocity(blw)
+
else:
raise RuntimeError('Incorrect number of dimensions', self.cfg['nDim'])
--
GitLab