Changeset 11191

Timestamp:

01/24/12 14:01:23 (13 years ago)

Author:

seroussi

Message:

finally finished to fix the NR !

Location:

issm/trunk-jpl/test

Files:

• Archives/Archive1108.mat (modified) ( previous)
• NightlyRun/test1108.m (modified) (2 diffs)

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

@@ -6 +6 @@
 
 for i=1:length(L_list),
-        L=3*L_list{i};
-        nx=60; %numberof nodes in x direction
-        ny=60;
+        L=L_list{i};
+        nx=30; %numberof nodes in x direction
+        ny=30;
         md=model;
         md=squaremesh(md,L,L,nx,ny);
         md=setmask(md,'',''); %ice sheet test
         md=parameterize(md,'../Par/ISMIPD.par');
-        md=extrude(md,6,1);
+        md=extrude(md,10,1);
 
-        md=setflowequation(md,'stokes','all');
+        md=setflowequation(md,'pattyn','all');
 
         %We need one grd on dirichlet: the 4 corners are set to zero
-        %md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
-        %md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
-        %md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+        md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+        md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
         
-        %pos=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
-        %md.diagnostic.spcvx(pos)=0;
-        %md.diagnostic.spcvy(pos)=0;
-        %md.diagnostic.spcvz(pos)=0;
+        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)));
+        md.diagnostic.spcvx(pos)=0;
+        md.diagnostic.spcvy(pos)=0;
+        md.diagnostic.spcvz(pos)=0;
 
-        %%Create MPCs to have periodic boundary conditions
-        %posx=find(md.mesh.x==0);
-        %posx2=find(md.mesh.x==max(md.mesh.x));
+        %Create MPCs to have periodic boundary conditions
+        posx=find(md.mesh.x==0);
+        posx2=find(md.mesh.x==max(md.mesh.x));
 
-        %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
-        %posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
+        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
+        posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
 
-        %md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
+        md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
 
         %Compute the diagnostic
         md.cluster=generic('name',oshostname(),'np',8);
+        md.verbose=verbose('convergence',true);
+        md=solve(md,DiagnosticSolutionEnum);
+        md.diagnostic.reltol=NaN;
+        md.diagnostic.abstol=NaN;
+        md.diagnostic.vertex_pairing=[];
+        %We need one grd on dirichlet: the 4 corners are set to zero
+        md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+        md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+        md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+        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
+        md.diagnostic.spcvx(pos)=md.results.DiagnosticSolution.Vx(pos);
+        md.diagnostic.spcvy(pos)=md.results.DiagnosticSolution.Vy(pos);
+        md=setflowequation(md,'stokes','all');
         md=solve(md,DiagnosticSolutionEnum);
 
@@ -57 +70 @@
         'Vx80km','Vy80km','Vz80km',...
         'Vx160km','Vy160km','Vz160km'
-}
+};
 field_tolerances={...
         1e-11,1e-11,1e-11,...