function [P, u, v, rho, NodeP, Nodeu, Nodev, Tsw, Ssw] = TestCase(problem, app)
% 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)
%       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 takinto 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 Deterministic
%   14: Double Gyre Test Case Stochastic
%   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;
app.NoPlot = false;

%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');
        
        if (~isfield(app, 'Nx') || (app.Nx < 14)), app.Nx = 20; end;
        if (~isfield(app, 'Ny') || (app.Nx < 14)), app.Ny = 30; end;
        if ~isfield(app, 'nu'), app.nu = 1; end;
        
        Heated_Plate_demo;
        
    case 2 % Tracer_Advection: Test case for advection QUICK scheme
        display('Tracer advection demo using QUICK scheme');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 1; end;
        if ~isfield(app, 'dt'), app.dt = app.T/1000; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 50; end;
        
        Tracer_Advection_demo;
        
    case 3 % Tracer_Advection: Test case for advection UPWIND scheme
        display('tracer advection demo using UPWIND');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 1; end;
        if ~isfield(app, 'dt'), app.dt = app.T/1000; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 50; end;
        
        Tracer_Advection_UW_demo;
        
    case 4 % Stokes equation: Analytical test case
        display('Analytical Stokes test case');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 2*pi; end;
        if ~isfield(app, 'dt'), app.dt = app.T/100; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 10; end;
        if ~isfield(app, 'nu'), app.nu = 1; end;
        
        SetupScript = 'StokesForcingSetup';
        PlotScript = 'StokesForcingPlots';
        
        solver = 1;
        
    case 5 % Stokes Poiseuille Flow test case
        display('Poiseuille Flow');
        
        if ~isfield(app, 'Nx'), app.Nx = 100; end;
        if ~isfield(app, 'Ny'), app.Ny = 50; end;
        if ~isfield(app, 'T'), app.T = 0.1; end;
        if ~isfield(app, 'dt'), app.dt = app.T/200; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 20; end;
        if ~isfield(app, 'nu'), app.nu = 1; end;
        
        SetupScript = 'StokesPoiseuilleSetup';
        PlotScript = 'StokesPoiseuillePlots';
        
        solver = 1;
        
    case 6 % Lock exchange experiment
        display('Lock Exchange');
        
        if ~isfield(app, 'Nx'), app.Nx = 50; end;
        if ~isfield(app, 'Ny'), app.Ny = 50; end;
        if ~isfield(app, 'T'), app.T = 10; end;
        if ~isfield(app, 'dt'), app.dt = 5e-3; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 50; end;
        if ~isfield(app, 'nu'), app.nu = 1e-4; end;
        if ~isfield(app, 'kappa'), app.kappa = 1e-4; end;
        if ~isfield(app, 'g'), app.g = 1; end;
        
        SetupScript = 'LockExchangeSetup';
        PlotScript = 'LockExchangePlots';
        
        solver = 1;
        
    case 6.1 % Lock exchange experiment Stochastic 1
        display('Stochastic Lock Exchange 1');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 2; end;
        if ~isfield(app, 'dt'), app.dt = 1e-3; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 100; end;
        if ~isfield(app, 'nu'), app.nu = 1e-4; end;
        if ~isfield(app, 'kappa'), app.kappa = 1e-4; end;
        if ~isfield(app, 'g'), app.g = 1; end;
        
        if ~isfield(app, 'S'), app.S = 2; end;
        if ~isfield(app, 'MC'), app.MC = 1000; end;
        
        SetupScript = 'LockExchangeSetup_DO';
        PlotScript = 'DORhoAdaptGenericPlots';
        
        solver = 3;
        
    case 6.2 % Lock exchange experiment Stochastic 2
        display('Stochastic Lock Exchange 2');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 2; end;
        if ~isfield(app, 'dt'), app.dt = 1e-3; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 100; end;
        if ~isfield(app, 'nu'), app.nu = 1e-4; end;
        if ~isfield(app, 'nu2'), app.nu2 = 1e-4; end;
        if ~isfield(app, 'kappa'), app.kappa = 1e-4; end;
        if ~isfield(app, 'kappa2'), app.kappa2 = 1e-4; end;
        if ~isfield(app, 'g'), app.g = 1; end;
        
        if ~isfield(app, 'S'), app.S = 2; end;
        if ~isfield(app, 'MC'), app.MC = 1000; end;
        
        SetupScript = 'LockExchangeSetup_DO2';
        PlotScript = 'LockExchangePlots_DO2';
        
        solver = 3;
        
    case 7 % Sudden Expansion
        display('Sudden Expansion (Look familiar?)');
        
        if ~isfield(app, 'Nx'), app.Nx = 80; end;
        if ~isfield(app, 'Ny'), app.Ny = 60; end;
        if ~isfield(app, 'T'), app.T = 50; end;
        if ~isfield(app, 'dt'), app.dt = 5/200; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 200; end;
        if ~isfield(app, 'nu'), app.nu = 1e-3; end;
        
        SetupScript = 'SuddenExpansionSetup';
        PlotScript = 'SuddenExpansionPlots';
        
        solver = 1;
        
    case 8 % Buoyant sudden expansion
        display('Buoyant Sudden Expansion');
        
        if ~isfield(app, 'Nx'), app.Nx = 80; end;
        if ~isfield(app, 'Ny'), app.Ny = 60; end;
        if ~isfield(app, 'T'), app.T = 50; end;
        if ~isfield(app, 'dt'), app.dt = 5/200; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 100; end;
        if ~isfield(app, 'nu'), app.nu = 1e-3; end;
        if ~isfield(app, 'kappa'), app.kappa = 1e-5; end;
        if ~isfield(app, 'g'), app.g = 1; end;
        
        SetupScript = 'Sudden_Expansion_Buoyant_Setup';
        PlotScript = 'Sudden_Expansion_Buoyant_Plot';
        
        solver = 1;
        
    case 9 % Buoyant bubble
        display('Buoyant Bubble');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 150; end;
        if ~isfield(app, 'T'), app.T = 1000; end;
        if ~isfield(app, 'dt'), app.dt = 1; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 50; end;
        if ~isfield(app, 'nu'), app.nu = 1e-4; end;
        if ~isfield(app, 'kappa'), app.kappa = 1e-4; end;
        if ~isfield(app, 'g'), app.g = 1; end;
        
        SetupScript = 'BuoyantBubbleSetup';
        PlotScript = 'BuoyantBubblePlots';
        
        solver = 1;
        
    case 10 % Lid Driven cavity flow
        display('Lid Driven Cavity Flow');
        display('For convergence test, run: LidDrivenCavityConvergenceTest.m')
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 100; end;
        if ~isfield(app, 'dt'), app.dt = 1/60; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 300; end;
        if ~isfield(app, 'nu'), app.nu = 1e-3; end;
        
        SetupScript = 'LidDrivenCavitySetup';
        PlotScript = 'LidDrivenCavityPlots';
        
        solver = 1;
        
    case 11 % Stochastic Lid Driven cavity flow
        display('Stochastic Lid Driven Cavity Flow');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 5; end;
        if ~isfield(app, 'dt'), app.dt = 1/200; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 40; end;
        if ~isfield(app, 'nu'), app.nu = 1e-3; end;
        
        if ~isfield(app, 'S'), app.S = 5; end;
        if ~isfield(app, 'var'), app.var = 1e-1; end;
        
        SetupScript = 'LidDrivenCavitySetup_DO';
        PlotScript = 'LidDrivenCavity_DO_Plots';
        
        solver = 2;
        
    case 12 % Adaptive Stochastic Lid Driven cavity flow
        display('Adaptive Stochastic Lid Driven Cavity Flow');
        display('Note, initialization is different than case 11');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 5; end;
        if ~isfield(app, 'dt'), app.dt = 1/200; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 40; end;
        if ~isfield(app, 'nu'), app.nu = 1e-3; end;
        
        if ~isfield(app, 'S'), app.S = [1, 4]; end;
        if ~isfield(app, 'var'), app.var = 1e-1; end;
        if ~isfield(app, 'DOAdapt'), app.DOAdapt = app.PlotIntrvl/2; end;
        
        SetupScript = 'LidDrivenCavitySetup';
        PlotScript = 'DOAdaptGenericPlots';
        solver = 2;
        
    case 13 % Double Gyre test case
        display('Deterministic Double Gyre Flow');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 1; end;
        if ~isfield(app, 'dt'), app.dt = 1e-4; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = app.T/app.dt/20; end;
        if ~isfield(app, 'nu'), app.nu = 1e-2; end;
        
        SetupScript = 'DoubleGyreSetup';
        PlotScript = 'DoubleGyrePlots';
        
        solver = 1;
        
    case 14 % Stochastic Double Gyre test case
        display('Stocastic Double Gyre Flow');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 1; end;
        if ~isfield(app, 'dt'), app.dt = 1e-4; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = app.T/app.dt/20; end;
        if ~isfield(app, 'nu'), app.nu = 1e-2; end;
        
        if ~isfield(app, 'S'), app.S = [1, 5]; end;
        if ~isfield(app, 'DOAdapt'), app.DOAdapt = 20; end;
        
        SetupScript = 'DoubleGyreSetup';
        PlotScript = 'DOAdaptGenericPlots';
        
        solver = 2;
        
    case 16 % PathPlanning with Analytical Flowfield
        display('PathPlanning with Analytical Flowfield');
        
        if ~isfield(app, 'Nx'), app.Nx = 100; end;
        if ~isfield(app, 'Ny'), app.Ny = 100; end;
        if ~isfield(app, 'T'), app.T = 200; end;
        if ~isfield(app, 'dt'), app.dt = 1e-3; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 100; end;
        if ~isfield(app, 'nu'), app.nu = 1e-2; end;
        
        
        app.V = @(x,y,t)  1.5*sqrt((x-0.4).^2 + (y-0.4).^2).*cos(PhaseMat(complex(x-0.4,y-0.4)));
        app.U = @(x,y,t)  -1.5*sqrt((x-0.4).^2 + (y-0.4).^2).*sin(PhaseMat(complex(x-0.4,y-0.4)));
        app.pathin.p = [0.4 0.4];
        app.pathin.v = [1 1]*1;
        app.Targets = [0.8 0.8];
        app.pathin.mag=sqrt(sum(app.pathin.v.^2));
        app.pathin.v = [0, app.pathin.mag];
        app.reinitLS = 1000;
        %app.reinitLS = floor(PlotIntrvl/30*100);
        
        SetupScript = 'PathPlanSetup';
        PlotScript = '';
        
        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');
        
        if ~isfield(app, 'Nx'), app.Nx = 30; end;
        if ~isfield(app, 'Ny'), app.Ny = 30; end;
        if ~isfield(app, 'T'), app.T = 20; end;
        if ~isfield(app, 'dt'), app.dt = 1/100; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 100; end;
        if ~isfield(app, 'nu'), app.nu = 1e-2; end;
        
        if ~isfield(app, 'g'), app.g = 1; end;
        if ~isfield(app, 'kappaT'), app.kappaT = 1e-2; end;
        if ~isfield(app, 'kappaS'), app.kappaS = 1e-2; end;
        
        SetupScript = 'ThermohalineSetup';
        PlotScript = 'ThermohalinePlots';
        
        solver = 5;
        
    case 19
        display('Flow past cylinder');
        
        if ~isfield(app, 'Nx'), app.Nx = 171; end;
        if ~isfield(app, 'Ny'), app.Ny = 21; end;
        if ~isfield(app, 'T'), app.T = 50; end;
        if ~isfield(app, 'dt'), app.dt = 0.01; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 100; end;
        if ~isfield(app, 'nu'), app.nu = 1e-2; end;
        
        SetupScript = 'CylinderCrossFlowSetup';
        PlotScript = 'SuddenExpansionPlots';
        
        solver = 1;
        
    case 20 % Rayleigh-Bernard Convection
        display('Rayleigh-Bernard Convection');
        
        if ~isfield(app, 'Nx'), app.Nx = 64; end;
        if ~isfield(app, 'Ny'), app.Ny = 64; end;
        if ~isfield(app, 'T'), app.T = 0.5; end;
        if ~isfield(app, 'dt'), app.dt = 1e-3; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 20; end;
        if ~isfield(app, 'nu'), app.nu = 0.7; end;
        
        if ~isfield(app, 'g'), app.g = -app.nu*1.5e4; end;
        if ~isfield(app, 'kappa'), app.kappa = 1; end;
        
        app.SetupScript = 'RayleighBenardSetup';
        app.PlotScript = 'RayleighBenardPlots';
        
        SetupScript = app.SetupScript;
        PlotScript = app.PlotScript;
        
        solver = 1;
        
    case 21 % Stochastic Rayleigh-Bernard Convection
        display('Stochastic Rayleigh-Bernard Convection');
        
        if ~isfield(app, 'Nx'), app.Nx = 64; end;
        if ~isfield(app, 'Ny'), app.Ny = 64; end;
        if ~isfield(app, 'T'), app.T = 1; end;
        if ~isfield(app, 'dt'), app.dt = 1e-3; end;
        if ~isfield(app, 'PlotIntrvl'), app.PlotIntrvl = 50; end;
        if ~isfield(app, 'nu'), app.nu = 0.7; end;
        
        if ~isfield(app, 'g'), app.g = -app.nu*1.5e4; end;
        if ~isfield(app, 'kappa'), app.kappa = 1; end;
        
        if ~isfield(app, 'S'), app.S = 3; end;
        if ~isfield(app, 'IP'), app.IP = [1,1,100]; end;
        
        SetupScript = 'RayleighBenardSetup_venturi_karniadakis_test';
        PlotScript = 'DORhoAdaptGenericPlots';
        % PlotScript = 'RayleighBenardPlots';
        
        solver = 3;
        
    otherwise
        error('No known problem selected');
