%clear all, clc, clf, close all
function [] = CylinderRunScript(Case)
%Function that submits case specific runs for the random inlet BCs flow
%over a square cylinder test case. Edit this file to add new cases.
%Currently, the cases are:
%Case = 1: A quick-test case -- just to check if the simulation runs. This
%          is low resolution with few time-steps, so the accuracy will be 
%          poor.
%       2: The standard run. Has reasonable resolution etc. A good starting
%          point.
%       3: A suggested modification to the case 2.

%The base or mean Reynolds number
Re = [105, 55, 43, 25, 73];

%The span of Reynolds numbers, i.e. the Re range for the calculation. 
%Re_span is an integer such that 0 < Re_span <= 1. 
%The range of RE will then be:
%  min(Re) = RE * (1 - Re_span)
%  max(Re) = RE * (1 + Re_span)
Re_span = [95./105., 45./55., 5./43., 15./25., 27./73.]; 

%The number of stochastic Monte-Carlo samples
MC = [1e6, 1e6, 1e6, 1e6, 1e6];

%The number of spatial points for 1 unit. That is, if the domain is 3x16,
%there will be 3N x 15N points in the domain. 
N =  [  15,  15,  15, 15, 15];

%The final non-dimensional time. 65 seems to give steady behaviour, but
%have gone as far as 186
T =  [60,  75,  75, 75, 75];

%The number of 'helper' DO modes. The total number of DO modes will be S+1,
%where the '1' comes from the 1 mode that introduces the stochastic
%boundary condition
S =  [  9,  9,  9,  9, 9];

%Multiplier for the timestep. The timestep size is automatically set for
%stability. To improve the accuracy, smaller timesteps could be used. That
%is, DT_MUL < 1. For DT_MUL > 1, the solution is likely to be unstable.
DT_MUL = [0.5, 1, 1, 1, 1]; 

%%
%Setup generics
d0=fileparts([pwd,filesep]);
d0=d0(1:strfind(d0,'trunk')+4);

%This script sets up the mesh, and boundary conditions. 
sscript = 'CylinderSetup_num';
%Use this setup script instead to RESTART a simulation
%sscript = 'CylinderRestart';

%This script controls the plotting of the outputs
pscript = 'CylinderPlots';
% pscript = '';


%% Assign the actual case-specific values
app.N = N(Case);
app.T = T(Case);
app.nu = 1/Re(Case);
Re = Re(Case);
app.Re_span = Re_span(Case);
app.S = S(Case);
app.MC = MC(Case);
app.var = 10^-16;
app.PlotScript = pscript;
app.DT_MUL = DT_MUL(Case);

%% Set the directories where images and .mat files will be saved
app.saveimg = ['randombcs/BConlyCylinder_RE',num2str(Re),'_Case',num2str(Case)];
app.savemat = ['randombcs/BConlyCylinder_RE',num2str(Re),'_Case',num2str(Case)];

%% Solve the equations
[P u ui v vi rho rhoi YYt NodeP Nodeu Nodev] = ...
    NavierStokesInc_DO_rhobcs(app,sscript,pscript);

%% Output some of the final plots
figure(2)
print('-dpng', '-r200', sprintf('%s/Save/%s/yyt_evol', d0, app.saveimg))
print('-deps', sprintf('%s/Save/%s/yyt_evol', d0, app.saveimg))
    
    
   

