%% Bottom Gravity Current Test Case Setup Script

% This test case is based on the paper "Dynamics of Two-Dimensional Turbulent Bottom
% Gravity Currents" by Ozgokmen and Chassignet in 2002 published in
% Journal of Physical Oceanography.

% Date: Dec. 2013
% Author: Jing Lin

%% Set up data and images saving directory

if ~isfield(app,'fullpathsgiven'), app.fullpathsgiven=0; end;

if ~app.fullpathsgiven
   d0 = fileparts([pwd, filesep]);
   d0 = d0(1 : strfind(d0,'trunk')+4);
else
   d0 = '';
end;

if ~isfield(app, 'savemat'), app.savemat = 'Save/mat/Bottom_Gravity_Current'; end;
if ~exist(fullfile(d0,app.savemat), 'dir')
    mkdir(fullfile(d0,app.savemat));
end;
if ~isfield(app, 'NoPlot'), app.NoPlot = false; end;
if ~app.NoPlot
    if ~isfield(app, 'saveimg'), app.saveimg = 'Save/img/Bottom_Gravity_Current'; end;
    if ~exist(fullfile(d0,app.saveimg), 'dir')
        mkdir(fullfile(d0,app.saveimg));
    end;
end;

%keyboard

%% Set up primitive parameters of the bottom gravity current.

if isfield (app, 'pass_params')
   A = app.pass_params.A; % Horizontal size of solution domain (unit: m)
   B = app.pass_params.B; % Vertical size of solution domain (unit: m)
   A1 = app.pass_params.A1; % Length of the initial plateau (unit: m)
   B1 = app.pass_params.B1; % Height of the initial plateau (unit: m)
   A2 = app.pass_params.A2; % Length of salinity forcing region (unit: m)

   h = app.pass_params.h; % Length scale of thickness of grivity current (unit: m)
%    (The paper does not directly provide this value.)
   rhoinitflag = app.pass_params.rhoinitflag;
else
   A = 10e3; % Horizontal size of solution domain (unit: m)
   B = 1e3; % Vertical size of solution domain (unit: m)
   A1 = 1e3; % Length of the initial plateau (unit: m)
   B1 = 0.8e3; % Height of the initial plateau (unit: m)
   A2 = 0.6e3; % Length of salinity forcing region (unit: m)

   h = 500; % Length scale of thickness of grivity current (unit: m)
%    (The paper does not directly provide this value.)
   rhoinitflag = 'cosine';
end;

% Maximum variation of salinity (unit: psu(nondimensional))
if isfield(app, 'DelS'), DelS = app.DelS; else DelS = 1.5; end;

% Slope angle of the inclined plane (unit: degree, not radian!)
if isfield(app, 'theta'), theta = app.theta; else theta = 2.5; end;

nu_x = 7; % Eddy viscosity in horizontal direction (unit: m^2/s)
r = 1e-5; % Ratio of vertical to horizontal diffusivity (nu_z / nu_x)
Pr = 7; % Turbulent Prandtl number (nu_x / K_x)
k = 1/21600; % Salinity forcing relaxation frequency (unit: Hz)

Ra = 506 * DelS; % Rayleigh number
kappa = k*h^2/nu_x; % Nondimensional relaxation frequency


%% Set up physical properties and forcing terms

% Define momentum forcing terms
% 'XY' is a 2-D array.
% XY(:, 1) contains x coordinates, XY(:, 2) contains y coordinates.
% 't' is just the current time (initial time is t = 0).
% Fu and Fv are required functions that need to be set for the solver
Fu = @(t, XY) 0;
Fv = @(t, XY) 0;

if ~isfield(app, 'nu'), app.nu = 1; end;
if ~isfield(app, 'nu2'), app.nu2 = r; end;
if ~isfield(app, 'kappa'), app.kappa = 1/Pr; end;
if ~isfield(app, 'kappa2'), app.kappa2 = r/Pr; end;
if ~isfield(app, 'g'), app.g = Ra; end;


%% Set up mesh and solution domain

% Set size of domain
a = A/h;
b = B/h;

a1 = A1/h;
b1 = B1/h;
a2 = A2/h;

tan_theta = tan(theta*pi/180);

if ~isfield(app, 'Nx'), app.Nx = 2000; end;
if ~isfield(app, 'Ny'), app.Ny = 200; end;

% Set number of timesteps.
if ~isfield(app, 'T'), app.T = 16200*nu_x/h^2; end;
if ~isfield(app, 'dt'), app.dt = 0.1*nu_x/h^2; end;
Nt = ceil(app.T/app.dt);
if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = round((app.T/app.dt)/270); end;

% Mesh the domain
% This is a standard call and should not require much modification.
% [0, a] gives the min/max x coordinates, and [0, b] gives the min/max y coordinates.
[XP, YP, XU, YU, XV, YV, dx, dy] = meshdomain(app, [0, a], [0, b]);

%% Set up the geometry mask

