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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   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)==hutterenum) + (md.elements_type(:,1)==macayealenum)  + (md.elements_type(:,1)==pattynenum)))
        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)==stokesenum) + (md.elements_type(:,2)==noneenum)))
        disp(['Types of elements have not been set properly, run setelementstype first'])
        bool=0;return;
end
if strcmpi(md.type,'2d'),
        if (ismember(pattynenum,md.elements_type(:,1)) |  ismember(stokesenum,md.elements_type(:,2))),
                disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
                bool=0;return;
        end
end

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

%FIELDS > 0 
fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
        'rho_ice','rho_water','B','elementoniceshelf','thickness','g','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(eval(['md.' char(fields(i))])),
                if find((eval(['md.' char(fields(i))]))<0),
                        disp(['model ' md.name ' has a <0 value in field ' char(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','segmentonneumann_diag','drag_type',...
        'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx',...
        'sparsity','deltaH','DeltaH','timeacc','timedec'};
for i=1:length(fields),
        if ~isempty(eval(['md.' char(fields(i))])),
                if find((eval(['md.' char(fields(i))]))==0),
                        disp(['model ' md.name ' has a =0 value in field ' char(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(eval(['md.' char(fields(i))]),1)~=md.numberofelements),
                disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']);
                bool=0; return;
        end
end

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

%THICKNESS = SURFACE - BED
if find((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 ~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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  SOLUTION CHECKS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%DIAGNOSTIC
if strcmpi(solutiontype,'diagnostic')

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

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

        %DIRICHLET VALUES
        if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
                disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
                bool=0;return;
        end

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

%PROGNOSTIC
if strcmp(solutiontype,'prognostic'),
        %old testing
        %       if isempty(find(md.gridondirichlet_diag)),
        %               disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!'])
        %               bool=0;return;
        %       end
        %       if isempty(find(md.gridondirichlet_prog)),
        %               disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!'])
        %               bool=0;return;
        %       end

        %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

%THERMAL TRANSIENT
if strcmp(solutiontype,'thermaltransient')

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

%THERMAL STEADY AND THERMAL TRANSIENT
if strcmpi(solutiontype,'thermalsteady') |  strcmp(solutiontype,'thermaltransient'),

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

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

%THERMAL TRANSIENT AND TRANSIENT
if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'transient'),

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

        %INITIAL TEMPERATURE
        if isstruct(md.temperature),
                disp(['The initial temperature should be empty or a list but not a structure'])
                bool=0;return;
        end
end

%PARAMETERS
if strcmp(solutiontype,'parameters')

        %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)==hutterenum); 
                disp(['The model has Hutter''s elements. Impossible to compute parameters']);
                bool=0;return;
        end
end

%CONTROL
if strcmpi(solutiontype,'control'),

        %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)
                disp('maxiter, optscal and fit 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 find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
                        disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
                        bool=0; return;
                end
        end
end

%QMU
if strcmpi(solutiontype,'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

%MESH
if strcmpi(solutiontype,'mesh'),
        %this solution is a little special. It should come right after the md=model;  operation. So a lot less checks!

        bool=1;
        return;
end

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

%CIELO
if strcmpi(package,'cielo'),

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

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

        %FIELD ~= 0
        fields={'time','sparsity'};
        for i=1:length(fields),
                if ~isempty(eval(['md.' char(fields(i))])),
                        if find((eval(['md.' char(fields(i))]))==0),
                                disp(['model ' md.name ' has a =0 value in field ' char(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

        %RIFTS
        if md.numrifts, 
                if ~strcmpi(md.cluster,'none')
                        if isempty(findstr('aijmumps',md.solverstring)),
                                disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']);
                                bool=0;return;
                        end
                end
        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

        %NP
        if ~strcmpi(md.cluster,'none'),
                if md.np==0,
                        disp(['model ' md.name ' has a =0 value in field np !']);
                elseif md.np<0,
                        disp(['model ' md.name ' has a negative value in field np !']);
                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
end

%No problems, just return 1;
bool=1;
return;