Skip to content
Snippets Groups Projects
Commit 46d33784 authored by Thomée Corentin's avatar Thomée Corentin
Browse files

[WIP] Implementation of UPM vii interface

parent d97807b8
No related branches found
No related tags found
No related merge requests found
Pipeline #20512 failed
Hide whitespace changes
Inline Side-by-side
blast/interfaces/UPM/UPMInterface.py 0 → 100644
+ 230
− 0
View file @ 46d33784
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2024 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.
# UPM interface
# Corentin Thomee
import numpy as np
from fwk.coloring import ccolors
from blast.interfaces.DSolversInterface import SolversInterface
class UPMInterface(SolversInterface):
def __init__(self, UPMCfg, vSolver, _cfg):
"""
Initialize the UPM solver and interface
"""
try:
from modules.UPM.upm.api.core import initUPM
self.solver = initUPM(UPMCfg)
except:
print(ccolors.ANSI_RED, 'Could not load UPM. Make sure it is installed.', ccolors.ANSI_RESET)
raise RuntimeError('Import error')
SolversInterface.__init__(self, vSolver, _cfg)
def run(self):
self.solver.solve()
return 1 # TODO: implémenter des retours dans la fonction solve en fonction du succès (echec solver linéaire, ITMAX, ...) !!
def writeCp(self, sfx=''):
self.solver.saveCp('Cp'+sfx+'.dat') # Directly use the UPM function
def save(self, sfx='inviscid'):
# TODO, fonction pour restart ?
self.solver.save(self.argOut['wrt'], sfx)
def getAoA(self):
return self.solver.AoA
def getMinf(self):
return 0 # Incompressible solver
def getCl(self):
return self.solver.cl
def getCd(self):
return self.solver.cd
def getCm(self):
# TODO (Shouldn't be necessary, but not hard to code in the UPM solver. Required to add a reference [in cfg probably *or harcoded @ (0,0)*], then do some maths around it)
return self.solver.Cm
def getVerbose(self):
# TODO Not necessary !
return self.solver.verbose
def reset(self):
self.solver.init() # Should be fine
def imposeInvBC(self, couplIter):
# TODO
""" Extract the inviscid paramaters at the edge of the boundary layer
and use it as a boundary condition in the viscous solver.
Params
------
- couplIter (int): Coupling iteration.
"""
# Retreive parameters.
for n in range(2):
for iRow, row in enumerate(self.iBnd[n].nodesCoord[:,3]):
row=int(row)
for iDim in range(3):
self.iBnd[n].V[iRow,iDim] = self.solver.U[row][iDim]
self.iBnd[n].M[iRow] = self.solver.M[row]
self.iBnd[n].Rho[iRow] = self.solver.rho[row]
if self.cfg['nDim']==3:
self.iBnd[n].V[:,[1,2]] = self.iBnd[n].V[:,[2,1]]
self.setViscousSolver(couplIter)
def imposeBlwVel(self):
# TODO
""" Extract the solution of the viscous calculation (blowing velocity)
and use it as a boundary condition in the inviscid solver.
"""
self.getBlowingBoundary()
# Impose blowing velocities.
for n in range(len(self.iBnd)):
if self.cfg['nDim'] == 2:
self.iBnd[n].connectBlowingVel()
for i, blw in enumerate(self.iBnd[n].blowingVel):
self.blw[n].setU(i, blw)
def setMesh(self):
if self.cfg['nDim'] == 2:
self.getAirfoil()
elif self.cfg['nDim'] == 3:
raise RuntimeError('UPMInterface::3D not implemented in UPM')
else:
raise RuntimeError('UPMInterface::Could not resolve how to set\
the mesh. Either ''nDim'' is incorrect or ''msh'' was not set for\
3D computations')
def getAirfoil(self):
"""Create data structure for information transfer.
"""
# TODO: WIP
# Aifoil boundary
N_nodes = self.solver.N+1
N_elem = self.solver.N
self.iBnd[0].initStructures(N_nodes) # Initialize the structure for the airfoil boundary
N1 = np.zeros(N_nodes, dtype=int) # Node number
# Index in boundary.nodes
connectListNodes = np.zeros(N_nodes, dtype=int)
# Index in boundary.elems
connectListElems = np.zeros((N_elem), dtype=int)
data = np.zeros((N_nodes, 4))
for i in range(N_nodes):
data[i,0] = i
data[i,1] = self.solver.x[i]
data[i,2] = self.solver.y[i]
data[i,3] = 0 # No z coordinate in 2D
eData = np.zeros((N_elem,3), dtype=int)
elemData = np.zeros((N_elem,4))
# Table containing the element and its nodes
for i in range(N_elem):
eData[i,0] = i
eData[i,1] = i
eData[i,2] = i+1
elemData[i,0] = self.solver.x[i] + self.solver.x[i+1]
elemData[i,1] = self.solver.y[i] + self.solver.y[i+1]
elemData[i,2] = 0 # No z coordinate in 2D
elemData[i,3] = i
elemData[i, :3] /= 2
# Find the stagnation point
idxStag = np.argmin(data[:,1])
globStag = int(data[idxStag,0]) # Position of the stagnation point
# Find TE nodes.
# Assuming suction side element is above pressure side element in the
# y direction.
# TE nodes
idxTE = np.where(data[:,1] == np.max(data[:,1]))[0]
# TE element 1
idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], \
eData[:,2] == data[idxTE[0],0]))[0]
# TE element 2
idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], \
eData[:,2] == data[idxTE[1],0]))[0]
# y coordinates of TE element 1 CG
ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] \
+ data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2])
# y coordinates of TE element 2 CG
ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] \
+ data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2])
if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
upperGlobTE = int(data[idxTE[0],0])
lowerGlobTE = int(data[idxTE[1],0])
else:
upperGlobTE = int(data[idxTE[1],0])
lowerGlobTE = int(data[idxTE[0],0])
# Connectivity
connectListElems[0] = np.where(eData[:,2] == upperGlobTE)[0]
# number of the stag node elems.nodes -> globStag
N1[0] = eData[connectListElems[0],1]
#connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
connectListNodes[0] = np.where(data[:,0] == upperGlobTE)[0]
i = 1
upperStag = 0
lowerTE = 0
# Sort the suction part
while upperStag == 0:
N1[i] = eData[connectListElems[i-1],1] # Second node of the element
# Index of the first node of the next element in elems.nodes
connectListElems[i] = np.where(eData[:,2] == N1[i])[0]
# Index of the node in boundary.nodes
connectListNodes[i] = np.where(data[:,0] == N1[i])[0]
if eData[connectListElems[i],1] == globStag:
upperStag = 1
i += 1
self.idxRt = i+1
# Sort the pressure side
connectListElems[i] = np.where(eData[:,2] == globStag)[0]
connectListNodes[i] = np.where(data[:,0] == globStag)[0]
N1[i] = eData[connectListElems[i],2]
while lowerTE == 0:
# First node of the element
N1[i+1] = eData[connectListElems[i],1]
# Index of the second node of the next element in elems.nodes
connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0]
# Index of the node in boundary.nodes
connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0]
if eData[connectListElems[i+1],1] == lowerGlobTE:
lowerTE = 1
i += 1
connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
data[:,0] = data[connectListNodes,0]
data[:,1] = data[connectListNodes,1]
data[:,2] = data[connectListNodes,2]
data[:,3] = data[connectListNodes,3]
data[:,4] = data[connectListNodes,4]
self.iBnd[0].nodesCoord = np.column_stack((data[:,1], data[:,2],\
data[:,3], data[:,4]))
self.iBnd[0].setConnectList(connectListNodes, connectListElems)
# Wake: TODO !!!!
blast/models/upm/n0012.dat 0 → 100644
+ 101
− 0
View file @ 46d33784
1 1.66533453693773e-17
0.999013364214136 -0.000140835441469506
0.996057350657239 -0.000561762174624333
0.991143625364344 -0.00125808467901582
0.984291580564316 -0.00222211850309532
0.975528258147577 -0.00344339291590901
0.964888242944126 -0.0049089163617143
0.95241352623301 -0.0066034880139773
0.938153340021932 -0.00851003605407132
0.922163962751008 -0.0106099620050805
0.904508497187474 -0.0128834706327159
0.885256621387895 -0.0153098665821892
0.864484313710706 -0.0178678019484713
0.842273552964344 -0.0205354631781272
0.818711994874345 -0.0232906907908307
0.793892626146236 -0.0261110310394699
0.767913397489498 -0.0289737244086614
0.740876837050858 -0.0318556413749577
0.712889645782536 -0.0347331807208502
0.684062276342339 -0.0375821495487363
0.654508497187474 -0.0403776466819014
0.624344943582427 -0.0430939721514063
0.593690657292862 -0.0457045848329603
0.562666616782152 -0.048182128006622
0.531395259764657 -0.0504985387651411
0.5 -0.0526252520005716
0.468604740235343 -0.0545335034546376
0.437333383217848 -0.056194729404379
0.406309342707138 -0.0575810534030492
0.375655056417573 -0.0586658435667832
0.345491502812526 -0.0594243176484838
0.315937723657661 -0.0598341679994961
0.287110354217464 -0.0598761748565317
0.259123162949142 -0.0595347744929392
0.232086602510502 -0.0587985488245103
0.206107373853763 -0.0576606051295196
0.181288005125655 -0.0561188185772977
0.157726447035656 -0.054175916084712
0.135515686289294 -0.0518393873480343
0.114743378612105 -0.0491212173440891
0.0954915028125263 -0.0460374436990048
0.0778360372489925 -0.0426075515760506
0.0618466599780682 -0.0388537276075468
0.0475864737669902 -0.0348000023735976
0.0351117570558743 -0.0304713175395326
0.0244717418524232 -0.0258925586035047
0.0157084194356845 -0.0210875969711737
0.00885637463565564 -0.0160783855767767
0.00394264934276106 -0.0108841504478057
0.000986635785864221 -0.00552071653494481
0 0
0.000986635785864221 0.00552071653494481
0.00394264934276106 0.0108841504478057
0.00885637463565569 0.0160783855767768
0.0157084194356845 0.0210875969711737
0.0244717418524232 0.0258925586035047
0.0351117570558744 0.0304713175395326
0.0475864737669903 0.0348000023735976
0.0618466599780683 0.0388537276075468
0.0778360372489925 0.0426075515760506
0.0954915028125264 0.0460374436990048
0.114743378612105 0.0491212173440891
0.135515686289294 0.0518393873480343
0.157726447035656 0.054175916084712
0.181288005125655 0.0561188185772977
0.206107373853763 0.0576606051295196
0.232086602510502 0.0587985488245103
0.259123162949142 0.0595347744929392
0.287110354217464 0.0598761748565317
0.315937723657661 0.0598341679994961
0.345491502812526 0.0594243176484838
0.375655056417573 0.0586658435667832
0.406309342707138 0.0575810534030492
0.437333383217848 0.056194729404379
0.468604740235343 0.0545335034546376
0.5 0.0526252520005716
0.531395259764657 0.0504985387651411
0.562666616782152 0.048182128006622
0.593690657292863 0.0457045848329603
0.624344943582427 0.0430939721514063
0.654508497187474 0.0403776466819014
0.684062276342339 0.0375821495487363
0.712889645782536 0.0347331807208502
0.740876837050858 0.0318556413749577
0.767913397489498 0.0289737244086614
0.793892626146237 0.0261110310394698
0.818711994874345 0.0232906907908307
0.842273552964344 0.0205354631781272
0.864484313710706 0.0178678019484713
0.885256621387895 0.0153098665821892
0.904508497187474 0.0128834706327159
0.922163962751008 0.0106099620050805
0.938153340021932 0.00851003605407132
0.95241352623301 0.0066034880139773
0.964888242944126 0.0049089163617143
0.975528258147577 0.00344339291590901
0.984291580564316 0.00222211850309532
0.991143625364344 0.00125808467901582
0.996057350657239 0.000561762174624333
0.999013364214136 0.000140835441469506
1 -1.66533453693773e-17
\ No newline at end of file
blast/tests/UPM/naca0012_2D.py 0 → 100644
+ 168
− 0
View file @ 46d33784
#!/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 2022
# Test the blast implementation. The test case is a compressible attached transitional flow.
# Tested functionalities;
# - Time integration.
# - Two time-marching techniques (agressive and safe).
# - Transition routines.
# - Forced transition.
# - Compressible flow routines.
# - Laminar and turbulent flow.
#
# Untested functionalities.
# - Separation routines.
# - Multiple failure case at one iteration.
# - Response to unconverged inviscid solver.
# Imports.
import blast.utils as viscUtils
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 {
'chord': 1,
'max_time_steps': 50,
'max_inner_it': 20,
'inner_tolerance': 1e-4,
'cl_tolerance': 1e-5,
'it_solver_tolerance': 1e-4,
'inner_relaxation': 1,
'dt': 10,
'aoa': 1,
'naca':False,
'filename': os.path.dirname(os.path.abspath(__file__)) + '/../../models/upm/n0012.dat'
}
def cfgBlast(verb):
return {
'Re' : 1e7, # Freestream Reynolds number
'Minf' : 0.2, # Freestream Mach number (used for the computation of the time step only)
'CFL0' : 1, # Inital CFL number of the calculation
'Verb': verb, # Verbosity level of the solver
'couplIter': 25, # Maximum number of iterations
'couplTol' : 1e-4, # Tolerance of the VII methodology
'iterPrint': 5, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
'sections' : [0],
'xtrF' : [None, 0.4],# Forced transition location
'interpolator' : 'Matching',
'nDim' : 2
}
def main():
# Timer.
tms = fwk.Timers()
tms['total'].start()
args = parseargs()
icfg = cfgInviscid(args.k, args.verb)
vcfg = cfgBlast(args.verb)
tms['pre'].start()
coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg, iSolver='UPM')
tms['pre'].stop()
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
aeroCoeffs = coupler.run()
tms['solver'].stop()
# 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, iSolverAPI.getMinf(), iSolverAPI.getAoA()*180/math.pi, iSolverAPI.getCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm()))
# Write results to file.
vSolution = viscUtils.getSolution(vSolver)
vSolution['Cdt_int'] = vSolution['Cdf'] + iSolverAPI.getCd()
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
print('CPU statistics')
print(tms)
print('SOLVERS statistics')
print(coupler.tms)
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('Cl', iSolverAPI.getCl(), 0.230, 5e-2)) # XFOIL 0.2325
tests.add(CTest('Cd wake', vSolution['Cdt_w'], 0.0056, 1e-3, forceabs=True)) # XFOIL 0.00531
tests.add(CTest('Cd integral', vSolution['Cdt_int'], 0.0059, 1e-3, forceabs=True)) # XFOIL 0.00531
tests.add(CTest('Cdf', vSolution['Cdf'], 0.00465, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
tests.add(CTest('xtr_top', vSolution['xtrT'], 0.20, 0.03, forceabs=True)) # XFOIL 0.1877
tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
tests.add(CTest('Iterations', len(aeroCoeffs), 17, 0, forceabs=True))
tests.run()
# Show results
if not args.nogui:
iCp = viscUtils.read('Cp_inviscid.dat')
vCp = viscUtils.read('Cp_viscous.dat')
xfoilCp = viscUtils.read('../../blast/models/references/cpXfoil_n0012.dat', delim=None, skip = 1)
xfoilCpInv = viscUtils.read('../../blast/models/references/cpXfoilInv_n0012.dat', delim=None, skip = 1)
plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
np.column_stack((xfoilCp[:,0], xfoilCp[:,1])),
np.column_stack((iCp[:,0], iCp[:,3])),
np.column_stack((xfoilCpInv[:,0], xfoilCpInv[:,1]))],
'labels': ['Blast (VII)', 'XFOIL VII', 'DART (inviscid)', 'XFOIL (inviscid)'],
'lw': [3, 3, 2, 2],
'color': ['firebrick', 'darkblue', 'firebrick', 'darkblue'],
'ls': ['-', '-', '--', '--'],
'reverse': True,
'xlim':[0, 1],
'yreverse': True,
'legend': True,
'xlabel': '$x/c$',
'ylabel': '$c_p$'
}
viscUtils.plot(plotcp)
xfoilBl = viscUtils.read('../../blast/models/references/blXfoil_n0012.dat', delim=None, skip = 1)
plotcf = {'curves': [np.column_stack((vSolution['x'], vSolution['Cf'])),
np.column_stack((xfoilBl[:,1], xfoilBl[:,6]))],
'labels': ['Blast', 'XFOIL'],
'lw': [3, 3],
'color': ['firebrick', 'darkblue'],
'ls': ['-', '-'],
'reverse': True,
'xlim':[0, 1],
'legend': True,
'xlabel': '$x/c$',
'ylabel': '$c_f$'
}
viscUtils.plot(plotcf)
# eof
print('')
if __name__ == "__main__":
main()
blast/utils.py
+ 2
− 0
View file @ 46d33784
......@@ -73,6 +73,8 @@ def initBlast(iconfig, vconfig, iSolver='DART'):
# Inviscid solver
if iSolver == 'DART':
from blast.interfaces.dart.DartInterface import DartInterface as interface
elif iSolver == 'UPM':
from blast.interfaces.UPM.UPMInterface import UPMInterface as interface
else:
print(ccolors.ANSI_RED + 'Solver '+iSolver+' currently not implemented' + ccolors.ANSI_RESET)
raise RuntimeError
......
modules/CMakeLists.txt
+ 2
− 1
View file @ 46d33784
ADD_SUBDIRECTORY( dartflo )
\ No newline at end of file
ADD_SUBDIRECTORY( dartflo )
ADD_SUBDIRECTORY( UPM )
\ No newline at end of file
UPM @ c342cc54
Compare eaa1ffe2...c342cc54
Subproject commit eaa1ffe224a3c8baff05ee56849d9462353acfdb
Subproject commit c342cc54b5625120f9ca958eed362b32aca9060f
dartflo @ 3a33458d
Compare 52a4059c...3a33458d
Subproject commit 52a4059c97365dd23e7dba39fb4623623be5769d
Subproject commit 3a33458dd16375c87ca86a702be84340a4ccb39e
run.py
+ 4
− 0
View file @ 46d33784
......@@ -24,12 +24,16 @@ def main():
thisdir = os.path.split(os.path.abspath(__file__))[0]
fwkdir = os.path.abspath(os.path.join(thisdir, 'modules', 'dartflo', 'ext', 'amfe'))
dartflodir = os.path.abspath(os.path.join(thisdir, 'modules', 'dartflo'))
UPMdir = os.path.abspath(os.path.join(thisdir, 'modules', 'UPM'))
if not os.path.isdir(fwkdir):
raise Exception('blaster/modules/dartflo/amfe not found!\n')
if not os.path.isdir(dartflodir):
raise Exception('blaster/modules/dartflo not found!\n')
if not os.path.isdir(UPMdir):
raise Exception('blaster/modules/UPM not found!\n')
sys.path.append(fwkdir)
sys.path.append(dartflodir)
sys.path.append(UPMdir)
# adds "." to the pythonpath
sys.path.append(thisdir)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment