% FGS forward solver default parameters overloading from sources
% P = fgsp(P)
%
% .tolF               Stopping tolerance for root finding algorithm
% .kmax               Maximum iterations for root finding algorithm
% .mkryl              Maximum krylov subspace dimension
% .algoF              Algorithm for root finding (jfnk,Newton,Picard);
% .epsilon_d          Step size to approximate Jacobian numerically
% .usepreconditioner  If ~=0, use preconditioner for JFNK
% .relax              Apply relaxation on solution
% .show_info_term     Show convergenve info during iterations.
% .robust_on          Allow subiteration with reduced steps in case of NaN in the operator
% .nrobust            Max number of subiterations for robust decrease of dx
% .Prec               Preconditionig matrix
% .workers            Number of workers for parallel computation
% .insrc              Data input source
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function [P] = fgsp(P,varargin)

D = meqp(); % meq default parameters

%% Default Parameters
% NL solver parameters 
D.mkryl = 150; % 
D.nrestarts = 5; % number of maximum randomized restarts of the solver
D.tolF = 1e-10; % default value of tol (minimum tolerance max(abs(F(x))))
D.kmax = 50; % default value of kmax (maximum number of Newton steps)
D.kinner = 10; % maximum inner iterations in Picard algorithm stabilizer  (only for algoNL='Picard')
D.kpic = 50;   % maximum Picard loop iterations (only for algoNL='Picard')
D.algoNL = 'all-nl'; % Plasma convergence method: 'all-nl': entire plasma state is unknown of nonlinear system. 'Picard': solve plasma with inner Picard loop.
D.algoF = 'jfnk'; % 'jfnk' 'Newton' Specify the algorithm to be used for nonlinear solver
D.algoGMRES = 'mgsaya'; % Modified Gram-Schmidt with Givens Rotations
D.epsilon_d = 1e-8; % default values of epsilon_d (finite step size used to approximate the Jacobian)
D.usepreconditioner = 0; % Preconditioning matrix
D.relax = 0; % Apply relaxation on solution
D.agcon = {'Ip','bp','qA'}; % default constraints (see meqagcon for details)
D.LYall = true; % return all time slices even if not converged
D.tolPicard = 1e-12; % Picard loop convergence tolerance, used for algoNL='Picard'

% Initialization source
D.insrc = 'fbt'; % 'fbt' 'liu'
D.shot = NaN;

%% Overwrites parameters from input and set defaults
 for k = 1:2:length(varargin), P.(varargin{k}) = varargin{k+1}; end
 for k = fieldnames(D)'
  if ~isfield(P,k{1}), P.(k{1}) = D.(k{1}); end
 end

%% order
P = orderfields(P);

%% Checks
assert(~P.usepreconditioner || ~strcmp(P.algoNL,'Picard'),'usepreconditioner must be false if algoNL=''Picard''')
end
