classdef fge_CDE_test < meq_test
  % Tests of FGE with current diffusion equation
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbosity = 0; % Make plot 1 to make the plot

    tok = 'ana';
    shot = 1;
    tstart = 0;
    tend = 0.01;
    dt = 5e-4;
    nu = 30; %
    wcompa = 0; % no regularization in compensation currents
    
    L, Lss, LX, LXss;
  end
  
  properties (ClassSetupParameter)
   cdetype = {'cde0Dss', 'OhmTor', 'OhmTor_rigid'};
  end
  
  methods (TestClassSetup)
    function setup_LLX(testCase, cdetype)
      t = testCase.tstart:testCase.dt:testCase.tend;
      
      [testCase.L, testCase.LX] = fge(testCase.tok,...
        testCase.shot,t,...
        'ssinit',false,...
        'debug',testCase.verbosity,... % debug ouput level of simulation
        'selu','e','nu',testCase.nu,...
        'izgrid',true,...
         'wcompa',testCase.wcompa);
      
      [testCase.Lss, testCase.LXss]  = fge(testCase.tok,...
        testCase.shot,t,...
        'ssinit',false,...
        'debug',testCase.verbosity,... % debug ouput level of simulation
        'selu','e','nu',testCase.nu,...
        'izgrid',true,...
        'cde',cdetype,...
        'wcompa',testCase.wcompa);
    end
  end
  
  methods (Test, TestTags = {'fge'})
    
    function test_initial_conditions(testCase)
      L = testCase.Lss; LX = meqxk(testCase.LXss,1);

      %% LX with steady-state initial condition
      L.P.ssinit = true;
      LX2 = fgex(L.P.tokamak,LX.t,L,LX);

      %% Same but with focus on LY region only for field compensation
      L.P.compregion = 'Iy';
      LX3 = fgex(L.P.tokamak,LX.t,L,LX);
      
      if testCase.verbosity
        clf;
        meqcompare(L,LX,LX2,LX3); 
        % inner flux surfaces are not touched, while vessel currents added.
      end
      
      % check that new equilibrium Iy is close to previous one as far as plasma is concerned.
      FAB = LX.FA-LX.FB;
      testCase.verifyLessThan(norm((LX.Iy-LX2.Iy)/LX.Ip),1e-4,'Iy mismatch');
      testCase.verifyLessThan(norm((LX.Iy-LX3.Iy)/LX.Ip),1e-4,'Iy mismatch');
      testCase.verifyLessThan(norm(LX.FA-LX2.FA)/FAB,1e-3,'FA mismatch');
      testCase.verifyLessThan(norm(LX.FA-LX3.FA)/FAB,1e-3,'FA mismatch');
      testCase.verifyLessThan(norm(LX.FB-LX2.FB)/FAB,1e-3,'FB mismatch');
      testCase.verifyLessThan(norm(LX.FB-LX3.FB)/FAB,1e-3,'FB mismatch');
    end
    
    function test_steady_state(testCase)
      % Check that steady-state conditions are maintained in time-dependent
      % simulation
      L = testCase.Lss; LX = testCase.LXss;
      L.P.ssinit = true;
      LX = fgex(L.P.tokamak,LX.t,L,LX); % compute new initial condition
      
      % FGE run
      LY = fget(L,LX);
      
      % Check that equilibrium is maintained
      testCase.verifyEqual(LY.Ip,LX.Ip,'AbsTol',10); 
      testCase.verifyEqual(LY.Ia,LX.Ia,'AbsTol',10);
      testCase.verifyEqual(LY.Iu,repmat(LX.Iu,1,numel(LX.t)),'AbsTol',5);
      testCase.verifyEqual(LY.zA,repmat(LX.zA,1,numel(LX.t)),'AbsTol',1e-5);
      testCase.verifyEqual(LY.rA,repmat(LX.rA,1,numel(LX.t)),'AbsTol',1e-5);
      testCase.verifyEqual(LY.FB-LY.FA,repmat(LX.FB-LX.FA,1,numel(LX.t)),'AbsTol',1e-5);
    end
    
    function test_fge_CDE_full_vs_ss(testCase)
      % Test open-loop FGE with CDE with various constraints   
   
      testCase.assumeFail('Ignoring test for now until we find a good test for SS vs Full CDE')
      
      LX = testCase.LX; %#ok<*PROP>
      Lss = testCase.Lss;
      L = testCase.L;
      
      L.P.debug = 1;
      
      LX.Ip  = repmat(LX.Ip(1),  1, numel(LX.t));
      LX.rBt = repmat(LX.rBt(1), 1, numel(LX.t));
      LX.bp  = repmat(LX.bp(1),  1, numel(LX.t));
      LX.qA  = repmat(LX.qA(1),  1, numel(LX.t));
      
      LYss = fget(Lss,LX);
      LY   = fget(L,LX);
      
      testCase.assertTrue(all(LY.isconverged),'FGE did not converge');
      testCase.verifyEqual(LY.Ip(end),LYss.Ip(end),'RelTol',1e-2,'Ip(end) tolerance exceeded')
      testCase.verifyEqual(LY.Ip(1),LYss.Ip(1),'Ip(1) tolerance exceeded')
      testCase.verifyEqual(LY.zA,LYss.zA,'RelTol',1e-3,'zA tolerance exceeded')
      testCase.verifyEqual(LY.rA,LYss.rA,'RelTol',1e-3,'rA tolerance exceeded')
      
      if testCase.verbose
        clf;
        meqplotevo(L,LY,LYss,'legstr',{'full','SS'});
        drawnow;
      end
    end
    
  
  end
end