compute_dt.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of NS2DDV. %
% %
% Copyright(C) 2011-2018 C. Calgaro (caterina.calgaro@math.univ-lille1.fr) %
% E. Creusé (emmanuel.creuse@math.univ-lille1.fr) %
% T. Goudon (thierry.goudon@inria.fr) %
% A. Mouton (alexandre.mouton@math.univ-lille1.fr) %
% %
% NS2DDV is free software: you can redistribute it and/or modify it under the terms %
% of the GNU General Public License as published by the Free Software Foundation, %
% either version 3 of the License, or (at your option) any later version. %
% %
% NS2DDV is distributed in the hope that it will be useful, but WITHOUT ANY %
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A %
% PARTICULAR PURPOSE. See the GNU General Public License for more details. %
% %
% You should have received a copy of the GNU General Public License along with %
% NS2DDV. If not, see <http://www.gnu.org/licenses/>. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[dt] = compute_dt(t0, MESH, PARAMETERS)
% INPUT
% PARAMETERS List of simulation parameters
% MESH Mesh metadata
% t0 Initial time
%
% OUTPUT
% dt Adapted FE time step
switch PARAMETERS.CV_STUDY
case {'NONE', 'IN_REYNOLDS', 'IN_SPACETIME'}
% Compute the time step thanks to the formula dt = Cmax*hmax^alphamax + Cmin*hmin^alphamin
% hmax and hmin are related to the current mesh, Cmax and Cmin are fixed across the simulation set
dtref = PARAMETERS.FE.CMAX_STEP_TIME*MESH.hmax^PARAMETERS.FE.ALPHAMAX_STEP_TIME + ...
PARAMETERS.FE.CMIN_STEP_TIME*MESH.hmin^PARAMETERS.FE.ALPHAMIN_STEP_TIME;
case {'IN_TIME'}
% Compute the time step thanks to the formula dt = C*hmax^alpha
% hmax is related to the mesh across the simulation set, C varies in the simulation set
dtref = PARAMETERS.FE.C_STEP_TIME*MESH.hmax^PARAMETERS.FE.ALPHA_STEP_TIME;
case {'IN_SPACE'}
% Compute the time step thanks to the formula dt = min(Cmax*hmax^alphamax + Cmin*hmin^alphamin)
% hmax and hmin vary across the simulation set, Cmax and Cmin are fixed
% --> We must estimate each hmin and hmax
switch PARAMETERS.MESH.GENERATION
case {'NS2DDV'}
switch PARAMETERS.DOMAIN.GEOMETRY
case {'RECTANGLE'}
if (isfield(PARAMETERS.MESH, 'NBSEG_X_LIST') && isfield(PARAMETERS.MESH, 'NBSEG_Y_LIST'))
dttest = zeros(1,numel(PARAMETERS.MESH.NBSEG_X_LIST));
for i=1:numel(PARAMETERS.MESH.NBSEG_X_LIST)
hx = (PARAMETERS.DOMAIN.XMAX-PARAMETERS.DOMAIN.XMIN) ...
/ PARAMETERS.MESH.NBSEG_X_LIST{i};
hy = (PARAMETERS.DOMAIN.YMAX-PARAMETERS.DOMAIN.YMIN) ...
/ PARAMETERS.MESH.NBSEG_Y_LIST{i};
hmax = max(hx,hy);
hmin = min(hx,hy);
dttest(i) = PARAMETERS.FE.CMAX_STEP_TIME * ...
hmax^PARAMETERS.FE.ALPHAMAX_STEP_TIME + ...
PARAMETERS.FE.CMIN_STEP_TIME * ...
hmin^PARAMETERS.FE.ALPHAMIN_STEP_TIME;
end
dtref = min(dttest);
else
error('PARAMETERS.MESH.NBSEG_X_LIST and PARAMETERS.MESH.NBSEG_Y_LIST are required');
end
case {'CIRCLE'}
if isfield(PARAMETERS.MESH, 'NBSEG_C_LIST')
dttest = zeros(1,numel(PARAMETERS.MESH.NBSEG_C_LIST));
for i=1:numel(PARAMETERS.MESH.NBSEG_C_LIST)
hmax = PARAMETERS.DOMAIN.RADIUS/PARAMETERS.MESH.NBSEG_C_LIST{i};
dttest(i) = PARAMETERS.FE.CMAX_STEP_TIME * ...
hmax^PARAMETERS.FE.ALPHAMAX_STEP_TIME + ...
PARAMETERS.FE.CMIN_STEP_TIME * ...
hmax^PARAMETERS.FE.ALPHAMIN_STEP_TIME;
end
dtref = min(dttest);
else
error('PARAMETERS.MESH.NBSEG_C_LIST is required');
end
otherwise
error('Only RECTANGLE and CIRCLE geometries can be used for convergence and/or stability studies');
end
case {'PDET'}
if isfield(PARAMETERS.MESH, 'H0_LIST')
dttest = zeros(1,numel(PARAMETERS.MESH.H0_LIST));
for i=1:numel(PARAMETERS.MESH.H0_LIST)
hmax = PARAMETERS.MESH.H0_LIST{i};
dttest(i) = PARAMETERS.FE.CMAX_STEP_TIME * ...
hmax^PARAMETERS.FE.ALPHAMAX_STEP_TIME + ...
PARAMETERS.FE.CMIN_STEP_TIME * ...
hmax^PARAMETERS.FE.ALPHAMIN_STEP_TIME;
end
dtref = min(dttest);
else
error('PARAMETERS.MESH.H0_LIST is required');
end
case {'FROM_FILE'}
error('Cannot identify the hmax of a mesh file set with PARAMETERS.MESH.GENERATION = ''FROM_FILE''');
end
otherwise
error('Unknown value of PARAMETERS.CV_STUDY');
end
if ~strcmp(PARAMETERS.FE.FIT_STEP_TIME, 'NONE')
if strcmp(PARAMETERS.MODEL, 'NSDV')
if max(strcmp(PARAMETERS.TESTCASE, {'EXAC'}))
switch PARAMETERS.TIME_SPLITTING
case {1}
N = 1+floor(PARAMETERS.FINAL_TIME/dtref);
dt = (PARAMETERS.FINAL_TIME-t0)/N;
case {2}
N = 1+floor((PARAMETERS.FINAL_TIME-t0)/(2.*dtref));
N = N + mod(N,2);
dt = (PARAMETERS.FINAL_TIME-t0)/(N+1);
otherwise
error('Unknown time splitting');
end
else
switch PARAMETERS.TIME_SPLITTING
case {1}
N = 1+floor((PARAMETERS.FINAL_TIME-t0)/dtref-dtref);
dt = (-N+sqrt(N*N+4.*(PARAMETERS.FINAL_TIME-t0)))*0.5;
case {2}
N = 1+floor((PARAMETERS.FINAL_TIME-t0)/(2.*dtref)-dtref);
dt = (-2.*N+sqrt(4.*N*N+8.*(PARAMETERS.FINAL_TIME-t0)))*0.25;
otherwise
error('Unknown time splitting');
end
end
elseif strcmp(PARAMETERS.MODEL, 'NS')
if max(strcmp(PARAMETERS.TESTCASE, {'EXAC','NONA','EXACNEU'}))
N = 1+floor((PARAMETERS.FINAL_TIME-t0)/dtref);
dt = (PARAMETERS.FINAL_TIME-t0)/N;
else
N = 1+floor((PARAMETERS.FINAL_TIME-t0)/dtref-dtref);
dt = (-N+sqrt(N*N+4.*(PARAMETERS.FINAL_TIME-t0)))*0.5;
end
else
error('Unknown model');
end
elseif strcmp(PARAMETERS.FE.FIT_STEP_TIME, 'NO')
dt = dtref;
else
error('Unknown value for PARAMETERS.FE.FIT_STEP_TIME');
end
end
