Changeset 5539


Ignore:
Timestamp:
08/24/10 11:30:11 (15 years ago)
Author:
seroussi
Message:

test1107 : ismipd pattyn

Location:
issm/trunk/test/NightlyRun
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/test/NightlyRun/test1103.m

    r5445 r5539  
    3939        vy=PatchToVec(md.results.DiagnosticSolution.Vy);
    4040        vz=PatchToVec(md.results.DiagnosticSolution.Vz);
    41         results{i}=md.results.DiagnosticSolution;
    4241        results{i}=md.results.DiagnosticSolution;
    4342        minvx(i)=min(vx(end-md.numberofgrids2d+1:end));
  • issm/trunk/test/NightlyRun/test1105.m

    r5516 r5539  
    1717        md=parameterize(md,'../Par/ISMIPC.par');
    1818        md=extrude(md,10,1);
    19         md.verbose=4;
    2019
    2120        md=setelementstype(md,'pattyn','all');
  • issm/trunk/test/NightlyRun/test1107.m

    r5112 r5539  
    11%This test is a test from the ISMP-HOM Intercomparison project
    22%Pattyn and Payne 2006
     3printingflag=true;
    34
    45L_list={5000,10000,20000,40000,80000,160000};
    56results={};
     7minvx=[];
     8maxvx=[];
    69
    710for i=1:length(L_list),
    8         L=3*L_list{i};
    9         nx=60; %numberof nodes in x direction
    10         ny=60;
     11        L=L_list{i};
     12        nx=30; %numberof nodes in x direction
     13        ny=30;
    1114        md=model;
    1215        md=squaremesh(md,L,L,nx,ny);
    1316        md=geography(md,'',''); %ice sheet test
    1417        md=parameterize(md,'../Par/ISMIPD.par');
    15         md=extrude(md,6,1);
     18        md=extrude(md,10,1);
    1619
    1720        md=setelementstype(md,'pattyn','all');
    1821
    1922        %We need one grd on dirichlet: the 4 corners are set to zero
    20         %md.spcvelocity=zeros(md.numberofgrids,6);
    21         %pos=find((md.x==0 | md.x==max(md.x)) & (md.y==0 | md.y==max(md.y)));
    22         %md.spcvelocity(pos,1:3)=1;
     23        md.spcvelocity=zeros(md.numberofgrids,6);
    2324
    24         %%Create MPCs to have periodic boundary conditions
    25         %posx=find(md.x==0);
    26         %posx2=find(md.x==max(md.x));
     25        %Create MPCs to have periodic boundary conditions
     26        posx=find(md.x==0 & ~(md.y==0 & md.gridonbed) & ~(md.y==L & md.gridonbed));
     27        posx2=find(md.x==max(md.x) & ~(md.y==0 & md.gridonbed) & ~(md.y==L & md.gridonbed));
    2728
    28         %posy=find(md.y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same grids two times
    29         %posy2=find(md.y==max(md.y) & md.x~=0 & md.x~=max(md.x));
     29        posy=find(md.y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same grids two times
     30        posy2=find(md.y==max(md.y) & md.x~=0 & md.x~=max(md.x));
    3031
    31         %md.penalties=[posx,posx2;posy,posy2];
     32        md.penalties=[posx,posx2;posy,posy2];
     33
     34        %Add spc on the corners
     35        pos=find((md.x==0 | md.x==L) & (md.y==0 | md.y==L) & md.gridonbed);
     36        md.spcvelocity(:,2)=1;
     37        md.spcvelocity(pos,1:2)=1;
     38        if(L==5000),
     39                md.spcvelocity(pos,4)=16.0912;
     40        elseif(L==10000),
     41                md.spcvelocity(pos,4)=16.52;
     42        elseif(L==20000),
     43                md.spcvelocity(pos,4)=17.77;
     44        elseif(L==40000),
     45                md.spcvelocity(pos,4)=19.88;
     46        elseif(L==80000),
     47                md.spcvelocity(pos,4)=18.65;
     48        elseif(L==160000),
     49                md.spcvelocity(pos,4)=16.91;
     50        end
     51       
     52        %Spc the bed at zero for vz
     53        pos=find(md.gridonbed);
     54        md.spcvelocity(pos,3)=1;
    3255
    3356        %Compute the diagnostic
     
    3962        vz=PatchToVec(md.results.DiagnosticSolution.Vz);
    4063        results{i}=md.results.DiagnosticSolution;
     64        minvx(i)=min(vx(end-md.numberofgrids2d+1:end));
     65        maxvx(i)=max(vx(end-md.numberofgrids2d+1:end));
    4166
    42         plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.numlayers)
     67        %Now plot vx, vy, vz and vx on a cross section
     68        plotmodel(md,'data',vx,'layer#all',md.numlayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
     69        if printingflag,
     70                set(gcf,'Color','w')
     71                printmodel(['ismipdpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
     72                system(['mv ismipdpattynvx' num2str(L) '.png ' ISSM_DIR '/../website/doc_pdf/validation/Images/ISMIP/TestD ']);
     73        end
     74        plotmodel(md,'data',vz,'layer#all',md.numlayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
     75        if printingflag,
     76                set(gcf,'Color','w')
     77                printmodel(['ismipdpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
     78                system(['mv ismipdpattynvz' num2str(L) '.png ' ISSM_DIR '/../website/doc_pdf/validation/Images/ISMIP/TestD ']);
     79        end
     80
     81        if(L==5000),
     82                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.numlayers,...
     83                        'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',4)
     84        elseif(L==10000),
     85                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.numlayers,...
     86                        'resolution',[10 10],'ylim',[0 20],'xlim',[0 10000],'title','','xlabel','','figure',4)
     87        elseif(L==20000),
     88                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.numlayers,...
     89                        'resolution',[10 10],'ylim',[0 30],'xlim',[0 20000],'title','','xlabel','','figure',4)
     90        elseif(L==40000),
     91                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.numlayers,...
     92                        'resolution',[10 10],'ylim',[10 60],'xlim',[0 40000],'title','','xlabel','','figure',4)
     93        elseif(L==80000),
     94                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.numlayers,...
     95                        'resolution',[10 10],'ylim',[0 200],'xlim',[0 80000],'title','','xlabel','','figure',4)
     96        elseif(L==160000),
     97                plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.numlayers,...
     98                        'resolution',[10 10],'ylim',[0 400],'xlim',[0 160000],'title','','xlabel','','figure',4)
     99        end
     100        if printingflag,
     101                set(gcf,'Color','w')
     102                printmodel(['ismipdpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
     103                system(['mv ismipdpattynvxsec' num2str(L) '.png ' ISSM_DIR '/../website/doc_pdf/validation/Images/ISMIP/TestD ']);
     104        end
    43105end
    44106
     107%Now plot the min and max values of vx for each size of the square
     108plot([5 10 20 40 80 160],minvx);ylim([2 18])
     109if printingflag,
     110        set(gcf,'Color','w')
     111        printmodel('ismipdpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
     112        system(['mv ismipdpattynminvx.png ' ISSM_DIR '/../website/doc_pdf/validation/Images/ISMIP/TestD ']);
     113end
     114plot([5 10 20 40 80 160],maxvx);ylim([0 300])
     115if printingflag,
     116        set(gcf,'Color','w')
     117        printmodel('ismipdpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
     118        system(['mv ismipdpattynmaxvx.png ' ISSM_DIR '/../website/doc_pdf/validation/Images/ISMIP/TestD ']);
     119end
    45120%Fields and tolerances to track changes
    46121field_names     ={...
     
    51126        'Vx80km','Vy80km','Vz80km',...
    52127        'Vx160km','Vy160km','Vz160km'
    53 }
     128};
    54129field_tolerances={...
    55130        1e-13,1e-13,1e-13,...
Note: See TracChangeset for help on using the changeset viewer.