classdef fge_constraint_tests < meq_test
  % Tests FGE with imposed constraints on various quantities
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    verbose = 0;
    tok = 'ana';
    shot = 1;
    tstart = 0;
    tend = 0.005;
    dt = 5e-5;
    tol = 1e-9; % solver tolerance
    nu = 30; %
    L;
  end
  
  properties(TestParameter)
    constraints = {{'Ip','bp','qA'},{'Ip','bp','li'},{'Ip','bt','qA'},{'Ip','bt','li'}};
  end
  
  methods(TestClassSetup)
    function setup_LLX(testCase)
      
      [L] = fge(testCase.tok,...
        testCase.shot,testCase.tstart,...
        'debug',testCase.verbose,... % debug ouput level of simulation
        'tolF', testCase.tol,...
        'selu','e','nu',testCase.nu,...
        'izgrid',true);
      
      testCase.L = L; %#ok<*PROP>
    end
  end
  
  methods(Test,TestTags={'fge'})
    
    function test_fge_constraints(testCase,constraints)
      t = testCase.tstart:testCase.dt:testCase.tend;

      L = testCase.L;   %#ok<*PROPLC>
      L.P.agcon = constraints;
      L = fgec(L.P,L.G); % re-consolidate L
      
      %% generate LX
      LX = fgex(testCase.tok,t,L);
      
      %% perturb constrained quantities
      for iC=1:L.nC
        field = L.agconc{iC,3};
        LX.(field) = LX.(field) + ...
          LX.(field)(1)*0.001*sin(2*pi*300*LX.t);
      end
      L = fgel(L,meqxk(LX,1));
      LY = fget(L,LX); % time stepper
      
      if testCase.verbose
        figure(1); clf;
        for iC=1:L.nC
          subplot(L.nC,1,iC)
          field = L.agconc{iC,3};
          if isfield(LX,field)
            plot(LX.t,LX.(field),'ko'); hold on;
            plot(LY.t,LY.(field),'b.--')
            title(field);
          end
        end
        figure(2); clf; meqplotevo(L,LX,LY);
        drawnow
      end
      
      %% Check that the constraint on any desired scalar quantities are indeed satisfied
      for ii=1:numel(L.P.agcon)
        myfield = L.P.agcon{ii};
        if isfield(LX,myfield)
          target = LX.(myfield);
          if size(target,2)==1
            target = repmat(target,numel(LY.t));
          end
          % specific scalings
          switch myfield
            case 'Ip', s = 1e6;
            case 'bt', s = 1e-3;
            otherwise, s = 1;
          end
          testCase.verifyEqual(LY.(myfield),target,'AbsTol',s*testCase.tol,...
            sprintf('Did not satisfy constraint on %s',myfield));
        end
      end
    end
  end
end
