3. Solving the flow
The third step consists in creating the solver and solving the flow. The solver is currently restricted to steady flow computations. Two types of solvers are available: direct, to compute the flow, and adjoint, to compute the flow sensitivities.
3.1 Using the direct solver
First, the linear (inner) solver that will solve the linear set of equations at each nonlinear (outer) step must be chosen. Five linear solvers are available: MUMPS, Intel's MKL PARDISO, Intel's MKL DSS, Eigen's ILU GMRES and Eigen's SparseLU. Note that MUMPS, DSS and PARDISO are only interfaces, and therefore require external packages (MUMPS and MKL) to be installed, while GMRES and SparseLU are implemented in Eigen and thus require no additional packages. To achieve the best performance, PARDISO or GMRES should be used. DSS is an alternative interface of PARDISO and can also be used. MUMPS will also yield good performance but is less efficient than PARDISO or GMRES. SparseLU is the default direct solver and should only be used for small 2D cases.
The linear solvers can be imported direclty from the tbox
module as,
linsol = tbox.Pardiso() # or tbox.Dss(), or tbox.Mumps(), or tbox.SparseLu()
If GMRES is used, some parameters must be defined,
linsol = tbox.Gmres(fill, drop, rst, tol, mit)
where fill
is the fill-in factor and drop
is the drop tolerance of the ILU preconditioner, and where rst
is the number of restar, tol
is the tolerance and mit
is the maximum number of iterations of the GMRES solver. To obtain best performance for flow solutions, use fill=2
, drop=1e-6
, rst=50
, tol=1e-5
and mit=200
.
The second step consists in choosing the nonlinear solution algorithm: Picard or Newton. The Newton solver is implemented for transonic flows and will yield the 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, rTol=1e-6, aTol=1e-8, mIt=25, nthrds=1, vrb=1)
and the Newton solver, with the dart::Newton class as,
solver = dart.Newton(pbl, linsol, rTol=1e-6, aTol=1e-8, mIt=25, nthrds=1, vrb=1)
where rTol
is the tolerance on the relative residual, aTol
is the tolerance on the absolute residual, mIt
is the maximum number of (outer) iterations, nthrds
is the number of threads (shared memory implementation) used to assembles the finite element matrices during the solver execution, and vrb
controls the level of verbosity. A level of 0 will not display any information about the solution process. A verbosity of 1 will show the first and the last iteration, and indicate the solver convergence status. A verbosity of 2 will also display all the intermediate iterations. Higher verbosity levels will display other information such internal timers or matrix sizes, mainly for testing purposes. Note that the solver will stop as soon as the first stopping criterion is met (relative or absolute tolerance, or maximum number of iterations).
A supplementary under-relaxation parameter can also be set for the Picard solver,
solver.relax = .7 # default value
Several supplementary parameters can also be set for the Newton solver,
solver.lsTol = 1e-6 # default values
solver.maxLsIt = 10
solver.avThrsh = 1e-2
where lsTol
and maxLsIt
are the tolerance on the residuals 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.
Once the solver has been initialized, it can be executed with,
status = solver.run()
where status
is an integer indicating the convergence of the sovler (0 = converged, 1 = maximum number of iterations reached, 2 = diverged).
3.3 Using the adjoint solver
An discrete adjoint method can be used if the user needs to compute the sensitivities of the flow functionals with respect to some design variables. The solver is created and executed using the dart::Adjoint C++ class as,
adjoint = dart.Adjoint(solver, morpher)
adjoint.run()
where solver
MUST be a Newton solver and morpher
is the mesh morpher. Note that, if an iterative solver such as GMRES is used as inner solver for the flow or mesh morphing solutions, its tolerance will be automatically lowered to 1e-12.
The adjoint solver can currently compute the (total) gradients of the lift and drag coefficients with respect to the angle of attack and to the motion of surface mesh node coordinates. All the gradients are computed using discretized analytical formula.