From aab1225e930bf221188c690324893ee8facf12dc Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Thu, 27 Feb 2025 18:39:57 +0100
Subject: [PATCH] (feat) Modify solver interface for new data structures and
new restart capability
---
blast/interfaces/blSolversInterface.py | 355 ++++++++++++++++---------
blast/tests/t_restart_2D.py | 168 ++++++++++++
2 files changed, 403 insertions(+), 120 deletions(-)
create mode 100644 blast/tests/t_restart_2D.py
diff --git a/blast/interfaces/blSolversInterface.py b/blast/interfaces/blSolversInterface.py
index 5f3c4cf..ced73ad 100644
--- a/blast/interfaces/blSolversInterface.py
+++ b/blast/interfaces/blSolversInterface.py
@@ -22,7 +22,9 @@ import numpy as np
import blast
import tbox
-from blast.interfaces.blDataStructure import Group
+# from blast.interfaces.blDataStructure import Group
+from blast.interfaces.blDataInviscid import InviscidData
+from blast.interfaces.blDataViscous import ViscousData
from scipy.interpolate import griddata
from scipy.interpolate import interp1d
@@ -90,17 +92,21 @@ class SolversInterface:
print(f'-Creating sections using {self.cfg["genSections"]}', end='...')
self._getSections()
print('done')
-
- # External quantities
+
self.deltaStarExt = None
- self.xxExt = None
- self.vtExt = None
+ self.xxExt = None
+ self.vtExt = None
- self.vSolver = None
- # Flags
self.vinit = False
+ # Flags
+ self.vSolver = None
self.hasvsol = False
+ # Check for restart file
+ if "restart_solution" in self.cfg and self.cfg["restart_solution"] == True:
+ print('Loading restart file', end='...')
+ self._load_restart(self.cfg["restart_inputName"], self.vBnd)
+ print('done')
self._printSetup()
########################
@@ -125,33 +131,41 @@ class SolversInterface:
couplIter: int
Current iteration of the coupling loop.
"""
- if couplIter != 0:
- for ibody, body in enumerate(self.vSolver.bodies):
- for isec, section in enumerate(body.sections):
- for iReg in range(len(self.vBnd[ibody][isec])):
- pystruct = self.vBnd[ibody][isec][iReg]
-
- if 'vWake' in pystruct.name or pystruct.type == 'fuselage':
- wakeId = section.getRegionIdByName("wake")
- if pystruct.blowingVel.shape[0] != len(section.regions[wakeId].blowingVelocity):
- raise RuntimeError('Blowing velocity size mismatch')
- pystruct.blowingVel = np.asarray(section.regions[wakeId].blowingVelocity)
- elif 'vAirfoil' in pystruct.name:
- upperId = section.getRegionIdByName("upper")
- lowerId = section.getRegionIdByName("lower")
- if pystruct.blowingVel.shape[0] != len(np.concatenate((np.flip(np.asarray(section.regions[upperId].blowingVelocity)),
- np.asarray(section.regions[lowerId].blowingVelocity)))):
- raise RuntimeError('Blowing velocity size mismatch')
- pystruct.blowingVel = np.concatenate((np.flip(np.asarray(section.regions[upperId].blowingVelocity)),
- np.asarray(section.regions[lowerId].blowingVelocity)))
+ if couplIter == 0:
+ return
+ for ibody, body in enumerate(self.vSolver.bodies):
+ for isec, section in enumerate(body.sections):
+ for iReg, pystruct in enumerate(self.vBnd[ibody][isec]):
- for i, reg in enumerate(section.regions):
- loc = np.zeros(reg.getnNodes())
- for inod in range(reg.getnNodes()):
- loc[inod] = reg.nodes[inod].xi
- self.xxExt[ibody][isec][i] = np.asarray(loc)
- self.deltaStarExt[ibody][isec][i] = np.asarray(reg.deltaStar)
- self.vtExt[ibody][isec][i] = np.asarray(reg.vt)
+ name, ptype = pystruct.getName(), pystruct.getType()
+ blowingVelocity = pystruct.getBlowingVelocity()
+
+ if 'vWake' in name or ptype == 'fuselage':
+ wakeId = section.getRegionIdByName("wake")
+ wakeVel = np.asarray(section.regions[wakeId].blowingVelocity)
+
+ if blowingVelocity.shape[0] != wakeVel.shape[0]:
+ raise RuntimeError('Blowing velocity size mismatch')
+
+ pystruct.setBlowingVelocity(wakeVel)
+
+ elif 'vAirfoil' in name:
+ upperId, lowerId = section.getRegionIdByName("upper"), section.getRegionIdByName("lower")
+ upperVel = np.asarray(section.regions[upperId].blowingVelocity)
+ lowerVel = np.asarray(section.regions[lowerId].blowingVelocity)
+
+ combinedVel = np.concatenate((upperVel[::-1], lowerVel))
+
+ if blowingVelocity.shape[0] != combinedVel.shape[0]:
+ raise RuntimeError('Blowing velocity size mismatch')
+
+ pystruct.setBlowingVelocity(combinedVel)
+
+ for i, reg in enumerate(section.regions):
+ nodes_xi = np.array([node.xi for node in reg.nodes], dtype=np.float64)
+ self.xxExt[ibody][isec][i] = nodes_xi
+ self.deltaStarExt[ibody][isec][i] = np.array(reg.deltaStar, dtype=np.float64)
+ self.vtExt[ibody][isec][i] = np.array(reg.vt, dtype=np.float64)
def interpolate(self, dir:str):
"""Interpolate solution between inviscid and viscous meshes.
@@ -162,34 +176,46 @@ class SolversInterface:
'i2v': inviscid to viscous.
'v2i': viscous to inviscid.
"""
- if dir == 'i2v':
- self.interpolator.inviscidToViscous(self.iBnd, self.vBnd)
- elif dir == 'v2i':
- self.interpolator.viscousToInviscid(self.iBnd, self.vBnd)
- else:
- raise RuntimeError('Invalid direction', dir)
+ interpolation_methods = {
+ "i2v": self.interpolator.inviscidToViscous,
+ "v2i": self.interpolator.viscousToInviscid
+ }
+ try:
+ interpolation_methods[dir](self.iBnd, self.vBnd)
+ except KeyError:
+ raise RuntimeError(f"Invalid direction: {dir}")
def setViscousSolver(self):
"""Sets the actual conditions in vBnd to the viscous solver.
"""
+ region_map = {
+ "wake": 1,
+ "upper": 0,
+ "lower": 0,
+ "fuselage": 0
+ }
+
for ibody, body in enumerate(self.vSolver.bodies):
for isec, section in enumerate(body.sections):
for i, reg in enumerate(section.regions):
- if reg.getName() in ['wake']:
- iReg = 1
- elif reg.getName() in ['upper', 'lower', 'fuselage']:
- iReg = 0
- else:
- raise RuntimeError('Invalid region', reg.getName())
+ reg_name = reg.getName()
+
+ if reg_name not in region_map:
+ raise RuntimeError(f"Invalid region: {reg_name}")
+
+ iReg = region_map[reg_name]
pystruct = self.vBnd[ibody][isec][iReg]
- Mach = pystruct.getMach(reg.getName())
- Rho = pystruct.getRho(reg.getName())
- Vt = pystruct.getVt(reg.getName())
- for inod in range(reg.getnNodes()):
- reg.Me[inod] = Mach[inod]
- reg.rhoe[inod] = Rho[inod]
- reg.vt[inod] = Vt[inod]
+ Mach, Rho, Vt = (
+ pystruct.getMach(reg_name),
+ pystruct.getDensity(reg_name),
+ pystruct.getVt(reg_name),
+ )
+
+ for inod, (mach, rho, vt) in enumerate(zip(Mach, Rho, Vt)):
+ reg.Me[inod] = mach
+ reg.rhoe[inod] = rho
+ reg.vt[inod] = vt
reg.deltaStarExt[inod] = self.deltaStarExt[ibody][isec][i][inod]
reg.nodes[inod].xiExt = self.xxExt[ibody][isec][i][inod]
reg.vtExt[inod] = self.vtExt[ibody][isec][i][inod]
@@ -210,9 +236,8 @@ class SolversInterface:
if len(body.sections) != len(self.vBnd[ibody]):
raise RuntimeError(f'Expected {len(self.vBnd[ibody])} sections in viscous solver, got {len(body.sections)}')
for isec, section in enumerate(body.sections):
-
- # Check that if self.vBnd[ibody][isec][0].type is 'wing', section mus have upper and lower
- if self.vBnd[ibody][isec][0].type == 'wing':
+ # Check that if self.vBnd[ibody][isec][0].getType() is 'wing', section mus have upper and lower
+ if self.vBnd[ibody][isec][0].getType() == 'wing':
hasupper = False
haslower = False
for reg in section.regions:
@@ -221,10 +246,10 @@ class SolversInterface:
elif reg.getName() == "lower":
haslower = True
if not hasupper or not haslower:
- raise RuntimeError(f'Expected upper and lower sections for type {self.vBnd[ibody][isec][0].type}')
+ raise RuntimeError(f'Expected upper and lower sections for type {self.vBnd[ibody][isec][0].getType()}')
# Check that the fuselage has no wake
- elif self.vBnd[ibody][isec][0].type == 'fuselage':
+ elif self.vBnd[ibody][isec][0].getType() == 'fuselage':
for reg in section.regions:
if reg.getName() == "wake":
raise RuntimeError('Fuselage should not have wake')
@@ -237,7 +262,7 @@ class SolversInterface:
bodyhaswake = body.hasWake()
pystructhaswake = False
for iReg, reg in enumerate(self.vBnd[ibody][isec]):
- if reg.type == 'wake':
+ if reg.getType() == 'wake':
pystructhaswake = True
break
if bodyhaswake != pystructhaswake:
@@ -255,12 +280,36 @@ class SolversInterface:
# Set mesh
self._setViscousMesh()
- self.xxExt = [[[np.zeros(0) for iReg in range(len(self.vSolver.bodies[ibody].sections[isec].regions))]\
- for isec in range(len(self.vSolver.bodies[ibody].sections))] for ibody in range(self.getnBodies())]
- self.deltaStarExt = [[[np.zeros(0) for iReg in range(len(self.vSolver.bodies[ibody].sections[isec].regions))]\
- for isec in range(len(self.vBnd[ibody]))] for ibody in range(self.getnBodies())]
- self.vtExt = [[[np.zeros(0) for iReg in range(len(self.vSolver.bodies[ibody].sections[isec].regions))]\
- for isec in range(len(self.vSolver.bodies[ibody].sections))] for ibody in range(self.getnBodies())]
+ # Initialize interaction law quantities
+ if not self.vinit:
+ self._initializeInteractionLaw()
+ print('Viscous solver initialized')
+ return 0
+
+ def _initializeInteractionLaw(self):
+ self.deltaStarExt = []
+ self.xxExt = []
+ self.vtExt = []
+ for ibody, body in enumerate(self.vBnd):
+ self.deltaStarExt.append([])
+ self.xxExt.append([])
+ self.vtExt.append([])
+ for isec, sec in enumerate(body):
+ self.deltaStarExt[ibody].append([])
+ self.xxExt[ibody].append([])
+ self.vtExt[ibody].append([])
+ for i, vreg in enumerate(sec):
+ if vreg.getType() == 'fuselage' or vreg.getType() == 'wake':
+ self.deltaStarExt[ibody][isec].append([]) # only one side
+ self.xxExt[ibody][isec].append([])
+ self.vtExt[ibody][isec].append([])
+ elif vreg.getType() == 'wing':
+ self.deltaStarExt[ibody][isec].append([]) # upper and lower
+ self.deltaStarExt[ibody][isec].append([])
+ self.xxExt[ibody][isec].append([])
+ self.xxExt[ibody][isec].append([])
+ self.vtExt[ibody][isec].append([])
+ self.vtExt[ibody][isec].append([])
# External variables
if not self.vinit:
@@ -279,9 +328,7 @@ class SolversInterface:
self.xxExt[ibody][isec][i] = loc
self.deltaStarExt[ibody][isec][i] = np.zeros(reg.getnNodes())
self.vtExt[ibody][isec][i] = np.zeros(reg.getnNodes())
- self.vinit = True
- print('Viscous solver initialized')
- return 0
+ self.vinit = True
def _setViscousMesh(self):
"""Set the mesh from vBnd to the viscous solver.
@@ -313,16 +360,16 @@ class SolversInterface:
pystruct = self.vBnd[ibody][isec][iReg]
# Create node
nodes = blast.std_vector_blNode()
- for inod, nod in enumerate(pystruct.getNodesCoord(reg.getName())):
+ for inod, nod in enumerate(pystruct.getNodes(reg.getName())):
pos = tbox.Vector3d(nod[xIdx], nod[yIdx], nod[zIdx])
- row = int(pystruct.getNodesRow(reg.getName())[inod]) if self.cfg['interpolator'] == 'Matching' else inod
+ row = int(pystruct.getNodeConnectList(reg.getName())[inod]) if self.cfg['interpolator'] == 'Matching' else inod
nodes.push_back(blast.blNode(inod, pos, row))
# Set nodes
reg.setMesh(nodes)
# Set elements connectivity
connectElems = tbox.std_vector_int()
if self.cfg['interpolator'] == 'Matching':
- for val in pystruct.getConnectList(reg.getName()):
+ for val in pystruct.getElmsConnectList(reg.getName()):
connectElems.push_back(int(val))
else:
for i in range(reg.getnElms()):
@@ -359,28 +406,23 @@ class SolversInterface:
# Update the mesh in data structures
for ibody in range(self.getnBodies()):
for ireg, reg in enumerate(self.iBnd[ibody]):
- for i, row in enumerate(reg.connectListNodes):
+ nodes = np.zeros((meshList[ibody][ireg].shape[0], 3), dtype=float)
+ elems = np.zeros((meshList[ibody][ireg].shape[0]-1, 3), dtype=float)
+ for i, row in enumerate(reg.getConnectListNodes()):
for idim in range(self.getnDim()):
- reg.nodesCoord[i, idim] = meshList[ibody][ireg][row, idim]
-
- # Recompute elements data
- for ibody in range(self.getnBodies()):
- for ireg, reg in enumerate(self.iBnd[ibody]):
- for i in range(len(reg.nodesCoord[:,0])-1):
+ nodes[i, idim] = meshList[ibody][ireg][row, idim]
+ for i in range(reg.getnElms()-1):
for idim in range(self.getnDim()):
- reg.elemsCoord[i, idim] = (reg.nodesCoord[i, idim] + reg.nodesCoord[i+1, idim])/2
+ elems[i, idim] = (nodes[i, idim] + nodes[i+1, idim])/2
+ reg.updateMesh(nodes, elems)
if self.cfg['interpolator'] == 'Matching':
for ibody in range(self.getnBodies()):
for iSec in range(len(self.vBnd)):
for iReg, vreg in enumerate(self.vBnd[ibody][iSec]):
ireg = self.iBnd[ibody][iReg]
- for i, row in enumerate(ireg.connectListNodes):
- for idim in range(self.getnDim()):
- vreg.nodesCoord[i, idim] = ireg.nodesCoord[i, idim]
- for i, row in enumerate(ireg.connectListElems):
- for idim in range(self.getnDim()):
- vreg.elemsCoord[i, idim] = ireg.elemsCoord[i, idim]
+ vreg.updateMesh(ireg.getNodes('all'), ireg.getElms('all'))
+
elif self.cfg['interpolator'] == 'Rbf':
upNodes, lwNodes, upElms, lwElms = self.getUpperLowerSides()
for ibody in range(self.getnBodies()):
@@ -479,8 +521,8 @@ class SolversInterface:
if not self.iBnd[ibody][0].sidesIdentified:
raise RuntimeError(f'Cannot generate sections in auto for body {bodyName}: sides not identified.')
- upper_nodes = self.iBnd[ibody][0].getUpperNodes().copy()
- lower_nodes = self.iBnd[ibody][0].getLowerNodes().copy()
+ upper_nodes = self.iBnd[ibody][0].getNodes('upper_fromLE')
+ lower_nodes = self.iBnd[ibody][0].getNodes('lower_fromLE')
# Remove points higher than span
if ndim == 3:
@@ -557,8 +599,9 @@ class SolversInterface:
# Wake
wake = np.zeros((0,3))
- for i in range(self.iBnd[ibody][1].nPoints):
- wake = np.row_stack((wake, [self.iBnd[ibody][1].nodesCoord[i,0], self.iBnd[ibody][1].nodesCoord[i,1], self.iBnd[ibody][1].nodesCoord[i,2]]))
+ wakeNodes = self.iBnd[ibody][1].getNodes('all')
+ for i in range(self.iBnd[ibody][1].getnNodes()):
+ wake = np.row_stack((wake, [wakeNodes[i,0], wakeNodes[i,1], wakeNodes[i,2]]))
interpolated_wake = []
@@ -602,12 +645,9 @@ class SolversInterface:
elems_coords_wake[isec][i, :] = (sec[i, :] + sec[i+1, :]) / 2
for iy in range(len(sections_final)):
- self.vBnd[ibody][iy][0].initStructures(sections_final[iy].shape[0], elems_coords_wing[iy].shape[0])
- self.vBnd[ibody][iy][0].nodesCoord = sections_final[iy]
- self.vBnd[ibody][iy][0].elemsCoord = elems_coords_wing[iy]
- self.vBnd[ibody][iy][1].initStructures(wake_final[iy].shape[0], elems_coords_wake[iy].shape[0])
- self.vBnd[ibody][iy][1].nodesCoord = wake_final[iy]
- self.vBnd[ibody][iy][1].elemsCoord = elems_coords_wake[iy]
+ self.vBnd[ibody][iy][0].setMesh(sections_final[iy], elems_coords_wing[iy])
+ if len(self.vBnd[ibody][iy]) > 1:
+ self.vBnd[ibody][iy][1].setMesh(wake_final[iy], elems_coords_wake[iy])
def _getSectionsFromVTK(self):
@@ -669,15 +709,12 @@ class SolversInterface:
elems_coords_wing = np.empty((section.shape[0]-1, section.shape[1]), dtype=float, order='C')
for ielm in range(section.shape[0]-1):
elems_coords_wing[ielm, :] = (section[ielm, :] + section[ielm+1, :]) / 2
-
- self.vBnd[ibody][isec][0].initStructures(section.shape[0], elems_coords_wing.shape[0])
- self.vBnd[ibody][isec][0].nodesCoord = section
- self.vBnd[ibody][isec][0].elemsCoord = elems_coords_wing
+ self.vBnd[ibody][isec][0].setMesh(section, elems_coords_wing)
# Wake
hasWake = False
for ibnd in self.iBnd[ibody]:
- if ibnd.type == 'wake':
+ if ibnd.getType() == 'wake':
hasWake = True
break
@@ -696,9 +733,7 @@ class SolversInterface:
for ielm in range(wake_final.shape[0]-1):
elems_coords_wake[ielm, :] = (wake_final[ielm, :] + wake_final[ielm+1, :]) / 2
- self.vBnd[ibody][isec][1].initStructures(wake_final.shape[0], elems_coords_wake.shape[0])
- self.vBnd[ibody][isec][1].nodesCoord = wake_final
- self.vBnd[ibody][isec][1].elemsCoord = elems_coords_wake
+ self.vBnd[ibody][isec][1].setMesh(wake_final, elems_coords_wake)
def _getSectionsFromFile(self):
"""Obtain the nodes' location and row and cg of all elements.
@@ -780,10 +815,8 @@ class SolversInterface:
for i in range(len(cgSec)):
cgSec[i,:3] = (nodesSec[i,:3] + nodesSec[i+1,:3])/2
cgSec[i, 3] = i
-
- self.vBnd[ibody][iy][iReg].initStructures(nodesSec.shape[0], cgSec.shape[0])
- self.vBnd[ibody][iy][iReg].nodesCoord = nodesSec
- self.vBnd[ibody][iy][iReg].elemsCoord = cgSec
+
+ self.vBnd[ibody][iy][iReg].setMesh(nodesSec[:,:3], cgSec[:,:3])
def _checkInput(self, cfg:dict)->dict:
"""Check the input configuration.
@@ -840,9 +873,9 @@ class SolversInterface:
else:
iname = 'iWing'
vname = 'vAirfoil'
- iBnd[ibody].append(Group(iname, typ, self.getnDim()))
+ iBnd[ibody].append(InviscidData(iname, typ, self.getnDim()))
for isec, s in enumerate(self.cfg['sections'][bodyName]):
- vBnd[ibody][isec].append(Group(vname, typ, self.getnDim()))
+ vBnd[ibody][isec].append(ViscousData(vname, typ, self.getnDim()))
return iBnd, vBnd
def _getInterpolator(self)->object:
@@ -860,11 +893,13 @@ class SolversInterface:
for iSec in range(len(self.vBnd)):
for iReg, vreg in enumerate(self.vBnd[ibody][iSec]):
ireg = self.iBnd[ibody][iReg]
- vreg.initStructures(ireg.nPoints, ireg.nElems)
- vreg.nodesCoord = ireg.nodesCoord.copy()
- vreg.elemsCoord = ireg.elemsCoord.copy()
- vreg.connectListNodes = ireg.connectListNodes.copy()
- vreg.connectListElems = ireg.connectListElems.copy()
+ nodes = ireg.getNodes('all').copy()
+ elems = ireg.getElms('all').copy()
+ vreg.setMesh(nodes, elems)
+ vreg.setConnectLists(ireg._connectNodes.copy(), ireg._connectElems.copy())
+ vreg.setRows(ireg.getNodeRows(), ireg.getElmRows())
+ # vreg.connectListNodes = ireg._connectNodes.copy()
+ # vreg.connectListElems = ireg._connectElems.copy()
elif self.cfg['interpolator'] == 'Rbf':
from blast.interfaces.interpolators.blRbfInterpolator import RbfInterpolator as interp
@@ -910,20 +945,24 @@ class SolversInterface:
# Get all named regions
allregs = []
for vreg in self.vBnd[ibody][0]:
- allregs.append(vreg.type)
+ allregs.append(vreg.getType())
# Compute average and max AR of viscous cells
avgAR = 0
nNodes = 0
if self.getnDim() == 3:
point_cpt = 0
for isec in range(1,len(self.vBnd[ibody])):
- if self.vBnd[ibody][isec][0].nodesCoord.shape[0] != self.vBnd[ibody][isec-1][0].nodesCoord.shape[0]:
- # Sections do not have the same number of elements
+ wing_section = self.vBnd[ibody][isec][0]
+ wing_section_prev = self.vBnd[ibody][isec-1][0]
+ if wing_section.getnNodes() != wing_section_prev.getnNodes():
+ # Sections do not have the same number of nodes
point_cpt = 1
continue
- for ipoint in range(1,self.vBnd[ibody][isec][0].nodesCoord.shape[0]):
- dx = abs((self.vBnd[ibody][isec][0].nodesCoord[ipoint,0] - self.vBnd[ibody][isec][0].nodesCoord[ipoint-1,0]))
- dy = abs((self.vBnd[ibody][isec][0].nodesCoord[ipoint,1] - self.vBnd[ibody][isec-1][0].nodesCoord[ipoint,1]))
+ nodes_sec = wing_section.getNodes('all')
+ nodes_prev = wing_section_prev.getNodes('all')
+ for ipoint in range(1,wing_section.getnNodes()):
+ dx = abs((nodes_sec[ipoint,0] - nodes_sec[ipoint-1,0]))
+ dy = abs((nodes_sec[ipoint,1] - nodes_prev[ipoint,1]))
# AR = b^2 / S
avgAR += dy**2 / (dx * dy)
point_cpt += 1
@@ -931,10 +970,86 @@ class SolversInterface:
avgAR /= point_cpt
elif self.getnDim() == 2:
for isec in range(len(self.vBnd[ibody])):
- for ipoint in range(1,self.vBnd[ibody][isec][0].nodesCoord.shape[0]):
+ for ipoint in range(1,self.vBnd[ibody][isec][0].getnNodes()):
nNodes += 1
print(f' - {bodyName}: {len(self.vBnd[ibody])} section(s) with sides {allregs},\n\t{nNodes} surface nodes, avg AR {avgAR:.2f}')
print('')
+
+ def save_restart(self, cd, sfx=''):
+ import h5py
+
+ if 'restart_outputName' in self.cfg:
+ filename = self.cfg['restart_outputName']
+ else:
+ filename = 'blaster_restart'
+
+ if self.getVerbose() > 0:
+ print(f'Saving restart file {filename+sfx}.h5', end='...')
+ with h5py.File(f'{filename+sfx}.h5', 'w') as f:
+ f.attrs["cd"] = cd
+ f.attrs["nDim"] = self.getnDim()
+ for ibody in range(len(self.vBnd)):
+ body = f.create_group(f"body{ibody}")
+ for isec in range(len(self.vBnd[ibody])):
+ section = body.create_group(f"section{isec}")
+ for ireg in range(len(self.vBnd[ibody][isec])):
+ pyreg = self.vBnd[ibody][isec][ireg]
+ region = section.create_group(f"region{ireg}")
+ region.attrs["name"] = pyreg.getName()
+ region.attrs["type"] = pyreg.getType()
+ if pyreg.getStagPoint() is None:
+ region.attrs["stagnation_index"] = -1
+ else:
+ region.attrs["stagnation_index"] = pyreg.getStagPoint()
+ region.create_dataset("nodes", data=pyreg.getNodes('all'))
+ region.create_dataset("elems", data=pyreg.getElms('all'))
+ region.create_dataset("blowing", data=pyreg.getBlowingVelocity('all'))
+ for ireg in range(len(self.deltaStarExt[ibody][isec])):
+ region = section.create_group(f"blInteraction{ireg}")
+ region.create_dataset("deltaStarExt", data=self.deltaStarExt[ibody][isec][ireg])
+ region.create_dataset("xxExt", data=self.xxExt[ibody][isec][ireg])
+ region.create_dataset("vtExt", data=self.vtExt[ibody][isec][ireg])
+ if self.getVerbose() > 0:
+ print('done')
+
+ def _load_restart(self, restart_file, vBnd):
+ # Override viscous initialization flag to avoid reinitializing the interaction laws
+ self.vinit = True
+ self._initializeInteractionLaw()
+
+ import h5py
+ with h5py.File(restart_file+'.h5', 'r') as f:
+ self.cdPrev = f.attrs["cd"]
+ if f.attrs["nDim"] != self.getnDim():
+ raise RuntimeError("Dimension mismatch in restart file")
+ for ibody in range(len(vBnd)):
+ body = f[f"body{ibody}"]
+ for isec in range(len(vBnd[ibody])):
+ section = body[f"section{isec}"]
+ # The regions are defined on wing and wake (no upper/lower separation)
+ for ireg in range(len(vBnd[ibody][isec])):
+ pyreg = vBnd[ibody][isec][ireg]
+ region = section[f"region{ireg}"]
+ if pyreg.getName() != region.attrs["name"]:
+ raise RuntimeError(f"Name mismatch in restart file for body {ibody}, section {isec}, region {ireg}, expected {pyreg.getName()}, got {region['name'][()]}")
+ if pyreg.getType() != region.attrs["type"]:
+ raise RuntimeError(f"Type mismatch in restart file for body {ibody}, section {isec}, region {ireg}, expected {pyreg.getType()}, got {region['type'][()]}")
+ if region.attrs["stagnation_index"] == -1:
+ stag = None # Wake
+ else:
+ stag = region.attrs["stagnation_index"]
+ pyreg.setStagPoint(stag)
+ nodes = region["nodes"][:]
+ elems = region["elems"][:]
+ pyreg.setMesh(nodes, elems)
+ pyreg.setBlowingVelocity(region["blowing"])
+ # Coupling parameters are defined on upper/lower sides and wake
+
+ for iReg in range(len(self.deltaStarExt[ibody][isec])):
+ blregion = section[f"blInteraction{iReg}"]
+ self.deltaStarExt[ibody][isec][iReg] = blregion["deltaStarExt"][:]
+ self.vtExt[ibody][isec][iReg] = blregion["vtExt"][:]
+ self.xxExt[ibody][isec][iReg] = blregion["xxExt"][:]
# Abstract methods for the interface
def run()->int:
diff --git a/blast/tests/t_restart_2D.py b/blast/tests/t_restart_2D.py
new file mode 100644
index 0000000..d4df788
--- /dev/null
+++ b/blast/tests/t_restart_2D.py
@@ -0,0 +1,168 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2022 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+# @author Paul Dechamps
+# @date 2025
+# Test the restart file
+
+# Imports.
+
+import blast.blUtils as vutils
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+ import os.path
+ # Parameters
+ return {
+ # Options
+ 'scenario' : 'aerodynamic',
+ 'task' : 'analysis',
+ 'Threads' : nthrds, # number of threads
+ 'Verb' : verb, # verbosity
+ # Model (geometry or mesh)
+ 'File' : os.path.dirname(os.path.abspath(__file__)) + '/../models/dart/n0012.geo', # Input file containing the model
+ 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'growthRatio': 1.3, 'msTe' : 0.01, 'msLe' : 0.01}, # parameters for input file model
+ 'Dim' : 2, # problem dimension
+ 'Format' : 'gmsh', # save format (vtk or gmsh)
+ # Markers
+ 'Fluid' : 'field', # name of physical group containing the fluid
+ 'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+ 'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+ 'Wake' : 'wake', # LIST of names of physical group containing the wake
+ 'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
+ 'Te' : 'te', # LIST of names of physical group containing the trailing edge
+ 'Upstream' : 'upstream',
+ # Freestream
+ 'M_inf' : 0.2, # freestream Mach number
+ 'AoA' : 1., # freestream angle of attack
+ # Geometry
+ 'S_ref' : 1., # reference surface length
+ 'c_ref' : 1., # reference chord length
+ 'x_ref' : .25, # reference point for moment computation (x)
+ 'y_ref' : 0.0, # reference point for moment computation (y)
+ 'z_ref' : 0.0, # reference point for moment computation (z)
+ # Numerical
+ 'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+ 'G_fill' : 2, # fill-in factor for GMRES preconditioner
+ 'G_tol' : 1e-5, # tolerance for GMRES
+ 'G_restart' : 50, # restart for GMRES
+ 'Rel_tol' : 1e-6, # relative tolerance on solver residual
+ 'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+ 'Max_it' : 20 # solver maximum number of iterations
+ }
+
+def cfgBlast(verb):
+ return {
+ 'Re' : 1e7, # Freestream Reynolds number
+ 'CFL0' : 1, # Inital CFL number of the calculation
+ 'couplIter': 50, # Maximum number of iterations
+ 'couplTol' : 1e-4, # Tolerance of the VII methodology
+ 'iterPrint': 1, # int, number of iterations between outputs
+ 'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
+ 'xtrF' : [None, 0.4], # Forced transition locations
+ 'interpolator' : 'Matching',# Interpolator for the coupling
+
+ # Restart file
+ 'restart_solution' : False, # Restart solution from file
+ 'restart_inputName': None, # Name of the restart file to read
+ 'restart_outputName': 'naca0012_restart', # Name of the restart file saved on disk
+ 'iterSaveRestart' : 1, # Number of iterations between restart file save
+ 'saveCouplingIters': True,
+ }
+
+def main():
+ import os.path
+
+ # Iteration to restart
+ iterRestart = 10
+
+ # Timer.
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ args = parseargs()
+ icfg = cfgInviscid(args.k, args.verb)
+ vcfg = cfgBlast(args.verb)
+
+ tms['pre'].start()
+ coupler, isol, vsol = vutils.initBlast(icfg, vcfg)
+ tms['pre'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ aeroCoeffs_ref = coupler.run()
+ tms['solver'].stop()
+
+ del isol, vsol, coupler
+
+ vcfg["restart_solution"] = True
+ vcfg["restart_inputName"] = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))\
+ + '/workspace/blast_tests_t_restart/'\
+ + vcfg["restart_outputName"] + f"_{iterRestart}"
+ vcfg["iterSaveRestart"] = float('inf')
+
+
+ coupler, isol, vsol = vutils.initBlast(icfg, vcfg)
+ aeroCoeff_restart = coupler.run()
+
+ # Save in vtk format
+ vutils.save(vsol.bodies[0].sections)
+
+ # Display results.
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' Re M alpha Cl Cd Cdp Cdf Cm')
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(vcfg['Re']/1e6, isol.getMinf(), isol.getAoA()*180/math.pi, isol.getCl(), vsol.Cdt, vsol.Cdp, vsol.Cdf, isol.getCm()))
+ print('')
+ total_iters_ref = len(aeroCoeffs_ref['Cl'])
+ total_iters_restart = len(aeroCoeff_restart['Cl'])
+ print(f"Number of iterations without restart: {total_iters_ref}")
+ print(f'Restarted at iteration {iterRestart}')
+ print(f"Number of iterations with restart: {total_iters_restart}")
+ tms['total'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+ print('SOLVERS statistics')
+ print(coupler.tms)
+
+ bl_solution = vutils.getSolution(vsol.bodies[0].sections, write=False)[0]
+
+ final_res_ref = np.log10((aeroCoeffs_ref['Cd'][-1] - aeroCoeffs_ref['Cd'][-2]) / aeroCoeffs_ref['Cd'][-1])
+ final_res_restart = np.log10((aeroCoeff_restart['Cd'][-1] - aeroCoeff_restart['Cd'][-2]) / aeroCoeff_restart['Cd'][-1])
+
+ # Test solution
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ tests.add(CTest('Cl compare', aeroCoeff_restart['Cl'][-1], aeroCoeffs_ref['Cl'][-1], 1e-8))
+ tests.add(CTest('Cd compare', aeroCoeff_restart['Cd'][-1], aeroCoeffs_ref['Cd'][-1], 1e-8))
+ tests.add(CTest('nIters', total_iters_restart, total_iters_ref-iterRestart-1, 0, forceabs=True))
+ tests.add(CTest('Final residual Cd (log)', final_res_restart, final_res_ref, 1e-6))
+ tests.run()
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
--
GitLab