Index: /issm/workshop/14_JKS/Front.exp =================================================================== --- /issm/workshop/14_JKS/Front.exp (revision 11038) +++ /issm/workshop/14_JKS/Front.exp (revision 11038) @@ -0,0 +1,15 @@ +## Name:Front +## Icon:0 +# Points Count Value +9 1.000000 +# X pos Y pos +-443691.7204476296 -2241744.0087129781 +-444020.9148401438 -2244124.0631026439 +-443856.3176438867 -2245205.9060070375 +-441551.9568962874 -2244773.1688452801 +-440399.7765224879 -2243474.9573600078 +-439905.9849337166 -2240662.1658085845 +-441387.3597000304 -2238065.7428380400 +-443856.3176438867 -2237633.0056762826 +-443691.7204476296 -2241744.0087129781 + Index: /issm/workshop/14_JKS/Jak.par =================================================================== --- /issm/workshop/14_JKS/Jak.par (revision 11038) +++ /issm/workshop/14_JKS/Jak.par (revision 11038) @@ -0,0 +1,107 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%Greenland parameters%%%%%%%%%% + + %Parameters for Jakobshavn tutorial + + modeldatapath = ['~/issmjpl/projects/ModelData/SeaRISE/Greenland5km_v1.2']; + + %Name and hemisphere + md.mesh.hemisphere='n'; + md.miscellaneous.name='SeaRISEgreenland'; + + %some temporary parameters + di=md.materials.rho_ice/md.materials.rho_water; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Some hardcoded parameters for this deck%%%%%%%%%% + + ncdata=[modeldatapath '/Greenland_5km_dev1.2.nc']; + x1=double(ncread(ncdata,'x1')); + y1=double(ncread(ncdata,'y1')); + + surface=double(ncread(ncdata,'usrf')); + md.geometry.surface=InterpFromGridToMesh(x1,y1,surface',md.mesh.x,md.mesh.y,0); + + bed=double(ncread(ncdata,'topg')); + md.geometry.bed=InterpFromGridToMesh(x1,y1,bed',md.mesh.x,md.mesh.y,0); + md.geometry.thickness=md.geometry.surface-md.geometry.bed; + + pos0=find(md.geometry.thickness<=0); + md.geometry.thickness(pos0)=1; + + md.geometry.surface=md.geometry.thickness+md.geometry.bed; + + disp(' reading velocities '); + velx=double(ncread(ncdata,'surfvelx')); + vely=double(ncread(ncdata,'surfvely')); + + md.inversion.vx_obs=InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0); + md.inversion.vy_obs=InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0); + md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2); + + md.initialization.vx=md.inversion.vx_obs; + md.initialization.vy=md.inversion.vy_obs; + md.initialization.vz=zeros(md.mesh.numberofvertices,1); + md.initialization.vel=md.inversion.vel_obs; + + disp(' friction coefficient'); + md.friction.coefficient=50*ones(md.mesh.numberofvertices,1); + md.friction.q=ones(md.mesh.numberofelements,1); + md.friction.p=ones(md.mesh.numberofelements,1); + + disp(' loading temperature '); + temp=double(ncread(ncdata,'surftemp')); + md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp',md.mesh.x,md.mesh.y,0)+273.15; + + disp(' creating flow law paramters'); + md.materials.rheology_B=paterson(md.initialization.temperature); + pos=find(md.materials.rheology_B<=0); + md.materials.rheology_B(pos)=-md.materials.rheology_B(pos); + md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + + disp(' creating smb'); + smb=double(ncread(ncdata,'smb')); + smb=InterpFromGridToMesh(x1,y1,smb',md.mesh.x,md.mesh.y,0); + md.surfaceforcings.mass_balance=[smb*1000/md.materials.rho_ice;1]; + + %Deal with boundary conditions: + %dynamic: + md.diagnostic.icefront=[md.mesh.segments 1*ones(size(md.mesh.segments,1),1)]; + + %set one node to velocity 0,so as to not get a singular problem. We take the node + %on 0 velocity part of Greenland. + resolution=20000; x_sing=4.6*10^5; y_sing=-2.22*10^6; + pos=find((md.mesh.x>(x_sing-resolution)) & (md.mesh.x<(x_sing+resolution)) & (md.mesh.y>(y_sing-resolution)) & (md.mesh.y<(y_sing+resolution))); + + 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(pos)=0; + md.diagnostic.spcvy(pos)=0; + md.diagnostic.spcvz(pos)=0; + md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6); + + %set all nodes with 0 thickness to observed velocity + %md.diagnostic.spcvelocity=zeros(md.mesh.numberofvertices,6); + pos=find(md.inversion.vel_obs==0); + md.diagnostic.spcvx(pos)=0; + md.diagnostic.spcvy(pos)=0; + md.diagnostic.spcvz(pos)=0; + + %freeze water + pos=find(md.mask.vertexonwater); + md.diagnostic.spcvx(pos)=0; + md.diagnostic.spcvy(pos)=0; + md.diagnostic.spcvz(pos)=0; + + %prognostic + md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1); + + disp(' thermal model'); + md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1); + md.thermal.spctemperature=[md.initialization.temperature;1]; %impose observed temperature on surface + + md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.mesh.z); + + gflux=double(ncread(ncdata,'bheatflx')); + md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,gflux',md.mesh.x,md.mesh.y,0); + + md=SetMarineIceSheetBC(md,'./Front.exp'); Index: /issm/workshop/14_JKS/WeakB.exp =================================================================== --- /issm/workshop/14_JKS/WeakB.exp (revision 11038) +++ /issm/workshop/14_JKS/WeakB.exp (revision 11038) @@ -0,0 +1,60 @@ +## Name:WeakB +## Icon:0 +# Points Count Value +24 1.000000 +# X pos Y pos +-441022.7144520371 -2238435.7116546673 +-441639.8422268176 -2239053.0142011684 +-439171.3311276956 -2240699.1543251718 +-437319.9478033541 -2241727.9919026736 +-435879.9829955329 -2241933.7594181742 +-432382.9256051101 -2242756.8294801759 +-430942.9607972889 -2242962.5969956759 +-429708.7052477279 -2243991.4345731782 +-428680.1589564271 -2245020.2721506800 +-427651.6126651262 -2246049.1097281822 +-425800.2293407847 -2247077.9473056840 +-424565.9737912237 -2247695.2498521851 +-419341.0001274022 -2251604.8326466926 +-416049.6519952395 -2253868.2753171972 +-416049.6519952395 -2252222.1351931938 +-418312.4538361013 -2249341.3899761885 +-423331.7182416627 -2245843.3422126817 +-425594.5200825246 -2245226.0396661805 +-426623.0663738254 -2244402.9696041788 +-427445.9034068661 -2243991.4345731782 +-428885.8682146872 -2242756.8294801759 +-431148.6700555491 -2241522.2243871735 +-436497.1107703134 -2239670.3167476696 +-441022.7144520371 -2238435.7116546673 + +## Name:WeakB +## Icon:0 +# Points Count Value +24 1.000000 +# X pos Y pos +-442668.3885181184 -2245226.0396661805 +-440405.5866772566 -2245843.3422126817 +-438759.9126111752 -2246049.1097281822 +-437525.6570616142 -2246460.6447591828 +-435879.9829955329 -2246872.1797901839 +-434028.5996711914 -2248106.7848831862 +-431354.3793138092 -2249547.1574916886 +-430120.1237642482 -2249341.3899761885 +-429297.2867312076 -2250164.4600381902 +-427857.3219233864 -2249958.6925226897 +-426623.0663738254 -2250164.4600381902 +-424771.6830494839 -2250575.9950691909 +-423537.4274999229 -2250781.7625846914 +-422714.5904668822 -2251399.0651311926 +-415226.8149621988 -2259629.7657512082 +-419135.2908691420 -2259218.2307202076 +-422838.0575178250 -2255720.1829567007 +-424977.3923077441 -2252016.3676776937 +-426623.0663738254 -2251810.6001621932 +-428268.7404399067 -2251399.0651311926 +-429708.7052477279 -2251399.0651311926 +-431560.0885720694 -2251399.0651311926 +-443902.6440676794 -2246460.6447591828 +-442668.3885181184 -2245226.0396661805 + Index: /issm/workshop/14_JKS/runme.m =================================================================== --- /issm/workshop/14_JKS/runme.m (revision 11038) +++ /issm/workshop/14_JKS/runme.m (revision 11038) @@ -0,0 +1,103 @@ +steps=[1 2 3 4]; + +%Hard coded parameters +cluster=generic('name',oshostname,'np',2); +ncdata='./Greenland_5km_dev1.2.nc'; + +%Run Steps +if any(steps==1) + + disp(' Step 1: Mesh creation'); + md=triangle(model,'Jks3_trimmed.exp',2000); + + %Get observed velocity field on mesh nodes + x1=double(ncread(ncdata,'x1')); + y1=double(ncread(ncdata,'y1')); + velx=double(ncread(ncdata,'surfvelx')); + vely=double(ncread(ncdata,'surfvely')); + vx=InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0); + vy=InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0); + vel=sqrt(vx.^2+vy.^2); + + %refine mesh using surface velocities as metric + md=bamg(md,'hmin',1200,'hmax',15000,'field',vel,'err',5); + + name='JKS.Mesh_generation'; + save(name,'md','-v7.3') +end +if any(steps==2) + + disp(' Step 2: Parameterization'); + md=loadmodel('JKS.Mesh_generation'); + + md=setmask(md,'',''); + md=parameterize(md,'Jak.par'); + + %zones of shear margin softening + weakb=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,'WeakB.exp','node',2); + pos=find(weakb); + md.materials.rheology_B(pos)=.3*md.materials.rheology_B(pos); + + name='JKS.Parameterization'; + save(name,'md','-v7.3') +end +if any(steps==3) + + disp(' Step 3: Control method friction'); + md=loadmodel('JKS.Parameterization'); + md=setflowequation(md,'macayeal','all'); + + %Control general + md.inversion.iscontrol=1; + md.inversion.nsteps=20; + md.inversion.step_threshold=0.999*ones(md.inversion.nsteps,1); + md.inversion.maxiter_per_step=15*ones(md.inversion.nsteps,1); + md.verbose=verbose('solution',true,'control',true); + + %Cost functions + md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1); + md.inversion.cost_functions(1:md.inversion.nsteps,1)=103; + md.inversion.cost_functions(floor(md.inversion.nsteps/2):md.inversion.nsteps,1)=101; + + %Controls + md.inversion.control_parameters={'FrictionCoefficient'}; + md.inversion.gradient_scaling(1:md.inversion.nsteps)=50; + md.inversion.gradient_scaling(5:end)=10; + md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1); + md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1); + + %Additional parameters + md.diagnostic.restol=0.01; + md.diagnostic.reltol=0.1; + md.diagnostic.abstol=NaN; + md.inversion.cost_function_threshold=NaN; + + %Go solve + md=setcluster(md,cluster); + md.solver=addoptions(md.solver,NoneAnalysisEnum,asmoptions); + md.solver=addoptions(md.solver,DiagnosticVertAnalysisEnum,jacobiasmoptions); + md=solve(md,DiagnosticSolutionEnum); + + name='JKS.Control_friction_macayeal'; + save(name,'md','-v7.3') +end +if any(steps==4) + + disp(' Plotting') + md=loadmodel('JKS.Control_friction_macayeal'); + + plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',... + 'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,... + 'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,... + 'FontSize#all',12,... + 'data',md.initialization.vel,'title','Observed velocity',... + 'data',md.results.DiagnosticSolution.Vel,'title','Modeled Velocity',... + 'colorbar#1','off','colorbar#2','on','colorbartitle#2','[m/yr]',... + 'caxis#1-2',[0,7000],... + 'data',md.geometry.bed,'title','Bed elevation',... + 'data',md.results.DiagnosticSolution.FrictionCoefficient,... + 'title','Friction Coefficient',... + 'colorbar#3','on','colorbartitle#3','[m]','colorbar#4','on'); + +end +