negative_f_correction.F90 File Reference

Functions/Subroutines
subroutine	negative_f_correction (f_new, vpar_all, vperp_all, vol_all, npar, nperp, mass, n0, u0, T0)
	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 Active-Set-Method 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 1E-7, 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!) More...

Function/Subroutine Documentation

subroutine negative_f_correction	(	real (kind=8), dimension(nperp,npar), intent(inout)	f_new,
		real (kind=8), dimension(npar), intent(in)	vpar_all,
		real (kind=8), dimension(nperp), intent(in)	vperp_all,
		real (kind=8), dimension(nperp), intent(in)	vol_all,
		integer	npar,
		integer	nperp,
		real (kind=8)	mass,
		real (kind=8)	n0,
		real (kind=8)	u0,
		real (kind=8)	T0
	)

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 Active-Set-Method 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 1E-7, 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)