Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • Paul.Dechamps/dartflo
1 result
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 418 additions and 877 deletions
dart/src/wWakeResidual.h
View file @ 5ba134b1
......@@ -41,4 +41,4 @@ public:
};
} // namespace dart
#endif //WWAKERESIDUAL_H
#endif // WWAKERESIDUAL_H
dart/tests/adjoint.py
View file @ 5ba134b1
......@@ -25,13 +25,13 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
def main():
# timer
......@@ -50,12 +50,12 @@ def main():
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.0075, 'msLe' : 0.0075}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
# create mesh deformation handler
morpher = floD.morpher(msh, dim, ['airfoil'])
morpher = floD.morpher(msh, dim, 'airfoil')
tms['msh'].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil')
tms['pre'].stop()
# solve problem
......@@ -94,9 +94,9 @@ def main():
# extract Cp and drag sensitivities
Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.4e', ',', 'x, y, cp', '')
dCd = floU.extract(bnd.groups[0].tag.elems, adjoint.dCdMsh)
tboxU.write(dCd, 'dCd_airfoil.dat', '%1.5e', ', ', 'x, y, z, dCd_dx, dCd_dy, dCd_dz', '')
tboxU.write(dCd, 'dCd_airfoil.dat', '%1.4e', ',', 'x, y, dCd_dx, dCd_dy, dCd_dz', '')
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' M alpha Cl Cd dCl_da dCd_da')
......@@ -114,17 +114,17 @@ def main():
if M_inf == 0. and alpha == 0*np.pi/180:
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 6.9, 1e-2))
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.0, 1e-3))
tests.add(CTest('dCl_dMsh', dClX, 55.439, 1e-3))
tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-3))
tests.add(CTest('dCl_dAoA (msh)', dClAoA, 6.9, 1e-2))
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.0, 1e-3))
tests.add(CTest('dCl_dMsh', dClX, 36.9, 1e-2))
tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-2))
tests.add(CTest('dCl_dAoA (msh)', dClAoA, adjoint.dClAoa, 1e-2))
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, adjoint.dCdAoa, 1e-3))
elif M_inf == 0.7 and alpha == 2*np.pi/180:
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 11.8, 1e-2)) # FD 11.835 (step = 1e-5)
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.127, 1e-2)) # 0.12692 (step = 1e-5)
tests.add(CTest('dCl_dMsh', dClX, 82.67, 1e-3))
tests.add(CTest('dCd_dMsh', dCdX, 0.908, 1e-3))
tests.add(CTest('dCl_dAoA (msh)', dClAoA, 11.8, 1e-2))
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.127, 1e-2))
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 12.0, 1e-2))
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.142, 1e-2))
tests.add(CTest('dCl_dMsh', dClX, 61.8, 1e-2))
tests.add(CTest('dCd_dMsh', dCdX, 0.815, 1e-2))
tests.add(CTest('dCl_dAoA (msh)', dClAoA, adjoint.dClAoa, 1e-2))
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, adjoint.dCdAoa, 1e-2))
else:
raise Exception('Test not defined for this flow')
tests.run()
......
dart/tests/bli.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 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.
## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
# Amaury Bilocq
#
# Test the viscous-inviscid interaction scheme
# Reference to the master's thesis: http://hdl.handle.net/2268/252195
# Reference test cases with Naca0012 (different from master's thesis):
# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
# -> msLe = 0.001, xtrTop = 0.061
# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
#
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import math
import dart.utils as floU
import dart.default as floD
import dart.viscous.solver as floVS
import dart.viscous.coupler as floVC
import tbox
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
def main():
# timer
tms = fwk.Timers()
tms['total'].start()
# define flow variables
Re = 1e7
alpha = 5*math.pi/180
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
M_inf = 0.5
c_ref = 1
dim = 2
tol = 1e-4 # tolerance for the VII (usually one drag count)
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
tms['pre'].stop()
# solve viscous problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
isolver = floD.newton(pbl)
vsolver = floVS.Solver(Re)
coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
coupler.run()
tms['solver'].stop()
# extract Cp
Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
# 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(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
# display timers
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
print('CPU statistics')
print(tms)
# visualize solution and plot results
floD.initViewer(pbl)
tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58
tests.add(CTest('Cd', vsolver.Cd, 0.0062, 1e-3, forceabs=True))
tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True))
tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.03, forceabs=True)) # Xfoil 0.056
tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.03, forceabs=True))
elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69
tests.add(CTest('Cd', vsolver.Cd, 0.0067, 1e-3, forceabs=True))
tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True))
tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.03, forceabs=True)) # Xfoil 0.038
tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.03, forceabs=True))
elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39
tests.add(CTest('Cd', vsolver.Cd, 0.011, 1e-3, forceabs=True))
tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 1e-3, forceabs=True))
tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.03, forceabs=True)) # Xfoil 0.008
tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.03, forceabs=True))
else:
raise Exception('Test not defined for this flow')
tests.run()
# eof
print('')
if __name__ == "__main__":
main()
dart/tests/cylinder.py
View file @ 5ba134b1
......@@ -25,11 +25,11 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import math
import dart.default as floD
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import math
def main():
# timer
......@@ -59,19 +59,19 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, l_ref, l_ref, 0., 0., 0., 'cylinder', te = 'te', dbc=True)
pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, l_ref, l_ref, 0., 0., 0., 'cylinder', dbc=True)
tms['pre'].stop()
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['picard'].start()
solver0 = floD.picard(pbl)
cnvrgd0 = solver0.run()
solver0.run()
solver0.save(gmshWriter)
tms['picard'].stop()
tms['newton'].start()
solver1 = floD.newton(pbl)
cnvrgd1 = solver1.run()
solver1.run()
solver1.save(gmshWriter)
tms['newton'].stop()
......@@ -126,8 +126,6 @@ def main():
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
if not cnvrgd0 or not cnvrgd1:
raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('Picard: iteration count', solver0.nIt, 12, 1, forceabs=True))
tests.add(CTest('Picard: Neumann B.C. mean error', errNeu0, 0., 1e-2))
......
dart/tests/cylinder2D5.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 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.
## Compute the (non)linear flow around extruded 2D cylinder, with various B.C.
# Adrien Crovato
#
# Test solver convergence and boundary conditions
#
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import math
import dart.default as floD
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
def main():
# timer
tms = fwk.Timers()
tms['total'].start()
# define flow variables
alpha = 3*math.pi/180
U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1
M_inf = 0.1
dim = 3
# define dimension and mesh size
dmt = 1.0 # cylinder diameter
lgt = 8.0 # channel length
hgt = 6.0 # channel height
wdt = 3.0 # channel width
l_ref = dmt # reference length
S_ref = dmt*wdt # reference area
fms = 0.5 # farfield mesh size
nms = 0.1 # nearfield mesh size
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms["msh"].start()
pars = {'dmt' : dmt, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
msh, gmshWriter = floD.mesh(dim, 'models/cylinder_2D5.geo', pars, ['field', 'cylinder', 'symmetry', 'downstream'])
tms["msh"].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'te', dbc=True)
tms['pre'].stop()
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
solver = floD.newton(pbl)
cnvrgd = solver.run()
solver.save(gmshWriter)
tms['solver'].stop()
# compute mean error in Body, Neumann boundary condition (works only for Tri3 elements)
gradn = []
for e in bnd.groups[0].tag.elems:
gradx = 0
grady = 0
gradz = 0
ux = e.nodes[1].pos[0] - e.nodes[0].pos[0]
uy = e.nodes[1].pos[1] - e.nodes[0].pos[1]
uz = e.nodes[1].pos[2] - e.nodes[0].pos[2]
vx = e.nodes[2].pos[0] - e.nodes[0].pos[0]
vy = e.nodes[2].pos[1] - e.nodes[0].pos[1]
vz = e.nodes[2].pos[2] - e.nodes[0].pos[2]
nx = uy*vz - uz*vy
ny = uz*vx - ux*vz
nz = ux*vy - uy*vx
for n in e.nodes:
gradx += solver.U[n.row][0] / e.nodes.size()
grady += solver.U[n.row][1] / e.nodes.size()
gradz += solver.U[n.row][2] / e.nodes.size()
gradn_ = gradx * nx + grady * ny + gradz * nz
gradn.append(gradn_)
errNeu = sum(gradn) / len(gradn)
# compute largest error in Dirichlet boundary condition
errDir = 0.
for n in dirichlet.nodes:
errCur = solver.phi[n.row] - (n.pos[0]*U_inf[0] + n.pos[1]*U_inf[1] + n.pos[2]*U_inf[2])
if (abs(errCur) > errDir):
errDir = abs(errCur)
# compute largest error in Wake boundary condition
errWak = 0.
for n in wake.nodMap:
errCur = solver.Cp[n.row] - solver.Cp[wake.nodMap[n].row]
if (abs(errCur) > errWak):
errWak = abs(errCur)
# display timers
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
print('CPU statistics')
print(tms)
# visualize solution
floD.initViewer(pbl)
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
if (not cnvrgd):
raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
tests.add(CTest('Neumann B.C. mean error', errNeu, 0., 1e-2))
tests.add(CTest('Dirichlet B.C. max error', errDir, 0., 1e-12))
tests.add(CTest('Wake/Kutta B.C. max error', errWak, 0., 1e-1))
tests.run()
# eof
print('')
if __name__ == "__main__":
main()
dart/tests/cylinder3.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 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.
## Compute the (non)linear flow around 3D cylinder, with various B.C.
# Adrien Crovato
#
# Test solver convergence and boundary conditions
#
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import math
import dart.utils as floU
import dart.default as floD
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
def main():
# timer
tms = fwk.Timers()
tms['total'].start()
# define flow variables
alpha = 3*math.pi/180
U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1
M_inf = 0.1
dim = 3
# define dimension and mesh size
dmt = 1.0 # cylinder diameter
spn = 2.0 # cylinder span
tpt = 5e-2 # cylinder tip thickness (ratio wrt span)
lgt = 8.0 # channel length
hgt = 6.0 # channel height
wdt = 3.0 # channel width
l_ref = dmt # reference length
S_ref = dmt*spn # reference area
fms = 0.5 # farfield mesh size
nms = 0.1 # nearfield mesh size
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
pars = {'dmt' : dmt, 'spn' : spn, 'tpt' : tpt, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
msh, gmshWriter = floD.mesh(dim, 'models/cylinder_3.geo', pars, ['field', 'cylinder', 'symmetry', 'downstream'], wktp = 'wakeTip')
tms['msh'].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, dirichlet, wake, bnd,_ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'teTip', dbc=True)
tms['pre'].stop()
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
solver = floD.newton(pbl)
cnvrgd = solver.run()
solver.save(gmshWriter)
tms['solver'].stop()
# compute mean error in Body, Neumann boundary condition (works only for Tri3 elements)
gradn = []
for e in bnd.groups[0].tag.elems:
gradx = 0
grady = 0
gradz = 0
ux = e.nodes[1].pos[0] - e.nodes[0].pos[0]
uy = e.nodes[1].pos[1] - e.nodes[0].pos[1]
uz = e.nodes[1].pos[2] - e.nodes[0].pos[2]
vx = e.nodes[2].pos[0] - e.nodes[0].pos[0]
vy = e.nodes[2].pos[1] - e.nodes[0].pos[1]
vz = e.nodes[2].pos[2] - e.nodes[0].pos[2]
nx = uy*vz - uz*vy
ny = uz*vx - ux*vz
nz = ux*vy - uy*vx
for n in e.nodes:
gradx += solver.U[n.row][0] / e.nodes.size()
grady += solver.U[n.row][1] / e.nodes.size()
gradz += solver.U[n.row][2] / e.nodes.size()
gradn_ = gradx * nx + grady * ny + gradz * nz
gradn.append(gradn_)
errNeu = sum(gradn) / len(gradn)
# compute largest error in Dirichlet boundary condition
errDir = 0.
for n in dirichlet.nodes:
errCur = solver.phi[n.row] - (n.pos[0]*U_inf[0] + n.pos[1]*U_inf[1] + n.pos[2]*U_inf[2])
if (abs(errCur) > errDir):
errDir = abs(errCur)
# compute largest error in Wake boundary condition
errWak = 0.
for n in wake.nodMap:
errCur = solver.Cp[n.row] - solver.Cp[wake.nodMap[n].row]
if (abs(errCur) > errWak):
errWak = abs(errCur)
# display timers
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
print('CPU statistics')
print(tms)
# visualize solution
floD.initViewer(pbl)
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
if (not cnvrgd):
raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('iteration count', solver.nIt, 4, 1, forceabs=True))
tests.add(CTest('Neumann B.C. mean error', errNeu, 0., 1e-2))
tests.add(CTest('Dirichlet B.C. max error', errDir, 0., 1e-12))
tests.add(CTest('Wake/Kutta B.C. max error', errWak, 0., 1e-1))
tests.run()
# eof
print('')
if __name__ == "__main__":
main()
dart/tests/lift.py
View file @ 5ba134b1
......@@ -25,13 +25,13 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import math
import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import math
def main():
# timer
......@@ -53,7 +53,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil')
tms['pre'].stop()
# solve problem
......@@ -66,7 +66,7 @@ def main():
# extract Cp
Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.4e', ',', 'x, y, cp', '')
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' M alpha Cl Cd Cm')
......@@ -80,23 +80,24 @@ def main():
# visualize solution and plot results
floD.initViewer(pbl)
tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), True)
floD.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
if M_inf == 0 and alpha == 3*math.pi/180:
tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.1, 1.5e-1))
tests.add(CTest('min(Cp)', min(Cp[:,-1]), -1.1, 1.5e-1))
tests.add(CTest('Cl', solver.Cl, 0.36, 5e-2))
tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
elif M_inf == 0.6 and alpha == 2*math.pi/180:
tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.03, 1e-1))
tests.add(CTest('min(Cp)', min(Cp[:,-1]), -1.03, 1e-1))
tests.add(CTest('Cl', solver.Cl, 0.315, 5e-2))
tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
elif M_inf == 0.7 and alpha == 2*math.pi/180:
tests.add(CTest('iteration count', solver.nIt, 12, 3, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.28, 5e-2))
tests.add(CTest('Cl', solver.Cl, 0.388, 5e-2))
tests.add(CTest('min(Cp)', min(Cp[:,-1]), -1.28, 5e-2))
tests.add(CTest('Cl', solver.Cl, 0.391, 5e-2))
tests.add(CTest('Cd', solver.Cd, 0.0013, 5e-4, forceabs=True))
tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2))
else:
raise Exception('Test not defined for this flow')
......
dart/tests/lift3.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 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.
## Compute lifting (linear or nonlinear) potential flow around a NACA 0012 rectangular wing
# Adrien Crovato
#
# Test the nonlinear shock-capturing capability and the Kutta condition
#
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import math
import dart.utils as floU
import dart.default as floD
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
def main():
# timer
tms = fwk.Timers()
tms['total'].start()
# define flow variables
alpha = 3*math.pi/180
M_inf = 0.3
dim = 3
# define dimension and mesh size
spn = 1.0 # wing span
lgt = 3.0 # channel half length
hgt = 3.0 # channel half height
wdt = 3.0 # channel width
S_ref = 1.*spn # reference area
c_ref = 1 # reference chord
fms = 1.0 # farfield mesh size
nms = 0.02 # nearfield mesh size
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
tms['msh'].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
tms['pre'].stop()
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
solver = floD.newton(pbl)
solver.run()
solver.save(gmshWriter)
tms['solver'].stop()
# display timers
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
print('CPU statistics')
print(tms)
# visualize solution
floD.initViewer(pbl)
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
if alpha == 3*math.pi/180 and M_inf == 0.3 and spn == 1.0:
tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
tests.add(CTest('CL', solver.Cl, 0.135, 5e-2))
tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
tests.add(CTest('CS', solver.Cs, 0.0121, 5e-2))
tests.add(CTest('CM', solver.Cm, -0.0278, 1e-2))
else:
raise Exception('Test not defined for this flow')
tests.run()
# eof
print('')
if __name__ == "__main__":
main()
dart/tests/meshDef.py → dart/tests/morpher.py
View file @ 5ba134b1
......@@ -25,13 +25,13 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
def main():
# timer
......@@ -55,7 +55,7 @@ def main():
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
# create mesh deformation handler
mshDef = floD.morpher(msh, dim, ['airfoil'])
morpher = floD.morpher(msh, dim, 'airfoil')
# mesh the reference geometry
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01, 'angle' : alfa, 'xRot' : xc}
msh_ref, gmshWriter_ref = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
......@@ -63,9 +63,9 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, c_ref, c_ref, xc+dx, dy, 0., 'airfoil', te = 'te')
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, c_ref, c_ref, xc+dx, dy, 0., 'airfoil')
# set the refrence problem and add fluid, initial/boundary conditions
pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, c_ref, c_ref, xc, 0., 0., 'airfoil', te = 'te')
pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, c_ref, c_ref, xc, 0., 0., 'airfoil')
tms['pre'].stop()
# solver problem for 0°AOA
......@@ -78,7 +78,7 @@ def main():
# deform the mesh (dummy rotation + translation)
print(ccolors.ANSI_BLUE + 'PyDeforming...' + ccolors.ANSI_RESET)
tms['deform'].start()
mshDef.savePos()
morpher.savePos()
rMat = np.matrix([[np.cos(alfa), np.sin(alfa)], [-np.sin(alfa), np.cos(alfa)]])
vNods = {}
for e in bnd.groups[0].tag.elems:
......@@ -91,7 +91,7 @@ def main():
vNods[no].pos[0] = new[0,0] + dx + xc
vNods[no].pos[1] = new[1,0] + dy
# deform the mesh
mshDef.deform()
morpher.deform()
gmshWriter.save("n0012_def")
tms['deform'].stop()
......@@ -112,9 +112,9 @@ def main():
# extract Cp
Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.4e', ',', 'x, y, cp', '')
Cp_ref = floU.extract(bnd_ref.groups[0].tag.elems, solver_ref.Cp)
tboxU.write(Cp, 'Cp_airfoil_ref.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
tboxU.write(Cp, 'Cp_airfoil_ref.dat', '%1.4e', ',', 'x, y, cp', '')
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' case M alpha Cl Cd Cm')
......@@ -128,15 +128,15 @@ def main():
print(tms)
# visualize solution and plot results
tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), True)
tboxU.plot(Cp_ref[:,0], Cp_ref[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver_ref.Cl, solver_ref.Cd, solver_ref.Cm, 4), True)
floD.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
floD.plot(Cp_ref[:,0], Cp_ref[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver_ref.Cl, solver_ref.Cd, solver_ref.Cm, 4), 'invert': True})
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('solver: iteration count', solver.nIt, 5, 3, forceabs=True)) # todo: tolerance should be reset to 1 as soon as gmsh4 is used
tests.add(CTest('solver_ref: iteration count', solver_ref.nIt, 5, 1, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,3]), min(Cp_ref[:,3]), 5e-2))
tests.add(CTest('min(Cp)', min(Cp[:,-1]), min(Cp_ref[:,-1]), 5e-2))
tests.add(CTest('Cl', solver.Cl, solver_ref.Cl, 5e-2))
tests.add(CTest('Cm', solver.Cm, solver_ref.Cm, 5e-2))
tests.run()
......
dart/tests/nonlift.py
View file @ 5ba134b1
......@@ -24,15 +24,14 @@
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
# The residual might not fully converge if this test is used with gmsh4
import math
import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import math
def main():
# timer
......@@ -48,13 +47,13 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.01, 'msTe' : 0.01, 'msLe' : 0.01}
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.0075, 'msLe' : 0.0075}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil')
tms['pre'].stop()
# solve problem
......@@ -67,7 +66,7 @@ def main():
# extract Cp
Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.4e', ',', 'x, y, cp', '')
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' M alpha Cl Cd Cm')
......@@ -81,18 +80,19 @@ def main():
# visualize solution and plot results
floD.initViewer(pbl)
tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), True)
floD.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
if M_inf == 0:
tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.405, 5e-2))
tests.add(CTest('min(Cp)', min(Cp[:,-1]), -0.405, 5e-2))
elif M_inf == 0.7:
tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.63, 5e-2))
tests.add(CTest('min(Cp)', min(Cp[:,-1]), -0.63, 5e-2))
elif M_inf == 0.8:
tests.add(CTest('iteration count', solver.nIt, 15, 3, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.89, 5e-2))
tests.add(CTest('iteration count', solver.nIt, 13, 3, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,-1]), -0.89, 5e-2))
tests.add(CTest('Cd', solver.Cd, 0.0058, 5e-4, forceabs=True))
else:
raise Exception('Test not defined for this flow')
tests.add(CTest('Cl', solver.Cl, 0., 1e-2))
......
dart/tests/adjoint3.py → dart/tests/rae_25.py
View file @ 5ba134b1
......@@ -16,22 +16,20 @@
# limitations under the License.
## Compute adjoint solution of lifting (linear or nonlinear) potential flow around a NACA 0012
## Compute the nonlinear flow around extruded 2D rae2822
# Adrien Crovato
#
# Test the adjoint solver
# Test transonic capabilities and kutta condition
#
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
def main():
# timer
......@@ -39,31 +37,38 @@ def main():
tms['total'].start()
# define flow variables
alpha = 3*np.pi/180
M_inf = 0.3
c_ref = 1
span = 1
alpha = 1*np.pi/180
M_inf = 0.72
dim = 3
# define dimension and mesh size
spn = 0.2 # wing span
lgt = 6.0 # channel length
hgt = 6.0 # channel height
wdt = 3.0 # channel width
c_ref = 1.0 # reference length
S_ref = spn # reference area
fms = 1.0 # farfield mesh size
nms = 0.01 # nearfield mesh size
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
pars = {'spn' : span, 'lgt' : 3, 'wdt' : 3, 'hgt' : 3, 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.}
msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
# create mesh deformation handler
morpher = floD.morpher(msh, dim, ['wing'])
tms['msh'].stop()
tms["msh"].start()
pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
msh, gmshWriter = floD.mesh(dim, 'models/rae_25.geo', pars, ['field', 'wing', 'symmetry', 'downstream'])
morpher = floD.morpher(msh, dim, 'wing')
tms["msh"].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, span*c_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0.25, 0., 0., 'wing')
tms['pre'].stop()
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
solver = floD.newton(pbl)
solver.run()
status = solver.run()
solver.save(gmshWriter)
tms['solver'].stop()
......@@ -83,7 +88,6 @@ def main():
dCdX += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]))
dClX = np.sqrt(dClX)
dCdX = np.sqrt(dCdX)
# recover gradient wrt to AoA from mesh gradient
drot = np.array([[-np.sin(alpha), 0, np.cos(alpha)], [0, 0, 0], [-np.cos(alpha), 0, -np.sin(alpha)]])
dClAoA = 0
......@@ -95,8 +99,8 @@ def main():
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' M alpha Cl Cd dCl_da dCd_da')
print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}'.format(M_inf, alpha*180/np.pi, solver.Cl, solver.Cd, adjoint.dClAoa, adjoint.dCdAoa))
print(' M alpha Cl Cd Cm')
print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alpha*180/np.pi, solver.Cl, solver.Cd, solver.Cm))
# display timers
tms['total'].stop()
......@@ -104,25 +108,25 @@ def main():
print('CPU statistics')
print(tms)
# visualize solution
floD.initViewer(pbl)
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
if status > 0:
raise Exception(ccolors.ANSI_RED + 'Solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
if M_inf == 0.3 and alpha == 3*np.pi/180 and span == 1:
tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 2.5, 5e-2)) # FD 2.526 (step 1e-5)
tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.146, 5e-2)) # FD 0.1460 (step 1e-5)
tests.add(CTest('dCl_dMsh', dClX, 2.016, 5e-2))
tests.add(CTest('dCd_dMsh', dCdX, 0.133, 5e-2))
tests.add(CTest('dCl/dAoA (msh)', dClAoA, 2.5, 5e-2))
tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.146, 5e-2))
elif M_inf == 0.7 and alpha == 5*np.pi/180 and span == 2:
tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 4.9, 5e-2))
tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.454, 5e-2))
tests.add(CTest('dCl_dMsh', dClX, 2.571, 5e-2))
tests.add(CTest('dCd_dMsh', dCdX, 0.204, 5e-2))
tests.add(CTest('dCl/dAoA (msh)', dClAoA, 4.9, 5e-2))
tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.454, 5e-2))
else:
raise Exception('Test not defined for this flow')
tests.add(CTest('iteration count', solver.nIt, 6, 1, forceabs=True))
tests.add(CTest('CL', solver.Cl, 0.59, 5e-2))
tests.add(CTest('CD', solver.Cd, 0.0017, 1e-2))
tests.add(CTest('CS', solver.Cs, 0.0000, 1e-3, forceabs=True))
tests.add(CTest('CM', solver.Cm, -0.111, 5e-2))
tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 12.9, 5e-2))
tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.045, 1e-2, forceabs=True))
tests.add(CTest('dCl_dMsh', dClX, 22.1, 5e-2))
tests.add(CTest('dCd_dMsh', dCdX, 1.1, 1e-1))
tests.add(CTest('dCl/dAoA (msh)', dClAoA, adjoint.dClAoa, 5e-2))
tests.add(CTest('dCd/dAoA (msh)', dCdAoA, adjoint.dCdAoa, 1e-2, forceabs=True))
tests.run()
# eof
......
dart/tests/meshDef3.py → dart/tests/rae_3.py
View file @ 5ba134b1
......@@ -16,21 +16,20 @@
# limitations under the License.
## Compute the flow around a NACA 0012 wing on a deformed mesh
## Compute the nonlinear flow around rectangular 3D rae2822
# Adrien Crovato
#
# Test the mesh deformation process
# Test the direct and adjoint solvers, and the mesh morpher in 3D
#
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
import numpy as np
import dart.utils as floU
import dart.default as floD
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
def main():
# timer
......@@ -38,50 +37,43 @@ def main():
tms['total'].start()
# define flow variables
alpha0 = 0 # must be zero for this testfile
M_inf = 0.3
alpha = 0
M_inf = 0.8
dim = 3
# define dimension and mesh size
spn = 1.0 # wing span
lgt = 3.0 # channel half length
hgt = 3.0 # channel half height
lgt = 6.0 # channel length
hgt = 6.0 # channel height
wdt = 3.0 # channel width
S_ref = 1.*spn # reference area
c_ref = 1 # reference chord
c_ref = 1.0 # reference length
S_ref = spn # reference area
fms = 1.0 # farfield mesh size
nms = 0.05 # nearfield mesh size
nms = 0.02 # nearfield mesh size
# parameters for mesh deformation
alfa = 3*np.pi/180
alfa = 1*np.pi/180
xc = 0.25
dz_max = 0.2*spn
# mesh an airfoil
# mesh the geometry and create mesh deformation handler
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
tms['msh'].stop()
# create mesh deformation handler
mshDef = floD.morpher(msh, dim, ['wing'])
tms["msh"].start()
pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
msh, gmshWriter = floD.mesh(dim, 'models/rae_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
morpher = floD.morpher(msh, dim, 'wing')
tms["msh"].stop()
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip')
solver = floD.newton(pbl) # create the solver now so element memory is up to date
tms['pre'].stop()
# solver problem for 0°AOA
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver0'].start()
solver = floD.newton(pbl)
solver.run()
tms['solver0'].stop()
# deform the mesh (dummy rotation + bending)
print(ccolors.ANSI_BLUE + 'PyDeforming...' + ccolors.ANSI_RESET)
tms['deform'].start()
mshDef.savePos()
morpher.savePos()
rMat = np.matrix([[np.cos(alfa), np.sin(alfa)], [-np.sin(alfa), np.cos(alfa)]])
vNods = {}
for e in bnd.groups[0].tag.elems:
......@@ -94,21 +86,46 @@ def main():
vNods[no].pos[0] = new[0,0]
vNods[no].pos[2] = new[1,0] + dz_max/(spn*spn)*(vNods[no].pos[1]*vNods[no].pos[1])
# deform the mesh
mshDef.deform()
morpher.deform()
gmshWriter.save(msh.name+"_def")
tms['deform'].stop()
# solve problem for alfa AOA
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver1'].start()
solver.run()
tms['solver'].start()
status = solver.run()
solver.save(gmshWriter)
tms['solver1'].stop()
tms['solver'].stop()
# solve adjoint problem
print(ccolors.ANSI_BLUE + 'PyGradient...' + ccolors.ANSI_RESET)
tms['adjoint'].start()
adjoint = floD.adjoint(solver, morpher)
adjoint.run()
adjoint.save(gmshWriter)
tms['adjoint'].stop()
# compute norm of gradient wrt to mesh
dClX = 0
dCdX = 0
for n in msh.nodes:
dClX += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]).dot(np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]))
dCdX += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]))
dClX = np.sqrt(dClX)
dCdX = np.sqrt(dCdX)
# recover gradient wrt to AoA from mesh gradient
drot = np.array([[-np.sin(alfa), 0, np.cos(alfa)], [0, 0, 0], [-np.cos(alfa), 0, -np.sin(alfa)]])
dClAoA = 0
dCdAoA = 0
for n in bnd.nodes:
dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1], n.pos[2]]))
dClAoA += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]).dot(dx)
dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(dx)
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' case M alpha Cl Cd Cm')
print(' true {0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alfa*180/np.pi, solver.Cl, solver.Cd, solver.Cm))
print(' M alpha Cl Cd dCl_da dCd_da')
print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}'.format(M_inf, alfa*180/np.pi, solver.Cl, solver.Cd, adjoint.dClAoa, adjoint.dCdAoa))
# display timers
tms['total'].stop()
......@@ -116,15 +133,25 @@ def main():
print('CPU statistics')
print(tms)
# visualize solution
floD.initViewer(pbl)
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
if status > 0:
raise Exception(ccolors.ANSI_RED + 'Solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
if alfa == 3*np.pi/180 and M_inf == 0.3 and spn == 1.0:
tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True))
tests.add(CTest('CL', solver.Cl, 0.135, 5e-1))
tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035
tests.add(CTest('CS', solver.Cs, -0.011, 5e-2))
tests.add(CTest('CM', solver.Cm, 0.0061, 1e-2))
tests.add(CTest('iteration count', solver.nIt, 7, 1, forceabs=True))
tests.add(CTest('CL', solver.Cl, 0.21, 5e-2))
tests.add(CTest('CD', solver.Cd, 0.008, 1e-2))
tests.add(CTest('CS', solver.Cs, -0.020, 5e-2))
tests.add(CTest('CM', solver.Cm, -0.080, 5e-2))
tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 3.0, 5e-2))
tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.2, 5e-2))
tests.add(CTest('dCl_dMsh', dClX, 2.2, 5e-2))
tests.add(CTest('dCd_dMsh', dCdX, 0.2, 5e-2))
tests.add(CTest('dCl/dAoA (msh)', dClAoA, adjoint.dClAoa, 5e-2))
tests.add(CTest('dCd/dAoA (msh)', dCdAoA, adjoint.dCdAoa, 5e-2))
tests.run()
# eof
......
dart/utils.py
View file @ 5ba134b1
......@@ -19,7 +19,7 @@
## Utilities for dart
# Adrien Crovato
def computeAeroCoef(_x, _y, _Cp, _alpha):
def compute_aero_coeff(_x, _y, _Cp, _alpha):
"""Compute 2D aerodynamic coefficients
The coefficients are computed from Cp averaged at nodes, which might leads to innacurate results
For accurate results, use coefficients from solver or body objects
......@@ -65,7 +65,7 @@ def convex_sort(vNodes, vData):
angle = np.zeros((vNodes.size()))
i = 0
while i < len(data[:,0]):
angle[i] = math.atan2(data[i,1]-cg[1],data[i,0]-cg[0])
angle[i] = np.arctan2(data[i,1]-cg[1], data[i,0]-cg[0])
if angle[i] < 0:
angle[i] = 2 * np.pi + angle[i]
i+=1
......@@ -97,21 +97,20 @@ def extract(vElems, vData):
else:
raise RuntimeError('dart.utils.extract: unrecognized format for vData!')
# store data (not efficient, but OK since only meant for small 2D cases)
data = np.zeros((vElems.size()+1,3+size))
data = np.zeros((vElems.size()+1,2+size))
i = 0
while i < vElems.size()+1:
data[i,0] = vElems[i%vElems.size()].nodes[0].pos[0]
data[i,1] = vElems[i%vElems.size()].nodes[0].pos[1]
data[i,2] = vElems[i%vElems.size()].nodes[0].pos[2]
if size == 1:
data[i,3] = vData[vElems[i%vElems.size()].nodes[0].row]
data[i,2] = vData[vElems[i%vElems.size()].nodes[0].row]
else:
for j in range(size):
data[i,3+j] = vData[vElems[i%vElems.size()].nodes[0].row][j]
data[i,2+j] = vData[vElems[i%vElems.size()].nodes[0].row][j]
i += 1
return data
def writeSlices(mshName, ys, wId, n = 0):
def write_slices(mshName, ys, wId, sfx=''):
"""Write slice data for each ys coordinate along the span (only works with VTK)
Parameters
......@@ -122,8 +121,8 @@ def writeSlices(mshName, ys, wId, n = 0):
y-coordinates of slices
wId : int
index of physical tag to slice (see gmsh)
n : int (optional, > 0)
index of mesh/solution file (will be appended to mesh file name)
sfx : string (optional)
suffix that will be appended to mesh/solution file name
"""
import numpy as np
try:
......@@ -134,20 +133,18 @@ def writeSlices(mshName, ys, wId, n = 0):
import tboxVtk.reader as vtkR
import tboxVtk.cutter as vtkC
reader = vtkR.Reader()
if n == 0:
reader.open(mshName)
else:
reader.open(mshName+'_'+str(n))
reader.open(mshName+sfx)
cutter = vtkC.Cutter(reader.reader.GetOutputPort())
for i in range(0, len(ys)):
pts, elems, vals = cutter.extract(cutter.cut(wId, [0., ys[i], 0.], [0., 1., 0.]), 2, ['Cp'])
x_c = np.zeros((pts.shape[0],1))
x_c[:,0] = (pts[:,0] - min(pts[:,0])) / (max(pts[:,0]) - min(pts[:,0]))
data = np.hstack((pts, x_c))
data = np.hstack((data, vals['Cp']))
x_c = np.zeros((pts.shape[0], 2))
ile = np.argmin(pts[:,0]) # LE index
crd = max(pts[:,0]) - min(pts[:,0]) # chord
x_c[:,0] = (pts[:,0] - pts[ile,0]) / crd
x_c[:,1] = (pts[:,2] - pts[ile,2]) / crd
data = np.hstack((x_c, vals['Cp']))
tU.sort(elems, data)
data = np.vstack((data, data[0,:]))
if n == 0:
tU.write(data, 'slice_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
else:
tU.write(data, 'slice_'+str(n)+'_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
hdr = f'y = {ys[i]}, c = {crd}\n'
hdr += '{:>9s}, {:>10s}, {:>10s}'.format('x/c', 'z/c', 'cp')
np.savetxt(mshName+sfx+'_slice_'+str(i)+'.dat', data, fmt='%+1.4e', delimiter=',', header=hdr)
dart/validation/agard.py
View file @ 5ba134b1
......@@ -19,13 +19,13 @@
## Compute flow around the Agard 445 wing at 1 degrees AOA and Mach 0.80
# Adrien Crovato
import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tbxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
try:
import tboxVtk
......@@ -56,7 +56,7 @@ def main():
rlems = 0.0028 # root leading edge mesh size
rtems = 0.0056 # root trailing edge mesh size
tlems = 0.0018 # tip leading mesh size
ttems = 0.0045 # tip trailing mesh size
ttems = 0.0036 # tip trailing mesh size
fms = 1.0 # farfield mesh size
# mesh the wing
......@@ -69,7 +69,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip')
tms['pre'].stop()
# solve problem
......@@ -83,9 +83,9 @@ def main():
# post process
tms['post'].start()
if canPost:
floU.writeSlices(msh.name,[0.01, 0.37, 0.75],5)
data = tbxU.read('slice_1.dat')
tbxU.plot(data[:,3], data[:,4], 'x', 'Cp', 'Pressure distribution at MAC', True)
floU.write_slices(msh.name,[0.01, 0.37, 0.75], 5)
data = np.loadtxt(f'{msh.name}_slice_1.dat', delimiter=',', skiprows=2)
floD.plot(data[:,0], data[:,-1], {'xlabel': 'x/c', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
tms['post'].stop()
# display timers
......@@ -97,8 +97,8 @@ def main():
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('CL', solver.Cl, 0.062, 1e-1))
tests.add(CTest('CD', solver.Cd, 0.0006, 1e-1))
tests.add(CTest('CL', solver.Cl, 0.061, 1e-1))
tests.add(CTest('CD', solver.Cd, 0.0007, 1e-1))
tests.add(CTest('CS', solver.Cs, 0.0017, 1e-1))
tests.add(CTest('CM', solver.Cm, -0.064, 1e-1))
tests.run()
......
dart/validation/crm.py 0 → 100644
View file @ 5ba134b1
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 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.
## Compute the flow around the NASA CRM (Common Research Model) at CL = 0.5 and Mach 0.85
# CAD created by Guillem Battle I Capa
# Valid mesh obtained with gmsh 4.8.4
try:
import tboxVtk
import vtk
hasVTK = True
except:
hasVTK = False
def cfg():
from fwk.wutils import parseargs
import os.path
# Mesh size
wrms = 0.121 # wing root trailing edge mesh size
wyms = 0.073 # wing root trailing edge mesh size
wtms = 0.028 # wing root trailing edge mesh size
trms = 0.055 # wing root trailing edge mesh size
ttms = 0.023 # wing root trailing edge mesh size
fums = 0.628 # fuselage mesh size
fams = 36.8 # farfield mesh size
# Parse
args = parseargs()
# Parameters
return {
# Options
'Threads' : args.k, # number of threads
'Verb' : args.verb + 2, # verbosity
# Model (geometry or mesh)
'File' : os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/crm.geo'), # Input file containing the model
'Pars' : {'msWingR' : wrms, 'msWingY' : wyms, 'msWingT' : wtms,
'msTailR' : trms, 'msTailT' : ttms,
'msFus' : fums, 'msFar' : fams}, # parameters for input file model
'Dim' : 3, # problem dimension (2 or 3)
'Format' : 'vtk', # 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)
'Wings' : ['wing', 'tail'], # LIST of names of physical groups containing the lifting surface boundary (first element will be the body of interest for aerostructural and optimization)
'Wakes' : ['wakeW', 'wakeT'], # LIST of names of physical group containing the wake
'WakeTips' : ['wakeTipW', 'wakeTipT'], # LIST of names of physical group containing the edge of the wake
'WakeExs' : ['wakeExW', 'wakeExT'], # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.)
'Tes' : ['teW', 'teT'], # LIST of names of physical group containing the trailing edge
'Fuselage' : 'fuselage', # name of physical group containing the fuselage boundaries
'Symmetry' : 'symmetry', # name of physical group containing the symmetry boundaries
# Freestream
'M_inf' : 0.85, # freestream Mach number
'AoA' : 1.2, # freestream angle of attack [deg] (optional, default=0)
'AoS' : 0.0, # freestream angle of sideslip [deg] (optional, default=0)
'Q_inf' : 0.0, # freesteam dynamic pressure (only required for aerostructural computations)
# Geometry
'S_ref' : 383.69/2, # reference surface length
'c_ref' : 7.005, # reference chord length
'x_ref' : 33.667, # x-coordinate of reference point for moment computation
'y_ref' : 11.906, # y-coordinate of reference point for moment computation
'z_ref' : 4.520, # z-coordinate of reference point for moment computation
# Numerical
'LSolver' : 'GMRES', # inner solver (PARDISO, MUMPS or GMRES)
'G_fill' : 3, # fill-in factor for GMRES preconditioner (optional, default=2)
'G_tol' : 1e-5, # tolerance for GMRES (optional, default=1e-5)
'G_restart' : 50, # restart for GMRES (optional, default=50)
'Rel_tol' : 1e-6, # relative tolerance on solver residual
'Abs_tol' : 1e-8, # absolute tolerance on solver residual
'Max_it' : 25 # maximum number of iterations for nonlinear solver
}
def main():
from dart.api.core import init_dart
from dart import utils
from fwk import Timers
from fwk.testing import CTest, CTests
import numpy as np
import matplotlib.pyplot as plt
# pre
tms = Timers()
tms['pre'].start()
_dart = init_dart(cfg(), scenario='aerodynamic', task='analysis')
msh, wrt, sol = (_dart.get(key) for key in ['msh', 'wrt', 'sol'])
tms['pre'].stop()
# solve
tms['sol'].start()
sol.run()
sol.save(wrt)
tms['sol'].stop()
# post
tms['pos'].start()
if hasVTK:
utils.write_slices(msh.name, [3.81973, 5.8765, 8.2271, 11.753, 14.6913, 17.6295, 19.0986, 21.4492, 24.9751, 27.91338], 12)
for i in range(10):
data = np.loadtxt(f'{msh.name}_slice_{i}.dat', delimiter=',', skiprows=2) # read
plt.plot(data[:,0], data[:,-1], lw=2) # plot results
font = {'family': 'serif', 'color': 'darkred', 'weight': 'normal', 'size': 16}
plt.xlabel('$x/c$', fontdict=font)
plt.ylabel('$C_p$', fontdict=font)
plt.gca().invert_yaxis()
plt.show()
tms['pos'].stop()
# display timers
print('CPU statistics')
print(tms)
# check results
tests = CTests()
tests.add(CTest('CL', sol.Cl, 0.5105, 1e-1))
tests.add(CTest('CD', sol.Cd, 0.0143, 1e-1))
tests.add(CTest('CM', sol.Cm, -0.1214, 1e-1))
tests.run()
# eof
print('')
if __name__ == '__main__':
main()
dart/validation/onera.py
View file @ 5ba134b1
......@@ -19,13 +19,13 @@
## Compute flow around the Onera M6 wing at 3 degrees AOA and Mach 0.84
# Adrien Crovato
import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tbxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
try:
import tboxVtk
......@@ -51,7 +51,7 @@ def main():
wdt = 3.0 # channel width
S_ref = 0.7528 # reference area
c_ref = 0.64 # reference chord (MAC)
fms = 0.5 # farfield mesh size
fms = 0.6 # farfield mesh size
rlems = 0.004 # root leading edge mesh size
rtems = 0.008 # root trailing edge mesh size
tlems = 0.002 # tip leading mesh size
......@@ -67,7 +67,7 @@ def main():
# set the problem and add fluid, initial/boundary conditions
tms['pre'].start()
pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip')
tms['pre'].stop()
# solve problem
......@@ -81,9 +81,9 @@ def main():
# post process
tms['post'].start()
if canPost:
floU.writeSlices(msh.name,[0.01, 0.239, 0.526, 0.777, 0.957, 1.076, 1.136, 1.184],5)
data = tbxU.read('slice_2.dat')
tbxU.plot(data[:,3], data[:,4], 'x', 'Cp', 'Pressure distribution at MAC', True)
floU.write_slices(msh.name,[0.01, 0.239, 0.526, 0.777, 0.957, 1.076, 1.136, 1.184], 5)
data = np.loadtxt(f'{msh.name}_slice_2.dat', delimiter=',', skiprows=2)
floD.plot(data[:,0], data[:,-1], {'xlabel': 'x/c', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
tms['post'].stop()
# display timers
......@@ -97,8 +97,8 @@ def main():
tests = CTests()
tests.add(CTest('CL', solver.Cl, 0.29, 1e-1))
tests.add(CTest('CD', solver.Cd, 0.011, 1e-1))
tests.add(CTest('CS', solver.Cs, 0.014, 1e-1))
tests.add(CTest('CM', solver.Cm, -0.212, 1e-1))
tests.add(CTest('CS', solver.Cs, 0.016, 1e-1))
tests.add(CTest('CM', solver.Cm, -0.218, 1e-1))
tests.run()
# eof
......
dart/viscous/__init__.py deleted 100644 → 0
View file @ 7b79ac85
# -*- coding: utf-8 -*-
dart/viscous/airfoil.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 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.
## Airfoil around which the boundary layer is computed
# Amaury Bilocq
from dart.viscous.boundary import Boundary
import numpy as np
class Airfoil(Boundary):
def __init__(self, _boundary):
Boundary.__init__(self, _boundary)
self.name = 'airfoil' # type of boundary
self.T0 = 0 # initial condition for the momentum thickness
self.H0 = 0
self.n0 = 0
self.Ct0 = 0
self.TEnd = [0,0]
self.HEnd = [0,0]
self.CtEnd = [0,0]
self.nEnd = [0,0]
def initialConditions(self, Re, dv):
if dv > 0:
self.T0 = np.sqrt(0.075/(Re*dv))
else:
self.T0 = 1e-6
self.H0 = 2.23 # One parameter family Falkner Skan
self.n0 = 0
self.Ct0 = 0
return self.T0, self.H0, self.n0, self.Ct0
def connectList(self):
''' Sort the value read by the viscous solver/ Create list of connectivity
'''
N1 = np.zeros(self.nN, dtype=int) # Node number
connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
data = np.zeros((self.boundary.nodes.size(), 10))
i = 0
for n in self.boundary.nodes:
data[i,0] = n.no
data[i,1] = n.pos[0]
data[i,2] = n.pos[1]
data[i,3] = n.pos[2]
data[i,4] = self.v[i,0]
data[i,5] = self.v[i,1]
data[i,6] = self.Me[i]
data[i,7] = self.rhoe[i]
data[i,8] = self.deltaStar[i]
data[i,9] = self.xx[i]
i += 1
# Table containing the element and its nodes
eData = np.zeros((self.nE,3), dtype=int)
for i in range(0, self.nE):
eData[i,0] = self.boundary.tag.elems[i].no
eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
# Find the stagnation point
idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
# Find TE nodes (assuming suction side element is above pressure side element in the y direction)
idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes
idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
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 1 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]) # y coordinates of TE element 2 CG
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[:,1] == globStag)[0]
N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
i = 1
upperTE = 0
lowerTE = 0
# Sort the suction part
while upperTE == 0:
N1[i] = eData[connectListElems[i-1],2] # Second node of the element
connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
if eData[connectListElems[i],2] == upperGlobTE:
upperTE = 1
i += 1
# Sort the pressure side
connectListElems[i] = np.where(eData[:,2] == globStag)[0]
connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
N1[i] = eData[connectListElems[i],2]
while lowerTE == 0:
N1[i+1] = eData[connectListElems[i],1] # First node of the element
connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
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]
data[:,5] = data[connectListNodes,5]
data[:,6] = data[connectListNodes,6]
data[:,7] = data[connectListNodes,7]
data[:,8] = data[connectListNodes,8]
data[:,9] = data[connectListNodes,9]
# Separated upper and lower part
data = np.delete(data,0,1)
uData = data[0:np.argmax(data[:,0])+1]
lData = data[np.argmax(data[:,0])+1:None]
lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
return connectListNodes, connectListElems,[uData, lData]
dart/viscous/boundary.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 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.
## Base class representing a physical boundary
# Amaury Bilocq
import numpy as np
class Boundary:
def __init__(self, _boundary):
self.boundary = _boundary # boundary of interest
self.nN = len(self.boundary.nodes) # number of nodes
self.nE = len(self.boundary.tag.elems) # number of elements
self.v = np.zeros((self.nN, 3)) # velocity at edge of BL
self.u = np.zeros(self.nE) # blowing velocity
self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer
self.rhoe = np.zeros(self.nN) # density at the edge of the boundary layer
self.connectListnodes = np.zeros(self.nN) # connectivity for nodes
self.connectListElems = np.zeros(self.nE) # connectivity for elements
self.deltaStar = np.zeros(self.nN) # displacement thickness
self.xx = np.zeros(self.nN) # coordinates in the reference of the physical surface. Used to store the values for the interaction law
dart/viscous/coupler.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 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.
## Viscous-inviscid coupler (quasi-simultaneous coupling)
# Amaury Bilocq
import numpy as np
from fwk.coloring import ccolors
from dart.viscous.airfoil import Airfoil
from dart.viscous.wake import Wake
class Coupler:
def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer):
self.isolver =_isolver # inviscid solver
self.vsolver = _vsolver # viscous solver
self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
self.tol = _tol # tolerance of the VII
self.writer = _writer
def run(self):
''' Perform coupling
'''
# initialize loop
it = 0
converged = 0 # temp
CdOld = self.vsolver.Cd
while converged == 0:
print(ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET)
# run inviscid solver
self.isolver.run()
print('--- Viscous solver parameters ---')
print('Tolerance (drag count):', self.tol * 1000)
print('--- Viscous problem definition ---')
print('Reynolds number:', self.vsolver.Re)
print('')
for n in range(0, len(self.group)):
print('Computing for', self.group[n].name, '...', end = ' ')
for i in range(0, len(self.group[n].boundary.nodes)):
self.group[n].v[i,0] = self.isolver.U[self.group[n].boundary.nodes[i].row][0]
self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row][1]
self.group[n].v[i,2] = self.isolver.U[self.group[n].boundary.nodes[i].row][2]
self.group[n].Me[i] = self.isolver.M[self.group[n].boundary.nodes[i].row]
self.group[n].rhoe[i] = self.isolver.rho[self.group[n].boundary.nodes[i].row]
# run viscous solver
self.vsolver.run(self.group[n])
if self.group[n].name == 'airfoil':
self.group[n+1].T0 = self.group[n].TEnd[0]+self.group[n].TEnd[1]
self.group[n+1].H0 = (self.group[n].HEnd[0]*self.group[n].TEnd[0]+self.group[n].HEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
self.group[n+1].Ct0 = (self.group[n].CtEnd[0]*self.group[n].TEnd[0]+self.group[n].CtEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
# get blowing velocity from viscous solver and update inviscid solver
for i in range(0, self.group[n].nE):
self.group[n].boundary.setU(i, self.group[n].u[i])
print('done!')
print(' Iter Cd xtr_top xtr_bot Error')
print('{0:8d} {1:12.5f} {2:12.5f} {3:12.5f} {4:12.5f}'.format(it, self.vsolver.Cd, self.vsolver.xtr[0], self.vsolver.xtr[1], abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd))
print('')
# Converged or not
if abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd < self.tol:
converged = 1
else:
CdOld = self.vsolver.Cd
it += 1
self.vsolver.it += 1
# save results
self.isolver.save(self.writer)
self.vsolver.writeFile()
print('\n')