Changeset 681


Ignore:
Timestamp:
06/01/09 16:41:17 (16 years ago)
Author:
Mathieu Morlighem
Message:

new results and model checks

Location:
issm/trunk/src/m/classes/public
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/classes/public/ismodelselfconsistent.m

    r632 r681  
    6767        'tolx','np','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
    6868for i=1:length(fields),
    69         if ~isempty(eval(['md.' char(fields(i))])),
    70                 if find(isnan(eval(['md.' char(fields(i))]))),
    71                         disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
     69        if ~isempty(eval(['md.' fields{i}])),
     70                if any(isnan(eval(['md.' fields{i}]))),
     71                        disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
    7272                        bool=0; return;
    7373                end
    … …  
    8080        'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
    8181for i=1:length(fields),
    82         if ~isempty(eval(['md.' char(fields(i))])),
    83                 if find((eval(['md.' char(fields(i))]))<0),
    84                         disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
     82        if ~isempty(eval(['md.' fields{i}])),
     83                if any((eval(['md.' fields{i}]))<0),
     84                        disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
    8585                        bool=0; return;
    8686                end
    … …  
    9797        'sparsity','deltaH','DeltaH','timeacc','timedec'};
    9898for i=1:length(fields),
    99         if ~isempty(eval(['md.' char(fields(i))])),
    100                 if find((eval(['md.' char(fields(i))]))==0),
    101                         disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
     99        if ~isempty(eval(['md.' fields{i}])),
     100                if any((eval(['md.' fields{i}]))==0),
     101                        disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
    102102                        bool=0; return;
    103103                end
    … …  
    108108fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
    109109for i=1:size(fields,2),
    110         if (size(eval(['md.' char(fields(i))]),1)~=md.numberofelements),
    111                 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']);
     110        if (size(eval(['md.' fields{i}]),1)~=md.numberofelements),
     111                disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofelements) '!']);
    112112                bool=0; return;
    113113        end
    … …  
    117117fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
    118118for i=1:length(fields),
    119         if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
    120                 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
     119        if length(eval(['md.' fields{i}]))~=md.numberofgrids,
     120                disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
    121121                bool=0; return;
    122122        end
    … …  
    124124
    125125%THICKNESS = SURFACE - BED
    126 if find((md.thickness-md.surface+md.bed)>tolerance),
     126if any((md.thickness-md.surface+md.bed)>tolerance),
    127127        disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
    128128        bool=0; return;
    … …  
    179179
    180180        %DIRICHLET IF THICKNESS <= 0
    181         if ~isempty(find(md.thickness<=0)),
     181        if any(md.thickness<=0),
    182182                pos=find(md.thickness<=0);
    183183                if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
    … …  
    249249        fields={'dt','ndt'};
    250250        for i=1:length(fields),
    251                 if find((eval(['md.' char(fields(i))]))<0),
    252                         disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
     251                if any((eval(['md.' fields{i}]))<0),
     252                        disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
    253253                        bool=0; return;
    254254                end
    … …  
    320320        fields={'vx_obs','vy_obs'};
    321321        for i=1:length(fields),
    322                 if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
    323                         disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
     322                if any(length(eval(['md.' fields{i}]))~=md.numberofgrids),
     323                        disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
    324324                        bool=0; return;
    325325                end
    … …  
    362362        fields={'sparsity'};
    363363        for i=1:length(fields),
    364                 if ~isempty(eval(['md.' char(fields(i))])),
    365                         if find(isnan(eval(['md.' char(fields(i))]))),
    366                                 disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
     364                if ~isempty(eval(['md.' fields{i}])),
     365                        if any(isnan(eval(['md.' fields{i}]))),
     366                                disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
    367367                                bool=0; return;
    368368                        end
    … …  
    373373        fields={'sparsity'};
    374374        for i=1:length(fields),
    375                 if ~isempty(eval(['md.' char(fields(i))])),
    376                         if find((eval(['md.' char(fields(i))]))<0),
    377                                 disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
     375                if ~isempty(eval(['md.' fields{i}])),
     376                        if any((eval(['md.' fields{i}]))<0),
     377                                disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
    378378                                bool=0; return;
    379379                        end
    … …  
    384384        fields={'sparsity'};
    385385        for i=1:length(fields),
    386                 if ~isempty(eval(['md.' char(fields(i))])),
    387                         if find((eval(['md.' char(fields(i))]))==0),
    388                                 disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
     386                if ~isempty(eval(['md.' fields{i}])),
     387                        if any((eval(['md.' fields{i}]))==0),
     388                                disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
    389389                                bool=0; return;
    390390                        end
    … …  
    435435                fields={'time','np'};
    436436                for i=1:length(fields),
    437                         if ~isempty(eval(['md.' char(fields(i))])),
    438                                 if find(isnan(eval(['md.' char(fields(i))]))),
    439                                         disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
     437                        if ~isempty(eval(['md.' fields{i}])),
     438                                if any(isnan(eval(['md.' fields{i}]))),
     439                                        disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
    440440                                        bool=0; return;
    441441                                end
    … …  
    446446                fields={'time','np'};
    447447                for i=1:length(fields),
    448                         if ~isempty(eval(['md.' char(fields(i))])),
    449                                 if find((eval(['md.' char(fields(i))]))<0),
    450                                         disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
     448                        if ~isempty(eval(['md.' fields{i}])),
     449                                if any((eval(['md.' fields{i}]))<0),
     450                                        disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
    451451                                        bool=0; return;
    452452                                end
    … …  
    457457                fields={'time','np'};
    458458                for i=1:length(fields),
    459                         if ~isempty(eval(['md.' char(fields(i))])),
    460                                 if find((eval(['md.' char(fields(i))]))==0),
    461                                         disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
     459                        if ~isempty(eval(['md.' fields{i}])),
     460                                if any((eval(['md.' fields{i}]))==0),
     461                                        disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
    462462                                        bool=0; return;
    463463                                end

  • issm/trunk/src/m/classes/public/isresultconsistent.m

    r465 r681  
    1 function bool=isresultconsistent(md)
     1function bool=isresultconsistent(md,solution)
    22%ISRESULTCONSISTENT: check that results are consistent
    33%
    … …  
    66%
    77%   Usage:
    8 %      bool=IsModelSelfConsistent(model)
     8%      bool=IsModelSelfConsistent(model,solution)
    99
    1010%tolerance we use in litmus tests for the consistency of the model
    1111tolerance=10^-2;
    1212
    13 if (nargin~=1  )
    14         IsResultConsistentUsage();
     13%some checks
     14if (nargin~=2 | nargout~=1)
     15        help isresultconsistent
    1516        error(' ');
    1617end
    17 %Check velocities i(no penalties for 2d meshes)
    18 if strcmp(md.analysis_type,'diagnostic'),
    19 
    20         if strcmpi(md.type,'3d'),
    21                 if ~isnan(md.penalties),
    22                         for i=1:size(md.penalties,1)
    23                                 for n=2:size(md.penalties,2)
    24                                         relativex=(md.vx(md.penalties(i,1))-md.vx(md.penalties(i,n)))./(md.vx(md.penalties(i,1)));
    25                                         if md.vx(md.penalties(i,1))==md.vx(md.penalties(i,n)), relativex=0; end
    26                                         relativey=(md.vy(md.penalties(i,1))-md.vy(md.penalties(i,n)))./(md.vy(md.penalties(i,1)));
    27                                         if md.vy(md.penalties(i,1))==md.vy(md.penalties(i,n)), relativey=0; end
    28                                         if abs(relativex)>tolerance | abs(relativey)>tolerance
    29                                                 disp(['''diagnostic'' result not consistent: increase penalty_offset (discontinuous horizontal velocity)']);
    30                                                 bool=0; return;
    31                                         end
    32 
    33                                 end
    34                         end
    35                 end
    36         end
    37 end
    38 
    39 if strcmp(md.analysis_type,'thermal')
    40 
    41         %check melting
    42         if any(abs(md.melting(md.numberofgrids2d+1:end))>tolerance)
    43                 disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']);
     18
     19%initialize output as 1 -> success
     20bool=1;
     21
     22%DIAGNOSTIC
     23if strcmp(solution,'diagnostic'),
     24
     25        if strcmpi(md.type,'3d')
     26                fields1={'results.diagnostic.vx','results.diagnostic.vy','results.diagnostic.vz','results.diagnostic.vel','results.diagnostic.pressure'};
     27                fields2={'results.diagnostic.vel','results.diagnostic.pressure'};
     28        else
     29                fields1={'results.diagnostic.vx','results.diagnostic.vy','results.diagnostic.vel'};
     30                fields2={'results.diagnostic.vel'};
     31        end
     32
     33        %check size
     34        if ~testsize(md,fields1,md.numberofgrids),
     35                bool=0; return
     36        end
     37
     38        %no NAN
     39        if ~testnan(md,fields1),
     40                bool=0; return
     41        end
     42
     43        %check value is real
     44        if ~testreal(md,fields1),
     45                bool=0; return
     46        end
     47
     48        %check value>=0
     49        if ~testpositive(md,fields2),
     50                bool=0; return
     51        end
     52
     53end
     54
     55%CONTROL
     56if strcmp(solution,'control'),
     57
     58        fields1={'results.control.vx','results.control.vy','results.control.vel','results.control.parameter'};
     59        fields2={'results.control.vel','results.control.J'};
     60
     61        %check size
     62        if ~testsize(md,fields1,md.numberofgrids),
     63                bool=0; return
     64        end
     65
     66        %no NAN
     67        if ~testnan(md,fields1),
     68                bool=0; return
     69        end
     70
     71        %check value is real
     72        if ~testreal(md,fields1),
     73                bool=0; return
     74        end
     75
     76        %check value>=0
     77        if ~testpositive(md,fields2),
     78                bool=0; return
     79        end
     80
     81        %check inversed parameter
     82        if any(md.control.parameter<md.mincontrolconstraint | md.control.parameter>md.maxcontrolconstraint)
     83                disp(['''control'' result not consistent: inverse parameter out of range [' md.mincontrolconstraint ' ' md.maxcontrolconstraint ']']);
    4484                bool=0; return;
    4585        end
    … …  
    4787end
    4888
    49 if strcmp(md.analysis_type,'thermaltransient')
    50 
    51         for iter=1:length(md.thermaltransient_results)
    52 
    53         %check melting
    54         if any(abs(md.thermaltransient_results(iter).melting(md.numberofgrids2d+1:end))>tolerance)
    55                 disp(['''thermaltransient'' result not consistent: increase penalty_melting (negative melting)']);
    56                 bool=0; return;
    57         end
    58 
    59         end
    60 end
    61 
    62 
    63 if (strcmp(md.analysis_type,'transient') & strcmpi(md.type,'3d')),
    64         for iter=1:length(md.transient_results)
    65 
    66                 %check melting
    67                 if any(abs(md.transient_results(iter).melting(md.numberofgrids2d+1:end))>tolerance)
    68                         disp(['''thermaltransient'' result not consistent: increase penalty_melting (negative melting)']);
    69                         bool=0; return;
    70                 end
    71 
    72 
    73                 %check velocities
    74                 if ~isnan(md.penalties),
    75                         for i=1:size(md.penalties,1)
    76                                 for n=2:size(md.penalties,2)
    77                                         relativex=(md.transient_results(iter).vx(md.penalties(i,1))-md.transient_results(iter).vx(md.penalties(i,n)))./(md.transient_results(iter).vx(md.penalties(i,1)));
    78                                         if md.transient_results(iter).vx(md.penalties(i,1))==md.transient_results(iter).vx(md.penalties(i,n)), relativex=0; end
    79                                         relativey=(md.transient_results(iter).vy(md.penalties(i,1))-md.transient_results(iter).vy(md.penalties(i,n)))./(md.transient_results(iter).vy(md.penalties(i,1)));
    80                                         if md.transient_results(iter).vy(md.penalties(i,1))==md.transient_results(iter).vy(md.penalties(i,n)), relativey=0; end
    81                                         if abs(relativex)>tolerance | abs(relativey)>tolerance
    82                                                 disp(['''transient'' result not consistent: increase penalty_offset (discontinuous horizontal velocity)']);
    83                                                 bool=0; return;
    84                                         end
    85 
    86                                 end
    87                         end
    88                 end
    89 
    90         end
    91 
    92 end
    93 
    94 
    95 %No problems, just return 1;
    96 bool=1;
    97 return;
    98        
    99 function IsResultConsistentUsage()
    100 disp(' ');
    101 disp('   IsResultConsistent usage: flag=IsResultConsistent(model)');
    102 disp('         where model is an instance of the @model class, and flag a boolean');
     89%THERMAL
     90if strcmp(solution,'thermal')
     91
     92        for iter=1:length(md.results.thermal)
     93
     94                fields1={['results.thermal(' num2str(iter) ').temperature'],['results.thermal(' num2str(iter) ').melting']};
     95                fields2={['results.thermal(' num2str(iter) ').temperature']};
     96
     97                %check size
     98                if ~testsize(md,fields1,md.numberofgrids),
     99                        bool=0; return
     100                end
     101
     102                %no NAN
     103                if ~testnan(md,fields1,md.numberofgrids),
     104                        bool=0; return
     105                end
     106
     107                %check value is real
     108                if ~testreal(md,fields1,md.numberofgrids),
     109                        bool=0; return
     110                end
     111
     112                %check value>=0
     113                if ~testpositive(md,fields2,md.numberofgrids),
     114                        bool=0; return
     115                end
     116
     117                %check melting (<=0 via penalties)
     118                if any(abs(md.results.thermal(iter).melting(md.numberofgrids2d+1:end))>tolerance)
     119                        disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']);
     120                        bool=0; return;
     121                end
     122        end
     123
     124end
     125
     126if strcmp(solution,'transient')
     127
     128        for iter=1:length(md.results.transient)
     129
     130                if strcmpi(md.type,'3d'),
     131                        fields1={['results.transient(' num2str(iter) ').vx'],['results.transient(' num2str(iter) ').vy'],...
     132                                ['results.transient(' num2str(iter) ').vz'],['results.transient(' num2str(iter) ').vel'],...
     133                                ['results.transient(' num2str(iter) ').bed'],['results.transient(' num2str(iter) ').surface'],...
     134                                ['results.transient(' num2str(iter) ').thickness'],...
     135                                ['results.transient(' num2str(iter) ').temperature'],['results.transient(' num2str(iter) ').melting']};
     136                        fields2={['results.transient(' num2str(iter) ').vel'],['results.transient(' num2str(iter) ').thickness'],...
     137                                ['results.transient(' num2str(iter) ').temperature']};
     138                else
     139                        fields1={['results.transient(' num2str(iter) ').vx'],['results.transient(' num2str(iter) ').vy'],...
     140                                ['results.transient(' num2str(iter) ').vz'],['results.transient(' num2str(iter) ').vel'],...
     141                                ['results.transient(' num2str(iter) ').bed'],['results.transient(' num2str(iter) ').surface'],...
     142                                ['results.transient(' num2str(iter) ').thickness']};
     143                        fields2={['results.transient(' num2str(iter) ').vel'],['results.transient(' num2str(iter) ').thickness']};
     144                end
     145
     146                %check size
     147                if ~testsize(md,fields1,md.numberofgrids),
     148                        bool=0; return
     149                end
     150
     151                %no NAN
     152                if ~testnan(md,fields1,md.numberofgrids),
     153                        bool=0; return
     154                end
     155
     156                %check value is real
     157                if ~testreal(md,fields1,md.numberofgrids),
     158                        bool=0; return
     159                end
     160
     161                %check value>=0
     162                if ~testpositive(md,fields2,md.numberofgrids),
     163                        bool=0; return
     164                end
     165
     166                %check melting (<=0 via penalties)
     167                if any(abs(md.results.transient(iter).melting(md.numberofgrids2d+1:end))>tolerance)
     168                        disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']);
     169                        bool=0; return;
     170                end
     171        end
     172end
     173end % end function
     174
     175function bool=testsize(md,fields,sizefield)
     176        %TESTSIZE - test size of a field
     177        bool=1;
     178        for i=1:length(fields),
     179                if length(eval(['md.' fields{i}]))~=sizefield
     180                        disp(['results not consistent: field ' fields{i} ' size should be ' num2str(md.numberofgrids)]);
     181                        bool=0; return;
     182                end
     183        end
     184end
     185
     186function bool=testnan(md,fields,nanfield)
     187        %TESTNAN - test nan values of a field
     188        bool=1;
     189        for i=1:length(fields),
     190                if any(isnan(eval(['md.' fields{i}]))),
     191                        disp(['results not consistent: NaN values in field ' fields{i}]);
     192                        bool=0; return;
     193                end
     194        end
     195end
     196
     197function bool=testreal(md,fields,realfield)
     198        %TESTREAL - test real values of a field
     199        bool=1;
     200        for i=1:length(fields),
     201                if any(eval(['~isreal(md.' fields{i} ')'])),
     202                        disp(['results not consistent: complex values in field ' fields{i}]);
     203                        bool=0; return;
     204                end
     205        end
     206end
     207
     208function bool=testpositive(md,fields,positivefield)
     209        %TESTPOSITIVE - test positive values of a field
     210        bool=1;
     211        for i=1:length(fields),
     212                if any(eval(['md.' fields{i} '<0'])),
     213                        disp(['results not consistent: negative values in field ' fields{i}]);
     214                        bool=0; return;
     215                end
     216        end
     217end

  • issm/trunk/src/m/classes/public/loadresultsfromdisk.m

    r673 r681  
    1616%Check result is consistent
    1717disp(sprintf('%s\n','checking result consistency'));
    18 if ~isresultconsistent(md),
     18if ~isresultconsistent(md,md.analysis_type),
    1919        %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
    2020        disp('!! results not consistent correct the model !!')

  • issm/trunk/src/m/classes/public/solve.m

    r629 r681  
    8282%Check result is consistent
    8383displaystring(md.debug,'%s\n','checking result consistency');
    84 if ~isresultconsistent(md),
     84if ~isresultconsistent(md,options.analysis_type),
    8585        disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
    8686end
Note: See TracChangeset for help on using the changeset viewer.