function [P u v rho NodeP Nodeu Nodev Tsw Ssw] = TestCase(problem, app, solver)
% function [P u v rho NodeP Nodeu Nodev] = TestCase(problem, app, solver))
%
% This is a control function to setup and run various test cases that solve
% Stokes and Navier stokes problems.
%
% INPUTS:
%       Problem:    Test case number (see below)
%       Solver:     See below for Solver choices
%       app:	    Application Data structure. TestCase automatically assigns 
%                   default values. Can contain entries such as:
%          app.Nx:           Number of control volumes in x direction (Optional)
%          app.Ny:           Number of control volumes in y direction (Optional)
%          app.T:            Final time (Optional)
%          app.dt:           Timestep size (Optional)
%          app.PlotIntrvl:   Number of timesteps till solution is plotted
%                             (Optional)
%          app.nu:           Viscosity used in momentum equations (Optional)
%                             app.nu = 1/Re
%          app.kappa:        Diffusivity used in Density equation (Optional)
%          app.savemat:      String indicating where to save mat files for 
%                             every PlotIntrvl timesteps. Files are saved in
%                             <rootdir>/Save/app.savemat, or in
%                             <rootdir>/Save/mat/ when app.savemat is empty
%                             (Optional)
%          app.saveimg:      Same as app.savemat except images are saved
%                             instead.
%          app.Stokes        If true, uses a stokes solver instead of 
%                             Navier-Stokes.
%          app.S:            Number of DO modes (Optional). 
%                              If length(S)==1 - non-adaptive.
%                              If length(S)==2 - adaptive S = [Smin, Smax]
%          app.DOAdapt:      Number of timesteps between adaptivly changing the
%                            number of DO mods (Optional).
%          app.var           Variance of largest mode at initialization
%                             (Optional) for DO solvers
%          app.addvar        If smallest variance mode has this variance,
%                              a new mode is added for adaptive DO solvers    
%                             (Optional) 
%          app.subbvar       If smallest variance mode has this variance, and
%                              second smallest variance is less than twice
%                              this variance, then the smallest variance mode
%                              is subtraced for adaptive DO solvers (Optional) 
%          app.pathin      Structure containing information about the
%                           path-planning section (Optional)
%          app.pathin.p:     Starting positions of assets p = [x;y] (Optional)
%          app.pathin.mag:   Magnitude of velocity that vehicle has (Optional)
%          app.pathin.v:     Prescribed vector velocities of assets relative to the
%                             asset coordinate system. That is, the coordinate
%                             system is fixed in space, not taking into account the
%                             flowfield velocity. The prescribed vector velocity of 
%                             asset i is: v(:,1:3,i) = [vx;vy;t]. (Optional)
%
%   If optional inputs are not specified or left empty [], default values
%   are used.
%
% OUTPUTS:
%       P:          Final Pressure
%       u:          Final U-component of velocity
%       v:          Final V-component of velocity
%       rho:        Final density
%       NodeP:  Pressure node-numbering matrix (useful for plotting)
%       Nodeu:  U-velocity node-numbering matrix (useful for plotting)
%       Nodev:  V-velocity node-numbering matrix (useful for plotting)
%       Tsw:    Temperature of sea-water (For solver 5)
%       Ssw:    Salinity of sea-water (For solver 5)
%
% PROBLEM CHOICES:
%   1: Heated plate demo: Demonstrates solution using Diffusion operator on
%       non-trivial geometry with Dirichlet and Neumann boundary conditions
%   2: Scalar advection: Demonstrates advection of a scalar quantity over a
%       periodic domain. This simulation uses the QUICK scheme
%   3: Scalar advection: Same as (2) but uses the UPWIND scheme
%   4: Analytical test case for Stokes Problem. Solution is:
%         u = pi*sin(time).*(sin(2*pi*y)).*(sin(pi*x)).^2);
%         v =-pi*sin(time).*(sin(2*pi*x)).*(sin(pi*y)).^2);
%         p = sin(time).* cos(  pi*x ).* sin(pi*y) ;
%   5: Poiseuille flow test case
%   6: Lock exchange experiment (buoyancy driven flow)
%   6.1: Lock exchange experiment with stochastically perturbed interface
%   6.2: Lock exchange experiment with completely stochastic initial
%        conditions
%   7: Sudden expansion test case (u velocity = 1, nu_default = 0.01,
%       Re=100)
%   8: Sudden expansion test case with warm inlet (buoyancy effects)
%   9:Warm bubble surrounded by cool air
%   10:Lid Driven Cavity. Lid is driven at velocity u = 1 in 1x1 domain
%   11:Stochastic Lid Driven Cavity. Lid is driven at velocity u = 1 in 1x1
%       domain
%   12: Adaptive Stochastic Lid Driven Cavity. Lid is driven at velocity 
%       u = 1 in 1x1 domain
%   13: Double Gyre Test Case Stochastic
%   14: Double Gyre Test Case Deterministic
%   15: Double Gyre Test Case Stochastic with PathPlanning
%   16: PathPlanning with analytically-provided flow-field (1)
%   17: PathPlanning with analytically-provided "Highway" flow-field (2)
%   18: Thermohaline circulation example
%   19: Flow past a cylinder
%   20: Deterministic Rayleigh-Bernard Convection
%   21: Stochastic Rayleigh-Bernard Convection
%
% SOLVER CHOICES:
%   1: Fluids Solver
%   2: Stochastic Fluids Solver
%   3: Stochastic Fluids Solver with Pathplanning (non-boussinesq)
%   4: Pathplanning Levelset Solver With analytical Flowfields
%   5: Boussinesq equations, Temperature, Salinity and a non-linear
%      equation of state
%
% Authors: Matt Ueckermann, Pierre Lermusiaux, Pat Haley, Themis Sapsis,
%          Konur Yigit, Theo Touvet. Written for MIT course 2.29.

