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
Showing with 83 additions and 143 deletions
......@@ -52,4 +52,4 @@ private:
} // namespace dart
#endif //WPICARD_H
#endif // WPICARD_H
......@@ -43,4 +43,4 @@ public:
};
} // namespace dart
#endif //WPOTENTIALRESIDUAL_H
#endif // WPOTENTIALRESIDUAL_H
......@@ -95,4 +95,4 @@ private:
} // namespace dart
#endif //WPROBLEM_H
#endif // WPROBLEM_H
......@@ -36,6 +36,7 @@
#include <tbb/spin_mutex.h>
#include <iomanip>
#include <fstream>
using namespace tbox;
using namespace dart;
......@@ -64,9 +65,9 @@ Solver::Solver(std::shared_ptr<Problem> _pbl,
<< "** _ /_/ /_ ___ | __ _/_ / **\n"
<< "** /_____/ /_/ |_/_/ |_| /_/ **\n"
<< "*******************************************************************************\n"
<< "** Hi! My name is DART v1.2.0-22.10 **\n"
<< "** Hi! My name is DART v1.3.0-24.12 **\n"
<< "** Adrien Crovato **\n"
<< "** ULiege 2018-2022 **\n"
<< "** ULiege 2018-2024 **\n"
<< "*******************************************************************************\n"
<< std::endl;
......@@ -170,7 +171,7 @@ void Solver::save(MshExport *mshWriter, std::string const &suffix)
<< pbl->alpha * 180 / 3.14159 << "deg AoA, "
<< pbl->beta * 180 / 3.14159 << "deg AoS)"
<< std::endl;
// setup results
// save results
Results results;
results.scalars_at_nodes["phi"] = &phi;
results.scalars_at_nodes["varPhi"] = &vPhi;
......@@ -179,10 +180,30 @@ void Solver::save(MshExport *mshWriter, std::string const &suffix)
results.scalars_at_nodes["rho"] = &rho;
results.scalars_at_nodes["Mach"] = &M;
results.scalars_at_nodes["Cp"] = &Cp;
// save (whole mesh and bodies)
mshWriter->save(pbl->msh->name + suffix, results);
for (auto bnd : pbl->bodies)
bnd->save(pbl->msh->name + '_' + bnd->groups[0]->tag->name + suffix, results);
// save aerodynamic loads
std::cout << "writing file: aeroloads" << suffix << ".dat... " << std::flush;
std::ofstream outfile;
outfile.open("aeroloads" + suffix + ".dat");
// total
outfile << "# Total - " << pbl->msh->name << std::endl;
outfile << std::fixed << std::setprecision(12)
<< "Cl = " << Cl << "\n"
<< "Cd = " << Cd << "\n"
<< "Cs = " << Cs << "\n"
<< "Cm = " << Cm << "\n";
// breakdown
for (auto b : pbl->bodies)
{
outfile << "# Body - " << b->groups[0]->tag->name << " (" << b->groups[0]->tag->elems.size() << " elements)" << std::endl;
outfile << std::fixed << std::setprecision(12)
<< "Cl = " << b->Cl << "\n"
<< "Cd = " << b->Cd << "\n"
<< "Cs = " << b->Cs << "\n"
<< "Cm = " << b->Cm << "\n";
}
outfile.close();
std::cout << "done." << std::endl;
}
/**
......
......@@ -80,10 +80,12 @@ public:
virtual Status run();
void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
int getRow(size_t i) const { return rows[i]; };
protected:
void computeFlow();
void computeLoad();
};
} // namespace dart
#endif //WSOLVER_H
#endif // WSOLVER_H
......@@ -112,10 +112,10 @@ void Wake::connectNodes()
// Create the node map
for (auto nUp : nodesUp)
{
//std::cout << "MASTER: " << nUp->pos(0) << ',' << nUp->pos(1) << ',' << nUp->pos(2) << std::endl;
// std::cout << "MASTER: " << nUp->pos(0) << ',' << nUp->pos(1) << ',' << nUp->pos(2) << std::endl;
for (auto nLw : nodesLw)
{
//std::cout << " SLAVE: " << nLw->pos(0) << ',' << nLw->pos(1) << ',' << nLw->pos(2) << std::endl;
// std::cout << " SLAVE: " << nLw->pos(0) << ',' << nLw->pos(1) << ',' << nLw->pos(2) << std::endl;
if ((nUp->pos - nLw->pos).norm() <= 1e-14)
{
nodMap[nUp] = nLw;
......
......@@ -59,4 +59,4 @@ private:
} // namespace dart
#endif //WWAKE_H
#endif // WWAKE_H
......@@ -67,4 +67,4 @@ private:
} // namespace dart
#endif //WWAKEELEMENT_H
#endif // WWAKEELEMENT_H
......@@ -41,4 +41,4 @@ public:
};
} // namespace dart
#endif //WWAKERESIDUAL_H
#endif // WWAKERESIDUAL_H
......@@ -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')
......
......@@ -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,22 +80,22 @@ def main():
# visualize solution and plot results
floD.initViewer(pbl)
tboxU.plot(Cp[:,0], Cp[:,3], {'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[:,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('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))
......
......@@ -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], {'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})
tboxU.plot(Cp_ref[:,0], Cp_ref[:,3], {'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})
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()
......
......@@ -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,18 +80,18 @@ def main():
# visualize solution and plot results
floD.initViewer(pbl)
tboxU.plot(Cp[:,0], Cp[:,3], {'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[:,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, 13, 3, forceabs=True))
tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.89, 5e-2))
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')
......
......@@ -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, sfx=''):
def write_slices(mshName, ys, wId, sfx=''):
"""Write slice data for each ys coordinate along the span (only works with VTK)
Parameters
......@@ -138,10 +137,14 @@ def writeSlices(mshName, ys, wId, 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,:]))
tU.write(data, mshName+sfx+'_slice_'+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)
......@@ -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('agard445_slice_1.dat')
tbxU.plot(data[:,3], data[:,4], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': 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
......
......@@ -84,7 +84,7 @@ def cfg():
}
def main():
from dart.api.core import initDart
from dart.api.core import init_dart
from dart import utils
from fwk import Timers
from fwk.testing import CTest, CTests
......@@ -93,10 +93,8 @@ def main():
# pre
tms = Timers()
tms['pre'].start()
_dart = initDart(cfg(), scenario='aerodynamic', task='analysis')
msh = _dart['msh']
sol = _dart['sol']
wrt = _dart['wrt']
_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()
......@@ -106,10 +104,10 @@ def main():
# post
tms['pos'].start()
if hasVTK:
utils.writeSlices(msh.name, [3.81973, 5.8765, 8.2271, 11.753, 14.6913, 17.6295, 19.0986, 21.4492, 24.9751, 27.91338], 12)
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=1) # read
plt.plot(data[:,3], data[:,4], lw=2) # plot results
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)
......
......@@ -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('oneraM6_slice_2.dat')
tbxU.plot(data[:,3], data[:,4], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': 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
......
Subproject commit 7b5a21b0c5fc3885a11e98ebc2e5acbd19bc55bf
Subproject commit c11ff6ee646765a19f9d43496d354e939a84b408
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# This script runs "clang-format" recursively on all C/C++ files
import sys
import os
import fnmatch
import re
import subprocess
def all_files(root,
patterns='*',
skips='*.git*;*build*',
single_level=False,
yield_folders=False):
# self.checkPath(root)
patterns = patterns.split(';')
skips = skips.split(';')
for path, subdirs, files in os.walk(root):
if yield_folders:
files.extend(subdirs)
files.sort()
for name in files:
for pattern in patterns:
if fnmatch.fnmatch(name, pattern):
fullname = os.path.join(path, name)
ok = True
for skip in skips:
if fnmatch.fnmatch(os.path.relpath(fullname, start=root), skip):
ok = False
if ok:
yield fullname
break
if single_level:
break
def main():
# loop over all files and format them
encs = {}
for f in all_files(os.getcwd(), patterns='*.cpp;*.c;*.h;*.hpp'):
# print(f)
cmd = ['clang-format-10', "-style=file", "-i", f]
retcode = subprocess.call(cmd)
if retcode != 0:
print(f'ERROR: retcode = {retcode}')
break
if __name__ == "__main__":
# print('running format_code.py...')
main()
# 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.
# Add source dir
ADD_SUBDIRECTORY( src )
ADD_SUBDIRECTORY( _src )
# Add test dir
MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests)
# Add to install
INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py
${CMAKE_CURRENT_LIST_DIR}/coupler.py
${CMAKE_CURRENT_LIST_DIR}/utils.py
DESTINATION vii)
INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api
${CMAKE_CURRENT_LIST_DIR}/interfaces
DESTINATION vii)