Skip to content
Snippets Groups Projects

use

Last edited by acrovato

Use

DART is designed to be used in two main fashions: starting from a scratch python file or through automated python scripts. The first method allows the user to customize his use of DART, while the second is destined to users who wants to use standard computing procedures, such as computing a wing polar or the angle of attack to reach a specified lift coefficient. Both methods follow the same flowchart and examples are given under tests (first method), and under scripts and cases (second method).
Note the tests under tests uses default methods. These methods are used to reduce to amount of copy-pasted code and are not meant to be general. They should not be used except for creating new tests, similar to those already existing.

Starting from scratch

This section describes the first method. A short description of the DART C++ classes can be found in C++ classes.

Plan

  1. Pre-processing
    1.1 Meshing
    1.2 Defining the problem
    1.2.1 Defining the fluid domain
    1.2.2 Adding the initial and boundary conditions
    1.2.3 Adding the Kutta condition
    1.2.4 Adding the boundaries of interest of the fluid domain
  2. Solving
    2.1 Initializing the solver
    2.2 Running the solver
  3. Post-processing
    3.1 Saving the results
    3.2 Computing derived quantities

0. Import
To use DART through python, some classes are compiled through swig. The resulting python classes (modules) need to be imported in the current environment. math or numpy contain useful mathematics constant/functions. Since math is a default module of python while numpy is not, this example will use math. fwk contains some useful functionalities, such as the timers. tbox contains the basic features of an FE environment. Finally, dart contains the features related to the solver.

import math
import fwk
import tbox as tbox
import tbox.gmsh as gmsh
import dart

1. Pre-processing
Before solving the problem, the user must first define it.

1.1 Meshing
The first step consists in creating a geometry and meshing it. Currently, the mesh can be given directly or the geometry file can be given with options. In both cases, the only supported format is the native gmsh format (.msh or .geo).

msh = gmsh.MeshLoader('path/to/mshfile','caller/path').execute()

or

pars = {key : value, ...}
msh = gmsh.MeshLoader('path/to/geofile','caller/path').execute(**pars)

where the 'caller/path' and 'path/to/gmshfile' will be concatenated to retrieve the gmsh file. In the second scenario, pars is a python dictionary containing a the parameters to be passed to gmsh. Note that the key is ALWAYS a string and that, if the first character is a - (dash), the key will be considered as a gmsh option instead. This returns the msh (tbox::MshData C++ class) containing the mesh data structure that can be further used.
For potential flows to be lifting, the potential value must be multivalued along a cut inside the computational domain. For aerodynamics applications, this cut is referred to as wake and extends from the trailing edge of any lifting configurations to the downstream boundary of the domain. In order to assign a discontinuous value to the potential on the wake, all nodes belonging to this surface must be duplicated. This is done by using the tbox::MshCrack C++ class,

mshCrck = tbox.MshCrack(msh, dim)
mshCrck.setCrack('wake')
mshCrck.addBoundaries(['boundary', ...])

where dim is the dimension of the problem (2 or 3), 'wake' is the name of the physical group of the mesh containing the wake surface, and 'boundary', ... are the name of the physical groups of the mesh having elements touching the wake surface. These groups typically include the fluid domain, the airfoil, wing or and fuselage, and the downstream and symmetry boundaries. Two methods are implemented to detect if the elements touching the wake are above or below it. Either the physical groups 'boundary' and 'boundary_' are present in the mesh, and the algorithm will consider that the elements of 'boundary' are above and that those of 'boundary_' are below, or only the group 'boundary' is present, and the algorithm will try to detect automatically were the element is. This is done by comparing the orientation of the normal of the closest wake element to the vector linking the center of gravity of the closest wake element to the considered element. The wake normals must be pointing in the positive (y in 2D and z in 3D) direction. Note that, in the first case, the group 'boundary_' will be merged within the group 'boundary'.
For a 3D flow, the potential at the wake tip (i.e., the edge of the wake downstream of the wingtip) will be continuous. This is enforced by not duplicating the nodes in this line,

