classdef (SharedTestFixtures={mexm_fixture}) gszr_test < meq_test
  % test GS solver in vacuum vs Green's functions
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  properties
    L;
    nz = 32; %  64  128  32  64  128};
    nr = 31; %  61  121  46  90  180};
    rel_tol     = 1.1e-4;
    rel_tol_dom = 5e-2;

    verbosity = 0;
  end
  
  properties(TestParameter)
    method = {'gszrmex','gszrmexm'};
  end
  
  properties (MethodSetupParameter)
    machine = {'ana'};
    shot    = {    1};
  end
  
  methods(TestMethodSetup, ParameterCombination='sequential')
    function gsrzTestSetup(testCase,machine,shot)
      testCase.L = liu(machine,shot,1,'nz',testCase.nz,'nr',testCase.nr);
    end
  end
  
  methods(Test, TestTags = {'Unit'})
    function test_gsrz_accuracy(testCase,method)
      % Tests accuracy of gsrz Poisson solver.
      % Compares results of gsrz with pre-computed Green's functions for
      % each coil
      % For coils inside the computational grid, the method is inaccurate
      % since the currents are assigned to discrete grid points.
      
      L = testCase.L; %#ok<*PROPLC>
      
      % Extract Green's functions stored in L.G.Mxa (active coils to grid points) for grid boundary points
      Mba = L.G.Mxa(~L.lxy,:);
      
      % Compute the maximum error of the Poisson solver for the vacuum flux produced by a unit current in each active coil
      eF = zeros(L.G.na,1);
      
      inDomain = false(L.G.na,1); inDomain(intersect(L.kew,1:L.G.na)) = true;
      for k = 1:L.G.na % iterate over coils
        Fb = Mba(:,k); % Boundary flux produced by coil k
        
        Iy = zeros(L.nzy,L.nry); % currents inside computational domain
        
        if inDomain(k)
          Iy(L.kyw) = L.Tww(:,k==L.kew); % Add conductor currents on computational grid
        end
        
        Fx = feval(method,Fb,Iy,L.cx,L.cq,L.cr,L.cs,L.ci,L.co,0);
        eFx = (Fx - reshape(L.G.Mxa(:,k),size(Fx)))./norm(Fx);
        eF(k,1) = max(abs(eFx(:)));
       
        if eF(k,1)>testCase.rel_tol && testCase.verbosity
          [rr,zz] = meshgrid(L.rx, L.zx);
          idx = find(L.G.Twa(:,k));

          figure
          subplot(121)
          contourf(rr,zz,eFx,21); colorbar; axis equal
          title('error')
          
          subplot(122)
          plot(L.G.rw(idx),L.G.zw(idx),'.m',...
            L.G.rl,L.G.zl,'k')
          hold on, grid on, axis equal
          
          contour(rr,zz,Fx,21,'b'); hold on;
          contour(rr,zz,reshape(L.G.Mxa(:,k),size(Fx)),21,'r--')
          title(['flux, coil: ' int2str(k)])
        end       
      end
      
      if any(~inDomain)
        testCase.verifyLessThan(eF(~inDomain), testCase.rel_tol,'Poisson solver accuracy failure');
      end
      if any(inDomain)
        % allow more error here since Green's function will be more precise
        testCase.verifyLessThan(eF(inDomain), testCase.rel_tol_dom,'Poisson solver accuracy failure');
      end
     
      if testCase.verbosity
        figure
        txt = sprintf('%dx%d',L.nzx,L.nrx);
        bar(eF(1:end-1,:)), set(gca,'yscale','log'), legend(txt,'location','best')
        title('Flux error')
      end
    end
  end
  
  
  
end
