classdef fge_convergence_tests < meq_test
  % Tests for Time Evolving solver
  % Check convergence of dt of final zA position
  % Check the convergence increasing number of vessel modes
  % Check done for Euler Explicit with preconditioner
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties(TestParameter)
    algoNL = {'Picard','all-nl'};
  end
  
  properties
    tok = 'ana';
    dosave = 0;
    shot = 2;
    tstart = 0;
    tinterval= 5e-4; % time interval length
    dt = [ 5e-5, 1e-5, 5e-6 2e-6];
    prec = 1;
    solvetol = 1e-8;
    nu = 50;
    selu = 'e';
    verbose = 0;
    mkryl = 20;
    iterq = 20;
  end
  
  % --------------------------------------------------------------------
  methods(Test,TestTags={'fge-convergence'})
    
    function test_algoNL_methods(testCase)
      myshot = testCase.shot;
      t = testCase.dt(1)*[0:100-1];
      % run for both methods and compare result
      for ii=1:numel(testCase.algoNL)
        mymethod = testCase.algoNL{ii};
        [L,~,LY(ii)] = fge('ana',myshot,t,'algoNL',mymethod,...
          'usepreconditioner',~isequal(mymethod,'Picard'),...
          'iterq',testCase.iterq); %#ok<AGROW>
        if ii>1
          Ip0 = L.P.b0/L.P.r0/(4e-7*pi);
          testCase.verifyEqual(LY(1).Fx,LY(ii).Fx,'AbsTol',100*L.P.tolF)
          testCase.verifyEqual(LY(1).Iy,LY(ii).Iy,'AbsTol',L.P.tolF*Ip0)
          testCase.verifyEqual(LY(1).bp,LY(ii).bp,'AbsTol',100*L.P.tolF)
          testCase.verifyEqual(LY(1).li,LY(ii).li,'AbsTol',100*L.P.tolF)
          testCase.verifyEqual(LY(1).Ip,LY(ii).Ip,'AbsTol',10*L.P.tolF*Ip0)
        end
      end
    end
    
    function dt_convergence(testCase,algoNL)
      
      % Run simulation at different dt
      dt_list = testCase.dt;
      zAlist = zeros(size(dt_list));
      
      [L] = fge(testCase.tok, testCase.shot,0,...
          'insrc','liu',...
          'tolF',testCase.solvetol,...
          'selu',testCase.selu,...
          'nu',testCase.nu,...
          'debug',testCase.verbose,...
          'algoNL',algoNL,...
          'usepreconditioner',~isequal(algoNL,'Picard'));
      
      LX0 = fgex(testCase.tok,testCase.tstart,L);
      if ~isequal(algoNL,'Picard')
        [L] = fgel(L,LX0); % linearize
      end
      for ii = 1:numel(dt_list)
        t = testCase.tstart + (0:dt_list(ii):testCase.tinterval);
        LX = fgex(testCase.tok,t,L,LX0);
        LX.Va(2,:) = LX.Va(2,:) + 20; % some step input
        LY = fget(L,LX);
        if testCase.verbose>1
          if ii==1; clf; end
          plot(LY.t,LY.zA,'displayname',sprintf('dt=%3.2e',dt_list(ii))); hold on;
          legend('show');
          drawnow;
        end
        zAlist(ii) = LY.zA(end);
      end
      
      % Compute residual with respect to smaller dt solution
      zAresidual = abs(zAlist(1:end-1)-zAlist(end));
      
      % Check that the total residual is decaying 
      testCase.verifyTrue(issorted(zAresidual, 'descend'));
     
      % Check convergence is linear in loglog plot
      logtt=log(dt_list(end-1:-1:1));
      logzA = log(zAresidual(end:-1:1));

      coeff_fit = polyfit(logtt,logzA,1);
      if testCase.verbose>1
        clf;
        plot(logtt,logzA,'x');  hold on;
        plot(logtt, coeff_fit(1)*logtt+ coeff_fit(2),'r--');
        legend('data','fit','location','best');
        ylabel('log error'); xlabel('log tt');
        title(sprintf('slope=%3.2f',coeff_fit(1)),'Interpreter','None')
        drawnow
      end
      testCase.verifyTrue( abs(coeff_fit(1)-1)<0.4,...
        'convergence is not linear in dt');
    end
    
    function eigenmode_convergence(testCase,algoNL)
      t = testCase.tstart + (0:2e-5:5e-4);
      % Test convergence for eigenmode vessel description
      RelTol = sqrt(testCase.solvetol)*10;
      %% basic case with full vessel
      
      [L] = fge(testCase.tok,testCase.shot,t,...
        'insrc','liu',...
        'debug',testCase.verbose);
      
      LX = fgex(testCase.tok,t,L);
      L  = fgel(L,meqxk(LX,1));
      LY = fget(L,LX);
      
      %% eigenmode case with all eigenmodes
      nu_ = L.G.nv; %#ok<*PROP>
      [~,~,LYe] = fge(testCase.tok,testCase.shot,t,...
        'insrc','liu',...
        'tolF',testCase.solvetol,...
        'selu','e','nu',nu_,....
        'debug',testCase.verbose);
      
      % verify equal
      testCase.verifyEqual(LY.Fx,LYe.Fx,...
        'RelTol',RelTol)
      
      if testCase.verbose
        clf;
        plot(LYe.t,LYe.zA,'displayname',sprintf('nu = %d',nu_));
        hold on; drawnow;
      end
      
      for ieigs = {50,30,10}
        nu_ = ieigs{1}; 
        fprintf('Run case with selu = ''e'', nu = %d\n',nu_)
        % try fewer eigenvalues
        [~,~,LYe] = fge(testCase.tok,testCase.shot,t,...
          'insrc','liu',...
          'tolF',testCase.solvetol,...
          'selu','e','nu',nu_,...
          'debug',testCase.verbose); 
        if testCase.verbose
          plot(LYe.t,LYe.zA,'displayname',sprintf('nu = %d',nu_)); drawnow;
        end
      end
      
      if testCase.verbose
        title('zA for various #eigenmodes');
        xlabel('time'); ylabel('zA [m]');
        legend show; legend(gca,'location','best')
      end
    end

    function gmres_accuracy(testCase)
      dt_list = testCase.dt;

      % setup FGE: use no preconditioner and large mkryl so that GMRES
      % performance matters most
      [L,LX0] = fge(testCase.tok, testCase.shot,0,...
        'insrc','liu',...
        'mkryl',20,...
        'usepreconditioner',0,...
        'tolF',testCase.solvetol,...
        'debug',0);

      for ii = numel(dt_list)  % last dt_list should be enough
        fprintf('[ Test %g ]\n\n', dt_list(ii));
        t = testCase.tstart + (0:dt_list(ii):testCase.tinterval);

        LX = fgex(testCase.tok,t,L,LX0);
        [L] = fgel(L,LX0); % pre-compute linearization for re-use
        
        algolist = {'cgs','mgsgiv','mgsaya'};
        telaps = zeros(numel(algolist),1); % init
        for ialgo = 1:numel(algolist) % loop over algorithms to try
          algoGMRES = algolist{ialgo};
          
          % run cgs
          L.P.algoGMRES = algoGMRES;
          start_time=tic;
          LY = fget(L,LX);
          telaps(ialgo) = toc(start_time);
          
          if ialgo==1
            LY0 = LY;
          else
            tol = sqrt(eps);
            testCase.verifyEqual(LY.rB, LY0.rB, 'AbsTol', tol);
            testCase.verifyEqual(LY.zB, LY0.zB, 'AbsTol', tol);
            testCase.verifyEqual(LY.FB, LY0.FB, 'AbsTol', tol);
            testCase.verifyEqual(LY.zA, LY0.zA, 'AbsTol', tol);
            testCase.verifyEqual(LY.FA, LY0.FA, 'AbsTol', tol);
            testCase.verifyEqual(LY.Fx, LY0.Fx, 'AbsTol', tol);
          end
        end
        
        fprintf('-- Comparitive times: CGS: %f[s], MGSGIV: %f[s] MGSAYA: %f[s] [%.2fx, %.2fx]\n',...
          telaps(1), telaps(2), telaps(3), telaps(1)/telaps(2), telaps(1)/telaps(3));
      end
    end
  end
end