mshCrck.setExcluded('wakeTip')

where 'wakeTip' is the name of the physical group of the mesh containing the wake free edge. The utility can now be executed, and directly deleted to save memory,

mshCrck.run()
del mshCrck

Since the execution will duplicate the nodes on the wake surface and modify the elements touching the wake, the existing mesh file must be overwritten. This is accomplished with the tbox::GmshExport class,

tbox.GmshExport(msh).save(msh.name)

1.2 Defining the problem
The next step is to define a problem which will manage the fluid domain, the boundary conditions, etc.

pbl = dart.Problem(msh, dim, alpha, beta M_inf, S_ref, l_ref, x_ref, y_ref, z_ref)

where alpha and beta are the angles of attack and sideslip of the oncoming flow, M_inf is the freestream Mach number, S_ref is the reference area, l_ref is the reference length, and x_ref, y_ref and z_ref are the coordinates of the reference point. The reference quantities are used to compute the aerodynamic loads coefficients.

1.2.1 Defining the fluid domain
If the flow is considered incompressible, M_inf must be set to 0, and the fluid domain is set with help of the C++ classes dart::Medium, dart::F0El and dart::F0Ps as,

pbl.set(dart.Medium(msh, 'field', dart.F0ElRhoL(), dart.F0ElMachL(), dart.F0ElCpL(), dart.F0PsPhiInf(dim, alpha, beta)))

where 'field' is the name of the physical group of the mesh containing the element belonging to the fluid domain. If the flow is compressible, M_inf must be different from 0, and the fluid domain is set as,

pbl.set(dart.Medium(msh, 'field', dart.F0ElRho(M_inf, M_crit), dart.F0ElMach(M_inf), dart.F0ElCp(M_inf), dart.F0PsPhiInf(dim, alpha, beta)))

where M_crit is the square of the critical Mach number above which the density will be clamped to avoid reaching near-void values. Note that, in both cases, beta can be omitted if it is zero.

1.2.2 Adding the initial and boundary conditions
The intial condition is added with the dart::Initial class as,

pbl.add(dart.Initial(msh, 'field', dart.F0PsPhiInf(dim, alpha, beta)))

Dirichlet boundary conditions can be added with the dart::Dirichlet class as,

pbl.add(dart.Dirichlet(msh, 'dirichlet', dart.F0PsPhiInf(dim, alpha, beta)))

where 'dirichlet' is the name of the physical group of the mesh containing the elements belonging to the boundary onto which the Dirichlet condition will be enforced.
Freestream (Neumann) boundary conditions can be added with the dart::Freestream class as,

pbl.add(dart.Freestream(msh, 'freestream', dart.F1ElVi(dim, alpha, beta)))

where 'freestream' is the name of the physical group of the mesh containing the elements belonging to the boundary onto which the Freestream condition will be enforced. In all the cases, beta can be omitted if it is zero. Note that the downstream boundary of the domain MUST ALWAYS be assigned a Freestream boundary condition since it will intersect a wake. Note that Freestream are usually preferred to Dirichlet boundary conditions for flow because they allow to use a smaller domain. In 2 dimensions, farfield boundaries are typically placed 5 chord lengths away from the body if Freestream boundary conditions are used, while they should be placed 20 chord lengths away if Dirichlet conditons are used. However, using only Freestream might cause the matrix of the linear system to be singular. It is there HIGHLY ADVISED to use a Dirichlet boundary condition on the upstream boundary. Note that zero mass flux boundary conditions are naturally enforced by the FEM. The slip boundary conditions hence do not require to be explicitly added.

1.2.3 Adding the Kutta condition
A Kutta condition must be added to allow lifting configurations to generate aerodynamic loads. For a 2D flow, both the dart::Wake and dart::Kutta classes should be used as,

