Changeset 23898

Timestamp:

04/26/19 21:15:26 (6 years ago)

Author:

seroussi

Message:

CHG: changing ice part of the test

File:

• issm/trunk-jpl/test/NightlyRun/test4003.m (modified) (10 diffs)

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

@@ -4 +4 @@
 %Script control parameters
 steps=[1 2 3 4 5 6 7 8 9 10 11 12];
+steps=[ 7 8 9 10];
 final_time=1/365;
 
@@ -161 +162 @@
 if perform(org,'CreateMesh'),
         
-        loaddata(org,'Parameters');
-        loaddata(org,'Bathymetry');
-        loaddata(org,'IceSheetGeometry');
-
         %create model: 
         md=model();
         
         %Grab lat,long from MITgcm: 
-        lat=lat(:);
-        long=long(:);
+        long=readbin('run/XG.data',[3 200]);
+        long=[long long(:,end)]; long=[long; -105.0*ones(1,size(long,2))];
+        lat=readbin('run/YG.data',[3 200]);
+        lat=[lat -73.8832*ones(size(lat,1),1)]; lat=[lat; lat(end,:)];
 
         %project lat,long: 
-        %[x,y]=ll2xy(lat,long,-1);
-        x=long;
-        y=lat;
-
+        [x,y]=ll2xy(lat(:),long(:),-1);
+
+        Nx=size(lat,1); Ny=size(lat,2)
         index=[];
         %  C  D
@@ -193 +191 @@
         %fill mesh and model: 
         md=meshconvert(md,index,x,y);
-        md.mesh.lat=lat;
-        md.mesh.long=long;
+        md.mesh.lat=lat(:);
+        md.mesh.long=long(:);
 
         savemodel(org,md);
@@ -203 +201 @@
 if perform(org,'MeshGeometry'),
         
-        loaddata(org,'Parameters');
         loaddata(org,'CreateMesh');
-        loaddata(org,'Bathymetry');
-        loaddata(org,'IceSheetGeometry');
 
         %transfer to vertices: 
-        bathymetry=bathymetry(:);
-        iceshelf_mask=iceshelf_mask(:);
-        ice_mask=ice_mask(:);
-        thickness=thickness(:);
-        draft=draft(:);
+        bathymetry=readbin('run/bathy.box',[3 200],1,'real*8');
+        bathymetry=[bathymetry bathymetry(:,end)]; bathymetry=[bathymetry(1,:); bathymetry];
+        iceshelf_mask=-1*ones(size(bathymetry));
+        ice_mask=readbin('run/hmask3.box',[3 200],1,'real*8');
+        ice_mask=[ice_mask ice_mask(:,end)]; ice_mask=[ice_mask(1,:); ice_mask];
+        thickness=readbin('run/h0.bin',[3 200],1,'real*8');
+        thickness=[thickness thickness(:,end)]; thickness=[thickness; thickness(end,:)];
 
         %start filling some of the fields 
-        md.geometry.bed=bathymetry;
-        md.geometry.thickness=thickness;
-        md.geometry.base=md.geometry.bed; 
-        pos=find(iceshelf_mask); md.geometry.base(pos)=draft(pos);
+        md.geometry.bed=bathymetry(:);
+        md.geometry.thickness=thickness(:);
+        md.geometry.base=-917/1028*md.geometry.thickness; 
         md.geometry.surface=md.geometry.base+md.geometry.thickness;
 
         %nothing passes icefront: 
-        pos=find(~ice_mask);
+        pos=find((~ice_mask(:) & ice_mask(:)~=0) | thickness(:)==0);
         md.geometry.thickness(pos)=1;
         md.geometry.surface(pos)=(1-di)*md.geometry.thickness(pos);
@@ -229 +225 @@
 
         %level sets: 
