source: issm/trunk/src/m/classes/public/ismodelselfconsistent.m@ 2211

function bool=ismodelselfconsistent(md),
%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
%
%   Usage:
%      bool=ismodelselfconsistent(md),

%tolerance we use in litmus tests for the consistency of the model
tolerance=10^-12;
if (nargin~=1  )
        help ismodelselfconsistent
        error('ismodelselfconsistent error message: wrong usage');
end

%initialize output
bool=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   COMMON CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%COUNTER
if md.counter<3,
        disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']);
        bool=0;return;
end

%NAME
if isempty(md.name),
        disp(['model is not correctly configured: missing name!']);
        bool=0;return;
end

%MESH
if md.numberofelements<=0,
        disp(['model ' md.name ' does not have any elements!']);
        bool=0; return;
end
if md.numberofgrids<=0,
        disp(['model ' md.name ' does not have any grids!']);
        bool=0; return;
end

%ELEMENTSTYPE
if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
        disp(['Types of elements have not been set properly, run setelementstype first'])
        bool=0;return;
end
if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==HutterFormulationEnum) + (md.elements_type(:,1)==MacAyealFormulationEnum)  + (md.elements_type(:,1)==PattynFormulationEnum)))
        disp(['Types of elements have not been set properly, run setelementstype first'])
        bool=0;return;
end
if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==StokesFormulationEnum) + (md.elements_type(:,2)==NoneFormulationEnum)))
        disp(['Types of elements have not been set properly, run setelementstype first'])
        bool=0;return;
end
if strcmpi(md.type,'2d'),
        if (ismember(PattynFormulationEnum,md.elements_type(:,1)) |  ismember(StokesFormulationEnum,md.elements_type(:,2))),
                disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
                bool=0;return;
        end
end
if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
        disp(['no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
        bool=0;return;
end
if (md.analysis_type==DiagnosticAnalysisEnum & any(ismember(MacAyealFormulationEnum,md.elements_type(:,1)) & ismember(PattynFormulationEnum,md.elements_type(:,1))))
        disp(['coupling MacAyeal/Pattyn not implemented yet']);
        bool=0;return;
end
if (md.isstokes & md.analysis_type==TransientAnalysisEnum());
        disp(['Stokes transient not implemented yet']);
        bool=0;return;
end

%ICEFRONT
if strcmpi(md.type,'2d'),
        if size(md.pressureload,2)~=3 ,
                disp(['pressureload in ' md.name ' shoud have three columns']);
                bool=0;return;
        end
elseif strcmpi(md.type,'3d'),
        if size(md.pressureload,2)~=5 ,
                disp(['pressureload in ' md.name ' shoud have five columns']);
                bool=0;return;
        end
end

%NO NAN
fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',...
        'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
        'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
for i=1:length(fields),
        if ~isempty(md.(fields{i})),
                if any(isnan(md.(fields{i}))),
                        disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
                        bool=0; return;
                end
        end
end

%FIELDS >= 0 
fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',...
        'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
        'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
for i=1:length(fields),
        if ~isempty(md.(fields{i})),
                if any(md.(fields{i})<0),
                        disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
                        bool=0; return;
                end
        end
end
if any(md.p<=0),
        disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
        bool=0; return;
end

%FIELDS ~=0
fields={'numberofelements','numberofgrids','elements','drag_type',...
        'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
        'sparsity','deltaH','DeltaH','timeacc','timedec'};
for i=1:length(fields),
        if ~isempty(md.(fields{i})),
                if any(md.(fields{i})==0),
                        disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
                        bool=0; return;
                end
        end
end

%SIZE NUMBEROFELEMENTS
fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
for i=1:size(fields,2),
        if (size(md.(fields{i}),1)~=md.numberofelements),
                disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofelements) '!']);
                bool=0; return;
        end
end

