use_custom_solver.md

+## 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. Four linear solvers are available: MUMPS, Intel's MKL PARDISO, Eigen's GMRES and Eigen's SparseLU. Note that MUMPS 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 external packages. TO achieve the best performance, PARDISO or GMRES should be used. MUMPS will also yield good performance. 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,
+```python
+linsol = tbox.Gmres() # or tbox.Pardiso(), or tbox.Mumps(), or tbox.SparseLu()
+```
+If you are unsure whether MKL or MUMPS have been compiled, but that you still want to use these solvers, a python module checking if the interface has been compiled and initializing the tbox::Pardiso and tbox::Mumps classes can first be imported,
+```python
+from tbox.solvers import LinearSolver # or from dartflo.tbox.solvers import LinearSolver
+linsol = LinearSolver().pardiso() # or LinearSolver().mumps()
+```
+If GMRES is used, some parameters must be set to obtain best performance,
+```python
+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 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,
+```python
+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,
+```python
+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 amount of information to display. 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,
+```python
+solver.relax = .7 # default value
+```
+Several supplementary parameters can also be set for the Newton solver,
+```python
+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,
+```python
+solver.run()
+```
+
+**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,
+```python
+adjoint = dart.Adjoint(solver, morpher, linsol)
+adjoint.run()
+```
+where `solver` MUST be a Newton solver, `morpher` is the mesh morpher and `linsol` MUST EITHER be a direct solver OR have a small tolerance (typically, 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.