pbl.add(dart.Wake(msh, ['wake', 'wake_', 'field']))
pbl.add(dart.Kutta(msh, ['te', 'wake_', 'body', 'field']))

where 'wake_' is the physical group of the mesh containing the lower wake elements automatically created by duplicating the wake nodes by tbox::MshCrack, and 'te' is the name of the physical group of the mesh containing THE point element defining the trailing edge of 'body'. Note that the wake surface does not need to be aligned with trailing edge bissector for 2D flows because dart::Kutta adds a supplementary term on the elements touching the trailing edge that correctly computes the local flow direction.
For 3D flows, only the dart::Wake class is used,

wake = dart.Wake(msh, ['wake', 'wake_', 'field', 'teTip'])

where 'teTip' is the physical group of the mesh containing the trailing edge of the lifting configuration AND the free edge of the wake (no contributions will be added to these nodes). Since dart::Kutta cannot be used for 3D flows, the wake surface MUST be aligned with the trailing edge bissector in order to yield consistent results.

1.2.4 Adding the boundaries of interest of the fluid domain
A boundary of interest can be added by using the dart::Boundary class as,

pbl.add(dart.Boundary(msh, ['body', 'field']))

This will allow the computation of the aerodynamic load coefficients on this boundary, as well as further data manipulation, e.g. for fluid-structure interaction problems.

2. Solving
The solver is currently restricted to steady flow problems. Two options are implemented.

2.1 Initializing the solver
The first step consists in choosing the linear (inner) solver that will solve the linear set of equations at each nonlinear (outer) step. Four linear solvers are available: MUMPS, Intel's MKL Pardiso, GMRES and SparseLU. Note that MUMPS and Pardiso are only interfaces, and therefore require the external packages (MUMPS and MKL) to be installed, while GMRES and SparseLU are implemented in Eigen and thus require no external packages. For the best performance, Pardiso or GMRES should be used. MUMPS will also yield acceptable performance. SparseLU is the default direct solver and should only be used for small 2D problems. In order to use Pardiso (or MUMPS), a python module checking if the interface has been compiled must first be imported,

from tbox.solvers import LinearSolver
linsol = LinearSolver().Pardiso() # or LinearSolver().Mumps()

This will initialize the tbox::Pardiso and tbox::Mumps classes. On the other hand, GMRES (tbox::Gmres) and SparseLU (tbox::SparseLu) can be imported directly,

linsol = tbox.Gmres() # or tbox.SparseLu()

If GMRES is used, some parameters can be set to improve performance,

linsol.setFillFactor(2)
linsol.setRestart(50)
linsol.setTolerance(1e-5)

The second step consists in choosing the nonlinear solution algorithm: Picard or Newton. The Newton solver is implemented for transonic flows and will yield best performance, while the Picard solver works only with subcritical flows (subsonic only) and is mainly used for testing/debugging purposes.
The Picard solver can be initialized using the dart::Picard class as,

solver = dart.Picard(pbl, linsol)

and the Newton solver, with the dart::Newton class as,

solver = dart.Newton(pbl, linsol)

Both solvers share numerical parameters that will be initialized to default values when creating the solvers. These parameters are public and can be accessed through the python interface,

solver.nthreads = 1
solver.verbose = 0
solver.relTol = 1e-6
solver.absTol = 1e-8
solver.maxIt = 50

where nthreads is the number of threads (shared memory implementation) used to assembles the FE matrices during the solver execution, verbose controls whether to display more information, relTol is the tolerance on the relative residual, absTol is the tolerance on the absolute residual and maxIt is the solver maximum number of iterations. Note that the solver will stop as soon as the first stopping criterion is met (relative or absolute tolerance, or maximum iterations).
A supplementary under-relaxation parameter also needs to be set for the Picard solver,

solver.relax = .7

Several supplementary parameters are also needed for the Newton solver,

solver.lsTol = 1e-6
solver.maxLsIt = 10
solver.avThrsh = 1e-2

