Smoothing routine to correct negative values in the distribution function after collisions while conserving mass, momentum and energy as far as possible. First, the velocity space is divided into quadratic tiles (allowed tile sizes are determined automatically). If more than one tile size is allowed, the routine loops over all tile sizes starting from the smallest one until an exit condition is met (see below). Inside the loop over tile sizes –> The routine loops over all tiles and checks each tile for negative values. If no negative values are found, the tile is not modified. If negative values are found, this routine minimizes the sum (over the tile) of the quadratic deviation from a shifted Maxwellian (for the density, temperature and flow of the input distribution function) with the constraints of mass, momentum and density conservation, and positivity of f. This is done with quadratic programming using the ActiveSetMethod algorithm by Goldfarb/Idnani (1983). If all negative values have been removed by this optimization and the relative errors of mass, momentum and energy are smaller than 1E7, the loop over tile sizes is interrupted. If there are negative values left, the routine proceeds to the next tile size. In case there are still negative values left at the largest allowed tile size, negative values on a tile are set to 1% of the corresponding shifted Maxwellian and all positive values are shifted such as to conserve at least the total mass of the tile. (If this shift would make a positive value on a tile negative, the accuracy of mass conservation is limited!)
Minimization process Minimize sum_i[f_i  f_(M,i)]^2 on each tile with negative values f_i —> Minimize: 1/2 * x^T.Q.x  d^T.x where : A_1^T.x = b_1 (Equality constraints –> mass, momentum, energy) a_2^T.x >= b_2 (Inequality constraints –> positivity of f) Q is diagonal –> Q = 2*I d^T = +2*f_M A_1^T =  dV1 dV2 ...   vpar1*dV1 vpar2*dV2 ...   v_1^2*dV1 v_2^2*dV2 ...  A_2^T = I b_1^T = (n0,p0,E0) (–> mass, momentum, density) b_2^T = (0,0,...,0) (–> positivity)
Since Q is diagonal, it is trivial to factorize it into Q = R^T.R –> R^T = R = sqrt(2)*I Then we pass only R^1 = 1/sqrt(2)*I to the solver
 Parameters

[in,out]  f_new  Distribution function to be smoothed, real(8) 
[in]  vpar_all  Parallel velocity grid, real(8) 
[in]  vperp_all  Perp. velocity grid, real(8) 
[in]  vol_all  Volume elements (d^3v = 2pi v_perp dv_perp dv_par, real(8) 
[in]  npar  Par. vel. grid size, integer 
[in]  nperp  Perp. vel. grid size, integer 
[in]  mass  Species mass, real(8) 
[in]  n0  Density moment of f_new, real(8) 
[in]  u0  Momentum density (!!!) moment of f_new, real(8) 
[in]  T0  Temperature moment of f_new, real(8) 