1 clear all;close all;clc;
   2 
   3 %Grid
   4 Nx = 1024; Nz = 256;
   5 Lx = 4.0; Lz = 1.0;
   6 hx = Lx/Nx; hz = Lz/Nz;
   7 x = -Lx/2:hx:(Lx/2)-hx;   z = -Lz/2:hz:(Lz/2)-hz;
   8 [xx,zz]= meshgrid(x,z);
   9 
  10 % Model parameters
  11 Ra = 1000;
  12 var_lnk=1; %variance of log(natural) permeability
  13 corr_lenx=10*hx; %correlation length in x
  14 corr_lenz=10*hz; %correlation length in z
  15 
  16 % Initial condition
  17 c = zeros(Nz,Nx);
  18 c0 = 1; % Prescribed concentration at the top
  19 c(1,:) = c0; % Prescribed concentration at the top
  20 % c(1,:) = 1e-03*rand(1,Nx); % Small perturbation if needed
  21 % Visualize the initial concentration field
  22 % surf(xx,zz,c,'facecolor','interp','edgecolor','none','facelighting','phong'); view([0,0,1]); 
  23 % axis equal tight; colorbar
  24 
  25 % Generate multigaussian random permeability field
  26 [perm,var_lnk_actual]=random_perm(var_lnk,corr_lenx,corr_lenz,Nx,Nz,Lx,Lz);
  27 
  28 % Visualize log permeability field
  29 surf(xx,zz,log(perm),'facecolor','interp','edgecolor','none','facelighting','phong');view([0,0,1]); 
  30 axis equal tight off; colorbar; %caxis([2 7]);
  31 
  32 % Save the permeability field to use in future runs. 
  33 % Since random number generator is used in the function below, every call generates a new field with presricbed variance and correlation lengths.
  34 % dlmwrite(['perm_'  num2str(var_lnk) '_' num2str(corr_lenx/hx) '_' num2str(corr_lenz/hz) '.txt'],perm);
  35 save perm.mat perm
  36 
  37 t = 0; istep= 0; 
  38 dt = 0.0001;
  39 % dtadv = min(hz,hx);
  40 % dtdiff = Ra/2*dtadv.^2;
  41 % dt = 1/2*min(dtdiff,dtadv);
  42 
  43 while (t<1)
  44     istep= istep+1;
  45 	
  46     % Solve for pressure and velocities implicitly using TPFA/MPFA 
  47     [ux,uz]= solve_vel(c,perm,hx,hz);
  48     
  49 	% Solve for concentration using RK explicit time integration 
  50     for istage= 1:3
  51         Res= residual_conc(c,ux,uz,Ra,hx,hz,c0); % Residual of the concentration eqn.
  52         if istage==1
  53             c= c + dt*( (8/15)*Res );
  54             Res1= Res;
  55         elseif istage==2
  56             c= c + dt*( (5/12)*Res + (-17/60)*Res1 );
  57             Res1= Res;
  58         else
  59             c= c + dt*( (3/4)*Res + (-5/12)*Res1 );
  60         end
  61     end
  62     t = t+dt;
  63     dt = min(hx,hz)/(4*max(max(ux)));
  64 	
  65     %%% save output %%%%%%%%%%%%%%%%%%        
  66 	surf(xx,zz,c,'facecolor','interp','edgecolor','none','facelighting','phong'); view([0,0,1]); 
  67 	axis equal tight; colorbar;
  68 	drawnow    
  69 end