% This mask is user-created. Zeros are unmasked, and non-zero integers
% indicate a masked region. The domain sides should not be considered here.
% Each interior masked region could have a different number to allow
% different boundary conditions.
% Create mask of correct size. This is a standard call and should not
% require modification.
Mask = zeros(app.Ny, app.Nx);

% User-defined true-false (boolean) functions that define the mask.
% True if masked, false if not masked.
forcing = @(x, y) (x <= a2) & (y >= b1);
plateau = @(x, y) (x <= a1) & (y <= b1);
slope = @(x, y) (x > a1) & (y <= max(0, (a1-x)*tan_theta + b1));

% Assign masked region.
Mask(plateau(XP, YP)) = 1;
Mask(slope(XP, YP)) = 1;

% Logical array add_bnd signifies whether to add (add if =true) boundaries for the 
% [west, east, south, north] sides of the solution domain.
% This is almost a standard procedure except when periodic boundaries are used.
add_bnd = [true, true, true, true];

% Create preliminary node numbering matrices using function "NodePad".
% Note: This is a standard procedure and should not be modified.
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask, add_bnd);

% Plot the current boundary IDs to facilitate setting up boundary conditions.
if ~app.NoPlot
    plot_bc_ids(NodeP, Mask);
    fprintf('Paused at boundary id plot. Hit any key to continue.\n\n');
    pause;
end;


%% Set up boundary conditions (BCs)

% User assigns BC types to each boundary ID for u, v, rho in three cell arrays.

% This requires user input. Use the above plot as reference to identify BC IDs.
% e.g. Element 1 in cell array bcid2type_u represents the BC type of the boundary
% with ID 1 for velocity u.

% Note: Only the following four types of BCs can be set for velocities and density:
% Single-value Dirichlet BC   (='D')
% Multiple-value Dirichlet BC (='Dmv')
% Single-value Neumann BC     (='N')
% Multiple-value Neumann BC   (='Nmv')

% Warning: Open boundary conditions (='O') are not allowed for velocities and density!

bcid2type_u = {'D', 'N', 'N', 'D', 'N'};      
bcid2type_v = {'D', 'N', 'N', 'D', 'D'};
bcid2type_rho = {'N', 'N', 'N', 'N', 'N'};

% Automatically assign BC types to each boundary ID for q which is an intermediate variable
% for which we need to solve a Poisson equation in the projection methods
% used by the Navier-Stokes solver.
% The BC types of q depends on those of velocities and can be either 'N' or 'O'.

% Note: This is a standard procedure and should not be modified.

if ~isfield (app,'qBC')
   app.qBC='Default';
end;

Nbcs = size(bcid2type_u); Nbcs = Nbcs(2);
bcid2type_q = cell(1, Nbcs);
for i = 1 : Nbcs
    if strcmp(bcid2type_u(i), 'D') || strcmp(bcid2type_u(i), 'Dmv')...
            || strcmp(bcid2type_v(i), 'D') || strcmp(bcid2type_v(i), 'Dmv')
        bcid2type_q{i} = 'N';
    else
        bcid2type_q{i} = 'O';
    end;
end;
% Allow over-ride of automatically defined q BC
switch(app.qBC)
   case 'Dinlet'
      bcid2type_q{1} = 'D';
   case 'Dinout'
      bcid2type_q{1} = 'D';
      bcid2type_q{2} = 'D';
   case 'Dirichlet'
      for i=1:Nbcs
         bcid2type_q{i} = 'D';
      end;
end;

% Set up the final node numbering matrix using function "set_bcs".
% For each single-value BC ('D', 'N', 'O'), all its boundary cells share one ID,
% like in the preliminary node numbering matrices.
% For each multiple-value BC ('Dmv', 'Nmv'), now each boundary cell will be assigned
% a different ID.
% Besides, the IDs are arranged to satisfy: 
%               Dirichlet_BC_IDs < Open_BC_IDs < Neumann_BC_IDs

% Note: This is a standard procedure and should not be modified.

% Without density equation:
% [NodeP, nbcP, Nodeu, nbcu, Nodev, nbcv] = set_bcs(...
%     bcid2type_P, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev);

% With density equation:
[NodeP, nbcq, Nodeu, nbcu, Nodev, nbcv, NodeRho, nbcrho] = set_bcs(... 
    bcid2type_q, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev, bcid2type_rho, NodeP);
   

% Set preliminary BC values to zero since most BCs are homogeneous.
% Some of them can be modified to be non-homogeneous BCs later.

% Note: This is a standard procedure and should not be modified.

[bcsDq, bcsOq, bcsNq] = init_bcs(nbcq);
[bcsDu, bcsOu, bcsNu] = init_bcs(nbcu);
[bcsDv, bcsOv, bcsNv] = init_bcs(nbcv);
[bcsDrho, bcsOrho, bcsNrho] = init_bcs(nbcrho);

% User sets non-zero boundary conditions.

