classdef liu_Newton_vs_Picard_test < meq_test
 % Compare Liuqe solution with Picard iteration vs Newton wrapper on liut for singlet and 1 case of droplet with stable Picard
 %
 % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
 
 properties (TestParameter)
  tok      = struct('TCV_diverted','tcv','TCV_limited_NT','tcv','TCV_droplet_Picard_stable','tcv');
  shot     = struct('TCV_diverted',69185,'TCV_limited_NT',69179,'TCV_droplet_Picard_stable',68466);
  t        = struct('TCV_diverted',   1 ,'TCV_limited_NT',    1,'TCV_droplet_Picard_stable', 2.43);
  idoublet = struct('TCV_diverted',false,'TCV_limited_NT',false,'TCV_droplet_Picard_stable', true);
 end
 
 properties
  verbosity = 0;
  psichco = 1e-11; % High convergence request for Picard iteration
  tol = 5e-7; % This should converge to 1e-14 asking for high enough convergence for both Picard and Newton
  toldzA = 1e-2; % tollerance [m] zA LIUQE Picard -  zA LIUQE Newton with dz = 0;
  % The convergence test in Picard is on zFx - Fx , while in Newton on Iy - tildeI(Iy, Ia, Ie, Bm, Ff, ....)
  % The two threshold are therefore not directly comparable
 end
 
 methods(Test,TestTags={'TCV'},ParameterCombination='sequential')
  function convergence_check(testCase,tok,shot,t, idoublet)
   % Compare the LIUQE solution solved with Picard and with Newton with same parameters
   
   %% Picard
   inPar = {'debug',  testCase.verbosity, ...
     'idoublet', idoublet,...
     'stabz', 1, ... % Doublet default has no dz stabilization hence `stabz = 1` re-introduced here for comparison of doublet case
     'wdz', [0;0;inf], ... % Disable dz in the mantle (no current there yet)
     'iters', 0, ... % For direct comparison liut Picard vs Newton iters = 0 for liut Picard
     'itera', 500, 'psichco',testCase.psichco,... % Strict convergence requirements
     'iterfrz', 500,  'izgrid',true};
   if idoublet
    % Get default parameter multiple domains and overload specific requests
    inPar = liup_doublet(inPar{:});
   end
   
   [LliuP,~,LYliuP] = liu(tok, shot,t, inPar{:});
   
   %% Newton
   inPar = {'debug',  testCase.verbosity, 'izgrid',true, 'idoublet', idoublet,...
    'stabz', 1, ... % Doublet default has no dz stabilization hence `stabz = 1` re-introduced here for comparison of doublet case
    'wdz', [0;0;inf], ... % Disable dz in the mantle (no current there yet)
    'algoNL', 'jfnk'};
   if idoublet
    % Get default parameter multiple domains and overload specific requests
    inPar = liup_doublet(inPar{:});
   end
   
   [LliuN,~, LYliuN] = liu(tok, shot,t,inPar{:} );
   
   testCase.assertTrue(~isempty(LYliuP) ,'LIUQE Picard did not converge')
   testCase.assertTrue(~isempty(LYliuN) ,'LIUQE Newton did not converge')
   % All fitted quantities {dz, Ia, Ie, ag} must converge to get the same equilibrium
   testCase.verifyEqual(LYliuN.dz,LYliuP.dz,'RelTol',testCase.tol,sprintf('abs(LY_Picard.dz - LY_Newton.dz)> tol'));
   testCase.verifyEqual(LYliuN.Ia,LYliuP.Ia,'RelTol',testCase.tol,sprintf('abs(LY_Picard.Ia - LY_Newton.Ia)> tol'));
   testCase.verifyEqual(LYliuN.Iu,LYliuP.Iu,'RelTol',testCase.tol,sprintf('abs(LY_Picard.Iu - LY_Newton.Iu)> tol'));
   testCase.verifyEqual(LYliuN.ag,LYliuP.ag,'RelTol',testCase.tol,sprintf('abs(LY_Picard.ag - LY_Newton.ag)> tol'));
   
   if testCase.verbosity
    clf;  figure; meqplotfancy(LliuN,LYliuN); figure; meqplotfancy(LliuP,LYliuP);
   end
  end
  
  function no_dz_stabilization(testCase,tok,shot,t, idoublet)
   % Compare the solution of Picard with stabilization and Newton with no dz stabilization
   
   %% Picard
   inPar = {'debug',  testCase.verbosity, ...
    'idoublet', idoublet,...
    'stabz', 1, ... % Doublet default has no dz stabilization hence `stabz = 1` re-introduced here for comparison of doublet case
    'wdz', [0;0;inf], ... % Disable dz in the mantle (no current there yet)
    'iters', 0, ... % For direct comparison liut Picard vs Newton iters = 0 for liut Picard
    'itera', 500, 'psichco',testCase.psichco,... % Strict convergence requirements
    'iterfrz', 500,  'izgrid',true};
   if idoublet
    % Get default parameter multiple domains and overload specific requests
    inPar = liup_doublet(inPar{:});
   end
   [LliuP,~,LYliuP] = liu(tok, shot,t, inPar{:});
   
   %% Newton
   inPar = { 'debug',  testCase.verbosity, 'izgrid',true, 'stabz', 0,'idoublet', idoublet, 'algoNL', 'jfnk'};
   if idoublet
    % Get default parameter multiple domains and overload specific requests
    inPar = liup_doublet(inPar{:});
   end
   [LliuN,~, LYliuN] = liu(tok, shot,t, inPar{:});
   
   %% Check solution
   testCase.assertTrue(~isempty(LYliuP) ,'LIUQE Picard did not converge')
   testCase.assertTrue(~isempty(LYliuN) ,'LIUQE Newton did not converge')
   testCase.assertTrue( LYliuN.dz ==0 ,'LIUQE Newton LY.dz~=0 while dz=0 requested')
   % Check that the solution of LIUQE Picard zA ~ LIUQE Newton zA .
   % One does not expect LIUQE Picard zA - LIUQE Newton zA = dz, since LIUQE Picard is not a proper solution of the GS.
   % See F.Carpanese PhD Thesis 2020 section 4.5.1.
   testCase.verifyEqual(LYliuN.zA(1),LYliuP.zA(1),'AbsTol',testCase.toldzA ,...
    'Significant difference between solution LIUQE Picard and LIUQE Newton dz= 0')
   if testCase.verbosity
    clf;  figure; meqplotfancy(LliuN,LYliuN); figure; meqplotfancy(LliuP,LYliuP);
   end
  end
  
 end
end