if nargin<2
    app=[];
end

%Initialize ouputs
P=[]; u=[]; v=[]; rho=[]; Nodeu=[]; Nodev=[]; NodeP=[]; Tsw=[]; Ssw=[];

switch problem
    case 1 % Heated Plate: Test case for Laplace sovler
        display('Heated Plate demo')
        display('Note, min(Nx) = min(Ny) = 14')
        if exist('solver','var')
            if (nargin>2 || ~isempty(solver))
                display('Note, no solver can be selected for this testcase')
            end
        end
        %Set default values
	    Nx = 20; Ny = 30;
        if isfield(app,'Nx')
            Nx = max(app.Nx,14);
        end
        if isfield(app,'Ny')
            Ny = max(app.Ny,14);
        end        
        Heated_Plate_demo;

    case 2 % Tracer_Advection: Test case for advection QUICK scheme
        display('Tracer advection demo using QUICK scheme')
        if exist('solver','var')
            if (nargin>2 || ~isempty(solver))
                display('Note, no solver can be selected for this testcase')
            end
        end
        %Set default values
        Nx = 30; Ny = 30;
        if isfield(app,'Nx')
            Nx = app.Nx;
        end
        if isfield(app,'Ny')
            Ny = app.Ny;
        end        
        Tracer_Advection_demo;

    case 3 % Tracer_Advection: Test case for advection UPWIND scheme 
        display('tracer advection demo using UPWIND')
        if exist('solver','var')
            if (nargin>2 || ~isempty(solver))
                display('Note, no solver can be selected for this testcase')
            end
        end
        Nx = 30; Ny = 30;
        if isfield(app,'Nx')
            Nx = app.Nx;
        end
        if isfield(app,'Ny')
            Ny = app.Ny;
        end   
        Tracer_Advection_UW_demo;

    case 4 % Stokes equation: Analytical test case
        display('Analytical Stokes test case')
        %Set default values
        Nx = 30; Ny =30; T = 2*pi; dt = T/100; PlotIntrvl = 10; nu=1;
        Stokes = true;
        SetupScript = 'StokesForcingSetup';
        PlotScript = 'StokesForcingPlots';
        if nargin<3 || isempty(solver)
            solver = 1;
        end

    case 5 %Poiseuille Flow test case
        display('Poiseuille Flow')
        %Set default values
        Nx = 100; Ny = 80; T = 2.5; dt=2.5/100; PlotIntrvl = 10; nu = 0.24;
        SetupScript = 'StokesPoiseuilleSetup';
        PlotScript = 'StokesPoiseuillePlots';
        Stokes = true;
        if nargin<3 || isempty(solver)
            solver = 1;
        end
        if solver==10
            dt = 1/10000;
            PlotIntrvl = 100;
        end

    case 6 %Lock exchange experiment
        display('Lock Exchange')
        Nx = 30; Ny = 30; T = 2; dt = 1/1000; nu = 0.0001; kappa=0.0001;
        PlotIntrvl = 100;
        SetupScript = 'LockExchangeSetup';
        PlotScript = 'LockExchangePlots';
        if nargin<3 || isempty(solver)
            solver = 1;
        end
        
    case 6.1 %Lock exchange experiment Stochastic 1
        display('Lock Exchange')
        Nx = 30; Ny = 30; T = 2; dt = 1/1000; nu = 0.0001; kappa=0.0001;        
        S=2;
        PlotIntrvl = 100;
        SetupScript = 'LockExchangeSetup_DO';
        PlotScript = 'DORhoAdaptGenericPlots';
        app.PlotScript = PlotScript; 
        if nargin<3 || isempty(solver)
            solver=2;
        end
         
    case 6.2 %Lock exchange experiment Stochastic2
        display('Lock Exchange')
        Nx = 30; Ny = 30; T = 2; dt = 1/1000; nu = 0.0001; kappa=0.0001;
        S=1;
        PlotIntrvl = 100;
        SetupScript = 'LockExchangeSetup_DO2';
        PlotScript = 'LockExchangePlots_DO2';
        app.PlotScript = PlotScript;
        if nargin<3 || isempty(solver)
            solver = 2;
        end

    case 7 %Sudden Expansion
        display('Sudden Expansion (Look familiar?)')
        Nx = 80; Ny = 60; T = 50; dt   = 5/200; nu=0.001;  kappa=0.01;
        PlotIntrvl = 1000/15;
        SetupScript = 'SuddenExpansionSetup';
        PlotScript = 'SuddenExpansionPlots';
        if nargin<3 || isempty(solver)
            solver = 1;
        end

    case 8 %Buoyant sudden expansion
        display('Buoyant Sudden Expansion')
        Nx = 80; Ny = 60; T = 50; dt   = 5/200; nu=0.0015;  kappa=0.01;
        PlotIntrvl = 1000/15;
        SetupScript = 'SuddenExpansionBuoyantSetup';
        PlotScript = 'SuddenExpansionPlots';
        if nargin<3 || isempty(solver)
            solver = 1;
        end

    case 9 %Buoyant bubble
        display('Buoyant Bubble')
        Nx = 40; Ny = 60; T = 1000; dt=1;  nu=0.0001;  kappa=0.0001;
        PlotIntrvl = 50;
        SetupScript = 'BuoyantBubbleSetup';
        PlotScript = 'BuoyantBubblePlots';
        if nargin<3 || isempty(solver)
            solver = 1;
        end

    case 10 %Lid Driven cavity flow
        display('Lid Driven Cavity Flow')
        display('For convergence test, run: LidDrivenCavityConvergenceTest.m')
        Nx = 21; Ny = 21; T = 100; dt = 1/60; PlotIntrvl = 300; nu=0.001;
        SetupScript = 'LidDrivenCavitySetup';
        PlotScript = 'LidDrivenCavityPlots';
        if nargin<3 || isempty(solver)
            solver = 1;
        end

    case 11 %Stochastic Lid Driven cavity flow
        display('Stochastic Lid Driven Cavity Flow')
        Nx = 30; Ny = 30; T = 5; dt = 1/200; PlotIntrvl = 40; nu=0.001;
        app.var = 1e-1;
        S=5;
        SetupScript = 'LidDrivenCavitySetup_DO';
        PlotScript = 'LidDrivenCavity_DO_Plots';
        if nargin<3 || isempty(solver) 
            solver = 2;
        end
        if (solver ~=2)
            error('This case does not work for non-DO solvers. Try Case 10.')
        end

     case 12 %Adaptive Stochastic Lid Driven cavity flow
        display('Adaptive Stochastic Lid Driven Cavity Flow')
        display('Note, initialization is different than case 11')
        Nx = 40; Ny = 30; T = 5; dt = 1/200; PlotIntrvl = 40; nu=0.001;
        app.var = 1e-1;
        DOAdapt = PlotIntrvl/2;
        S=[1 4];
        SetupScript = 'LidDrivenCavitySetup';
        PlotScript = 'DOAdaptGenericPlots';
        if nargin<3 || isempty(solver)
            solver = 2;
        end
        if (solver ~= 2)
            error('This case does not work for non-DO solvers. Try Case 10.')
        end
    
    case 13 %Double Gyre test case
        display('Stocastic Double Gyre Flow')
        Nx = 30; Ny = 30;T = 1; dt = 1e-4; PlotIntrvl = T/dt/20; nu=0.01;
        S=[1 5];
        DOAdapt = 20;
        SetupScript = 'DoubleGyreSetup';
        PlotScript = 'DOAdaptGenericPlots';
        if nargin<3 || isempty(solver)
            solver = 2;
        end
        if (solver ~= 2)
            error('This case does not work for non-DO solvers. Try Case 14.')
        end

    case 14 %Double Gyre test case
        display('Deterministic Double Gyre Flow')
        Nx = 30; Ny = 30;T = 1; dt = 1e-4; PlotIntrvl = T/dt/20; nu=0.01;
        SetupScript = 'DoubleGyreSetup';
        PlotScript = 'DoubleGyrePlots'; 
        if nargin<3 || isempty(solver)
            solver = 1;
        end

     case 15 %Double Gyre test case with path-planning
        display('Stochastic Double Gyre Flow with path-planning')
        Nx = 40; Ny = 40; T = 0.15; dt = 1e-3/6; PlotIntrvl = 20; nu=0.0005;
        app.var = 10;
        S=[2];
        pathin.p = [0.25 0.45];
        pathin.v = [1 0.5]*2;
        pathin.mag=sqrt(sum(pathin.v.^2));
        pathin.v = [0 pathin.mag];

        %This is a new variable we introduce into app that is specific to our
        %plotting script. The Level Set require re-initialization at times, and
        %the following parameter determines how many plot cycles between 
        %re-initialization.
        app.reinitLS = 200;
        %Setupscript
        SetupScript = 'DoubleGyrePathPlanSetup';
        %This is so that the Plotscript isn't overwritten when the DOmodes are 
        %initialized in DOInit.m
        app.PlotScript = 'DoubleGyrePathPlanPlots'; 
        PlotScript = 'DoubleGyrePathPlanPlots'; 
        if nargin>3 && ~isempty(solver)
            if solver ~= 3
                error('Please choose solver 3.')
            end
        end
        solver = 3;

     case 16 %PathPlanning with Analytical Flowfield
        display('PathPlanning with Analytical Flowfield')
        Nx = 40; Ny = 40;T = 2; dt = 1e-3; PlotIntrvl = 100;
        app.reinitLS = floor(PlotIntrvl/30*100);
        U = @(x,y,t) cos(pi*t).*(y-0.5).*sin(pi*x).*sin(pi*y);
        V = @(x,y,t) -cos(pi*t).*(x-0.5).*sin(pi*x).*sin(pi*y);
        pathin.p = [0.55 0.38];
        pathin.v = [1 1]*0.1;
        pathin.mag=sqrt(sum(pathin.v.^2));
        pathin.v = [0 pathin.mag];
        SetupScript = 'PathPlanSetup';
        PlotScript = 'PathPlanPlot';  
         if nargin<3 || ~isempty(solver)
            display('Choosing solver 4.')
        end
        solver = 4;

      case 17 %PathPlanning with Highway Analytical Flowfield
        display('PathPlanning with Highway Analytical Flowfield')
        Highway;
        solver= 0;

      case 18 %Thermohaline circulation
        display('Temperature driven deterministic Thermohaline circulation')
        Nx = 30; Ny=30; T = 20; dt = 1/100; nu = 0.01;
        kappaT = 0.01; kappaS = 0.01; PlotIntrvl = 100;
        SetupScript = 'ThermohalineSetup';
        PlotScript = 'ThermohalinePlots'; 
        if nargin>2
            if solver~=5
                display('Using solver 5 instead.')
            end
        end 
        solver=5;
        
    case 19
        display('Flow past cylindrical column')
        Nx = ceil(171); Ny = ceil(21); T = 30; dt   = 1/75; nu=0.01;  kappa=0.00;
        PlotIntrvl = 1/dt/3;
        SetupScript = 'CylinderCrossFlowSetup';
        PlotScript = 'SuddenExpansionPlots';
        if nargin<3 || isempty(solver)
            solver = 1;
        end
        
    case 20 %Rayleigh-Bernard Convection
        display('Rayleigh-Bernard Convection')
        Nx = 64; Ny = 64; T = 1; dt = 10^-3; nu = 0.7; kappa = 1;
        app.PrRa= -nu*15000; 
        PlotIntrvl = 50;
        app.SetupScript = 'RayleighBenardSetup';
        app.PlotScript = 'RayleighBenardPlots';
        
        PlotScript =app.PlotScript; 
        SetupScript=app.SetupScript;

        if nargin<3 || isempty(solver)
            solver = 1;
        end    

    case 21 %Stochastic Rayleigh-Bernard Convection
        display('Stochastic Rayleigh-Bernard Convection')
        Nx = 64; Ny = 64; T = 1; dt = 10^-3; nu = 0.7; kappa = 1;
        app.PrRa=-nu*15000; 
        PlotIntrvl = 50;
        app.SetupScript = 'RayleighBenardSetup_venturi_karniadakis_test';
        app.PlotScript = 'DORhoAdaptGenericPlots';%'LockExchangePlots';
        app.IP=[1,1,100];
        app.S=3;
        
        PlotScript =app.PlotScript; 
        SetupScript=app.SetupScript;
         %         PlotScript = 'RayleighBenardPlots';
        if nargin<3 || isempty(solver)
            solver = 2;
        end    

      otherwise
            display('No known Problem selected')
    
