function EqnBIC = EqnBIC_2phase_StanfordPerm
%% EqnBIC: 2-phase flow with Stanford permeability on [0,5]x[0,1]
% James Liu, ColoState, Victor Ginting, UWYO; 2012/07--2018/12 

EqnBIC = struct('fxnK',@fxnK, 'fxnKs',@fxnKs, 'fxnf',@fxnf, ...
  'fxnpD',@fxnpD, 'fxnuN',@fxnuN, ...
  'fxnp',@fxnp, 'fxnpg',@fxnpg, ... 
  'fxnTotalMobility',@fxnTotalMobility, ...
  'fxnFractionalFlow',@fxnFractionalFlow);

% Darcy: Permeability as 2x2 SPD matrix 
function K = fxnK(pt)
  tmp = load('2phase_StanfordPerm_coeff100x100log.txt');
  PERMK = tmp(2:101,2:101);
  PERMK = exp(2.0*PERMK);
  % Temporarily 
  x = pt(:,1);  y = pt(:,2);
  I = 100 - floor(100*y);
  J = floor(100*x/5) + 1;
  NumPts = size(pt,1);
  K = zeros(NumPts,2,2);
  for k=1:NumPts
    K(k,1,1) = PERMK(I(k),J(k));
    K(k,2,2) = PERMK(I(k),J(k));
  end
end

% Darcy: Permeability as scalar 
function Ks = fxnKs(pt)
  tmp = load('2phase_StanfordPerm_coeff100x100log.txt');
  PERMK = tmp(2:101,2:101);
  PERMK = exp(2.0*PERMK);
  % Temporarily 
  x = pt(:,1);  y = pt(:,2);
  I = 100 - floor(100*y);
  J = floor(100*x/5) + 1;
  NumPts = size(pt,1);
  Ks = zeros(NumPts,1);
  for k=1:NumPts
    Ks(k) = PERMK(I(k),J(k));
  end
end

% Darcy: Source (the right-hand function in the Darcy equation)
function f = fxnf(pt)
  x = pt(:,1);  % y = pt(:,2);
  f = zeros(size(x));
end

% Darcy: Dirichlet boundary conditions 
function pD = fxnpD(pt) 
  x = pt(:,1);  % y = pt(:,2);
  bndry0 = (x==0);
  bndry1 = (x==5);
  pD = 1*bndry0 + 0*bndry1;
end

% Darcy: Neumann boundary conditions 
function uN = fxnuN(pt)
  % x = pt(:,1);  
  y = pt(:,2);
  bndry = (y==0 | y==1);
  uN = 0*bndry;
end

% Darcy: No known exact solution for pressure or set it to zero
function p = fxnp(pt)
  x = pt(:,1);  % y = pt(:,2);
  p = zeros(size(x));
end

% Darcy: No known exact solution for pressure gradient or set it to zero
function pg = fxnpg(pt) 
  % x = pt(:,1);  y = pt(:,2);
  pg = zeros(size(pt));
end

% 2-phase: Total Mobility: To be multiplied to permeability 
function tb = fxnTotalMobility(s)
  tb = s.*s + (1-s).*(1-s)/5;
end

% 2-phase: Fractional Flow: A part of the saturation advection equation 
function ff = fxnFractionalFlow(s)
  ff = (s.^2)./fxnTotalMobility(s);
end

end