Changeset 4883


Ignore:
Timestamp:
07/29/10 12:42:22 (15 years ago)
Author:
seroussi
Message:

Validation tests ISMIP-HOM

Location:
issm/trunk/test/Validation/ISMIP
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/test/Validation/ISMIP/TestA/Square.par

    r3768 r4883  
    2222disp('      boundary conditions for diagnostic model');
    2323%Create grid on boundary fist (because we cannot use mesh)
    24 md.gridonboundary=zeros(md.numberofgrids,1);
    25 pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
    26 md.gridonboundary(pos)=1;
    27 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    28 pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
    29 md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    30 md.segmentonneumann_diag=zeros(0,3);
    31 md.gridondirichlet_thermal=zeros(md.numberofgrids,1);
     24md=SetIceSheetBC(md);
     25
     26md.np=3;

  • issm/trunk/test/Validation/ISMIP/TestA/runme.m

    r899 r4883  
    1212md=extrude(md,10,1);
    1313
    14 %Create dirichlet on the bed on top of the borders
    15 pos=find(md.gridonbed);
    16 md.gridondirichlet_diag(pos)=1;
    17 
    1814md=setelementstype(md,'pattyn','all');
    1915%md=setelementstype(md,'pattyn','all','stokes','all');
    2016
    2117%Create dirichlet on the bed only
    22 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
     18md.spcvelocity=zeros(md.numberofgrids,6);
    2319pos=find(md.gridonbed);
    24 md.gridondirichlet_diag(pos)=1;
     20md.spcvelocity(pos,1:2)=1;
    2521
    2622%Create MPCs to have periodic boundary conditions
    27 posx=find(md.x==0);
    28 posx2=find(md.x==max(md.x));
    29 
    30 posy=find(md.y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same grids two times
    31 posy2=find(md.y==max(md.y) & md.x~=0 & md.x~=max(md.x));
    32 
    33 md.penalties=[posx,posx2;posy,posy2];
     23%posx=find(md.x==0);
     24%posx2=find(md.x==max(md.x));
     25%
     26%posy=find(md.y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same grids two times
     27%posy2=find(md.y==max(md.y) & md.x~=0 & md.x~=max(md.x));
     28%
     29%md.penalties=[posx,posx2;posy,posy2];
    3430
    3531%Compute the diagnostic
    36 md=solve(md,'analysis_type','diagnostic','package','ice');
     32md=solve(md,'analysis_type',DiagnosticSolutionEnum);
     33vx=zeros(md.numberofgrids,1);
     34vx(md.results.DiagnosticSolution.Vx.index)=md.results.DiagnosticSolution.Vx.value;
     35vy=zeros(md.numberofgrids,1);
     36vy(md.results.DiagnosticSolution.Vy.index)=md.results.DiagnosticSolution.Vy.value;
     37vz=zeros(md.numberofgrids,1);
     38vz(md.results.DiagnosticSolution.Vz.index)=md.results.DiagnosticSolution.Vz.value;
    3739
    38 plotmodel(md,'data',md.results.diagnostic.vx,'data',md.results.diagnostic.vy,'data',md.results.diagnostic.vz,'layer#all',md.numlayers,'contourlevels',{9.2,18.4,27.6,46,73.6,82.8},'contourlevels',{0,-1.2,-0.9,-0.6,-0.3,1.2},'contourlevels',{-5.2,-3.9,-1.3,0,1.3,3.9,5.2})
     40plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.numlayers,'contourlevels',{9.2,18.4,27.6,46,73.6,82.8},'contourlevels',{0,-1.2,-0.9,-0.6,-0.3,1.2},'contourlevels',{-5.2,-3.9,-1.3,0,1.3,3.9,5.2})

  • issm/trunk/test/Validation/ISMIP/TestB/Square.par

    r3768 r4883  
    2020md.rheology_n=3*ones(md.numberofelements,1);
    2121
    22 disp('      boundary conditions for diagnostic model: ');
    23 %Create grid on boundary fist (because wi can not use mesh)
    24 md.gridonboundary=zeros(md.numberofgrids,1);
    25 pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
    26 md.gridonboundary(pos)=1;
    27 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    28 pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
    29 md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    30 md.segmentonneumann_diag=zeros(0,3);
    31 md.gridondirichlet_thermal=zeros(md.numberofgrids,1);
     22disp('      boundary conditions for diagnostic model');
     23%Create grid on boundary fist (because we cannot use mesh)
     24md=SetIceSheetBC(md);
     25
     26md.np=3;

  • issm/trunk/test/Validation/ISMIP/TestB/runme.m

    r899 r4883  
    1212md=extrude(md,10,1);
    1313
    14 %Create dirichlet on the bed on top of the boundary
    15 pos=find(md.gridonbed);
    16 md.gridondirichlet_diag(pos)=1;
    17 
    1814%md=setelementstype(md,'pattyn','all','stokes','all');
    1915md=setelementstype(md,'pattyn','all');
    2016
    2117%Create dirichlet on the bed only
    22 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
     18md.spcvelocity=zeros(md.numberofgrids,6);
    2319pos=find(md.gridonbed);
    24 md.gridondirichlet_diag(pos)=1;
     20md.spcvelocity(pos,1:2)=1;
    2521
    26 %Create MPCs to have periodic boundary conditions
    27 posx=find(md.x==0);
    28 posx2=find(md.x==max(md.x));
    29 
    30 posy=find(md.y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same grids two times
    31 posy2=find(md.y==max(md.y) & md.x~=0 & md.x~=max(md.x));
    32 
    33 md.penalties=[posx,posx2;posy,posy2];
     22%%Create MPCs to have periodic boundary conditions
     23%posx=find(md.x==0);
     24%posx2=find(md.x==max(md.x));
     25%
     26%posy=find(md.y==0 & md.x~=0 & md.x~=max(md.x)); %Don't take the same grids two times
     27%posy2=find(md.y==max(md.y) & md.x~=0 & md.x~=max(md.x));
     28%
     29%md.penalties=[posx,posx2;posy,posy2];
    3430
    3531%Compute the diagnostic
    36 md=solve(md,'analysis_type','diagnostic','package','ice');
     32md=solve(md,'analysis_type',DiagnosticSolutionEnum);
     33vx=zeros(md.numberofgrids,1);
     34vx(md.results.DiagnosticSolution.Vx.index)=md.results.DiagnosticSolution.Vx.value;
     35vy=zeros(md.numberofgrids,1);
     36vy(md.results.DiagnosticSolution.Vy.index)=md.results.DiagnosticSolution.Vy.value;
     37vz=zeros(md.numberofgrids,1);
     38vz(md.results.DiagnosticSolution.Vz.index)=md.results.DiagnosticSolution.Vz.value;
    3739
    38 plotmodel(md,'data',md.results.diagnostic.vx,'data',md.results.diagnostic.vy,'data',md.results.diagnostic.vz,'layer#all',md.numlayers)
     40plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.numlayers,'contourlevels',{9.2,18.4,27.6,46,73.6,82.8},'contourlevels',{0,-1.2,-0.9,-0.6,-0.3,1.2},'contourlevels',{-5.2,-3.9,-1.3,0,1.3,3.9,5.2})

  • issm/trunk/test/Validation/ISMIP/TestC/Square.par

    r3768 r4883  
    2222disp('      boundary conditions for diagnostic model: ');
    2323%Create grid on boundary fist (because wi can not use mesh)
    24 md.gridonboundary=zeros(md.numberofgrids,1);
    25 pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
    26 md.gridonboundary(pos)=1;
    27 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    28 pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
    29 md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    30 md.segmentonneumann_diag=zeros(0,3);
    31 md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
     24md=SetIceSheetBC(md);
     25md.np=3;

  • issm/trunk/test/Validation/ISMIP/TestC/runme.m

    r899 r4883  
    1515%md=setelementstype(md,'pattyn','all','stokes','all');
    1616
    17 %Add dirichlet values on the corners (we need dirichlet somewhere
    18 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    19 pos=find((md.x==0 |md.x==L) & (md.y==0 | md.y==L));
    20 %pos=find((md.x==0 |md.x==L | md.y==0 | md.y==L)); %dirichlet on all the boundaries
    21 md.gridondirichlet_diag(pos)=1;
    22 
    23 %Create MPCs to have periodic boundary conditions
    24 posx=find(md.x==0);
    25 posx2=find(md.x==L);
    26 
    27 posy=find(md.y==0 & md.x~=0 & md.x~=L); %Don't take the same grids two times
    28 posy2=find(md.y==L & md.x~=0 & md.x~=L);
    29 
    30 md.penalties=[posx,posx2;posy,posy2];
     17%%Create MPCs to have periodic boundary conditions
     18%posx=find(md.x==0);
     19%posx2=find(md.x==L);
     20%
     21%posy=find(md.y==0 & md.x~=0 & md.x~=L); %Don't take the same grids two times
     22%posy2=find(md.y==L & md.x~=0 & md.x~=L);
     23%
     24%md.penalties=[posx,posx2;posy,posy2];
    3125
    3226%Compute the diagnostic
    33 md=solve(md,'analysis_type','diagnostic','package','ice');
     27md=solve(md,'analysis_type',DiagnosticSolutionEnum);
     28vx=zeros(md.numberofgrids,1);
     29vx(md.results.DiagnosticSolution.Vx.index)=md.results.DiagnosticSolution.Vx.value;
     30vy=zeros(md.numberofgrids,1);
     31vy(md.results.DiagnosticSolution.Vy.index)=md.results.DiagnosticSolution.Vy.value;
     32vz=zeros(md.numberofgrids,1);
     33vz(md.results.DiagnosticSolution.Vz.index)=md.results.DiagnosticSolution.Vz.value;
    3434
    35 plotmodel(md,'data',md.results.diagnostic.vx,'data',md.results.diagnostic.vy,'data',md.results.diagnostic.vz,'layer#all',md.numlayers,'contourlevels',{10.878,16.317,21.756,43.51,59.829},'contourlevels',{-2.421,-1.614,-0.807,3.228},'contourlevels',{-2.48,-1.86,0,0.62,1.24,1.86,2.48})
     35plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.numlayers,'contourlevels',{10.878,16.317,21.756,43.51,59.829},'contourlevels',{-2.421,-1.614,-0.807,3.228},'contourlevels',{-2.48,-1.86,0,0.62,1.24,1.86,2.48})

  • issm/trunk/test/Validation/ISMIP/TestD/Square.par

    r3768 r4883  
    2222disp('      boundary conditions for diagnostic model: ');
    2323%Create grid on boundary fist (because wi can not use mesh)
    24 md.gridonboundary=zeros(md.numberofgrids,1);
    25 pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
    26 md.gridonboundary(pos)=1;
    27 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    28 pos=find(md.gridonboundary);md.gridondirichlet_diag(pos)=1;
    29 md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    30 md.segmentonneumann_diag=zeros(0,3);
    31 md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
     24md=SetIceSheetBC(md);
     25
     26md.np=3;

  • issm/trunk/test/Validation/ISMIP/TestD/runme.m

    r899 r4883  
    1616%md=setelementstype(md,'pattyn','all','stokes','all');
    1717
    18 %No dirichlet in this model, just periodic boundary conditions
    19 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    20 pos=find((md.x==0 |md.x==80000 ) & (md.y==0 | md.y==80000));
    21 md.gridondirichlet_diag(pos)=1;
    22 
    2318%We need one grd on dirichlet (in the middle)
    2419%md.gridondirichlet_diag(33)=1;
    … …  
    3530
    3631%Compute the diagnostic
    37 md=solve(md,'analysis_type','diagnostic','package','ice');
     32md=solve(md,'analysis_type',DiagnosticSolutionEnum);
     33vx=zeros(md.numberofgrids,1);
     34vx(md.results.DiagnosticSolution.Vx.index)=md.results.DiagnosticSolution.Vx.value;
     35vy=zeros(md.numberofgrids,1);
     36vy(md.results.DiagnosticSolution.Vy.index)=md.results.DiagnosticSolution.Vy.value;
     37vz=zeros(md.numberofgrids,1);
     38vz(md.results.DiagnosticSolution.Vz.index)=md.results.DiagnosticSolution.Vz.value;
    3839
    39 plotmodel(md,'data',md.results.diagnostic.vx,'data',md.results.diagnostic.vy,'data',md.results.diagnostic.vz,'layer#all',md.numlayers)
     40plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.numlayers)

  • issm/trunk/test/Validation/ISMIP/TestE/Square.par

    r3768 r4883  
    3434disp('      boundary conditions for diagnostic model: ');
    3535%Create grid on boundary fist (because wi can not use mesh)
    36 md.gridonboundary=zeros(md.numberofgrids,1);
    37 pos=find(md.x==0 | md.y==0 | md.x==max(md.x) | md.y==max(md.y));
    38 md.gridonboundary(pos)=1;
    39 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
    40 pos=find(md.y==0 | md.y==max(md.y));md.gridondirichlet_diag(pos)=1;
    41 md.dirichletvalues_diag=zeros(md.numberofgrids,2);
    42 md.segmentonneumann_diag=zeros(0,3);
    43 md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
     36md=SetIceSheetBC(md);
     37
     38md.np=3;

  • issm/trunk/test/Validation/ISMIP/TestE/runme.m

    r899 r4883  
    1313md=extrude(md,10,1);
    1414
    15 md.gridondirichlet_diag(find(md.gridonbed))=1; %case 1: no slip
    16 %md.gridondirichlet_diag(find(md.gridonbed & (md.y<2200 | md.y> 2500)))=1; %case 2: slip between 2200 and 2500m
     15md=setelementstype(md,'pattyn','all');
     16%md=setelementstype(md,'pattyn','all','stokes','all');
    1717
    18 %md=setelementstype(md,'pattyn','all');
    19 md=setelementstype(md,'pattyn','all','stokes','all');
    20 
    21 md.gridonstokes(find(md.x==0 | md.x==Lx))=0;
    22 %Create MPCs to have periodic boundary conditions
    23 posx=find(md.x==0);
    24 posx2=find(md.x==max(md.x));
    25 
    26 md.penalties=[posx,posx2];
     18%md.gridonstokes(find(md.x==0 | md.x==Lx))=0;
     19%%Create MPCs to have periodic boundary conditions
     20%posx=find(md.x==0);
     21%posx2=find(md.x==max(md.x));
     22%
     23%md.penalties=[posx,posx2];
    2724
    2825%Compute the diagnostic
    29 md=solve(md,'analysis_type','diagnostic','package','ice');
     26md=solve(md,'analysis_type',DiagnosticSolutionEnum);
     27vx=zeros(md.numberofgrids,1);
     28vx(md.results.DiagnosticSolution.Vx.index)=md.results.DiagnosticSolution.Vx.value;
     29vy=zeros(md.numberofgrids,1);
     30vy(md.results.DiagnosticSolution.Vy.index)=md.results.DiagnosticSolution.Vy.value;
     31vz=zeros(md.numberofgrids,1);
     32vz(md.results.DiagnosticSolution.Vz.index)=md.results.DiagnosticSolution.Vz.value;
    3033
    31 plotmodel(md,'data',md.results.diagnostic.vx,'data',md.results.diagnostic.vy,'data',md.results.diagnostic.vz,'layer#all',md.numlayers)
     34plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.numlayers)

  • issm/trunk/test/Validation/ISMIP/TestF/Square.par

    r3768 r4883  
    3737md.connectivity=200;
    3838md.np=8;
    39 md.cluster='astrid';
     39md.cluster='wilkes';
    4040md.time=50;
    4141

  • issm/trunk/test/Validation/ISMIP/TestF/runme.m

    r2527 r4883  
    88md=model;
    99%md=squaremesh(md,L,L,nx,ny);
    10 md=mesh(md,'DomainOutline.exp',2500);
     10md=mesh(md,'DomainOutline.exp',5500);
    1111md=geography(md,'',''); %ice sheet test
    1212md=parameterize(md,'Square.par');
    13 md=extrude(md,10,1);
     13md=extrude(md,5,1);
    1414
    1515md=setelementstype(md,'pattyn','all');
    … …  
    3131
    3232%Compute the diagnostic
    33 md=solve(md,'analysis_type','transient');
     33md=solve(md,'analysis_type',Transient3DSolutionEnum);
     34vx=zeros(md.numberofgrids,1);
     35vx(md.results.Transient2DSolution(end).Vx.index)=md.results.Transient3DSolution.Vx(end).value;
     36vy=zeros(md.numberofgrids,1);
     37vy(md.results.Transient3DSolution.Vy(end).index)=md.results.Transient3DSolution.Vy(end).value;
     38vz=zeros(md.numberofgrids,1);
     39vz(md.results.Transient3DSolution.Vz(end).index)=md.results.Transient3DSolution(end).Vz.value;
    3440
    35 %plotmodel(md,'data',md.results.diagnostic.vx,'data',md.results.diagnostic.vy,'data',md.results.diagnostic.vz,'layer#all',md.numlayers,'contourlevels',{9.2,18.4,27.6,46,73.6,82.8},'contourlevels',{0,-1.2,-0.9,-0.6,-0.3,1.2},'contourlevels',{-5.2,-3.9,-1.3,0,1.3,3.9,5.2})
     41%plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.numlayers,'contourlevels',{9.2,18.4,27.6,46,73.6,82.8},'contourlevels',{0,-1.2,-0.9,-0.6,-0.3,1.2},'contourlevels',{-5.2,-3.9,-1.3,0,1.3,3.9,5.2})
Note: See TracChangeset for help on using the changeset viewer.