end;

app

%Call solvers
if (problem > 3)
    switch solver
        case 1
            % Standard Fluids Solver
            display('Fluids solver')
            [P, u, v, rho, NodeP, Nodeu, Nodev] = ...
                FluidsSolver2(app, SetupScript, PlotScript);
            
        case 2
            % Stochastic Fluids Solver without density
            display('Stochastic Fluids solver without density');
            [P, u, v, NodeP, Nodeu, Nodev] = ...
                FluidsSolver2_DO_NS(app, SetupScript, PlotScript);    
            
        case 3
            % Stochastic Fluids Solver with density
            display('Stochastic Fluids solver with density');
            [P, u, v, rho, NodeP, Nodeu, Nodev] = ...
                FluidsSolver2_DO(app, SetupScript, PlotScript);
            
        case 4
            % Pathplanning with analytical flow-field
            display('PathPlanning using Level sets and Analytical flowfield Solver')
            [pathout, NodeP, Nodeu, Nodev] = ...
                Solver_pathplan_mod(app, SetupScript, PlotScript);
            
            display('Performing Backtracking Calculations');
            Backtrack(app,pathout);
            
        case 5
            % Boussinesq with seawater temperature and salinity (old solver)
            display('Boussinesq with seawater temperature and salinity instead of density');
            [P, u, v, rho, NodeP, Nodeu, Nodev, Tsw, Ssw] = ...
                BoussinesqTS(app, SetupScript, PlotScript);
        otherwise
            error('No known solver selected');
    end;
end;
