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

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

%tolerance we use in litmus checks for the consistency of the model
tolerance=10^-12;

%check usage
if nargin~=1,
        help ismodelselfconsistent
        error('ismodelselfconsistent error message: wrong usage');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   TRANSIENT   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%check analysis
if (md.analysis_type==Transient2DSolutionEnum | md.analysis_type==Transient3DSolutionEnum),
        if md.dt<=0,
                error('model not consistent: field dt must be positive for a transient run')
        end

        %recursive call to ismodelselfconsistent
        if (md.dim==2),
                analysis=[DiagnosticSolutionEnum PrognosticSolutionEnum];
        else
                analysis=[DiagnosticSolutionEnum PrognosticSolutionEnum ThermalSolutionEnum];
        end

        for i=1:length(analysis),
                md.analysis_type=analysis(i);
                ismodelselfconsistent(md);
        end
end


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

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

%NAME
if isempty(md.name),
        error(['model is not correctly configured: missing name!']);
end

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

%ICEFRONT
if (md.dim==2),
        fields={'pressureload'};
        checksize(md,fields,[NaN 4]);
elseif md.dim==3,
        fields={'pressureload'};
        checksize(md,fields,[NaN 6]);
else
        error('type should be either ''2d'' or ''3d''');
end

%ELEMENTS
fields={'elements'};
if (md.dim==2),
        checksize(md,fields,[md.numberofelements 3]);
else
        checksize(md,fields,[md.numberofelements 6]);
end

