classdef solveF_test < meq_test
  % Tests for general F(x)=0 solver
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

  properties
    tol = 1e-8;
  end
  
  properties(TestParameter)
    algoF = {'jfnk','jfnk','jfnk','Newton'}
    algoGMRES = {'cgs','mgsgiv','mgsaya',''};
  end
  
  methods(Test, TestTags={'Unit'}, ParameterCombination='sequential')
    function test_quadratic(testCase,algoF,algoGMRES)
      
      xSol = [1 -1]'; x0=[0 0]';
      F = @(x,doplot) (x-xSol).^2;

      x = solveF(F,x0,...
        'robust_on',1, 'algoF', algoF,'algoGMRES',algoGMRES,...
        'mkryl', 200, 'Prec', 1, 'relax', 0, 'debug', 1,...
        'kmax', 150, 'tolF', testCase.tol);
      
      testCase.verifyTrue(norm(x-xSol)<sqrt(testCase.tol),...
        sprintf('Norm exceeds error tolerance %4.3e>%4.3e',norm(x-xSol),sqrt(testCase.tol)));
    end
    
    function test_linear(testCase,algoF,algoGMRES)
      A = rand(100); b=rand(100,1);
      x0 = rand(100,1);
      F = @(x,u) A*x-b;
      xSol = A\b;
      
      P = struct('robust_on',1, 'algoF', algoF,'algoGMRES',algoGMRES,...
        'mkryl', 200, 'Prec', 1, 'relax', 0, 'debug', 1,...
        'kmax', 150, 'tolF', testCase.tol');
      
      x = solveF(F,x0,P);
      
      testCase.verifyTrue(norm(x-xSol)<sqrt(testCase.tol),...
        sprintf('Norm exceeds error tolerance %4.3e>%4.3e',norm(x-xSol),sqrt(testCase.tol)));
      
    end

    function test_reduce_step(testCase,algoF,algoGMRES)
      % First step will inevitably imply a division by zero, so the
      % variational step must be decreased in order to converge to the
      % solution.
      epsilon_d = 0.5;
      xSol = 0.1; x0 = 0;
      F = @(x,doplot) (x-xSol)/(x-epsilon_d)/(x+epsilon_d);

      x = solveF(F,x0,...
        'robust_on',1, 'algoF', algoF,'algoGMRES',algoGMRES,...
        'mkryl', 200, 'Prec', 1, 'relax', 0, 'debug', 1,...
        'kmax', 150, 'tol', testCase.tol, 'epsilon_d', epsilon_d);
      
      testCase.verifyTrue(norm(x-xSol)<sqrt(testCase.tol),...
        sprintf('Norm exceeds error tolerance %4.3e>%4.3e',norm(x-xSol),sqrt(testCase.tol)));
    end
  end
  
end
