Changeset 123


Ignore:
Timestamp:
04/29/09 11:45:30 (16 years ago)
Author:
Mathieu Morlighem
Message:

reshaped runme.m in test 1

Location:
issm/trunk/test/Verification/IceSheetNoIceFront2d_1
Files:
7 added
4 deleted
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/test/Verification/IceSheetNoIceFront2d_1/runme.m

    r115 r123  
    1 % This file can be run to check that the current version of macayeal and issm models are giving
     1% This file can be run to check that the current version of issm is giving
    22% coherent results. This test deals with an icesheet without icefront for a 2d model. The geometry
    33% is square. Just run this file in Matlab, with a properly setup ISSM code.
    … …  
    66
    77% Errors  between archived results and the current version will get flagged if they are not within
    8 % a certain tolerance. The current tolerance is 10^-14. If you have good reasons to believe this
     8% a certain tolerance. The current tolerance is 10^-12. If you have good reasons to believe this
    99% tolerance should be lowered (for example, if you are running single precision compilers?), feel
    1010% free to tweak the tolerance variable.
    … …  
    1212%packages and solutions to be tested
    1313packages={'macayeal','ice','cielo_serial','cielo_parallel'};
     14solutions={'diagnostic','prognostic'};
    1415
    1516%Initialize log message for nightly runs.
    16 testname='IceSheetNoIceFront2d_1';
     17testname='IceSheetNoIceFrontM2d_16';
    1718tolerance=10^-12;
    1819
    … …  
    2122        package=packages{i};
    2223
    23         %initialize model
    24         md=model;
    25         md=mesh(md,'DomainOutline.exp',50000);
    26         md=geography(md,'','');
    27         md=parameterize(md,'Square.par');
    28         md=setelementstype(md,'macayeal','all');
     24        for j=1:length(solutions),
     25                solution=solutions{j};
    2926
    30         %Compute solution using requested package
    31         if strcmpi(package,'macayeal'),
    32 
    33                 %DIAGNOSTIC
    34                 md=solve(md,'diagnostic','macayeal');
    35                 %check result
    36                 load ArchiveMacAyealDiag
    37                 pos=find(ArchiveMacAyealDiag); error_vel=abs(norm((ArchiveMacAyealDiag(pos)-md.vel(pos))./ArchiveMacAyealDiag(pos),2));
    38                 if (error_vel>tolerance);
    39                         disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> ERROR (difference=' num2str(error_vel) ')'])
    40                 else
    41                         disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> SUCCESS (difference=' num2str(error_vel) ')'])
     27                %check package
     28                if ~(strcmpi(package,'macayeal') | strcmpi(package,'ice') | strcmpi(package,'cielo_serial') | strcmpi(package,'cielo_parallel'));
     29                        error(['package: ' package  ' in test: ' testname  ' not supported yet']);
     30                %check solution
     31                elseif ~(strcmpi(solution,'diagnostic') | strcmpi(solution,'prognostic'));
     32                        error(['solution: ' solution  ' in test: ' testname  ' not supported yet']);
     33                end
     34                %check solution requested
     35                if (~strcmpi(package,'ice') & strcmpi(solution,'prognostic')),
     36                        disp(sprintf(['\nsolution: ' solution  ', with package: ' package  ', in test: ' testname  ', not supported yet.\n']));
     37                        continue
    4238                end
    4339
    44         elseif strcmpi(package,'ice'),
    45 
    46                 %DIAGNOSTIC
    47                 md=solve(md,'diagnostic','ice');
    48                 %check result
    49                 load ArchiveIceDiag
    50                 pos=find(ArchiveIceDiag); error_vel=abs(norm((ArchiveIceDiag(pos)-md.vel(pos))./ArchiveIceDiag(pos),2));
    51                 if (error_vel>tolerance);
    52                         disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> ERROR (difference=' num2str(error_vel) ')'])
    53                 else
    54                         disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> SUCCESS (difference=' num2str(error_vel) ')'])
     40                %initialize model
     41                md=model;
     42                md=mesh(md,'DomainOutline.exp',50000);
     43                md=geography(md,'','');
     44                md=parameterize(md,'Square.par');
     45                md=setelementstype(md,'macayeal','all');
     46                if strcmpi(package,'cielo_parallel'), md.cluster='wilkes'; end
     47                if md.numberofgrids==388
     48                        load Velocities; md.vx=0.5*vx; md.vy=0.5*vy;
    5549                end
    5650
    57         elseif strcmpi(package,'cielo_serial'),
    58 
    59                 %DIAGNOSTIC
    60                 md=solve(md,'diagnostic','cielo');
    61                 %check result
    62                 load ArchiveCieloSerialDiag
    63                 pos=find(ArchiveCieloSerialDiag); error_vel=abs(norm((ArchiveCieloSerialDiag(pos)-md.vel(pos))./ArchiveCieloSerialDiag(pos),2));
    64                 if (error_vel>tolerance);
    65                         disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> ERROR (difference=' num2str(error_vel) ')'])
     51                %compute solution
     52                if strcmpi(package,'cielo_parallel') & strcmpi(solution,'diagnostic'),
     53                        md=solve(md,'diagnostic_horiz','cielo');
     54                elseif strcmpi(package,'cielo_serial'),
     55                        eval(['md=solve(md,''' solution ''',''cielo'');']);
    6656                else
    67                         disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> SUCCESS (difference=' num2str(error_vel) ')'])
     57                        eval(['md=solve(md,''' solution ''',''' package ''');']);
    6858                end
    6959
    70         elseif strcmpi(package,'cielo_parallel'),
    71 
    72                 %DIAGNOSTIC
    73                 md.cluster='wilkes';
    74                 md=solve(md,'diagnostic_horiz','cielo');
    75                 %check result
    76                 load ArchiveCieloParallelDiag
    77                 pos=find(ArchiveCieloParallelDiag); error_vel=abs(norm((ArchiveCieloParallelDiag(pos)-md.vel(pos))./ArchiveCieloParallelDiag(pos),2));
    78                 if (error_vel>tolerance);
    79                         disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> ERROR (difference=' num2str(error_vel) ')'])
    80                 else
    81                         disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> SUCCESS (difference=' num2str(error_vel) ')'])
     60                %compute field to be checked
     61                if strcmpi(solution,'diagnostic'),
     62                        fieldtest=md.vel;
     63                elseif strcmpi(solution,'prognostic'),
     64                        fieldtest=md.new_thickness';
    8265                end
    8366
    84         else
    85                 disp(sprintf(['\npackage ' package  ' not supported yet.\n']));
     67                %load archive
     68                eval(['load Archive' package solution ]);
     69
     70                %compare to archive
     71                eval(['pos=find(Archive' package solution ');']);
     72                eval(['error_diff=abs(norm((Archive' package solution '(pos)-fieldtest(pos))./Archive' package solution '(pos),2));']);
     73
     74                %disp test result
     75                if (error_diff>tolerance);
     76                        disp(['test: ' testname  ', solution: ' solution  ', package: ' package ' -> ERROR (difference=' num2str(error_diff) ' > ' num2str(tolerance) ')'])
     77                else
     78                        disp(['test: ' testname  ', solution: ' solution  ', package: ' package ' -> SUCCESS (difference=' num2str(error_diff) ' < ' num2str(tolerance) ')'])
     79                end
    8680        end
    8781end

  • issm/trunk/test/Verification/IceSheetNoIceFront2d_1/updatearchive.m

    r115 r123  
    88%packages and solutions to be tested
    99packages={'macayeal','ice','cielo_serial','cielo_parallel'};
     10solutions={'diagnostic','prognostic'};
    1011
    1112%go through the solutions requested
    … …  
    1314        package=packages{i};
    1415
    15         %initialize model
    16         md=model;
    17         md=mesh(md,'DomainOutline.exp',50000);
    18         md=geography(md,'','');
    19         md=parameterize(md,'Square.par');
    20         md=setelementstype(md,'macayeal','all');
     16        for j=1:length(solutions),
     17                solution=solutions{j};
    2118
    22         %Compute solution using requested package
    23         if strcmpi(package,'macayeal'),
     19                %check package
     20                if ~(strcmpi(package,'macayeal') | strcmpi(package,'ice') | strcmpi(package,'cielo_serial') | strcmpi(package,'cielo_parallel'));
     21                        error(['package: ' package  ' in test: ' testname  ' not supported yet']);
     22                        %check solution
     23                elseif ~(strcmpi(solution,'diagnostic') | strcmpi(solution,'prognostic'));
     24                        error(['solution: ' solution  ' in test: ' testname  ' not supported yet']);
     25                end
     26                %check solution requested
     27                if (~strcmpi(package,'ice') & strcmpi(solution,'prognostic')),
     28                        disp(sprintf(['\nsolution: ' solution  ', with package: ' package  ', in test: ' testname  ', not supported yet.\n']));
     29                        continue
     30                end
    2431
    25                 %DIAGNOSTIC
    26                 md=solve(md,'diagnostic','macayeal');
    27                 %save archive
    28                 ArchiveMacAyealDiag=md.vel;
    29                 save ArchiveMacAyealDiag ArchiveMacAyealDiag
     32                %initialize model
     33                md=model;
     34                md=mesh(md,'DomainOutline.exp',50000);
     35                md=geography(md,'','');
     36                md=parameterize(md,'Square.par');
     37                md=setelementstype(md,'macayeal','all');
     38                if strcmpi(package,'cielo_parallel'), md.cluster='wilkes'; end
     39                if md.numberofgrids==388
     40                        load Velocities; md.vx=0.5*vx; md.vy=0.5*vy;
     41                end
    3042
    31         elseif strcmpi(package,'ice'),
     43                %compute solution
     44                if strcmpi(package,'cielo_parallel') & strcmpi(solution,'diagnostic'),
     45                        md=solve(md,'diagnostic_horiz','cielo');
     46                elseif strcmpi(package,'cielo_serial'),
     47                        eval(['md=solve(md,''' solution ''',''cielo'');']);
     48                else
     49                        eval(['md=solve(md,''' solution ''',''' package ''');']);
     50                end
    3251
    33                 %DIAGNOSTIC
    34                 md=solve(md,'diagnostic','ice');
    35                 %save archive
    36                 ArchiveIceDiag=md.vel;
    37                 save ArchiveIceDiag ArchiveIceDiag
    38 
    39         elseif strcmpi(package,'cielo_serial'),
    40 
    41                 %DIAGNOSTIC
    42                 md=solve(md,'diagnostic','cielo');
    43                 %save archive
    44                 ArchiveCieloSerialDiag=md.vel;
    45                 save ArchiveCieloSerialDiag ArchiveCieloSerialDiag
    46 
    47         elseif strcmpi(package,'cielo_parallel'),
    48 
    49                 %DIAGNOSTIC
    50                 md.cluster='wilkes';
    51                 md=solve(md,'diagnostic_horiz','cielo');
    52                 %save archive
    53                 ArchiveCieloParallelDiag=md.vel;
    54                 save ArchiveCieloParallelDiag ArchiveCieloParallelDiag
    55 
    56         else
    57                 disp(sprintf(['\npackage ' package  ' not supported yet.\n']));
     52                %save new archive
     53                if strcmpi(solution,'diagnostic'),
     54                        fieldtest=md.vel;
     55                elseif strcmpi(solution,'prognostic'),
     56                        fieldtest=md.new_thickness';
     57                end
     58                eval(['Archive' package solution '=fieldtest;']);
     59                eval(['save Archive' package solution ' Archive' package solution]);
    5860        end
    5961end
Note: See TracChangeset for help on using the changeset viewer.