-        md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
-        md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
-
-        pos=find(ice_mask); md.mask.ice_levelset(pos)=-1; 
-        pos=find(~iceshelf_mask & ice_mask); md.mask.groundedice_levelset(pos)=1;
-
-        %identify edges of grounded ice: 
-        groundedice_levelset=md.mask.groundedice_levelset;
-        for i=1:md.mesh.numberofelements,
-                m=groundedice_levelset(md.mesh.elements(i,:));
-                if abs(sum(m))~=3,
-                        pos=find(m==1); md.mask.groundedice_levelset(md.mesh.elements(i,pos))=0;
-                end
-        end
-
-        %identify edges of ice: 
-        ice_levelset=md.mask.ice_levelset;
-        for i=1:md.mesh.numberofelements,
-                m=ice_levelset(md.mesh.elements(i,:));
-                if abs(sum(m))~=3,
-                        pos=find(m==-1); md.mask.ice_levelset(md.mesh.elements(i,pos))=0;
-                end
-        end
+        md.mask.groundedice_levelset=iceshelf_mask(:);
+        md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
+
+        pos=find(~ice_mask(:) & thickness(:)==0); md.mask.ice_levelset(pos)=1; 
 
         savemodel(org,md);
@@ -259 +236 @@
 if perform(org,'ParameterizeIce'),
         
-        loaddata(org,'Parameters');
-        loaddata(org,'CreateMesh');
         loaddata(org,'MeshGeometry');
 
         %miscellaneous
-        md.miscellaneous.name='test4002';
+        md.miscellaneous.name='test4003';
 
         %initial velocity: 
@@ -281 +256 @@
         md.initialization.temperature=(273.15-22)*ones(md.mesh.numberofvertices,1);
         md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base);
-        md.smb.mass_balance = [1*ones(md.mesh.numberofvertices,1); 1];
+        md.smb.mass_balance = [0*ones(md.mesh.numberofvertices,1); 1];
 
         %Flow law 
@@ -297 +272 @@
         md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1); 
 
-        %deal with water: 
-        pos=find(md.mask.ice_levelset>0); 
-        md.stressbalance.spcvx(pos)=0;
+        %get some flux at the ice divide: 
+        pos=find(md.mesh.lat==min(md.mesh.lat));
+        md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
+        md.stressbalance.spcvx(pos)=-800;
         md.stressbalance.spcvy(pos)=0;
-        md.stressbalance.spcvz(pos)=0;
-        md.masstransport.spcthickness(pos)=0;
-
-        %get some flux at the ice divide: 
-        pos=find(md.mesh.y==min(md.mesh.y));
-        md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
-        md.stressbalance.spcvx(pos)=0;
-        md.stressbalance.spcvy(pos)=1500;
 
         %deal with boundaries, excluding icefront: 
-        vertex_on_boundary=zeros(md.mesh.numberofvertices,1);
-        vertex_on_boundary(md.mesh.segments(:,1:2))=1;
-        pos=find(vertex_on_boundary & md.mask.groundedice_levelset<=0);
-        md.stressbalance.spcvx(pos)=0;
+        pos=find(md.mesh.long==min(md.mesh.long) | md.mesh.long==max(md.mesh.long));
+        md.stressbalance.spcvy(pos)=0;
+
+        point1=find(md.mesh.y==min(md.mesh.y)); point2=find(md.mesh.x==max(md.mesh.x));
+        costheta=(md.mesh.x(point2)-md.mesh.x(point1))/sqrt((md.mesh.x(point2)-md.mesh.x(point1)).^2+(md.mesh.y(point2)-md.mesh.y(point1)).^2);
+        sintheta=(md.mesh.y(point2)-md.mesh.y(point1))/sqrt((md.mesh.x(point2)-md.mesh.x(point1)).^2+(md.mesh.y(point2)-md.mesh.y(point1)).^2);
+        md.stressbalance.referential(:,1:3)=repmat([costheta,sintheta,0],md.mesh.numberofvertices,1);
+        md.stressbalance.referential(:,4:6)=repmat([-sintheta,costheta,0],md.mesh.numberofvertices,1);
 
         md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
@@ -335 +307 @@
         %timestepping: 
         md.timestepping.final_time=100;
-        md.timestepping.time_step=1;
+        md.timestepping.time_step=0.25;
         md.transient.isgroundingline=0;
         md.transient.isthermal=0;
@@ -346 +318 @@
 
         savemodel(org,md);
+
+        plotmodel(md,'data',md.results.TransientSolution(end).Vel,'data',md.results.TransientSolution(end).Thickness)
 end
 % }}}