restart_nsdv.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[] = restart_nsdv(A)
% INPUT
% A Raw contents of a backup file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1 - Extract material for resuming the code %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Reading the backup file... ');
[unm1, unm2, rhonm1, rhonm2, pnm1, pnm2, phinm1, phinm2, unm1correction, tnm1, dtnm1_FE, dtn_FE, tend, ...
num_iter_start, MESH, RENUMBERING, refinement_lvl, idbackup_start, isave_start, imesh_start, ...
READPARAMETERS] = extract_backup_nsdv(A);
READPARAMETERS.RESTART = true;
PARAMETERS_SET = complete_manual_setup(READPARAMETERS);
truncate_logfile_nsdv(num_iter_start, PARAMETERS_SET{refinement_lvl});
fprintf('done\n');
switch PARAMETERS_SET{1}.CV_STUDY
case {'IN_REYNOLDS', 'IN_SPACE', 'IN_TIME', 'IN_SPACETIME'}
temp_output = PARAMETERS_SET{1}.OUTPUT.DIRECTORY_NAME(1:numel(PARAMETERS_SET{1}.OUTPUT.DIRECTORY_NAME)-numel(PARAMETERS_SET{1}.MODEL)-numel(PARAMETERS_SET{1}.TESTCASE)-4);
case {'NONE'}
temp_output = PARAMETERS_SET{1}.OUTPUT.DIRECTORY_NAME;
otherwise
error('Unknown value of CV_STUDY');
end
for kr=refinement_lvl:numel(PARAMETERS_SET)
PARAMETERS = PARAMETERS_SET{kr};
if PARAMETERS.FE.STEP_TIME_MODIFIED
fprintf('Warning: the time step has been modified to take into account absorbing boundary conditions\nwith a perturbed reference outlet profile\n');
end
if ~strcmp(PARAMETERS.CV_STUDY, 'NONE')
fprintf('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n');
fprintf('Start the simulation %d\n', kr);
fprintf('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n');
else
fprintf('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n');
fprintf('Start the simulation\n');
fprintf('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2 - Start a parallel environment if it is asked %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(PARAMETERS.PARALLELIZATION.TOOLBOX, 'PCT')
delete(gcp('nocreate'));
switch PARAMETERS.PARALLELIZATION.SERVER
case {'ZEUS'}
pc = parcluster('local');
pc.JobStorageLocation = strcat('/home/',getenv('USER'),'/.matlab/local_cluster_jobs/R', version('-release'));
parpool(PARAMETERS.PARALLELIZATION.NBWORKERS);
case {'MATHCALC', 'LOCAL'}
parpool(PARAMETERS.PARALLELIZATION.NBWORKERS);
otherwise
error(sprintf('Cluster profile not defined.\n'));
end
end
mess = sprintf('Mesh generation method: %s\n', PARAMETERS.MESH.GENERATION);
switch PARAMETERS.MESH.GENERATION
case {'NS2DDV'}
if isfield(PARAMETERS.MESH, 'NBSEG_PERUNIT_X')
mess = strcat(mess, ...
sprintf('\nNumber of segments per unit in x: %d\n', PARAMETERS.MESH.NBSEG_PERUNIT_X));
end
if isfield(PARAMETERS.MESH, 'NBNODES_INBL_Y')
mess = strcat(mess, ...
sprintf('\nNumber of segments in boundary layer in y: %d\n', ...
PARAMETERS.MESH.NBNODES_INBL_Y));
end
if isfield(PARAMETERS.MESH, 'NBSEG_OUTBL_Y')
mess = strcat(mess, ...
sprintf('\nNumber of segments out of boundary layer per unit in y: %d\n', ...
PARAMETERS.MESH.NBSEG_OUTBL_Y));
end
if isfield(PARAMETERS.MESH, 'NBSEG_X')
mess = strcat(mess, ...
sprintf('\nNumber of segments in x: %d\n', PARAMETERS.MESH.NBSEG_X));
end
if isfield(PARAMETERS.MESH, 'NBSEG_Y')
mess = strcat(mess, ...
sprintf('\nNumber of segments in y: %d\n', PARAMETERS.MESH.NBSEG_Y));
end
case {'PDET'}
if isfield(PARAMETERS.MESH, 'H0')
mess = strcat(mess, ...
sprintf('\nh0: %d\n', PARAMETERS.MESH.H0));
end
case {'FROM_FILE'}
if isfield(PARAMETERS.MESH, 'SRC_FILES')
mess = strcat(mess, sprintf('\nSource files: {'));
for p=1:numel(PARAMETERS.MESH.SRC_FILES)
mess = strcat(mess, sprintf('%s ', PARAMETERS.MESH.SRC_FILES{p}));
end
end
mess = strcat(mess, sprintf('}\n'));
otherwise
error('Unknown mesh generation method');
end
mess = strcat(mess, sprintf('\nhmin = %f\n', MESH.hmin));
mess = strcat(mess, sprintf('\nhmax = %f\n', MESH.hmax));
mess = strcat(mess, sprintf('\nnumber of P1 nodes = %d\n', MESH.np1));
mess = strcat(mess, sprintf('\nnumber of P1 triangles = %d\n', MESH.nt1));
disp(mess);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3 - Build the mesh and the renumbering matrices %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (kr > refinement_lvl)
% The mesh that has been extracted from the backup file is now useless
if strcmp(PARAMETERS.MESH.ADAPTATIVE, 'NO')
[MESH, RENUMBERING, PARAMETERS.DOMAIN] = build_mesh(PARAMETERS);
else
error(sprintf('Mesh adaptativity is not implemented yet.\n'));
end
end
%if ~PARAMETERS.FE.PSEUDOTRACTION
if PARAMETERS.FE.FIXPRESSURE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4 - Discretize the pressure constrain %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Identification of zero pressure node... ');
pressure_constrain = pressure_point_nsdv(MESH, PARAMETERS);
if strcmp(PARAMETERS.MESH.ADAPTATIVE, 'NO')
PARAMETERS.FE.NFX_X = pressure_constrain.nfx_x;
PARAMETERS.FE.NFX_Y = pressure_constrain.nfx_y;
PARAMETERS.FE.NFX = pressure_constrain.nfx;
else
error(sprintf('Mesh adaptativity is not implemented yet.\n'));
end
fprintf('done\n');
else
pressure_constrain = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5 - Identify the boundaries %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Discretization of boundary conditions... ');
BC = identify_BC_nsdv(MESH, pressure_constrain, PARAMETERS);
fprintf('done\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6 - Compute the Finite Element time step %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(PARAMETERS.INIT_FILE, '')
t0 = 0.;
else
diags{1} = 'u';
diags{2} = 'p';
diags{3} = 'METHOD_U';
diags{4} = 'time';
diags{5} = 'rho';
res = read_2D(PARAMETERS.INIT_FILE, diags);
if strcmp(PARAMETERS.RESET_TIME, 'YES')
t0 = 0.;
else
t0 = res{4};
end
end
PARAMETERS.FE.STEP_TIME = compute_dt(t0, MESH, PARAMETERS);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 7 - Compute the stationary matrices %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(PARAMETERS.MESH.ADAPTATIVE, 'NO')
% Matrix of int(grad(phi_i).grad(phi_j))
% phi_i are the DOFs for the velocity
% The approximation of the vector part is given in PARAMETERS.FE.TYPE (velocity discretization)
% Stiffness matrix for velocity (increase of 1 the Gauss quadrature precision order)
fprintf('Computing %s laplacian matrix... ', PARAMETERS.FE.TYPE);
CONSTANT_MATRICES.MLAPL_U = assembling_mlapl_analytic(MESH, PARAMETERS.FE.TYPE, @idscalar, 1, ...
PARAMETERS);
fprintf('done\n');
% Divergence matrices for (velocity,pressure)
% Matrices of int((d phi_j/dx)*psi_i) and int((d phi_j/dx)*psi_i)
% phi_i are the DOFs for the velocity
% psi_i are the DOFs for the pressure
% The approximation of the scalar part will be P1 (pressure discretization)
% The approximation of the vector part is given in PARAMETERS.FE.TYPE (velocity discretization)
fprintf('Computing divergence matrices... ');
[CONSTANT_MATRICES.DIVU_XPART, CONSTANT_MATRICES.DIVU_YPART] = ...
assembling_div(MESH, 'P1', PARAMETERS.FE.TYPE, PARAMETERS);
fprintf('done\n');
% Gradient matrices for (velocity,pressure)
% Matrices of int((d phi_j/dx)*psi_i) and int((d phi_j/dx)*psi_i)
% phi_i are the DOFs for the pressure
% psi_i are the DOFs for the velocity
% The approximation of the scalar part will be P1 (pressure discretization)
% The approximation of the vector part is given in PARAMETERS.FE.TYPE (velocity discretization)
fprintf('Computing gradient matrices... ');
[CONSTANT_MATRICES.GRADP_XPART, CONSTANT_MATRICES.GRADP_YPART] = ...
assembling_grad(MESH, 'P1', PARAMETERS.FE.TYPE, PARAMETERS);
fprintf('done\n');
% Matrix of int(phi_i*phi_j)
% phi_i are the DOFs for the velocity
% The approximation of the vector part is given in PARAMETERS.FE.TYPE (velocity discretization)
% Mass matrix for velocity (increase of 1 the Gauss quadrature precision order)
fprintf('Computing %s mass matrix... ', PARAMETERS.FE.TYPE);
CONSTANT_MATRICES.MASS_U = assembling_mass_analytic(MESH, PARAMETERS.FE.TYPE, @idscalar, 1, ...
PARAMETERS);
fprintf('done\n');
% P1 laplacian matrix for density (increase of 1 the Gauss quadrature precision order)
fprintf('Computing P1 laplacian matrix... ');
CONSTANT_MATRICES.MLAPL_P = assembling_mlapl_analytic(MESH, 'P1', @idscalar, 1, PARAMETERS);
fprintf('done\n');
% P1 mass matrix for density (increase of 1 the Gauss quadrature precision order)
fprintf('Computing P1 mass matrix... ');
CONSTANT_MATRICES.MASS_P = assembling_mass_analytic(MESH, 'P1', @idscalar, 1, PARAMETERS);
fprintf('done\n');
% If a null average condition is considered for the pressure, add a modified version of
% divergence matrices or laplacian matrices
if PARAMETERS.FE.FIXPRESSURE
% If a null average condition is considered for the pressure, add a modified version of
% divergence matrices or laplacian matrices
if strcmp(PARAMETERS.FE.PRESSURE_CONSTRAIN, 'ZERO_AVERAGE')
switch PARAMETERS.FE.SCHEME
case {'BDF2'}
% Add a modified version of divergence matrices
[CONSTANT_MATRICES.DIVU_XPART_NULLAV, CONSTANT_MATRICES.DIVU_YPART_NULLAV, ...
CONSTANT_MATRICES.NULLAV_P] ...
= null_average_divmatrix_ns(CONSTANT_MATRICES.DIVU_XPART, ...
CONSTANT_MATRICES.DIVU_YPART, pressure_constrain, MESH);
case {'BDF2_PROJ'}
% Add a modified version of laplacian matrix for pressure
CONSTANT_MATRICES.MLAPL_P_NULLAV = ...
null_average_mlaplmatrix_ns(CONSTANT_MATRICES.MLAPL_P, ...
pressure_constrain, MESH);
if isfield(PARAMETERS.PHYSICAL, 'FORCE_INCOMP')
if (PARAMETERS.PHYSICAL.FORCE_INCOMP > 0)
[CONSTANT_MATRICES.DIVU_XPART_NULLAV, ...
CONSTANT_MATRICES.DIVU_YPART_NULLAV, CONSTANT_MATRICES.NULLAV_P] ...
= null_average_divmatrix_ns(CONSTANT_MATRICES.DIVU_XPART, ...
CONSTANT_MATRICES.DIVU_YPART, pressure_constrain, MESH);
end
end
case {'BDF2_PROJ_VAR', 'BDF2_PROJ_VAR_A'}
% Nothing to do
otherwise
error('Unknown value for PARAMETERS.FE.SCHEME');
end
end
end
else
error(sprintf('Mesh adaptativity is not implemented yet.\n'));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8 - Compute the initial states %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
prefixfile = strcat(PARAMETERS.OUTPUT.DIRECTORY_NAME, '/', PARAMETERS.OUTPUT.FILE_NAME);
if (kr == refinement_lvl)
% The code crashed during a time loop : we use the backup data to resume the code
% Nothing more to compute
% Start the visualization interface
if max(strcmp(PARAMETERS.TESTCASE, {'RTIN', 'DROP'}))
nbrows = 1;
nbcols = numel(PARAMETERS.OUTPUT.XDISPLAY_LIST);
else
nbrows = ceil(sqrt(numel(PARAMETERS.OUTPUT.XDISPLAY_LIST)));
nbcols = ceil(numel(PARAMETERS.OUTPUT.XDISPLAY_LIST)/nbrows);
end
HANDLES = start_figs(PARAMETERS.OUTPUT.XDISPLAY_LIST, PARAMETERS.OUTPUT.XDISPLAY, nbrows, nbcols);
fprintf('Beginning of time loop\n');
else
% The time loop that has been crashed is now finished : we start a new time loop
% --> Compute (rho^(n-2),u^(n-2),p^(n-2)) and (rho^(n-1),u^(n-1),p^(n-1)) for n = 2
fprintf('Computing the initial states... ');
[u0,p0,rho0] = set_initialization_nsdv(MESH, PARAMETERS);
PARAMETERS.PHYSICAL.CHI = min(rho0);
fprintf('done\n');
fprintf('Beginning of time loop\n');
[unm1, unm2, pnm1, pnm2, rhonm1, rhonm2, phinm1, phinm2, unm1correction, isave_start, ...
imesh_start, HANDLES, num_iter_start, tnm1, dtn_FE, tend] ...
= initialization_time_splitting(u0, p0, rho0, t0, BC, ...
CONSTANT_MATRICES, MESH, prefixfile, PARAMETERS);
end
%%%%%%%%%%%%%%%%%%%%%%%%%
% 9 - Run the time loop %
%%%%%%%%%%%%%%%%%%%%%%%%%
idbackup = run_time_splitting_nsdv(unm1, unm2, rhonm1, rhonm2, pnm1, pnm2, phinm1, phinm2, ...
unm1correction, tnm1, dtn_FE, tend, HANDLES, MESH, BC, pressure_constrain, RENUMBERING, ...
CONSTANT_MATRICES, prefixfile, kr, idbackup_start, isave_start, num_iter_start, ...
imesh_start, PARAMETERS);
idbackup_start = idbackup;
fprintf('End of time loop\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 10 - Stop a parallel environment if it is asked %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(PARAMETERS.PARALLELIZATION.TOOLBOX, 'PCT')
delete(gcp('nocreate'));
end
if (kr < numel(PARAMETERS_SET))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 11 - Modify the parameter list for preparing a time/space mesh refinement %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stop_figs(HANDLES);
end
end
for kr=1:numel(PARAMETERS_SET)
OUTPUTDIRS{kr} = PARAMETERS_SET{kr}.OUTPUT.DIRECTORY_NAME;
end
last_postprocessing_nsdv(OUTPUTDIRS, PARAMETERS.CV_STUDY, temp_output);
end