end

%Over-writes
if ~isfield(app,'Stokes')
    if exist('Stokes','var')
        app.Stokes = Stokes;
    end
end
if ~isfield(app,'Nx')
    if exist('Nx','var')
        app.Nx = Nx;
    end
end
if ~isfield(app,'Ny')
    if exist('Ny','var')
        app.Ny = Ny;
    end
end
if ~isfield(app,'T')
    if exist('T','var')
        app.T = T;
    end
end
if ~isfield(app,'dt')
    if exist('dt','var')
        app.dt = dt;
    end
end
if ~isfield(app,'PlotIntrvl')
    if exist('PlotIntrvl','var')
        app.PlotIntrvl = PlotIntrvl;
    end
end
if ~isfield(app,'nu')
    if exist('nu','var')
        app.nu = nu;
    end
end
if ~isfield(app,'kappa')
    if exist('kappa','var')
        app.kappa = kappa;
    end
end
if ~isfield(app,'g')
    if exist('g','var')
        app.g = g;
    end
end
if ~isfield(app,'kappaT')
    if exist('kappaT','var')
        app.kappaT = kappaT;
    end
end
if ~isfield(app,'kappaS')
    if exist('kappaS','var')
        app.kappaS = kappaS;
    end
end
if ~isfield(app,'S')
    if exist('S','var')
        app.S = S;
    end