%SIZE NUMBEROFGRIDS
fields={'x','y','z','B','drag','spcvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
for i=1:length(fields),
        if length(md.(fields{i}))~=md.numberofgrids,
                disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
                bool=0; return;
        end
end

%THICKNESS = SURFACE - BED
if any((md.thickness-md.surface+md.bed)>tolerance),
        disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
        bool=0; return;
end

%RIFTS
if md.numrifts,
        if ~strcmpi(md.type,'2d'),
                disp(['Models with rifts are only supported in 2d for now!']);
                bool=0;return;
        end
end
if ~isstruct(md.rifts),
        if ~isnan(md.rifts),
                if ~isempty(find(md.segmentmarkers>=2)),
                        %We have segments with rift markers, but no rift structure!
                        disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
                        bool=0; return;
                end
        end
end

%ARTIFICIAL DIFFUSIVITY
if ~isscalar(md.artificial_diffusivity),
        disp('artificial_diffusivity should be a scalar (1 or 0)');
        bool=0;return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  SOLUTION CHECKS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%QMU
if md.qmu_analysis,
        if md.qmu_params.evaluation_concurrency~=1,
                disp(['concurrency should be set to 1 when running dakota in library mode']);
                bool=0;return;
        end
        if ~isempty(md.part),
                if numel(md.part)~=md.numberofgrids,
                        disp(['user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
                        bool=0;return;
                end
                if find(md.part)>=md.numberofgrids,
                        disp(['user supplied partition should be indexed from 0 (c-convention)']);
                        bool=0;return;
                end
                if md.npart~=md.numberofgrids,
                        disp(['user supplied partition should have same size as md.npart']);
                        bool=0;return;
                end

        end
        if md.eps_rel>10^-3,
                disp(['for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
                bool=0;return;
        end
end

%DIAGNOSTIC
if md.analysis_type==DiagnosticAnalysisEnum,

        %CHECK THAT WE ARE NOT FULLY CONSTRAINED
        if strcmpi(md.type,'2d'),
                if isempty(find(~md.spcvelocity(:,1:2))),
                        disp(['model ' md.name ' is totally constrained, no need to solve!']);
                        bool=0;return;
                end
        end

        %HUTTER ON ICESHELF WARNING
        if any(md.elements_type(:,1)==HutterFormulationEnum & md.elementoniceshelf),
                disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
        end

        %SINGULAR
        if ~any(sum(md.spcvelocity(:,1:3),2)==3),
                disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
                bool=0;return;
        end

        %DIRICHLET IF THICKNESS <= 0
        if any(md.thickness<=0),
                pos=find(md.thickness<=0);
                if any(find(md.spcthickness(pos,1)==0)),
                        disp(['model ' md.name ' has some grids with 0 thickness']);
                        bool=0; return;
                end
        end
end

%PROGNOSTIC
if md.analysis_type==PrognosticAnalysisEnum,

        %VELOCITIES
        if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
                disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
                bool=0; return;
        end
end

%STEADYSTATE
if md.analysis_type==SteadystateAnalysisEnum,
        %NDT
        if md.dt~=0,
                disp(['For a steadystate computation, dt must be zero.']);
                bool=0;return;
        end
        %PRESSURE
        if isnans(md.pressure),
                disp(['For a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
                bool=0;return;
        end

        %eps: 
        if isnan(md.eps_rel),
                disp(['For a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
                bool=0;return;
        end

        %dim: 
        if strcmpi(md.type,'2d'),
                disp(['For a steadystate computation, model needs to be 3d']);
                bool=0;return;
        end
end

%THERMAL STEADY AND THERMAL TRANSIENT
if md.analysis_type==ThermalAnalysisEnum,

        %EXTRUSION
        if strcmp(md.type,'2d'),
                disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
                bool=0;return;
        end

        %VELOCITIES AND PRESSURE
        if (length(md.vx)~=md.numberofgrids | length(md.vy)~=md.numberofgrids | length(md.vz)~=md.numberofgrids),
                disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
                bool=0;return;
        end
        if (length(md.pressure)~=md.numberofgrids),
                disp(['pressure is required. Run ''diagnostic'' solution first!'])
                bool=0;return;
        end
end

%THERMAL TRANSIENT
if md.analysis_type==ThermalAnalysisEnum & md.dt~=0

        %DT and NDT
        fields={'dt','ndt'};
        for i=1:length(fields),
                if any(md.(fields{i})<0),
                        disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
                        bool=0; return;
                end
        end

        %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
        if (length(md.temperature)~=md.numberofgrids),
                disp(['An  initial temperature is needed for a transient thermal computation'])
                bool=0;return;
        end
        if (length(md.temperature)~=md.numberofgrids | length(md.accumulation)~=md.numberofgrids | length(md.melting)~=md.numberofgrids),
                disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
                bool=0;return;
        end
end

%PARAMETERS
if md.analysis_type==ParametersAnalysisEnum,

        %OUTPUT
        if ~iscell(md.parameteroutput)
                disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
                bool=0; return;
        end
        for i=1:length(md.parameteroutput)
                if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress')  & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
                                & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed')  & ~strcmpi(md.parameteroutput(i),'stress_surface')
                        disp(['one of the parameteroutput is not supported yet']);
                        bool=0; return;
                end
        end
        %VELOCITY
        if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
                disp(['velocities are required!']);
                bool=0;return;
        end

        %HUTTER
        if any(md.elements_type(:,1)==HutterFormulationEnum); 
                disp(['The model has Hutter''s elements. Impossible to compute parameters']);
                bool=0;return;
        end
end

%CONTROL
if md.control_analysis,

        %CONTROL TYPE
        if ~ischar(md.control_type),
                disp('control_type should be a string');
                bool=0;return;
        end

        %LENGTH CONTROL FIELDS
        if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps | length(md.cmjump)~=md.nsteps)
                disp('maxiter, optscal, fit and cmjump must have the length specified by nsteps')
                bool=0;return;
        end

        %FIT
        if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
                disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
                bool=0;return;
        end

        %OBSERVED VELOCITIES
        fields={'vx_obs','vy_obs'};
        for i=1:length(fields),
                if any(length(md.(fields{i}))~=md.numberofgrids),
                        disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
                        bool=0; return;
                end
        end

        %SINGULAR
        if ~any(sum(md.spcvelocity(:,1:3),2)==3),
                disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
                bool=0;return;
        end

        %DIRICHLET IF THICKNESS <= 0
        if any(md.thickness<=0),
                pos=find(md.thickness<=0);
                if any(find(md.spcthickness(pos,1)==0)),
                        disp(['model ' md.name ' has some grids with 0 thickness']);
                        bool=0; return;
                end
        end
end

%QMU
if strcmpi(md.analysis_type,'qmu'),
        if ~strcmpi(md.cluster,'none'),
                if md.waitonlock==0,
                        disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
                        bool=0;return;
                end
        end
end

%MESH2GRID
if md.analysis_type==Mesh2gridAnalysisEnum,
        if ~strcmpi(md.cluster,'none'),
                disp(['model is not correctly configured: mesh2grid not supported in parallel yet!']);
                bool=0;return;
        end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  PACKAGE CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%NAN VALUES
fields={'sparsity'};
for i=1:length(fields),
        if ~isempty(md.(fields{i})),
                if any(isnan(md.(fields{i}))),
                        disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
                        bool=0; return;
                end
        end
end

%FIELD > 0
fields={'sparsity'};
for i=1:length(fields),
        if ~isempty(md.(fields{i})),
                if any(md.(fields{i})<0),
                        disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
                        bool=0; return;
                end
        end
end

%FIELD ~= 0
fields={'sparsity'};
for i=1:length(fields),
        if ~isempty(md.(fields{i})),
                if any(md.(fields{i})==0),
                        disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
                        bool=0; return;
                end
        end
end

%SPARSITY BETWEEN 0 AND 1
if ( (md.sparsity<=0) | (md.sparsity>1)),
        disp(['model ' md.name ' sparsity should be inside the [0 1] range']);
        bool=0; return;
end

%CONNECTIVITY
if strcmpi(md.type,'2d'),
        if md.connectivity<9, 
                disp('connectivity should be at least 9 for 2d models');
                bool=0;return;
        end
end
if strcmpi(md.type,'3d'),
        if md.connectivity<24, 
                disp('connectivity should be at least 24 for 3d models');
                bool=0;return;
        end
end

%LOWMEM = 0 or 1
if ((md.lowmem ~= 1) & (md.lowmem~=0)),
        disp(['model ' md.name ' lowmem field should be 0 or 1']);
        bool=0; return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  PARALLEL CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~strcmpi(md.cluster,'none'),

        %NAN VALUES
        fields={'time','np'};
        for i=1:length(fields),
                if ~isempty(md.(fields{i})),
                        if any(isnan(md.(fields{i}))),
                                disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
                                bool=0; return;
                        end
                end
        end

        %FIELD > 0
        fields={'time','np'};
        for i=1:length(fields),
                if ~isempty(md.(fields{i})),
                        if any(md.(fields{i})<0),
                                disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
                                bool=0; return;
                        end
                end
        end

        %FIELD ~= 0
        fields={'time','np'};
        for i=1:length(fields),
                if ~isempty(md.(fields{i})),
                        if any(md.(fields{i})==0),
                                disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
                                bool=0; return;
                        end
                end
        end

end