% The coordinates of the boundary points on each grid can be obtained...
% by function "get_bc_coord".
% Each of the three XY arrays contains 2 columns with x coordinates in the first and 
% y coordinates in the second stored in order of BC IDs.
% The three XY arrays correspond to boundary points of type Dirichlet ('D' and 'Dmv'),
% open ('O'), and Neumann ('N' and 'Nmv').

% [XYD, XYO, XYN] = get_bc_coord('u', XU, YU, Nodeu, bcsDu, bcsOu, bcsNu);
% ids = (abs(XYD(:, 1)) < eps);
% bcsDu(ids) = 1*h/nu_x;

% [XYD, XYO, XYN] = get_bc_coord('P', XP, YP, NodeRho, bcsDrho, bcsOrho, bcsNrho);
% ids = (abs(XYD(:, 1)) < eps);
% bcsDrho(ids) = 0.5*(1 - cos(pi*(b-XYD(ids, 2)/(b-b1))));


%% Set up initial conditions

% User specifies the functions that provide initial values at given coordinates.

[rhoinit_choice,rhoinit_parms] = strtok (rhoinitflag);

switch rhoinit_choice
   case 'cosine'
      rho_init = @(X, Y) 0.5*(1 - cos(pi*(b-Y)/(b-b1))).*forcing(X,Y);
      P_init   = @(X, Y) (app.g*0.5)*((b-Y) - ((b-b1)/pi).*sin(pi*(b-Y)/(b-b1))).*forcing(X,Y);

   case 'linear'
      rho_init = @(X, Y) (b-Y)./(b-b1).*forcing(X,Y);
      P_init   = @(X, Y) app.g*((b-Y).^2)./(2.*(b-b1)).*forcing(X,Y);

   case 'exponential'
      rhodk = str2num(rhoinit_parms);
      rho_init = @(X, Y) (1-exp(rhodk.*(Y-b)))./(1-exp(rhodk.*(b1-b))).*forcing(X,Y);
      P_init   = @(X, Y) app.g*((b-Y)-(1-exp(rhodk.*(Y-b)))./rhodk)./(1-exp(rhodk.*(b1-b))).*forcing(X,Y);

   case 'cosine+exponential'
      rhodk = str2num(rhoinit_parms);
      rho_init = @(X, Y) 0.5.*( 0.5*(1 - cos(pi*(b-Y)/(b-b1))) + (1-exp(rhodk.*(Y-b)))./(1-exp(rhodk.*(b1-b))) ).*forcing(X,Y);
      P_init   = @(X, Y) (app.g*0.5)*( 0.5*((b-Y) - ((b-b1)/pi).*sin(pi*(b-Y)/(b-b1))) + ((b-Y)-(1-exp(rhodk.*(Y-b)))./rhodk)./(1-exp(rhodk.*(b1-b))) ).*forcing(X,Y);

   case 'cosine+linear'
      rho_init = @(X, Y) 0.5.*( 0.5*(1 - cos(pi*(b-Y)/(b-b1))) + (b-Y)./(b-b1) ).*forcing(X,Y);
      P_init   = @(X, Y) (app.g*0.5)*( 0.5*((b-Y) - ((b-b1)/pi).*sin(pi*(b-Y)/(b-b1))) + ((b-Y).^2)./(2.*(b-b1)) ).*forcing(X,Y);

   case 'linear+exponential'
      rhodk = str2num(rhoinit_parms);
      rho_init = @(X, Y) 0.5.*( (b-Y)./(b-b1) + (1-exp(rhodk.*(Y-b)))./(1-exp(rhodk.*(b1-b))) ).*forcing(X,Y);
      P_init   = @(X, Y) (app.g*0.5)*( ((b-Y).^2)./(2.*(b-b1)) + ((b-Y)-(1-exp(rhodk.*(Y-b)))./rhodk)./(1-exp(rhodk.*(b1-b))) ).*forcing(X,Y);
end;
u_init   = @(X, Y) zeros(size(X));
v_init   = @(X, Y) zeros(size(X));

% Initialize the state variables with values from both boundary and internal cells.

% Note: This is a standard procedure and should not be modified.

rho = [bcsDrho; bcsNrho; rho_init(XP(idsP), YP(idsP))];
P   = [bcsDq; bcsNq; P_init(XP(idsP), YP(idsP))];
u   = [bcsDu; bcsNu; u_init(XU(idsu), YU(idsu))];
v   = [bcsDv; bcsNv; v_init(XV(idsv), YV(idsv))];

%keyboard


%% Set up density forcing term

% Define density forcing term
if ~isfield (app,'rhofrc_type')
   app.rhofrc_type = 'fullvert';
end;
switch app.rhofrc_type
   case 'fullvert'
      app.Frho = @(t, XY, rho) -kappa*(rho - 1).*((XY(:,1) <= a1)&(XY(:,2) >= b1));
   case 'halfvert'
      app.Frho = @(t, XY, rho) -kappa*(rho - (XY(:,2)<=(b+b1)*0.5)).*((XY(:,1) <= a1)&(XY(:,2) >= b1));
end;


%% Save all the setup data in a mat-file.

save(fullfile(d0,app.savemat,'setup'));

