Changeset 5375

Timestamp:

08/18/10 13:58:29 (15 years ago)

Author:

seroussi

Message:

test eismint mass conservation

File:

• issm/trunk/test/NightlyRun/test1201.m (modified) (1 diff)

issm/trunk/test/NightlyRun/test1201.m

@@ -1 @@
-%This test is a test from the EISMINT for Ice shelves 
-% Vincent Rommelaere 1996
+%This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996
 
-%The goal is to test the prognostic model
-%md=bamg(model,'domain','DomainOutline.exp','hmax',4550);
-md=bamg(model,'domain','../Exp/SquareEISMINT.exp','hmax',3000);
-md=geography(md,'all','');
-md=parameterize(md,'../Par/SquareEISMINT.par');
-md=setelementstype(md,'macayeal','all');
+results={};
 
-disp('      initial velocity');
-md.vx=zeros(md.numberofgrids,1);
-md.vy=-400*ones(md.numberofgrids,1);
-md.vz=zeros(md.numberofgrids,1);
+for type=1:3;
+        %The goal is to test the prognostic model
+        md=bamg(model,'domain','../Exp/SquareEISMINT.exp','hmax',3000);
+        md=geography(md,'all','');
+        md=parameterize(md,'../Par/SquareEISMINT.par');
+        md=setelementstype(md,'macayeal','all');
 
-%analytical test
-md.spcvelocity(:,1:3)=1;
-md.spcvelocity(:,4)=0;
-md.spcvelocity(:,5)=-400;
-md.spcvelocity(1,1)=0;
-md.prognostic_DG=1;
+        disp('      initial velocity');
+        md.vx=zeros(md.numberofgrids,1);
+        md.vy=-400*ones(md.numberofgrids,1);
+        md.vz=zeros(md.numberofgrids,1);
 
-%Launch many prognostic as we don't really want a transient
-%FOR NOW:
-i=0;
-time=0;
-md.dummy=struct();
-pos=find(md.y>199999.9);
-connectivity=full(sparse(md.elements(:),1,1));
-while(time<500),
-        disp(['step ' num2str(i) '/' num2str(500/md.dt)]);
-        i=i+1;
-        time=time+md.dt;
-        md=solve(md,'analysis_type',PrognosticSolutionEnum);
-        thickness   =zeros(md.numberofgrids,1);
-        connectivity=sparse(md.results.PrognosticSolution.Thickness.index(:),1,1);
-        thickness   =sparse(md.results.PrognosticSolution.Thickness.index(:),1,md.results.PrognosticSolution.Thickness.value(:));
-        md.thickness=full(thickness./connectivity);
-        md.thickness(pos)=500+100*sin(2*pi*time/200)*ones(size(pos,1),1);
-        md.surface=md.bed+md.thickness;
-        md.dummy(i).thickness=md.thickness;
-        %plotmodel(md,'data',md.thickness,'sectionvalue','CrossLine.exp')
+        %analytical test
+        md.spcvelocity(:,1:3)=1;
+        md.spcvelocity(:,4)=0;
+        md.spcvelocity(:,5)=-400;
+        md.spcvelocity(1,1)=0;
+        if type==1,
+                md.artificial_diffusivity=1;
+                md.prognostic_DG=0;
+        elseif type==2,
+                md.artificial_diffusivity=0;
+                md.prognostic_DG=0;
+        elseif type==3,
+                md.prognostic_DG=1;
+        end
+
+        %Launch many prognostic as we don't really want a transient
+        %FOR NOW:
+        i=0;
+        time=0;
+        md.dummy=struct();
+        pos=find(md.y>199999.9);
+        connectivity=full(sparse(md.elements(:),1,1));
+        while(time<500),
+                disp(['step ' num2str(i) '/' num2str(500/md.dt)]);
+                i=i+1;
+                time=time+md.dt;
+                md=solve(md,'analysis_type',PrognosticSolutionEnum);
+                md.thickness=PatchToVec(md.results.PrognosticSolution.Thickness);
+                md.thickness(pos)=500+100*sin(2*pi*time/200)*ones(size(pos,1),1);
+                md.surface=md.bed+md.thickness;
+                md.dummy(i).thickness=md.thickness;
+        end
+        results{type}=md.thickness;
 end
 
 %plot results
-plotmodel(md,'data',md.thickness)
-%plotmodel(md,'data',[md.thickness 500+100*sin(2*pi/200*(500-md.y/400))],'sectionvalue','../Exp/CrossLineEISMINT.exp')
+[elements,x,y,z,s,h1]=SectionValues(md,results{1},'../Exp/CrossLineEISMINT.exp',100);
+[elements,x,y,z,s,h2]=SectionValues(md,results{2},'../Exp/CrossLineEISMINT.exp',100);
+[elements,x,y,z,s,h3]=SectionValues(md,results{3},'../Exp/CrossLineEISMINT.exp',100);
+[elements,x,y,z,s,hth]=SectionValues(md, 500+100*sin(2*pi/200*(500-md.y/400)),'../Exp/CrossLineEISMINT.exp',100);
+plot(s,h1,'r',s,h2,'b',s,h3,'g',s,hth,'k')
+legend('Art. diff.','No Art. diff.','Discontinue','Theorique')
+
+set(gcf,'Color','w')
+printmodel('eismintmasscthickness','png','margin','on','marginsize',25,'frame','off','resolution',2,'hardcopy','off');
 
 %Fields and tolerances to track changes
 field_names     ={ ...
-        'Thickness', ...
+        'ThicknessArtDigg','ThicknessNoArtDiff','ThicknessDG' ...
 };
 field_tolerances={...
-        1e-13, ...
+        1e-13, 1e-13, 1e-13...
 };
 field_values={
-        md.thickness, ...
+        results{1}, ...
+        results{2}, ...
+        results{3}, ...
 };