end
if ~isfield(app,'DOAdapt')
    if exist('DOAdapt','var')
        app.DOAdapt = DOAdapt;
    end
end
if ~isfield(app,'pathin')
    if exist('pathin','var')
        app.pathin = pathin;
    end
end
if ~isfield(app,'U')
    if exist('U')
        app.U = U;
    end
end
if ~isfield(app,'V')
    if exist('V')
        app.V = V;
    end
end
app
%Call solvers
if problem >3
    switch solver
        case 1 %Standard Fluids Solver
            display('Fluids solver')
            [P u v rho NodeP Nodeu Nodev] = ...
                FluidsSolver(app, SetupScript, PlotScript);

        case 2 %Standard Stochastic Fluids Solver
            display('Stochastic Fluids solver')            
            [P u v rho NodeP Nodeu Nodev] = NavierStokesInc_DO_rho(app,...
                SetupScript, PlotScript);     
            
        case 3 %Stochastic Fluids Solver with Pathplanning
            display('Stochastic Fluids solver with pathplanning')            
            [P u v rho NodeP Nodeu Nodev] = NavierStokesInc_DO_coriolis_pathplan(app,...
                SetupScript, PlotScript);              

        case 4 %Pathplanning with analytical flow-field
            display('PathPlanning using Level sets and Analytical flowfield Solver')
            [P u v rho] = Solver_pathplan_NB(app, SetupScript, PlotScript);
            
        case 5 %Boussinesq with seawater Temperature and salinity instead of density
            display('Boussinesq with seawater Temperature and salinity instead of density')
            [P u v rho NodeP Nodeu Nodev Tsw Ssw] = BoussinesqTS(app,...
                SetupScript, PlotScript);                                                    
        otherwise
            display('No known solver selected')
    end
end