%NO NAN
fields={'numberofelements','numberofgrids','x','y','z','drag_coefficient','drag_type','drag_p','drag_q',...
        'rho_ice','rho_water','rheology_B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
        'tolx','np','eps_res','max_nonlinear_iterations','exclusive','rheology_n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
checknan(md,fields);

%FIELDS >= 0 
fields={'numberofelements','numberofgrids','elements','drag_coefficient','drag_type','drag_p','drag_q',...
        'rho_ice','rho_water','rheology_B','elementoniceshelf','thickness','g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
        'sparsity','lowmem','rheology_n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
checkgreater(md,fields,0);

%FIELDS > 0
fields={'numberofelements','numberofgrids','elements','drag_type','drag_p',...
        'rho_ice','rho_water','rheology_B','thickness','g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter','tolx',...
        'sparsity','deltaH','DeltaH','timeacc','timedec'};
checkgreaterstrict(md,fields,0);

%SIZE NUMBEROFELEMENTS
fields={'drag_p','drag_q','elementoniceshelf','rheology_n','elementonbed'};
checksize(md,fields,[md.numberofelements 1]);

%SIZE NUMBEROFGRIDS
fields={'x','y','z','rheology_B','drag_coefficient','melting_rate','accumulation_rate','surface','thickness','bed','gridonbed','gridonsurface'};
checksize(md,fields,[md.numberofgrids 1]);

%SIZE 6
fields={'spcvelocity'};
checksize(md,fields,[md.numberofgrids 6]);

%THICKNESS = SURFACE - BED
if any((md.thickness-md.surface+md.bed)>tolerance),
        error(['model not consistent: model ' md.name ' violates the equality thickness=surface-bed!']);
end

%RIFTS
if md.numrifts,
        if ~(md.dim==2),
                error(['model not consistent: models with rifts are only supported in 2d for now!']);
        end
        if ~isstruct(md.rifts),
                error(['model not consistent: md.rifts should be a structure!']);
        end
        if ~isempty(find(md.segmentmarkers>=2)),
                %We have segments with rift markers, but no rift structure!
                error(['model not consistent: model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
        end
        %Check that rifts are filled with proper material
        checkvalues(md,{'rifts.fill'},[WaterEnum() AirEnum() IceEnum() MelangeEnum()]);
else
        if ~isnans(md.rifts),
                error(['model not consistent: md.rifts shoud be NaN since md.numrifts is 0!']);
        end
end

%ARTIFICIAL DIFFUSIVITY
if ~ismember(md.artificial_diffusivity,[0 1]),
        error('model not consistent: artificial_diffusivity should be a scalar (1 or 0)');
end

%PARAMETEROUTPUT
if md.numoutput~=length(md.parameteroutput),
        error('model not consistent: numoutput should be the same size as parameteroutput');
end

%CONNECTIVITY
if (md.dim==2),
        if md.connectivity<9, 
                error('model not consistent: connectivity should be at least 9 for 2d models');
        end
end
if md.dim==3,
        if md.connectivity<24, 
                error('model not consistent: connectivity should be at least 24 for 3d models');
        end
end

%LOWMEM = 0 or 1
if ((md.lowmem ~= 1) & (md.lowmem~=0)),
        error(['model not consistent: model ' md.name ' lowmem field should be 0 or 1']);
end

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

        %NAN VALUES
        fields={'time','np'};
        checknan(md,fields);

        %FIELD > 0
        fields={'time','np'};
        checkgreaterstrict(md,fields,0);

end

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

%QMU
if md.qmu_analysis,
        if md.qmu_params.evaluation_concurrency~=1,
                error(['model not consistent: concurrency should be set to 1 when running dakota in library mode']);
        end
        if ~isempty(md.part),
                if numel(md.part)~=md.numberofgrids,
                        error(['model not consistent: user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
                end
                if find(md.part)>=md.numberofgrids,
                        error(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
                end
                if min(md.part)~=0,
                        error(['model not consistent: partition vector not indexed from 0 on']);
                end
                if max(md.part)>=md.numberofgrids,
                        error(['model not consistent: partition vector cannot have maximum index larger than number of grids']);
                end
                if ~isempty(find(md.part<0)),
                        error(['model not consistent: partition vector cannot have values less than 0']);
                end
                if ~isempty(find(md.part>=md.npart)),
                        error(['model not consistent: partition vector cannot have values more than md.npart-1']);
                end
                if max(md.part)>=md.npart,
                        error(['model not consistent: for qmu analysis, partitioning vector cannot go over npart, number of partition areas']);
                end
        end
        if md.eps_rel>1.1*10^-5,
                error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
        end
end

%DIAGNOSTIC
if md.analysis_type==DiagnosticSolutionEnum,

        %CHECK THAT WE ARE NOT FULLY CONSTRAINED
        if isempty(find(~md.spcvelocity(:,1:2))),
                error(['model not consistent: model ' md.name ' is totally constrained horizontally, no need to solve!']);
        end
        if (md.dim==3),
                if isempty(find(~md.spcvelocity(:,3))),
                        error(['model not consistent: model ' md.name ' is totally constrained vertically, no need to solve!']);
                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:2),2)==2),
                error(['model not consistent: model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
        end

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

        %INITIAL VELOCITY
        if length(md.vx)==md.numberofgrids & length(md.vy)==md.numberofgrids,
                fields={'vx','vy'};
                checknan(md,fields);
        end
end

%PROGNOSTIC
if md.analysis_type==PrognosticSolutionEnum,

        %INITIAL VELOCITIES
        fields={'vx','vy'};
        checksize(md,fields,[md.numberofgrids 1]);
        checknan(md,fields);

        %CHECK THAT WE ARE NOT FULLY CONSTRAINED
        if (md.dim==2),
                if isempty(find(~md.spcthickness(:,1))),
                        error(['model not consistent: model ' md.name ' is totally constrained for prognostic, no need to solve!']);
                end
        end

end

%STEADYSTATE
if md.analysis_type==SteadystateSolutionEnum,

        %NDT
        if md.dt~=0,
                error(['model not consistent: for a steadystate computation, dt must be zero.']);
        end

        %PRESSURE
        if isnans(md.pressure),
                error(['model not consistent: for a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
        end

        %eps: 
        if isnan(md.eps_rel),
                error(['model not consistent: for a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
        end

        %dim: 
        if (md.dim==2),
                error(['model not consistent: for a steadystate computation, model needs to be 3d']);
        end
end

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

        %EXTRUSION
        if (md.dim==2),
                error(['model not consistent: for a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
        end

        %CHECK THAT WE ARE NOT FULLY CONSTRAINED
        if isempty(find(~md.spctemperature(:,1))),
                error(['model not consistent: model ' md.name ' is totally constrained for temperature, no need to solve!']);
        end

        %VELOCITIES AND PRESSURE
        fields={'vx','vy','vz','pressure'};
        checksize(md,fields,[md.numberofgrids 1]);
        checknan(md,fields);

end

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

        %DT and NDT
        fields={'dt','ndt'};
        checkgreaterstrict(md,fields,0);

        %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
        fields={'temperature','accumulation_rate','melting_rate'};
        checksize(md,fields,[md.numberofgrids 1]);
        checknan(md,fields);

        %INITIAL TEMPERATURE
        fields={'temperature','spctemperature(:,2)','observed_temperature'};
        checkgreater(md,fields,0)

end

%BALANCEDTHICKNESS
if md.analysis_type==BalancedthicknessSolutionEnum

        %VELOCITIES MELTING AND ACCUMULATION
        fields={'vx','vy','accumulation_rate','melting_rate'};
        checksize(md,fields,[md.numberofgrids 1]);
        checknan(md,fields);

        %SPC                             
        if any(md.spcthickness(find(md.gridonboundary))~=1),             
                error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']);                      
        end 
end

%BALANCEDTHICKNESS2
if md.analysis_type==Balancedthickness2SolutionEnum

        %VELOCITIES MELTING AND ACCUMULATION
        fields={'vx','vy','accumulation_rate','melting_rate','dhdt'};
        checksize(md,fields,[md.numberofgrids 1]);
        checknan(md,fields);

         %CHECK THE COMPATIBILITY OF THE EDGES
         fields={'edges'};
         checksize(md,fields,[NaN 4]);
         checknan(md,fields);
end

%BALANCEDVELOCITIES
if md.analysis_type==BalancedvelocitiesSolutionEnum

        %VELOCITIES MELTING AND ACCUMULATION
        fields={'vx','vy','accumulation_rate','melting_rate'};
        checksize(md,fields,[md.numberofgrids 1]);
        checknan(md,fields);

        %SPC
        if any(md.spcvelocity(find(md.gridonboundary),[1:2])~=1),
                error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcvelocity']);
        end
end

%CONTROL
if md.control_analysis,

        %CONTROL TYPE
        if ~ischar(md.control_type),
                error('model not consistent: control_type should be a string');
        end

        %LENGTH CONTROL FIELDS
        fields={'maxiter','optscal','fit','cm_jump'};
        checksize(md,fields,[md.nsteps 1]);

        %FIT
        checkvalues(md,{'fit'},[0 1 2 3 4]);

        %WEIGHTS
        fields={'weights'};
        checksize(md,fields,[md.numberofgrids 1]);
        checkgreater(md,fields,0);

        %OBSERVED VELOCITIES
        fields={'vx_obs','vy_obs'};
        checksize(md,fields,[md.numberofgrids 1]);
        checknan(md,fields);

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

        %parameters
        fields={'cm_noisedmp'};
        checknan(md,fields);
end

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

end %end function

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHECK FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function checklength(md,fields,fieldlength)
        %CHECKSIZE - check length of a field
        for i=1:length(fields),
                if length(eval(['md.' fields{i}]))~=fieldlength,
                        error(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]);
                end
        end
end

function checksize(md,fields,fieldsize)
        %CHECKSIZE - check size of a field
        for i=1:length(fields),
                if isnan(fieldsize(1)),
                        if (size(eval(['md.' fields{i}]),2)~=fieldsize(2)),
                                %LOG warnings:
                                if strcmpi(fields{i},'pressureload'),
                                        disp(' ');
                                        disp('   If pressureload does not have the right size (one column missing), this is due to a change in md.pressure.');
                                        disp('      To update your model, use the following line');
                                        disp('      md.pressureload=[md.pressureload WaterEnum*md.elementoniceshelf(md.pressureload(:,end))+AirEnum*md.elementonicesheet(md.pressureload(:,end))];'); 
                                        disp(' ');
                                end
                                error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']);
                        end
                elseif isnan(fieldsize(2)),
                        if (size(eval(['md.' fields{i}]),1)~=fieldsize(1)),
                                error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']);
                        end
                else
                        if ((size(eval(['md.' fields{i}]),1)~=fieldsize(1)) |  (size(eval(['md.' fields{i}]),2)~=fieldsize(2)))
                                %LOG warnings:
                                if strcmpi(fields{i},'weights'),
                                        disp(' ');
                                        disp('   ''weights'' is a new field that should be of length numberofgrids.');
                                        disp('      To update your model, use the following line');
                                        disp('      md.weights=ones(md.numberofgrids,1)');
                                        disp(' ');
                                end
                                error(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]);
                        end
                end
        end
end

function checknan(md,fields)
        %CHECKNAN - check nan values of a field
        for i=1:length(fields),
                if any(isnan(eval(['md.' fields{i}]))),
                        error(['model not consistent: NaN values in field ' fields{i}]);
                end
        end
end

function checkreal(md,fields)
        %CHECKREAL - check real values of a field
        for i=1:length(fields),
                if any(eval(['~isreal(md.' fields{i} ')'])),
                        error(['model not consistent: complex values in field ' fields{i}]);
                end
        end
end

function checkgreaterstrict(md,fields,lowerbound)
        %CHECKGREATERSTRICT - check values of a field
        for i=1:length(fields),
                if any(eval(['md.' fields{i} '<=' num2str(lowerbound) ])),
                        error(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]);
                end
        end
end

function checkgreater(md,fields,lowerbound)
        %CHECKGREATER - check values of a field
        for i=1:length(fields),
                if any(eval(['md.' fields{i} '<' num2str(lowerbound) ])),
                        error(['model not consistent: field ' fields{i} ' should have values above ' num2str(lowerbound)]);
                end
        end
end

function checklessstrict(md,fields,upperbound)
        %CHECKLESSSTRICT - check values of a field
        for i=1:length(fields),
                if any(eval(['md.' fields{i} '>=' num2str(upperbound) ])),
                        error(['model not consistent: field ' fields{i} ' should have values stricly below ' num2str(upperbound)]);
                end
        end
end

function checkless(md,fields,upperbound)
        %CHECKLESS - check values of a field
        for i=1:length(fields),
                if any(eval(['md.' fields{i} '>' num2str(upperbound) ])),
                        error(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]);
                end
        end
end

function checkvalues(md,fields,values)
        %CHECKVALUE - check that a field has a certain value
        for i=1:length(fields),
                if eval(['~ismember( md.' fields{i} ',values)']),
                        error(['model not consistent: field ' fields{i} ' should have values in ' num2str(values)]);
                end
        end
end