classdef call_fge_from_fgs_and_liu < meq_test
 % Test the following calls for singlet and droplet case
 % fgs from liu
 % fge from liu
 % fge from fgs
 % This indirectly check that meqxconvert and fgex works properly
 %
 % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
 
 properties (TestParameter)
  shot =     struct('TCV_std', 69194              ,'TCV_droplet_stable', 68466);
  tin  =     struct('TCV_std', 0.8                ,'TCV_droplet_stable', 2.43);
  idoublet = struct('TCV_std', false              ,'TCV_droplet_stable', true);
  bfp =      struct('TCV_std', [1 2]              ,'TCV_droplet_stable', [0 1]);
  agcon =    struct('TCV_std', {{'Ip','qA', 'bp'}},'TCV_droplet_stable', {{'Ip'}}); % Doublet only work with bfp = [0 1]
 end
 
 properties
  nt = 5;
  dtliu = 1e-5;
  verbosity = 0;
  tok = 'tcv';
  cde = 'cde0Dss';
  npq = 40; % meqxconvert does not handle difference in nQ between codes for the moment
  kmax = 200; % Max iteration for fgs to reach convergence
 end
 
 methods(Test,TestTags={'TCV'},ParameterCombination='sequential')
  
  function liu_from_mds_to_fgs_and_fge(testCase,shot, tin, idoublet, agcon, bfp )
   %% Compute equilibrium reconstruction
   t = tin + (0:testCase.dtliu:testCase.nt*testCase.dtliu);
   tend = t(end);
   if idoublet
    inPar = liup_doublet('idoublet',true,'algoNL','picard','itert',150,'psichco',1e-10,...
     'stabz',1,'wdz',[0;0;Inf]);
   else
    inPar = {};
   end
   [Lliu,~,LYliu] = liu('tcv',shot,t(1),inPar{:},'npq',testCase.npq,'bfp',bfp,'debug',testCase.verbosity);
   % Replicate liuqe structure in order to check if meqxconvert is handled properly
   LYliu = meqxk(LYliu,ones(1,numel(t)+1)); 
   strout = sprintf('LIUQE failed shot= %d', shot);
   testCase.assertTrue(~isempty(LYliu),strout)
   
   %% Get fgs and fge parameters
   % Common parameters
   Params = {'tcv',shot,t,...
    'bfct',Lliu.P.bfct, ...
    'bfp', bfp,...
    'agcon', agcon,...
    'izgrid',true,...
    'npq', Lliu.npq ,... % meqxconvert does not handle cases with different nQ  
    'idoublet', idoublet,...
    'debug', testCase.verbosity};
   
   % fgs parameters
   [Lfgs] = fgs(Params{:});
   
   % Add specific fge parameters
   [Lfge] = fge(Params{:},'cde',testCase.cde);
   
   %% fgs run from liu
   LXfgs = meqxconvert(Lliu,meqxk(LYliu,1),Lfgs);
   LYfgs = fgst(Lfgs,LXfgs,  'kmax', testCase.kmax);
   strout = sprintf('fgs from liu failed shot= %d', shot);
   testCase.assertTrue(all(LYfgs.isconverged),strout)
   
   % Compute one initial condition
   % This allows to compute the linerization and hence the preconditioner only once and speed up the test
   LXfge = fgex('tcv',t,Lfge,meqxconvert(Lliu,meqxk(LYliu,1),Lfge));
   Ltmp = fgel(Lfge, meqxk(LXfge,1)); % Store the linearization for preconditioner
   for ii= [true,false]
    % Test both stationary state and non stationary state initial condition for fge
    Lfge.P.ssinit = ii;
    Lfge = fgec(Lfge.P,Lfge.G);
    Lfge.lin = Ltmp.lin;
    %% fge run from liu
    LXfge = fgex('tcv',t,Lfge,meqxconvert(Lliu,meqxk(LYliu,1),Lfge));
    LYfge = fget(Lfge, LXfge);
    strout = sprintf('fge from liu failed shot= %d, ssinit = %d', shot, ii);
    testCase.assertTrue(all(LYfge.isconverged),strout)
    
    %% fge run from fgs
    LXfge = fgex('tcv',t,Lfge,meqxconvert(Lfgs,meqxk(LYfgs,1),Lfge));
    LYfge = fget(Lfge, LXfge);
    strout = sprintf('fge from fgs shot= %d, ssinit = %d', shot, ii);
    testCase.assertTrue(all(LYfge.isconverged),strout)
   end
   
  end
  
 end
end