% Vertical mixing from a steady state transverse line source
% solve 2D ADE equation: C(x,y) in a rectangular box for
% uniform u_x, e_y
% u_x*dC/dx=+e_y*d^2C/dy^2
% continuous source at (x0,y0)=(0,ymin+frac*(ymax-ymin))
%-------------------------------------------------
clear all
close all
global M w h;
global u_x e_y;
% problem data
M=10.; % [kg/s] mass discharged from point load
Q=0.5; % [m^3/s] flow rate
w=2.; % [m] channel width
h=1.; % [m] channel height
A=w*h; % [m^2] channel cross section
u_x=Q/A; % [m/s] flow velocity
e_y=0.01; % [m^2/s] turbulent diffusion in the vertical direction
% space discretizazion: x=[xmin:xmax]
% xmin > 0 to avoid divergence of analytical solution 
xmax=25.;
xmin=0.01*xmax;
% discrete points to be used for calculation
Nx=100;
% domain extension: y=[0:h] 
ymin=1/100*h; % free surface
ymax=99/100*h; % bottom height
% discrete points to be used for calculation
Ny=50;
% 2D grid dimensional
% linspace is used to define a row vector with Nx points between xmin e xmax
% the vector can be also defined as
% x=[xmin:dx_increment:xmax]
% coordinate system is centered at bottom left of box domain
x=linspace(xmin,xmax,Nx); % [m]
y=linspace(ymin,ymax,Ny); % [m]
% point load data
frac=0.65; % position of point load (from bottom/surface)
h0=frac*h; % [m] height above the bottom
% analytical solution for concentration field downstream the real source
yS=y-h0; % yS: y-coordinate with origin at source position 
CR=conc_unbounded(Nx,Ny,x,yS);
% figure 1: concentration field due to real source in unbounded domain
figure(1)
xlabel("downstream distance, [m]");
ylabel("vertical coordinate, [m]");
% isocontour plot:filled contours
contourf(x,y,CR);
% isocontour plot: line contours at specified contour levels
Clev=[1:2:20] % contour range [1:20], step 2
contour(x,y,CR,Clev);
colorbar;
% CR does not satisfy the boundary conditions at top/bottom of channel
% dC/dy <>0 at y=h and y=0 (flux across the top and bottom surface!)
% add virtual sources to account for the effect of boundaries
% virtual source 1: below the bottom
yV=y+h0; 
CV=conc_unbounded(Nx,Ny,x,yV);
% figure 2: concentration field due to real source and first virtual source
figure(2)
xlabel("downstream distance, [m]");
ylabel("vertical coordinate, [m]");
% effect of real and first virtual source
CRV=CR+CV;
contourf(x,y,CRV);
% isocontour plot: line contours at specified contour levels
Clev=[1:2:20]; % contour range [1:20], step 2
contour(x,y,CRV,Clev);
colorbar;
% CRV satisfies the boundary conditions at bottom of channel but not at top
% dC/dy=0 at y=0 and dC/dY<>0 at y=h (flux across the top!)
% add virtual source at 2h-h0 to account for the effect of top boundary
% virtual source 2: above the top
yV2=y-(2*h-h0);
CV2=conc_unbounded(Nx,Ny,x,yV2);
% figure 3: concentration field due to real source, first and second virtual source
figure(3)
xlabel("downstream distance, [m]");
ylabel("vertical coordinate, [m]");
% effect of real and first virtual source
CRV12=CR+CV+CV2;
contourf(x,y,CRV12);
% isocontour plot: line contours at specified contour levels
Clev=[1:2:20]; % contour range [1:20], step 2
contour(x,y,CRV12,Clev);
colorbar;
% CRV12 satisfies the boundary conditions at bottom/bottom of channel but not at large distances
% add additional virtual sources! Position of virtual sources can be defined as 
% V2=(0,2h-h0)
% V3=(0,-(2h-h0))
% V4=(0,4h-h0)
% V5=(0,4h+h0)
% V6=(0,-4h-h0)
% V7=(0,-4h+h0)
% or, in general (0,+-2nh+-h0) for n=1,... (n=number of reflections at boundaries)
%---------------------------------------------
% cicle to superpose the effect of virtual sources
Csum=0;  % concentration field to build up overall effect 
Csum=Csum+CR+CV; % real + first virtual source
for i=1:3;
yV1=y-(2*i*h-h0);
yV2=y+(2*i*h-h0);
yV3=y-(2*i*h+h0);
yV4=y+(2*i*h+h0);
CV1=conc_unbounded(Nx,Ny,x,yV1);
CV2=conc_unbounded(Nx,Ny,x,yV2);
CV3=conc_unbounded(Nx,Ny,x,yV3);
CV4=conc_unbounded(Nx,Ny,x,yV4);
Csum=Csum+CV1+CV2+CV3+CV4;
% isocontour plot: line contours at specified contour levels
figure(2)
Clev=[1:1:40]; % increase the number/range of contours [1:40], step 1
contour(x,y,Csum,Clev);
colorbar;
endfor
%--------------------------------------------------
% compact coding for the same implementation of superposition
%--------------------------------------------------------
%Csum=0;  % concentration field to build up overall effect 
%Csum=Csum+CR+CV; % real + first virtual source
%for i=1:3;
%for is=1:2;
%for it=1:2;
%yVit=y+(-1)^is*(2*i*h-(-1)^it*h0);
%CVit=conc_unbounded(Nx,Ny,x,yVit);
%Csum=Csum+CVit;
%% isocontour plot: line contours at specified contour levels
%figure(2)
%Clev=[1:1:40]; % increase the number/range of contours [1:40], step 1
%contour(x,y,Csum,Clev);
%colorbar;
%endfor
%---------------------------------------
% homogeneization index: Pm=[0:1]; 
% if Pm=1, concentration field perfectly mixed along the vertical
% if Pm=0, there are unmixed regions along the vertical
% the distance downstream the point load to achieve complete mixing is calculated 
% from the plot of Pm(x)=min(C(x,y))/max(C(x,y)) vs x
% Alternatively: concentration dis-homogeneity along the vertical can be given by
% Pm(x)=variance or std deviation of C(y) 
for i=1:Nx;
Cmin=min(Csum(:,i));
Cmax=max(Csum(:,i));
Pm(i)=Cmin/Cmax;
endfor
figure(5)
plot(x,Pm);


