%% run_2phase_WG_FVM_QuadriQ0RT0.m 
% Main program for a 2-phase flow model (Darcy Eqn. + Saturation Eqn.)
% -- Darcy eqn. (flow) solved by weak Gaerkin (Q0,Q0;RT[0]) finite elements
% -- Transport eqn. (saturation) solved by the finite volume method (FVM)
% James Liu, ColoState; Victor Ginting, UWYO; 2012/07--2018/12 

% clc;  
clear all;  close all;
format compact;  format short e;

%% Setting up the problem 
% For EqnBC_2phase_StanfordPerm: On the rectangular domain [0,5]x[0,1] 
xa = 0;  xb = 5;  yc = 0;  yd = 1;
BndryDescMat = [xa,yc,xb,yc, 0,-1;... 
                xb,yc,xb,yd, 1, 0;... 
                xb,yd,xa,yd, 0, 1;... 
                xa,yd,xa,yc,-1, 0];
EqnBIC = EqnBIC_2phase_StanfordPerm;
BndryCondType = [0;2;1;2;1];

%% Discretization: Mesh generation and preparation 
% For a rectangular or smooth quadrilateral mesh 
n = 100  
nx = n;  ny= n;  delta = 0;  status = 2;
QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbSmth( ...
  xa,xb,nx,yc,yd,ny,delta,status);
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich3(QuadriMesh,BndryDescMat);

%% Setting up quadratures 
GAUSSQUAD = SetGaussQuad(5,25,9);

%% Sampling permeability 
PermK = Darcy_SmplnPerm_QuadriMesh(EqnBIC.fxnK, QuadriMesh, GAUSSQUAD);
PermKs = PermK(:,1,1);

%% Presenting: heterogeneity ratio, permeability profile 
max(max(PermKs)) / min(min(PermKs))
show_QuadriMesh_ScaEm_ClrImg(QuadriMesh, PermKs, ...
  11, 'Permeability');

%% 2-phase flow: Simulation with operator splitting: WG + FVM 
% WG(Q0,Q0;RT[0]) for Darcy;  FVM for transport;  On same quadri.mesh 
% Setting an initial condition for saturation 
sat = zeros(nx+2,ny+2); 
% Setting boundary conds. for saturation: Inflow on the left boundary 
sat(2:ny+1,1) = 1.0; 
% Setting time steps: Coarse-scale for flow; Fine-scale for saturation 
finaltime = 10; 
cNstep = 20;  % Number of coarse time steps 
cDt = finaltime / cNstep;
fNstep = 50;  % Number of fine time steps  
fDt = cDt / fNstep;
fact = fDt / max(QuadriMesh.area);

%% Time-marching (Coarse-scale for flow; Fine-scale for saturation) 
tyme = 0.0; 
% Coarse-scale time-marching 
satur = zeros(cNstep,nx,ny);  % A 3-dim array for storing data 
for itc = 1 : cNstep
  itc
  sati = sat(2:nx+1,2:ny+1);
  LambdaPermK = reshape(EqnBIC.fxnTotalMobility(sati), nx*ny, 1);
  PermK(:,1,1) = LambdaPermK .* PermKs;
  PermK(:,2,2) = LambdaPermK .* PermKs;
  % 
  % Darcy: WG(Q0,Q0;RT[0])Quadri: Assembling and Solving 
  [sln] = Darcy_WG_QuadriQ0Q0RT0_AsmSlv( ...
    EqnBIC, BndryCondType, QuadriMesh, PermK, GAUSSQUAD);
  [NumerPresEm, NuemrVelCofRT0, NumerFlux] = ...
    Darcy_WG_QuadriQ0Q0RT0_PresVelFlux(QuadriMesh, PermK, sln);
  % 
  % Fine-scale time-marching 
  for itf = 1 : fNstep 
    fc = AdvecSatur_FVM_QuadriQ0(EqnBIC, nx, ny, NumerFlux, sat);
    sat(2:nx+1,2:ny+1) = sat(2:nx+1,2:ny+1) - fc * fact;
  end 
  % 
  % Marching forward 
  tyme = tyme + cDt;
  % 
%   % Output to data files for later plotting with GLE 
%   fname = sprintf('satforGLE_%0.6f.z', tyme);
%   fid = fopen(fname, 'wt');
%   fprintf(fid,...
%     '! nx %d ny %d xmin %0.6f xmax %0.6f ymin %0.6f ymax %0.6f \n',...
%     nx, ny, xa, xb, yc, yd);
%   for j=2:ny+1
%     for i=2:nx+1
%       fprintf(fid,'%0.6f ',sat(j,i));
%     end
%     fprintf(fid,'\n');
%   end
%   fclose(fid);
  % Saved to the 3-dim array 
  satur(itc,:,:) = sat(2:nx+1,2:ny+1);
  saturation = squeeze(satur(itc,:,:));
  saturation = saturation(:);
  % 
  % Presenting numerical and graphical results 
  show_QuadriMesh_ScaEm_ClrImg(QuadriMesh, saturation, ...
    31, 'Saturation by WG+FVM.Quadri(Q0,RT[0])');
end
save('2phase_satur_WG_FVM_QuadriQ0RT0_n100t20.mat', 'satur'); 

return;