where lsTol and maxLsIt are the tolerance on the residual and the maximum number of iterations for the line search procedure combined to the Newton algorithm, and avThrsh is the threshold on the relative residual below which the numerical viscosity (for shock capturing) will be decreased. Note that there are two line search algorithms available, Bank & Rose or quadratic (default is Bank & Rose), but the option CANNOT be chosen through python. The user must modify the method dart::Newton::run() and recompile the code in order to change the algorithm.

2.2 Running the solver
Once the solver has been initialized, it can be executed with,

solver.run()

3. Post-processing
After the case has been solved, the user can access the results generated by the solver and post-process them. The results can be saved to disk for further processing or directly used in python.

3.1 Saving the results
The results are stored in the solver object (dart::Solver and derived classes). Before they can be written to disk, a gmsh or vtk output utility must first be initialized. This is done with the tbox::GmshExport class as,

mshWriter = tbox.GmshExport(msh)

or, with the tboxVtk::VtkExport class as,

import tboxVtk
mshWriter = tboxVtk.VtkExport(msh)

Note that the module tboxVtk must be imported before the VTK export utility can be used.
The data can then be written to disk with,

solver.save(n, mshWriter)

where n is the iteration number when solving cases in sequence. If only one case is solved, n=0can be used. The save operation will trigger two processes. First, the flow variables (total and perturbation potential), the residual and the derived quantities (density, Mach and pressure coefficient) will be written at the mesh nodes in gmsh (.pos) or VTK (.vtu) format. Second, these variables will be written in ascii (.dat) format on each zone specified as a boundary of interest (with the dart::Boundary class).

3.2 Computing derived quantities
The data can also manipulated directly in python. For example, the pressure coefficient on an airfoil surface (previously defined as bnd = dart.Boundary(msh, ['wing', 'field'])) can be extracted and plotted as,

import dart.utils as dartU
import tbox.utils as tbxU
Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
tbxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'title', True)

where 'title' will be the plot title and the flag True is used to reverse the y-axis.
As another example, the user might also want to extract slices along the span of a wing. This can be done as,

floU.writeSlices(msh.name, [y_i, ...], id)

where y_i is the y coordinate of the i spanwise section and id is the number of the physical group (defined in gmsh) of the wing. The data will then be written in Selig order (i.e., from the trailing edge to the trailing edge, passing by the leading edge) in ascii format with the filename(s) 'slice_i.dat'. They can further be plotted with,

data = tbxU.read('slice_2.dat')
tbxU.plot(data[:,3], data[:,4], 'x', 'Cp', 'title', True)

This will generate a x/c against Cp plot.

Automated python script

This section describes the second method. A brief description of the DART python classes can be found in Python classes. At the present moment, two computing procedures are implemented: polar (to compute several angles of attack) and trim (to reach a specified lift coefficient).

Plan

  1. Defining
  2. Running
    2.1 Configuring the solver
    2.2 Running the script
    2.3 Displaying the results

1. Defining
The fist step consists in defining the parameters (as a python dictionary) that be used to initialize the solver and drive the computation. The mandatory parameters are briefly described in Python classes and examples are given in cases.

2. Running
After the parameters have been defined, the solver can be set up and run, and the results can be generated and post-processed.

2.1 Configuring the solver
The solver is automatically configured when the user creates the object handling a computing procedure. A call to,

import dart.scripts.polar as p
polar = p.Polar(getParam())

or to,

import dart.scripts.trim as t
trim = t.Trim(getParam())

will setup the polar or trim computation procedure, and will initialize the solver through the Config python class. Note that Config is not designed be used as standalone (it merely configures the solver), but should be called in the constructor of any computing procedure.

2.2 Running the script
The script can then be run with,

polar.run()

or,

trim.run()

2.3 Displaying the results
Finally, the results can be displayed with,

polar.disp()

or,

trim.disp()