%% Buoyant "Bubble" test case (Navier-Stokes with Boussinesq Approximation)

% Last modified by Jing Lin in Apr. 2014


%% 01 Specify problem-specific settings


%% 02 Set up general physical properties

if ~isfield(app, 'nu'), error('app.nu not specified'); end;
if ~isfield(app, 'kappa'), error('app.kappa not specified'); end;
if ~isfield(app, 'g'), error('app.g not specified'); end;


%% 03 Generate mesh and geometry mask

% Define size of solution domain
a = 1;
b = 5;

if ~isfield(app, 'Nx'), error('app.Nx not specified'); end;
if ~isfield(app, 'Ny'), error('app.Ny not specified'); 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);


%% 04 Set up boundary conditions (BCs)

% 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 (~isfield(app, 'NoPlot'))||(~app.NoPlot)
    plot_bc_ids(NodeP, idsP);
    fprintf('Paused at boundary id plot. Hit any key to continue.\n\n');
    pause;
end;

% 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', 'D', 'D', 'N'};
bcid2type_v = {'D', 'D', 'D', 'D'};
bcid2type_rho = {'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.

Nbcs = size(bcid2type_u, 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;

% 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, nbcq, Nodeu, nbcu, Nodev, nbcv] = set_bcs(...
%     bcid2type_q, 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.


%% 05 Specify time marching settings

% Set number of timesteps.
if ~isfield(app, 'T'), error('app.T not specified'); end;
if ~isfield(app, 'dt'), error('app.dt not specified'); end;
Nt = ceil(app.T/app.dt);


%% 06 Set up initial conditions

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

P_init   = @(X, Y) zeros(size(X));
u_init   = @(X, Y) zeros(size(X));
v_init   = @(X, Y) zeros(size(X));
rho_init = @(X, Y) -0.001*(sqrt((X-1/2).^2+(Y-1).^2) < 0.25) ;

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

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

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))];
rho = [bcsDrho; bcsNrho; rho_init(XP(idsP), YP(idsP))];


%% 07 Specify general solver options


%% 08 Set up user-defined 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 = @(time,XY) 0;
Fv = @(time,XY) 0;


%% 09 Specify monitoring, plotting and saving settings

if ~isfield(app, 'PlotIntrvl'), error('app.PlotIntrvl not specified'); end;
if ~isfield(app, 'NoPlot'), error('app.NoPlot not specified'); end;


%% 10 Set up other problem-specific setting     