%----------------------------------------------
% saving data: save into file "data.out" the given list of variables (x,y,Csum,Pm)
%----------------------------------------------
save data065.out x y Csum Pm
% clean all the working space (delete all variables)
clear all
% close all previous figures
close all 
% load data from files saved before to compare solutions
% load without list of variables loads ALL data with their original name 
% (see variable name, type and values in the workspace)
load data065.out
% assign a different name to variables if necessary
Pm025=Pm; % define a new vector Pm1 copying Pm values
Csum025=Csum;
% second case study
load data05.out
% assign a different name to variables if necessary
Pm050=Pm; % define a new vector Pm1 copying Pm values
Csum050=Csum;
% comparison between concentration fields
figure(1)
contour(x,y,Csum025,1:1:40,"linewidth",2); % increased line thickness
title("Point load at 0.25 h") % titolo figura
xlabel("downstream distance, [m]");
ylabel("vertical coordinate, [m]");
colorbar;
figure(2)
contour(x,y,Csum050,1:1:40,"linewidth",2); % increased line thickness
title("Point load at 0.50 h") % titolo figura
xlabel("downstream distance, [m]");
ylabel("vertical coordinate, [m]");
colorbar;
% comparison between downstream distance to complete vertical mixing
figure(3)
plot(x,Pm025,'b',"linewidth",2,x,Pm050,'r',"linewidth",2);
legend('h0=0.25 h','h0=0.50 h'); % labels for line plot (same order as in plot)
ylabel("Degree of mixing, [-]");
xlabel("downstream distance, [m]");



