Changeset 14134


Ignore:
Timestamp:
12/11/12 15:25:15 (12 years ago)
Author:
jschierm
Message:

NEW: Working python tests 1101-1108 (plus other cosmetic changes).

Location:
issm/trunk-jpl/test
Files:
14 added
23 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt

    r14107 r14134  
    2424propagate into the sixth digit of the thickness, which is enough to fail NR.
    2525
     26test1109,test1110
     27Ugly crashes in solver, but same behavior as Matlab.
     28

  • issm/trunk-jpl/test/NightlyRun/test1101.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%Pattyn and Payne 2006
    33printingflag=false;
    44
    5 L_list={5000,10000,20000,40000,80000,160000};
     5L_list={5000.,10000.,20000.,40000.,80000.,160000.};
    66results={};
    77minvx=[];
    … …  
    2626
    2727        pos=find(md.mesh.vertexonbed);
    28         md.diagnostic.spcvx(pos)=0;
    29         md.diagnostic.spcvy(pos)=0;
     28        md.diagnostic.spcvx(pos)=0.;
     29        md.diagnostic.spcvy(pos)=0.;
    3030
    3131        %Create MPCs to have periodic boundary conditions
    32         posx=find(md.mesh.x==0);
     32        posx=find(md.mesh.x==0.);
    3333        posx2=find(md.mesh.x==max(md.mesh.x));
    3434
    35         posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
    36         posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
     35        posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
     36        posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
    3737
    3838        md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
    … …  
    5555                set(gcf,'Color','w')
    5656                printmodel(['ismipapattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    57                 system(['mv ismipapattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     57                system(['mv ismipapattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    5858        end
    5959        plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
    … …  
    6161                set(gcf,'Color','w')
    6262                printmodel(['ismipapattynvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    63                 system(['mv ismipapattynvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     63                system(['mv ismipapattynvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    6464        end
    6565        plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
    … …  
    6767                set(gcf,'Color','w')
    6868                printmodel(['ismipapattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    69                 system(['mv ismipapattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     69                system(['mv ismipapattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    7070        end
    7171
    72         if(L==5000),
     72        if(L==5000.),
    7373                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
    7474                        'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','')
    75         elseif(L==10000),
     75        elseif(L==10000.),
    7676                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
    7777                        'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','')
    78         elseif(L==20000),
     78        elseif(L==20000.),
    7979                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
    8080                        'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','')
    81         elseif(L==40000),
     81        elseif(L==40000.),
    8282                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
    8383                        'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','')
    84         elseif(L==80000),
     84        elseif(L==80000.),
    8585                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
    8686                        'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','')
    87         elseif(L==160000),
     87        elseif(L==160000.),
    8888                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
    8989                        'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
    … …  
    9292                set(gcf,'Color','w')
    9393                printmodel(['ismipapattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    94                 system(['mv ismipapattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     94                system(['mv ismipapattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    9595        end
    9696end
    … …  
    101101        set(gcf,'Color','w')
    102102        printmodel('ismipapattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    103         system(['mv ismipapattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     103        system(['mv ismipapattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    104104end
    105105plot([5 10 20 40 80 160],maxvx);ylim([0 120]);xlim([0 160])
    … …  
    107107        set(gcf,'Color','w')
    108108        printmodel('ismipapattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    109         system(['mv ismipapattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     109        system(['mv ismipapattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    110110end
     111
    111112%Fields and tolerances to track changes
    112 field_names     ={ ...
     113field_names     ={...
    113114        'Vx5km','Vy5km','Vz5km',...
    114115        'Vx10km','Vy10km','Vz10km',...

  • issm/trunk-jpl/test/NightlyRun/test1102.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%Pattyn and Payne 2006
    33printingflag=false;
    44
    5 L_list={5000,10000,20000,40000,80000,160000};
     5L_list={5000.,10000.,20000.,40000.,80000.,160000.};
    66results={};
    77minvx=[];
    … …  
    1717
    1818%       %Find elements at the corner and extract model
    19 %       posnodes=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
     19%       posnodes=find((md.mesh.x==0. | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0. | md.mesh.y==max(md.mesh.y)));
    2020%       [a,b]=find(ismember(md.mesh.elements,posnodes));
    2121%       elements=ones(md.mesh.numberofelements,1);
    … …  
    2929        %Create dirichlet on the bed only
    3030        pos=find(md.mesh.vertexonbed);
    31         md.diagnostic.spcvx(pos)=0;
    32         md.diagnostic.spcvy(pos)=0;
    33         md.diagnostic.spcvz(pos)=0;
     31        md.diagnostic.spcvx(pos)=0.;
     32        md.diagnostic.spcvy(pos)=0.;
     33        md.diagnostic.spcvz(pos)=0.;
    3434
    35         %Create MPCs to have periodic boundary conditions
    36         %posx=find(md.mesh.x==0);
    37         %posx2=find(md.mesh.x==max(md.mesh.x));
    38         %posx=find(md.mesh.x==0 & md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
    39         %posx2=find(md.mesh.x==max(md.mesh.x) &  md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
     35%       %Create MPCs to have periodic boundary conditions
     36%       posx=find(md.mesh.x==0.);
     37%       posx2=find(md.mesh.x==max(md.mesh.x));
     38%       posx=find(md.mesh.x==0. & md.mesh.y~=0. & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
     39%       posx2=find(md.mesh.x==max(md.mesh.x) &  md.mesh.y~=0. & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
    4040
    41         %posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed); %Don't take the same nodes two times
    42         %posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed);
     41%       posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed); %Don't take the same nodes two times
     42%       posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed);
    4343
    44         %md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
     44%       md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
    4545
    4646        %Compute the diagnostic
    4747        md.diagnostic.abstol=NaN;
    4848        md.diagnostic.reltol=NaN;
    49         md.diagnostic.restol=1;
     49        md.diagnostic.restol=1.;
    5050        md.cluster=generic('name',oshostname(),'np',8);
    5151        md=solve(md,DiagnosticSolutionEnum());
    … …  
    6262        %Now plot vx, vy, vz and vx on a cross section
    6363        plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
    64         if printingflag, 
     64        if printingflag,
    6565                set(gcf,'Color','w')
    6666                printmodel(['ismipastokesvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    67                 system(['mv ismipastokesvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     67                system(['mv ismipastokesvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    6868        end
    6969        plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
    70         if printingflag, 
     70        if printingflag,
    7171                set(gcf,'Color','w')
    7272                printmodel(['ismipastokesvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    73                 system(['mv ismipastokesvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     73                system(['mv ismipastokesvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    7474        end
    7575        plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4)
    76         if printingflag, 
     76        if printingflag,
    7777                set(gcf,'Color','w')
    7878                printmodel(['ismipastokesvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    79                 system(['mv ismipastokesvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     79                system(['mv ismipastokesvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    8080        end
    8181
    82         if(L==5000),
     82        if(L==5000.),
    8383                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
    8484                        'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','')
    85         elseif(L==10000),
     85        elseif(L==10000.),
    8686                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
    8787                        'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','')
    88         elseif(L==20000),
     88        elseif(L==20000.),
    8989                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
    9090                        'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','')
    91         elseif(L==40000),
     91        elseif(L==40000.),
    9292                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
    9393                        'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','')
    94         elseif(L==80000),
     94        elseif(L==80000.),
    9595                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
    9696                        'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','')
    97         elseif(L==160000),
     97        elseif(L==160000.),
    9898                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
    9999                        'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
    100100        end
    101         if printingflag, 
     101        if printingflag,
    102102                set(gcf,'Color','w')
    103103                printmodel(['ismipastokesvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    104                 system(['mv ismipastokesvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     104                system(['mv ismipastokesvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    105105        end
    106106end
    … …  
    108108%Now plot the min and max values of vx for each size of the square
    109109plot([5 10 20 40 80 160],minvx);ylim([0 18])
    110 if printingflag, 
     110if printingflag,
    111111        set(gcf,'Color','w')
    112112        printmodel('ismipastokesminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    113         system(['mv ismipastokesminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     113        system(['mv ismipastokesminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    114114end
    115115plot([5 10 20 40 80 160],maxvx);ylim([0 120])
    116 if printingflag, 
     116if printingflag,
    117117        set(gcf,'Color','w')
    118118        printmodel('ismipastokesmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    119         system(['mv ismipastokesmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
     119        system(['mv ismipastokesmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
    120120end
     121
    121122%Fields and tolerances to track changes
    122123field_names     ={...
    … …  
    127128        'Vx80km','Vy80km','Vz80km',...
    128129        'Vx160km','Vy160km','Vz160km'
    129 }
     130};
    130131field_tolerances={...
    131132        1e-12,1e-12,1e-12,...

  • issm/trunk-jpl/test/NightlyRun/test1103.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%Pattyn and Payne 2006
    33printingflag=false;
    44
    5 L_list={5000,10000,20000,40000,80000,160000};
     5L_list={5000.,10000.,20000.,40000.,80000.,160000.};
    66results={};
     7minvx=[];
     8maxvx=[];
    79
    810for i=1:length(L_list),
    … …  
    2325        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
    2426        pos=find(md.mesh.vertexonbed);
    25         md.diagnostic.spcvx(pos)=0;
    26         md.diagnostic.spcvy(pos)=0;
     27        md.diagnostic.spcvx(pos)=0.;
     28        md.diagnostic.spcvy(pos)=0.;
    2729
    2830        %Create MPCs to have periodic boundary conditions
    29         posx=find(md.mesh.x==0);
     31        posx=find(md.mesh.x==0.);
    3032        posx2=find(md.mesh.x==max(md.mesh.x));
    3133
    32         posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
    33         posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
     34        posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
     35        posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
    3436
    3537        md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
    … …  
    4951        %Now plot vx, vy, vz and vx on a cross section
    5052        plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
    51         if printingflag, 
     53        if printingflag,
    5254                set(gcf,'Color','w')
    5355                printmodel(['ismipbpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    54                 system(['mv ismipbpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
     56                system(['mv ismipbpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
    5557        end
    5658        plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
    57         if printingflag, 
     59        if printingflag,
    5860                set(gcf,'Color','w')
    5961                printmodel(['ismipbpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    60                 system(['mv ismipbpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
     62                system(['mv ismipbpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
    6163        end
    6264
    63         if(L==5000),
     65        if(L==5000.),
    6466                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
    6567                        'resolution',[10 10],'ylim',[6 16],'xlim',[0 5000],'title','','xlabel','')
    66         elseif(L==10000),
     68        elseif(L==10000.),
    6769                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
    6870                        'resolution',[10 10],'ylim',[0 40],'xlim',[0 10000],'title','','xlabel','')
    69         elseif(L==20000),
     71        elseif(L==20000.),
    7072                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
    7173                        'resolution',[10 10],'ylim',[0 60],'xlim',[0 20000],'title','','xlabel','')
    72         elseif(L==40000),
     74        elseif(L==40000.),
    7375                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
    7476                        'resolution',[10 10],'ylim',[0 100],'xlim',[0 40000],'title','','xlabel','')
    75         elseif(L==80000),
     77        elseif(L==80000.),
    7678                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
    7779                        'resolution',[10 10],'ylim',[0 120],'xlim',[0 80000],'title','','xlabel','')
    78         elseif(L==160000),
     80        elseif(L==160000.),
    7981                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
    8082                        'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
    8183        end
    82         if printingflag, 
     84        if printingflag,
    8385                set(gcf,'Color','w')
    8486                printmodel(['ismipbpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    85                 system(['mv ismipbpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
     87                system(['mv ismipbpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
    8688        end
    8789end
    … …  
    8991%Now plot the min and max values of vx for each size of the square
    9092plot([5 10 20 40 80 160],minvx);ylim([0 14]);xlim([0 160])
    91 if printingflag, 
     93if printingflag,
    9294        set(gcf,'Color','w')
    9395        printmodel('ismipbpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    94         system(['mv ismipbpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
     96        system(['mv ismipbpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
    9597end
    9698plot([5 10 20 40 80 160],maxvx);ylim([0 120]);xlim([0 160])
    97 if printingflag, 
     99if printingflag,
    98100        set(gcf,'Color','w')
    99101        printmodel('ismipbpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    100         system(['mv ismipbpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
     102        system(['mv ismipbpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
    101103end
     104
    102105%Fields and tolerances to track changes
    103106field_names     ={...

  • issm/trunk-jpl/test/NightlyRun/test1104.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%Pattyn and Payne 2006
    33
    4 L_list={5000,10000,20000,40000,80000,160000};
     4L_list={5000.,10000.,20000.,40000.,80000.,160000.};
    55results={};
    66
    … …  
    2222
    2323        pos=find(md.mesh.vertexonbed);
    24         md.diagnostic.spcvx(pos)=0;
    25         md.diagnostic.spcvy(pos)=0;
     24        md.diagnostic.spcvx(pos)=0.;
     25        md.diagnostic.spcvy(pos)=0.;
    2626
    2727        %Create MPCs to have periodic boundary conditions
    28         posx=find(md.mesh.x==0);
     28        posx=find(md.mesh.x==0.);
    2929        posx2=find(md.mesh.x==max(md.mesh.x));
    3030
    31         posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
    32         posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
     31        posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
     32        posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
    3333
    3434        md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
    … …  
    3838        md.cluster=generic('name',oshostname(),'np',8);
    3939        md=solve(md,DiagnosticSolutionEnum());
    40         pos=find(md.mesh.x==0 | md.mesh.y==0 | md.mesh.x==max(md.mesh.x) | md.mesh.y==max(md.mesh.y));
     40        pos=find(md.mesh.x==0. | md.mesh.y==0. | md.mesh.x==max(md.mesh.x) | md.mesh.y==max(md.mesh.y));
    4141        md.diagnostic.spcvx(pos)=md.results.DiagnosticSolution.Vx(pos);
    4242        md.diagnostic.spcvy(pos)=md.results.DiagnosticSolution.Vy(pos);

  • issm/trunk-jpl/test/NightlyRun/test1105.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%Pattyn and Payne 2006
    33printingflag=false;
    44
    5 L_list={5000,10000,20000,40000,80000,160000};
     5L_list={5000.,10000.,20000.,40000.,80000.,160000.};
    66results={};
    77minvx=[];
    … …  
    99
    1010for i=1:length(L_list),
    11         L=L_list{i};  %in m (3 times the desired lenght for BC problems) 
     11        L=L_list{i};  %in m (3 times the desired length for BC problems) 
    1212        nx=30; %number of nodes in x direction
    1313        ny=30;
    … …  
    1818        md=extrude(md,10,1.);
    1919
    20         md=setflowequation(md,'pattyn','all'); 
     20        md=setflowequation(md,'pattyn','all');
    2121
    2222        %Create MPCs to have periodic boundary conditions
    … …  
    2525        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
    2626
    27         posx=find(md.mesh.x==0 & md.mesh.y~=0 & md.mesh.y~=L);
    28         posx2=find(md.mesh.x==L & md.mesh.y~=0 & md.mesh.y~=L);
     27        posx=find(md.mesh.x==0. & md.mesh.y~=0. & md.mesh.y~=L);
     28        posx2=find(md.mesh.x==L & md.mesh.y~=0. & md.mesh.y~=L);
    2929
    30         posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=L); %Don't take the same nodes two times
    31         posy2=find(md.mesh.y==L & md.mesh.x~=0 & md.mesh.x~=L);
     30        posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=L); %Don't take the same nodes two times
     31        posy2=find(md.mesh.y==L & md.mesh.x~=0. & md.mesh.x~=L);
    3232
    3333        md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
    3434
    3535        %Add spc on the corners
    36         pos=find((md.mesh.x==0 | md.mesh.x==L) & (md.mesh.y==0 | md.mesh.y==L) & md.mesh.vertexonbed);
    37         md.diagnostic.spcvx(pos)=0;
    38         md.diagnostic.spcvy(pos)=0;
    39         if(L==5000),
     36        pos=find((md.mesh.x==0. | md.mesh.x==L) & (md.mesh.y==0. | md.mesh.y==L) & md.mesh.vertexonbed);
     37        md.diagnostic.spcvx(pos)=0.;
     38        md.diagnostic.spcvy(pos)=0.;
     39        if(L==5000.),
    4040                md.diagnostic.spcvx(pos)=15.66;
    4141                md.diagnostic.spcvy(pos)=-0.1967;
    42         elseif(L==10000),
     42        elseif(L==10000.),
    4343                md.diagnostic.spcvx(pos)=16.04;
    4444                md.diagnostic.spcvy(pos)=-0.1977;
    45         elseif(L==20000),
     45        elseif(L==20000.),
    4646                md.diagnostic.spcvx(pos)=16.53;
    4747                md.diagnostic.spcvy(pos)=-1.27;
    48         elseif(L==40000),
     48        elseif(L==40000.),
    4949                md.diagnostic.spcvx(pos)=17.23;
    5050                md.diagnostic.spcvy(pos)=-3.17;
    51         elseif(L==80000),
     51        elseif(L==80000.),
    5252                md.diagnostic.spcvx(pos)=16.68;
    5353                md.diagnostic.spcvy(pos)=-2.69;
    54         elseif(L==160000),
     54        elseif(L==160000.),
    5555                md.diagnostic.spcvx(pos)=16.03;
    5656                md.diagnostic.spcvy(pos)=-1.27;
    … …  
    5959        %Spc the bed at zero for vz
    6060        pos=find(md.mesh.vertexonbed);
    61         md.diagnostic.spcvz(pos)=0;
     61        md.diagnostic.spcvz(pos)=0.;
    6262
    6363        %Compute the diagnostic
    … …  
    7575        %Now plot vx, vy, vz and vx on a cross section
    7676        plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
    77         if printingflag, 
     77        if printingflag,
    7878                set(gcf,'Color','w')
    7979                printmodel(['ismipcpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    80                 system(['mv ismipcpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
     80                system(['mv ismipcpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
    8181        end
    8282        plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
    83         if printingflag, 
     83        if printingflag,
    8484                set(gcf,'Color','w')
    8585                printmodel(['ismipcpattynvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    86                 system(['mv ismipcpattynvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
     86                system(['mv ismipcpattynvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
    8787        end
    8888        plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4)
    89         if printingflag, 
     89        if printingflag,
    9090                set(gcf,'Color','w')
    9191                printmodel(['ismipcpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    92                 system(['mv ismipcpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
     92                system(['mv ismipcpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
    9393        end
    9494
    95         if(L==5000),
     95        if(L==5000.),
    9696                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
    9797                        'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',5)
    98         elseif(L==10000),
     98        elseif(L==10000.),
    9999                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
    100100                        'resolution',[10 10],'ylim',[13 18],'xlim',[0 10000],'title','','xlabel','')
    101         elseif(L==20000),
     101        elseif(L==20000.),
    102102                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
    103103                        'resolution',[10 10],'ylim',[14 22],'xlim',[0 20000],'title','','xlabel','')
    104         elseif(L==40000),
     104        elseif(L==40000.),
    105105                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
    106106                        'resolution',[10 10],'ylim',[10 40],'xlim',[0 40000],'title','','xlabel','')
    107         elseif(L==80000),
     107        elseif(L==80000.),
    108108                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
    109109                        'resolution',[10 10],'ylim',[0 80],'xlim',[0 80000],'title','','xlabel','')
    110         elseif(L==160000),
     110        elseif(L==160000.),
    111111                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
    112112                        'resolution',[10 10],'ylim',[0 200],'xlim',[0 160000],'title','','xlabel','')
    113113        end
    114         if printingflag, 
     114        if printingflag,
    115115                set(gcf,'Color','w')
    116116                printmodel(['ismipcpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    117                 system(['mv ismipcpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
     117                system(['mv ismipcpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
    118118        end
    119119end
    … …  
    121121%Now plot the min and max values of vx for each size of the square
    122122plot([5 10 20 40 80 160],minvx);ylim([4 18]);xlim([0 160])
    123 if printingflag, 
     123if printingflag,
    124124        set(gcf,'Color','w')
    125125        printmodel('ismipcpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    126         system(['mv ismipcpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
     126        system(['mv ismipcpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
    127127end
    128128plot([5 10 20 40 80 160],maxvx);ylim([0 200]); xlim([0 160])
    129 if printingflag, 
     129if printingflag,
    130130        set(gcf,'Color','w')
    131131        printmodel('ismipcpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    132         system(['mv ismipcpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
     132        system(['mv ismipcpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
    133133end
     134
    134135%Fields and tolerances to track changes
    135136field_names     ={...

  • issm/trunk-jpl/test/NightlyRun/test1106.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%Pattyn and Payne 2006
    33
    4 L_list={5000,10000,20000,40000,80000,160000};
     4L_list={5000.,10000.,20000.,40000.,80000.,160000.};
    55results={};
    66
    … …  
    1010        md=setmask(md,'',''); %ice sheet test
    1111        md=parameterize(md,'../Par/ISMIPC.par');
    12         md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/L).*sin(md.mesh.y*2*pi/L)));
     12        md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/L).*sin(md.mesh.y*2.*pi/L)));
    1313        md=extrude(md,10,1.);
    1414
    1515        %Add spc on the borders
    16         pos=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x) | md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
    17         md.diagnostic.spcvx(pos)=0;
    18         md.diagnostic.spcvy(pos)=0;
    19         if(L==5000),
     16        pos=find(md.mesh.x==0. | md.mesh.x==max(md.mesh.x) | md.mesh.y==0. | md.mesh.y==max(md.mesh.y));
     17        md.diagnostic.spcvx(pos)=0.;
     18        md.diagnostic.spcvy(pos)=0.;
     19        if(L==5000.),
    2020                md.diagnostic.spcvx(pos)=15.66;
    2121                md.diagnostic.spcvy(pos)=-0.1967;
    22         elseif(L==10000),
     22        elseif(L==10000.),
    2323                md.diagnostic.spcvx(pos)=16.04;
    2424                md.diagnostic.spcvy(pos)=-0.1977;
    25         elseif(L==20000),
     25        elseif(L==20000.),
    2626                md.diagnostic.spcvx(pos)=16.53;
    2727                md.diagnostic.spcvy(pos)=-1.27;
    28         elseif(L==40000),
     28        elseif(L==40000.),
    2929                md.diagnostic.spcvx(pos)=17.23;
    3030                md.diagnostic.spcvy(pos)=-3.17;
    31         elseif(L==80000),
     31        elseif(L==80000.),
    3232                md.diagnostic.spcvx(pos)=16.68;
    3333                md.diagnostic.spcvy(pos)=-2.69;
    34         elseif(L==160000),
     34        elseif(L==160000.),
    3535                md.diagnostic.spcvx(pos)=16.03;
    3636                md.diagnostic.spcvy(pos)=-1.27;
    3737        end
    3838
    39         md=setflowequation(md,'stokes','all'); 
     39        md=setflowequation(md,'stokes','all');
    4040
    4141        %Compute the diagnostic

  • issm/trunk-jpl/test/NightlyRun/test1107.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%Pattyn and Payne 2006
    33printingflag=false;
    44
    5 L_list={5000,10000,20000,40000,80000,160000};
     5L_list={5000.,10000.,20000.,40000.,80000.,160000.};
    66results={};
    77minvx=[];
    … …  
    2626
    2727        %Create MPCs to have periodic boundary conditions
    28         posx=find(md.mesh.x==0 & ~(md.mesh.y==0 & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
    29         posx2=find(md.mesh.x==max(md.mesh.x) & ~(md.mesh.y==0 & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
     28        posx=find(md.mesh.x==0. & ~(md.mesh.y==0. & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
     29        posx2=find(md.mesh.x==max(md.mesh.x) & ~(md.mesh.y==0. & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
    3030
    31         posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
    32         posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
     31        posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
     32        posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
    3333
    3434        md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
    3535
    3636        %Add spc on the corners
    37         pos=find((md.mesh.x==0 | md.mesh.x==L) & (md.mesh.y==0 | md.mesh.y==L) & md.mesh.vertexonbed);
    38         md.diagnostic.spcvy(:)=0;
    39         md.diagnostic.spcvx(pos)=0;
    40         if(L==5000),
     37        pos=find((md.mesh.x==0. | md.mesh.x==L) & (md.mesh.y==0. | md.mesh.y==L) & md.mesh.vertexonbed);
     38        md.diagnostic.spcvy(:)=0.;
     39        md.diagnostic.spcvx(pos)=0.;
     40        if(L==5000.),
    4141                md.diagnostic.spcvx(pos)=16.0912;
    42         elseif(L==10000),
     42        elseif(L==10000.),
    4343                md.diagnostic.spcvx(pos)=16.52;
    44         elseif(L==20000),
     44        elseif(L==20000.),
    4545                md.diagnostic.spcvx(pos)=17.77;
    46         elseif(L==40000),
     46        elseif(L==40000.),
    4747                md.diagnostic.spcvx(pos)=19.88;
    48         elseif(L==80000),
     48        elseif(L==80000.),
    4949                md.diagnostic.spcvx(pos)=18.65;
    50         elseif(L==160000),
     50        elseif(L==160000.),
    5151                md.diagnostic.spcvx(pos)=16.91;
    5252        end
    … …  
    5454        %Spc the bed at zero for vz
    5555        pos=find(md.mesh.vertexonbed);
    56         md.diagnostic.spcvz(pos)=0;
     56        md.diagnostic.spcvz(pos)=0.;
    5757
    5858        %Compute the diagnostic
    … …  
    7070        %Now plot vx, vy, vz and vx on a cross section
    7171        plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
    72         if printingflag, 
     72        if printingflag,
    7373                set(gcf,'Color','w')
    7474                printmodel(['ismipdpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    75                 system(['mv ismipdpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
     75                system(['mv ismipdpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
    7676        end
    7777        plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
    78         if printingflag, 
     78        if printingflag,
    7979                set(gcf,'Color','w')
    8080                printmodel(['ismipdpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    81                 system(['mv ismipdpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
     81                system(['mv ismipdpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
    8282        end
    8383
    84         if(L==5000),
     84        if(L==5000.),
    8585                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
    8686                        'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',4)
    87         elseif(L==10000),
     87        elseif(L==10000.),
    8888                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
    8989                        'resolution',[10 10],'ylim',[0 20],'xlim',[0 10000],'title','','xlabel','','figure',4)
    90         elseif(L==20000),
     90        elseif(L==20000.),
    9191                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
    9292                        'resolution',[10 10],'ylim',[0 30],'xlim',[0 20000],'title','','xlabel','','figure',4)
    93         elseif(L==40000),
     93        elseif(L==40000.),
    9494                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
    9595                        'resolution',[10 10],'ylim',[10 60],'xlim',[0 40000],'title','','xlabel','','figure',4)
    96         elseif(L==80000),
     96        elseif(L==80000.),
    9797                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
    9898                        'resolution',[10 10],'ylim',[0 200],'xlim',[0 80000],'title','','xlabel','','figure',4)
    99         elseif(L==160000),
     99        elseif(L==160000.),
    100100                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
    101101                        'resolution',[10 10],'ylim',[0 400],'xlim',[0 160000],'title','','xlabel','','figure',4)
    102102        end
    103         if printingflag, 
     103        if printingflag,
    104104                set(gcf,'Color','w')
    105105                printmodel(['ismipdpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    106                 system(['mv ismipdpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
     106                system(['mv ismipdpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
    107107        end
    108108end
    … …  
    110110%Now plot the min and max values of vx for each size of the square
    111111plot([5 10 20 40 80 160],minvx);ylim([2 18]);xlim([0 160])
    112 if printingflag, 
     112if printingflag,
    113113        set(gcf,'Color','w')
    114114        printmodel('ismipdpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    115         system(['mv ismipdpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
     115        system(['mv ismipdpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
    116116end
    117117plot([5 10 20 40 80 160],maxvx);ylim([0 300]);xlim([0 160])
    118 if printingflag, 
     118if printingflag,
    119119        set(gcf,'Color','w')
    120120        printmodel('ismipdpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    121         system(['mv ismipdpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
     121        system(['mv ismipdpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
    122122end
     123
    123124%Fields and tolerances to track changes
    124125field_names     ={...

  • issm/trunk-jpl/test/NightlyRun/test1108.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%Pattyn and Payne 2006
    33
    4 L_list={5000,10000,20000,40000,80000,160000};
     4L_list={5000.,10000.,20000.,40000.,80000.,160000.};
    55results={};
    66
    … …  
    1717        md=setflowequation(md,'pattyn','all');
    1818
    19         %We need one grd on dirichlet: the 4 corners are set to zero
     19        %We need one grid on dirichlet: the 4 corners are set to zero
    2020        md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
    2121        md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
    2222        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
    2323       
    24         pos=find(md.mesh.vertexonbed & (md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
    25         md.diagnostic.spcvx(pos)=0;
    26         md.diagnostic.spcvy(pos)=0;
    27         md.diagnostic.spcvz(pos)=0;
     24        pos=find(md.mesh.vertexonbed & (md.mesh.x==0. | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0. | md.mesh.y==max(md.mesh.y)));
     25        md.diagnostic.spcvx(pos)=0.;
     26        md.diagnostic.spcvy(pos)=0.;
     27        md.diagnostic.spcvz(pos)=0.;
    2828
    2929        %Create MPCs to have periodic boundary conditions
    30         posx=find(md.mesh.x==0);
     30        posx=find(md.mesh.x==0.);
    3131        posx2=find(md.mesh.x==max(md.mesh.x));
    3232
    33         posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
    34         posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
     33        posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
     34        posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
    3535
    3636        md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
    … …  
    4747        md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
    4848        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
    49         pos=find(md.mesh.y==0 | md.mesh.x==0 | md.mesh.x==max(md.mesh.x) | md.mesh.y==max(md.mesh.y)); %Don't take the same nodes two times
     49        pos=find(md.mesh.y==0. | md.mesh.x==0. | md.mesh.x==max(md.mesh.x) | md.mesh.y==max(md.mesh.y)); %Don't take the same nodes two times
    5050        md.diagnostic.spcvx(pos)=md.results.DiagnosticSolution.Vx(pos);
    5151        md.diagnostic.spcvy(pos)=md.results.DiagnosticSolution.Vy(pos);

  • issm/trunk-jpl/test/NightlyRun/test1109.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%TestE
    33%Four tests to run: - Pattyn frozen
    … …  
    99
    1010for i=1:4,
    11         Lx=10; %in m
    12         Ly=5000; %in m
     11        Lx=10.; %in m
     12        Ly=5000.; %in m
    1313        nx=3; %number of nodes in x direction
    1414        ny=51;
    … …  
    2626
    2727        %Create MPCs to have periodic boundary conditions
    28         posx=find(md.mesh.x==0);
     28        posx=find(md.mesh.x==0.);
    2929        posx2=find(md.mesh.x==max(md.mesh.x));
    3030        md.diagnostic.vertex_pairing=[posx,posx2];
    … …  
    3535        md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
    3636        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
    37         md.diagnostic.spcvx(pos)=0;
    38         md.diagnostic.spcvy(pos)=0;
    39         md.diagnostic.spcvz(pos)=0;
     37        md.diagnostic.spcvx(pos)=0.;
     38        md.diagnostic.spcvy(pos)=0.;
     39        md.diagnostic.spcvz(pos)=0.;
    4040
    4141        %Remove the spc where there is some sliding (case 3 and 4):
    … …  
    5858        if i==1,
    5959                plotmodel(md,'data',vy,'ylim',[-10 80],'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIPE.exp','resolution',[10 10],'title','','xlabel','')
    60                 if printingflag, 
     60                if printingflag,
    6161                        set(gcf,'Color','w')
    6262                        printmodel('ismipepattynvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    63                         system(['mv ismipepattynvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE ']);
     63                        system(['mv ismipepattynvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE']);
    6464                end
    6565        elseif i==2,
    6666                plotmodel(md,'data',vy,'ylim',[-10 80],'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIPE.exp','resolution',[10 10],'title','','xlabel','')
    67                 if printingflag, 
     67                if printingflag,
    6868                        set(gcf,'Color','w')
    6969                        printmodel('ismipestokesvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    70                         system(['mv ismipestokesvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE ']);
     70                        system(['mv ismipestokesvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE']);
    7171                end
    7272        elseif i==3,
    7373                plotmodel(md,'data',vy,'ylim',[-50 200],'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIPE.exp','resolution',[10 10],'title','','xlabel','')
    74                 if printingflag, 
     74                if printingflag,
    7575                        set(gcf,'Color','w')
    7676                        printmodel('ismipepattynvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    77                         system(['mv ismipepattynvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE ']);
     77                        system(['mv ismipepattynvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE']);
    7878                end
    7979        elseif i==4,
    8080                plotmodel(md,'data',vy,'ylim',[-50 200],'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIPE.exp','resolution',[10 10],'title','','xlabel','')
    81                 if printingflag, 
     81                if printingflag,
    8282                        set(gcf,'Color','w')
    8383                        printmodel('ismipestokesvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    84                         system(['mv ismipestokesvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE ']);
     84                        system(['mv ismipestokesvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE']);
    8585                end
    8686        end
    8787end
     88
    8889%Fields and tolerances to track changes
    89 field_names     ={ ...
     90field_names     ={...
    9091        'VyPattynSliding','VzPattynSliding',...
    9192        'VxStokesSliding','VyStokesSliding','VzStokesSliding',...

  • issm/trunk-jpl/test/NightlyRun/test1110.m

    r13671 r14134  
    1 %This test is a test from the ISMP-HOM Intercomparison project
     1%This test is a test from the ISMP-HOM Intercomparison project.
    22%TestF
    33printingflag=false;
     4results={};
    45
    56for i=1:4,
    6         L=100000; %in m
     7        L=100000.; %in m
    78        nx=30; %numberof nodes in x direction
    89        ny=30;
    910        md=model();
    1011        md=squaremesh(md,L,L,nx,ny);
    11         %md=triangle(md,'../Exp/SquareISMIP.exp',5500.);
     12%       md=triangle(md,'../Exp/SquareISMIP.exp',5500.);
    1213        md=setmask(md,'',''); %ice sheet test
    1314        md=parameterize(md,'../Par/ISMIPF.par');
    … …  
    2627                %Create dirichlet on the bed if no slip
    2728                pos=find(md.mesh.vertexonbed);
    28                 md.diagnostic.spcvx(pos)=0;
    29                 md.diagnostic.spcvy(pos)=0;
    30                 md.diagnostic.spcvz(pos)=0;
     29                md.diagnostic.spcvx(pos)=0.;
     30                md.diagnostic.spcvy(pos)=0.;
     31                md.diagnostic.spcvz(pos)=0.;
    3132        else
    32                 pos=find(md.mesh.vertexonbed & (md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
    33                 md.diagnostic.spcvx(pos)=100; %because we need a dirichlet somewhere
    34                 md.diagnostic.spcvy(pos)=0;
    35                 md.diagnostic.spcvz(pos)=0;
     33                pos=find(md.mesh.vertexonbed & (md.mesh.x==0. | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0. | md.mesh.y==max(md.mesh.y)));
     34                md.diagnostic.spcvx(pos)=100.; %because we need a dirichlet somewhere
     35                md.diagnostic.spcvy(pos)=0.;
     36                md.diagnostic.spcvz(pos)=0.;
    3637        end
    3738        pos=find(~md.mesh.vertexonbed);
    38         md.thermal.spctemperature(pos)=255;
     39        md.thermal.spctemperature(pos)=255.;
    3940
    4041        %Create MPCs to have periodic boundary conditions
    41         posx=find(md.mesh.x==0);
     42        posx=find(md.mesh.x==0.);
    4243        posx2=find(md.mesh.x==max(md.mesh.x));
    4344
    44         posy=find(md.mesh.y==0);
     45        posy=find(md.mesh.y==0.);
    4546        posy2=find(md.mesh.y==max(md.mesh.y));
    4647
    … …  
    4849        md.prognostic.vertex_pairing=[posx,posx2;posy,posy2];
    4950
    50         md.timestepping.time_step=3;
    51         md.timestepping.final_time=300;
     51        md.timestepping.time_step=3.;
     52        md.timestepping.final_time=300.;
    5253        md.settings.output_frequency=50;
    5354        md.prognostic.stabilization=1;
    … …  
    6869                plotmodel(md,'data',(md.results.TransientSolution(end).Vx),'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIP100000.exp','title','','xlabel','','ylabel','Velocity (m/yr)','linewidth',3,'grid','on','unit','km','ylim',[185 200])
    6970        end
    70         if printingflag, 
     71        if printingflag,
    7172                set(gcf,'Color','w')
    7273                if i==1,
    7374                        printmodel('ismipfpattynvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    74                         system(['mv ismipfpattynvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
     75                        system(['mv ismipfpattynvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
    7576                elseif i==2,
    7677                        printmodel('ismipfpattynvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    77                         system(['mv ismipfpattynvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
     78                        system(['mv ismipfpattynvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
    7879                elseif i==3,
    7980                        printmodel('ismipfstokesvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    80                         system(['mv ismipfstokesvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
     81                        system(['mv ismipfstokesvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
    8182                elseif i==4,
    8283                        printmodel('ismipfstokesvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    83                         system(['mv ismipfstokesvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
     84                        system(['mv ismipfstokesvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
    8485                end
    8586        end
    8687
    8788        plotmodel(md,'data',(md.results.TransientSolution(end).Surface)-md.geometry.surface,'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIP100000.exp','title','','xlabel','','ylabel','Surface (m)','linewidth',3,'grid','on','unit','km','ylim',[-30 50])
    88         if printingflag, 
     89        if printingflag,
    8990                set(gcf,'Color','w')
    9091                if i==1,
    9192                        printmodel('ismipfpattyndeltasurfacefrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    92                         system(['mv ismipfpattyndeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
     93                        system(['mv ismipfpattyndeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
    9394                elseif i==2,
    9495                        printmodel('ismipfpattyndeltasurfacesliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    95                         system(['mv ismipfpattyndeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
     96                        system(['mv ismipfpattyndeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
    9697                elseif i==3,
    9798                        printmodel('ismipfstokesdeltasurfacefrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    98                         system(['mv ismipfstokesdeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
     99                        system(['mv ismipfstokesdeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
    99100                elseif i==4,
    100101                        printmodel('ismipfstokesdeltasurfacesliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
    101                         system(['mv ismipfstokesdeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
     102                        system(['mv ismipfstokesdeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
    102103                end
    103104        end
    … …  
    105106
    106107%Fields and tolerances to track changes
    107 field_names     ={ ...
     108field_names     ={...
    108109        'VxPattynFrozen','VyPattynFrozen','VzPattynFrozen','SurfacePattynFrozen',...
    109110        'VxPattynSliding','VyPattynSliding','VzPattynSliding','SurfacePattynSliding',...

  • issm/trunk-jpl/test/NightlyRun/test1201.m

    r14107 r14134  
    1 %This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996
     1%This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996.
    22printingflag=false;
    33

  • issm/trunk-jpl/test/NightlyRun/test1201.py

    r14107 r14134  
    1010
    1111"""
    12 This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996
     12This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996.
    1313"""
    1414

  • issm/trunk-jpl/test/NightlyRun/test1202.m

    r14107 r14134  
    1 %Test on the diagnostic model and the prognostic in 2d
     1%Test on the diagnostic model and the prognostic in 2d.
    22printingflag=false;
    33

  • issm/trunk-jpl/test/NightlyRun/test1205.m

    r13671 r14134  
    1 %The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
     1%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling.
    22printingflag=false;
    33
    44numlayers=10;
    5 resolution=30000;
     5resolution=30000.;
    66
    77%To begin with the numerical model
    88md=model();
    9 md=roundmesh(md,750000,resolution);
     9md=roundmesh(md,750000.,resolution);
    1010md=setmask(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
    1111md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
    … …  
    1313%Calculation of the analytical 2d velocity field
    1414constant=0.3;
    15 vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1;
    16 vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1;
    17 vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2));
     15vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1;
     16vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1;
     17vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2);
    1818
    1919%We extrude the model to have a 3d model
    … …  
    2323%Spc the nodes on the bed
    2424pos=find(md.mesh.vertexonbed);
    25 md.diagnostic.spcvx(pos)=0;
    26 md.diagnostic.spcvy(pos)=0;
    27 md.diagnostic.spcvz(pos)=0;
     25md.diagnostic.spcvx(pos)=0.;
     26md.diagnostic.spcvy(pos)=0.;
     27md.diagnostic.spcvz(pos)=0.;
    2828
    2929%Now we can solve the problem
    … …  
    3636vel=zeros(md.mesh.numberofvertices2d,1);
    3737
    38 node_vel=0;
    3938for i=1:md.mesh.numberofvertices2d
    40         for j=1:(md.mesh.numberoflayers-1)
    41                 node_vel=node_vel+1/(2*(md.mesh.numberoflayers-1))*(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+...
    42                         vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
     39        node_vel=0.;
     40        for j=1:md.mesh.numberoflayers-1
     41                node_vel=node_vel+1./(2.*(md.mesh.numberoflayers-1))*...
     42                        (sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
    4343                        sqrt(vx(i+(j-1)*md.mesh.numberofvertices2d,1).^2+vy(i+(j-1)*md.mesh.numberofvertices2d,1).^2));
    4444        end
    4545        vel(i,1)=node_vel;
    46         node_vel=0;
    4746end
    4847
    … …  
    8281caxis([0 100]);
    8382
    84 if printingflag, 
     83if printingflag,
    8584        set(gcf,'Color','w')
    8685        printmodel('hutterstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
    … …  
    9291        'Vx','Vy','Vel', ...
    9392};
    94 field_tolerances={...
     93field_tolerances={ ...
    9594        1e-13,1e-13,1e-13, ...
    9695};
    97 field_values={
     96field_values={ ...
    9897        vx,vy,vel, ...
    9998};

  • issm/trunk-jpl/test/NightlyRun/test1206.m

    r13671 r14134  
    1 %The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
     1%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling.
    22printingflag=false;
    33
    44numlayers=10;
    5 resolution=30000;
     5resolution=30000.;
    66
    77%To begin with the numerical model
    88md=model();
    9 md=roundmesh(md,750000,resolution);
     9md=roundmesh(md,750000.,resolution);
    1010md=setmask(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
    1111md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
    … …  
    1313%Calculation of the analytical 2d velocity field
    1414constant=0.3;
    15 vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1;
    16 vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1;
    17 vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2));
     15vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1;
     16vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1;
     17vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2);
    1818
    1919%We extrude the model to have a 3d model
    … …  
    2323%Spc the nodes on the bed
    2424pos=find(md.mesh.vertexonbed);
    25 md.diagnostic.spcvx(pos)=0;
    26 md.diagnostic.spcvy(pos)=0;
    27 md.diagnostic.spcvz(pos)=0;
     25md.diagnostic.spcvx(pos)=0.;
     26md.diagnostic.spcvy(pos)=0.;
     27md.diagnostic.spcvz(pos)=0.;
    2828
    2929%Now we can solve the problem
    … …  
    3636vel=zeros(md.mesh.numberofvertices2d,1);
    3737
    38 node_vel=0;
    3938for i=1:md.mesh.numberofvertices2d
    40         for j=1:(md.mesh.numberoflayers-1)
    41                 node_vel=node_vel+1/(2*(md.mesh.numberoflayers-1))*(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+...
    42                         vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
     39        node_vel=0.;
     40        for j=1:md.mesh.numberoflayers-1
     41                node_vel=node_vel+1./(2.*(md.mesh.numberoflayers-1))*...
     42                        (sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
    4343                        sqrt(vx(i+(j-1)*md.mesh.numberofvertices2d,1).^2+vy(i+(j-1)*md.mesh.numberofvertices2d,1).^2));
    4444        end
    4545        vel(i,1)=node_vel;
    46         node_vel=0;
    4746end
    4847
    … …  
    8180caxis([0 100]);
    8281
    83 if printingflag, 
     82if printingflag,
    8483        set(gcf,'Color','w')
    8584        printmodel('pattynstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
    … …  
    9190        'Vx','Vy','Vel', ...
    9291};
    93 field_tolerances={...
     92field_tolerances={ ...
    9493        1e-12,1e-12,1e-12, ...
    9594};
    96 field_values={
     95field_values={ ...
    9796        vx,vy,vel, ...
    9897};

  • issm/trunk-jpl/test/NightlyRun/test1207.m

    r13671 r14134  
    1 %The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
     1%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling.
    22printingflag=false;
    33
    44numlayers=10;
    5 resolution=30000;
     5resolution=30000.;
    66
    77%To begin with the numerical model
    88md=model();
    9 md=roundmesh(md,750000,resolution);
     9md=roundmesh(md,750000.,resolution);
    1010md=setmask(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
    1111md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
    … …  
    1313%Calculation of the analytical 2d velocity field
    1414constant=0.3;
    15 vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1;
    16 vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1;
    17 vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2));
     15vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1;
     16vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1;
     17vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2);
    1818
    1919%We extrude the model to have a 3d model
    … …  
    2323%Spc the nodes on the bed
    2424pos=find(md.mesh.vertexonbed);
    25 md.diagnostic.spcvx(pos)=0;
    26 md.diagnostic.spcvy(pos)=0;
    27 md.diagnostic.spcvz(pos)=0;
     25md.diagnostic.spcvx(pos)=0.;
     26md.diagnostic.spcvy(pos)=0.;
     27md.diagnostic.spcvz(pos)=0.;
    2828
    2929%Now we can solve the problem
    … …  
    3636vel=zeros(md.mesh.numberofvertices2d,1);
    3737
    38 node_vel=0;
    3938for i=1:md.mesh.numberofvertices2d
    40         for j=1:(md.mesh.numberoflayers-1)
    41                 node_vel=node_vel+1/(2*(md.mesh.numberoflayers-1))*(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+...
    42                         vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
     39        node_vel=0.;
     40        for j=1:md.mesh.numberoflayers-1
     41                node_vel=node_vel+1./(2.*(md.mesh.numberoflayers-1))*...
     42                        (sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
    4343                        sqrt(vx(i+(j-1)*md.mesh.numberofvertices2d,1).^2+vy(i+(j-1)*md.mesh.numberofvertices2d,1).^2));
    4444        end
    4545        vel(i,1)=node_vel;
    46         node_vel=0;
    4746end
    4847
    … …  
    8180caxis([0 100]);
    8281
    83 if printingflag, 
     82if printingflag,
    8483        set(gcf,'Color','w')
    8584        printmodel('stokesstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
    … …  
    9190        'Vx','Vy','Vel', ...
    9291};
    93 field_tolerances={...
     92field_tolerances={ ...
    9493        1e-12,1e-12,1e-12, ...
    9594};
    96 field_values={
     95field_values={ ...
    9796        vx,vy,vel, ...
    9897};

  • issm/trunk-jpl/test/Par/ISMIPA.par

    r13137 r14134  
    22
    33disp('      creating thickness');
    4 md.geometry.surface=-md.mesh.x*tan(0.5*pi/180);
    5 md.geometry.bed=md.geometry.surface-1000+500*sin(md.mesh.x*2*pi/max(md.mesh.x)).*sin(md.mesh.y*2*pi/max(md.mesh.x));
     4md.geometry.surface=-md.mesh.x*tan(0.5*pi/180.);
     5md.geometry.bed=md.geometry.surface-1000.+500.*sin(md.mesh.x*2.*pi/max(md.mesh.x)).*sin(md.mesh.y*2.*pi/max(md.mesh.x));
    66md.geometry.thickness=md.geometry.surface-md.geometry.bed;
    77
    88disp('      creating drag');
    9 md.friction.coefficient=200*ones(md.mesh.numberofvertices,1); %q=1.
     9md.friction.coefficient=200.*ones(md.mesh.numberofvertices,1); %q=1.
    1010%Take care of iceshelves: no basal drag
    1111pos=find(md.mask.elementonfloatingice);
    12 md.friction.coefficient(md.mesh.elements(pos,:))=0;
     12md.friction.coefficient(md.mesh.elements(pos,:))=0.;
    1313md.friction.p=ones(md.mesh.numberofelements,1);
    1414md.friction.q=ones(md.mesh.numberofelements,1);
    1515
    16 disp('      creating flow law paramter');
     16disp('      creating flow law parameter');
    1717md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    18 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     18md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
    1919
    2020disp('      boundary conditions for diagnostic model');
    21 %Create node on boundary fist (because we cannot use mesh)
     21%Create node on boundary first (because we cannot use mesh)
    2222md=SetIceSheetBC(md);

  • issm/trunk-jpl/test/Par/ISMIPB.par

    r13137 r14134  
    22
    33disp('      creating thickness');
    4 md.geometry.surface=-md.mesh.x*tan(0.5*pi/180);
    5 md.geometry.bed=md.geometry.surface-1000+500*sin(md.mesh.x*2*pi/max(md.mesh.x));
     4md.geometry.surface=-md.mesh.x*tan(0.5*pi/180.);
     5md.geometry.bed=md.geometry.surface-1000.+500.*sin(md.mesh.x*2.*pi/max(md.mesh.x));
    66md.geometry.thickness=md.geometry.surface-md.geometry.bed;
    77
    88disp('      creating drag');
    9 md.friction.coefficient=200*ones(md.mesh.numberofvertices,1); %q=1.
     9md.friction.coefficient=200.*ones(md.mesh.numberofvertices,1); %q=1.
    1010%Take care of iceshelves: no basal drag
    1111pos=find(md.mask.elementonfloatingice);
    12 md.friction.coefficient(md.mesh.elements(pos,:))=0;
     12md.friction.coefficient(md.mesh.elements(pos,:))=0.;
    1313md.friction.p=ones(md.mesh.numberofelements,1);
    1414md.friction.q=ones(md.mesh.numberofelements,1);
    1515
    16 disp('      creating flow law paramter');
     16disp('      creating flow law parameter');
    1717md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    18 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     18md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
    1919
    2020disp('      boundary conditions for diagnostic model');
    21 %Create node on boundary fist (because we cannot use mesh)
     21%Create node on boundary first (because we cannot use mesh)
    2222md=SetIceSheetBC(md);

  • issm/trunk-jpl/test/Par/ISMIPC.par

    r13137 r14134  
    22
    33disp('      creating thickness');
    4 md.geometry.surface=2000-md.mesh.x*tan(0.1*pi/180); %to have z>0
    5 md.geometry.bed=md.geometry.surface-1000;
     4md.geometry.surface=2000.-md.mesh.x*tan(0.1*pi/180.); %to have z>0
     5md.geometry.bed=md.geometry.surface-1000.;
    66md.geometry.thickness=md.geometry.surface-md.geometry.bed;
    77
    88disp('      creating drag');
    9 %md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x/2)).*sin(md.mesh.y*2*pi/max(md.mesh.x/2)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)));
    10 md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x)).*sin(md.mesh.y*2*pi/max(md.mesh.x))));
     9%md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x/2.)).*sin(md.mesh.y*2.*pi/max(md.mesh.x/2.)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)));
     10md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x)).*sin(md.mesh.y*2.*pi/max(md.mesh.x))));
    1111%Take care of iceshelves: no basal drag
    1212pos=find(md.mask.elementonfloatingice);
    13 md.friction.coefficient(md.mesh.elements(pos,:))=0;
     13md.friction.coefficient(md.mesh.elements(pos,:))=0.;
    1414md.friction.p=ones(md.mesh.numberofelements,1);
    1515md.friction.q=zeros(md.mesh.numberofelements,1);
    1616
    17 disp('      creating flow law paramter');
     17disp('      creating flow law parameter');
    1818md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    19 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     19md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
    2020
    21 disp('      boundary conditions for diagnostic model: ');
    22 %Create node on boundary fist (because wi can not use mesh)
     21disp('      boundary conditions for diagnostic model:');
     22%Create node on boundary first (because we can not use mesh)
    2323md=SetIceSheetBC(md);

  • issm/trunk-jpl/test/Par/ISMIPD.par

    r13137 r14134  
    22
    33disp('      creating thickness');
    4 md.geometry.surface=2000-md.mesh.x*tan(0.1*pi/180); %to have z>0
    5 md.geometry.bed=md.geometry.surface-1000;
     4md.geometry.surface=2000.-md.mesh.x*tan(0.1*pi/180.); %to have z>0
     5md.geometry.bed=md.geometry.surface-1000.;
    66md.geometry.thickness=md.geometry.surface-md.geometry.bed;
    77
    88disp('      creating drag');
    9 md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x))));
     9md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x))));
    1010%Take care of iceshelves: no basal drag
    1111pos=find(md.mask.elementonfloatingice);
    12 md.friction.coefficient(md.mesh.elements(pos,:))=0;
     12md.friction.coefficient(md.mesh.elements(pos,:))=0.;
    1313md.friction.p=ones(md.mesh.numberofelements,1);
    1414md.friction.q=zeros(md.mesh.numberofelements,1);
    1515
    16 disp('      creating flow law paramter');
     16disp('      creating flow law parameter');
    1717md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    18 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     18md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
    1919
    20 disp('      boundary conditions for diagnostic model: ');
    21 %Create node on boundary fist (because wi can not use mesh)
     20disp('      boundary conditions for diagnostic model:');
     21%Create node on boundary first (because we can not use mesh)
    2222md=SetIceSheetBC(md);

  • issm/trunk-jpl/test/Par/ISMIPE.par

    r13233 r14134  
    77for i=1:md.mesh.numberofvertices
    88        y=md.mesh.y(i);
    9         point1=floor(y/100)+1;
    10         point2=min(point1+1,51);
    11         coeff=(y-(point1-1)*100)/100;
    12         md.geometry.bed(i)=(1-coeff)*data(point1,2)+coeff*data(point2,2);
    13         md.geometry.surface(i)=(1-coeff)*data(point1,3)+coeff*data(point2,3);
     9        point1=floor(y/100.)+1.;
     10        point2=min(point1+1,51.);
     11        coeff=(y-(point1-1.)*100.)/100.;
     12        md.geometry.bed(i)=(1.-coeff)*data(point1,2)+coeff*data(point2,2);
     13        md.geometry.surface(i)=(1.-coeff)*data(point1,3)+coeff*data(point2,3);
    1414end
    1515md.geometry.thickness=md.geometry.surface-md.geometry.bed;
    … …  
    2222md.friction.q=ones(md.mesh.numberofelements,1);
    2323
    24 disp('      creating flow law paramter');
     24disp('      creating flow law parameter');
    2525md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    26 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     26md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
    2727
    28 disp('      boundary conditions for diagnostic model: ');
    29 %Create node on boundary fist (because wi can not use mesh)
     28disp('      boundary conditions for diagnostic model:');
     29%Create node on boundary first (because we can not use mesh)
    3030md=SetIceSheetBC(md);

  • issm/trunk-jpl/test/Par/ISMIPF.par

    r13137 r14134  
    33
    44disp('      creating thickness');
    5 md.geometry.surface=-md.mesh.x*tan(3*pi/180);
    6 %md.geometry.bed=md.geometry.surface-1000;
    7 md.geometry.bed=md.geometry.surface-1000+100*exp(-((md.mesh.x-max(md.mesh.x)/2).^2+(md.mesh.y-max(md.mesh.y)/2).^2)/(10000^2));
     5md.geometry.surface=-md.mesh.x*tan(3.*pi/180.);
     6%md.geometry.bed=md.geometry.surface-1000.;
     7md.geometry.bed=md.geometry.surface-1000.+100.*exp(-((md.mesh.x-max(md.mesh.x)/2.).^2+(md.mesh.y-max(md.mesh.y)/2.).^2)/(10000.^2));
    88md.geometry.thickness=md.geometry.surface-md.geometry.bed;
    99
    1010disp('      creating drag');
    11 md.friction.coefficient=sqrt(md.constants.yts/(2.140373*10^-7*1000))*ones(md.mesh.numberofvertices,1);
     11md.friction.coefficient=sqrt(md.constants.yts/(2.140373*10^-7*1000.))*ones(md.mesh.numberofvertices,1);
    1212md.friction.p=ones(md.mesh.numberofelements,1);
    1313md.friction.q=zeros(md.mesh.numberofelements,1);
    1414
    15 disp('      creating flow law paramter');
     15disp('      creating flow law parameter');
    1616md.materials.rheology_B=1.4734*10^14*ones(md.mesh.numberofvertices,1);
    17 md.materials.rheology_n=1*ones(md.mesh.numberofelements,1);
     17md.materials.rheology_n=1.*ones(md.mesh.numberofelements,1);
    1818md.materials.rheology_law='None';
    1919
    2020disp('      boundary conditions for diagnostic model');
    21 %Create node on boundary fist (because we cannot use mesh)
     21%Create node on boundary first (because we cannot use mesh)
    2222md=SetIceSheetBC(md);
    23 md.diagnostic.spcvx=100*ones(md.mesh.numberofvertices,1);
     23md.diagnostic.spcvx=100.*ones(md.mesh.numberofvertices,1);
    2424md.initialization.vx=zeros(md.mesh.numberofvertices,1);
    2525md.initialization.vy=zeros(md.mesh.numberofvertices,1);
    … …  
    2727md.initialization.vel=zeros(md.mesh.numberofvertices,1);
    2828md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
    29 md.initialization.temperature=255*ones(md.mesh.numberofvertices,1);
     29md.initialization.temperature=255.*ones(md.mesh.numberofvertices,1);
    3030pos=find(md.mesh.x==min(md.mesh.x) | md.mesh.x==max(md.mesh.x) | md.mesh.y==min(md.mesh.y) | md.mesh.y==max(md.mesh.y));
    3131md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
    … …  
    3333md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
    3434md.prognostic.spcthickness(pos)=md.geometry.thickness(pos);
    35 md.thermal.spctemperature=255*ones(md.mesh.numberofvertices,1);
     35md.thermal.spctemperature=255.*ones(md.mesh.numberofvertices,1);
    3636md.basalforcings.geothermalflux=0.4*ones(md.mesh.numberofvertices,1);
    3737
    … …  
    4040
    4141%Transient options
    42 md.timestepping.time_step=1;
    43 md.timestepping.final_time=10;
     42md.timestepping.time_step=1.;
     43md.timestepping.final_time=10.;
    4444md.prognostic.stabilization=1;
    4545md.thermal.stabilization=1;
Note: See TracChangeset for help on using the changeset viewer.