Changeset 4884

Timestamp:

07/29/10 12:49:25 (15 years ago)

Author:

seroussi

Message:

Validation test: Round Ice Sheet

Location:

issm/trunk/test/Validation/RoundIceSheet

Files:

• Circ.par (modified) (2 diffs)
• runme.m (modified) (2 diffs)

issm/trunk/test/Validation/RoundIceSheet/Circ.par

@@ -1 +1 @@
 %Ok, start defining model parameters here
+md.np=3;
+md.cluster='wilkes';
 
 disp('      creating thickness');
@@ -39 +41 @@
 disp('      boundary conditions for diagnostic model: ');
 %Build gridonicefront, array of boundary grids belonging to the icefront:
-gridonicefront=double(md.gridonboundary);
-
-md.gridondirichlet_diag=zeros(md.numberofgrids,1);
-pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
-md.dirichletvalues_diag=zeros(md.numberofgrids,2);
+md=SetMarineIceSheetBC(md,'Front.exp');
 
 radius=sqrt((md.x).*md.x+(md.y).*md.y);
-pos=find(radius==min(radius));md.gridondirichlet_diag(pos)=1;
+pos=find(radius==min(radius));
 md.x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center
 
-pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-md.segmentonneumann_diag=md.segments(pos,:);
-md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
-
-disp('      boundary conditions for thermal model: ');
-md.gridondirichlet_thermal=zeros(md.numberofgrids,1);
+md.spcvelocity(pos,1:3)=1;
+md.spcvelocity(pos,4:6)=0;

issm/trunk/test/Validation/RoundIceSheet/runme.m

@@ -3 +3 @@
 %This solution is valid for a circular contour
 
-%You can adjust the parameters numberofpoints and density to refine the mesh and check the precision of your solution. You can change the number of layers of the domain with the parameter layernumber.
+%You can adjust the parameters resolution(density) and the numberoflayers of the model to refine the mesh and check the precision of your solution. You can change the number of layers of the domain with the parameter layernumber.
+numlayers=10;
+resolution=40000;
 
 %To begin with the numerical model
 md=model;
 
-numberofgrids=20;
-density=190000;
-layernumber=4;
-
-%First we have to create the domain outline 
-contour=CreateDomainOutline(numberofgrids,750000);
-
-expgen('DomainOutline.exp',contour,1);
-md=mesh(md,'DomainOutline.exp',density)
-
-%move the closest node to the center
-[mini pos]=min(md.x.^2+md.y.^2);
-md.x(pos)=0;
-md.y(pos)=0;
-
-%Then we have to create the geography and the parameters of the mesh
+md=roundmesh(md,750000,resolution);
 md=geography(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
 
 %We now create a file Front.exp needed to parameterize but the analytical solution can't support an icefront so this file is empty
-front=CreateDomainOutline(numberofgrids,1000000);
-
+front=CreateDomainOutline(40,1000000);
 expgen('Front.exp',front,1);
-
 md=parameterize(md,'Circ.par');
 
@@ -40 +25 @@
 
 %We extrude the model to have a 3d model
-md=extrude(md,layernumber,3);
+md=extrude(md,numlayers,1);
 md=setelementstype(md,'Pattyn','all');
 
 %Spc the grids on the bed
 pos=find(md.gridonbed);
-md.gridondirichlet_diag(pos)=1;
+md.spcvelocity(pos,1:3)=1;
 
+md.verbose=10;
 %Now we can solve the problem 
-md=solve(md,'analysis_type','diagnostic','package','ice');
+md=solve(md,'analysis_type',DiagnosticSolutionEnum);
 
 %Calculate the depth averaged velocity field (2d):
 vel=zeros(md.numberofgrids2d,1);
 
+vx=zeros(md.numberofgrids,1);
+vx(md.results.DiagnosticSolution.Vx.index)=md.results.DiagnosticSolution.Vx.value;
+vy=zeros(md.numberofgrids,1);
+vy(md.results.DiagnosticSolution.Vy.index)=md.results.DiagnosticSolution.Vy.value;
+
 grid_vel=0;
 for i=1:md.numberofgrids2d
         for j=1:(md.numlayers-1)
-                grid_vel=grid_vel+1/(2*(md.numlayers-1))*(sqrt(md.results.diagnostic.vx(i+j*md.numberofgrids2d,1).^2+...
-                        md.results.diagnostic.vy(i+j*md.numberofgrids2d,1).^2)+...
-                        sqrt(md.results.diagnostic.vx(i+(j-1)*md.numberofgrids2d,1).^2+md.results.diagnostic.vy(i+(j-1)*md.numberofgrids2d,1).^2));
+                grid_vel=grid_vel+1/(2*(md.numlayers-1))*(sqrt(vx(i+j*md.numberofgrids2d,1).^2+...
+                        vy(i+j*md.numberofgrids2d,1).^2)+...
+                        sqrt(vx(i+(j-1)*md.numberofgrids2d,1).^2+vy(i+(j-1)*md.numberofgrids2d,1).^2));
         end
         vel(i,1)=grid_vel;