Index: /issm/trunk-jpl/examples/Jakobshavn/Domain.exp =================================================================== --- /issm/trunk-jpl/examples/Jakobshavn/Domain.exp (revision 11676) +++ /issm/trunk-jpl/examples/Jakobshavn/Domain.exp (revision 11676) @@ -0,0 +1,34 @@ +## Name:JksOutline +## Icon:0 +# Points Count Value +28 1.000000 +# X pos Y pos +-442320.3882147236 -2277080.6549030584 +-415649.3149820761 -2315476.9356749961 +-356259.6577274340 -2351834.4758754564 +-309585.2795717255 -2384114.5349311619 +-279812.9187561168 -2353193.6362587260 +-274850.8586214480 -2287953.9379578433 +-294854.1635456266 -2263149.2609987380 +-309585.2795753506 -2232228.3623235333 +-348041.2456311518 -2191113.7607873622 +-376728.1557929466 -2197909.5626934194 +-391614.3362011510 -2203685.9943137392 +-417355.0231571503 -2206404.3150757281 +-420122.3381705935 -2211830.0535394503 +-421889.2049176788 -2218344.8015625700 +-420705.1891998384 -2226200.5387845989 +-425774.4833068820 -2230306.5097973733 +-429242.2244052034 -2235353.9400299159 +-427025.7933000360 -2240741.0910965428 +-431045.5636035603 -2255158.9014255330 +-424708.0734753361 -2261343.7626088522 +-423137.8820362946 -2268466.9593191240 +-428054.2739333302 -2272589.0010017036 +-432362.3262468725 -2272393.9161962867 +-435068.8255516959 -2273140.2906106352 +-437641.3726878429 -2272612.1952726739 +-439529.7800752302 -2272929.1351353633 +-442475.4525940326 -2275381.7044264446 +-442320.3882147236 -2277080.6549030584 + Index: /issm/trunk-jpl/examples/Jakobshavn/Front.exp =================================================================== --- /issm/trunk-jpl/examples/Jakobshavn/Front.exp (revision 11676) +++ /issm/trunk-jpl/examples/Jakobshavn/Front.exp (revision 11676) @@ -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/trunk-jpl/examples/Jakobshavn/Jks.par =================================================================== --- /issm/trunk-jpl/examples/Jakobshavn/Jks.par (revision 11676) +++ /issm/trunk-jpl/examples/Jakobshavn/Jks.par (revision 11676) @@ -0,0 +1,63 @@ + +%Name and hemisphere +md.mesh.hemisphere='n'; +md.miscellaneous.name='Jakobshavn'; + +%Load SeaRISE dataset +disp(' Loading SeaRISE data from NetCDF'); +ncdata = './Greenland_5km_v1.1.nc'; +x1 = ncread(ncdata,'x1'); +y1 = ncread(ncdata,'y1'); +thk = ncread(ncdata,'thk')'; +topg = ncread(ncdata,'topg')'; +velx = ncread(ncdata,'surfvelx')'; +vely = ncread(ncdata,'surfvely')'; +temp = ncread(ncdata,'airtemp2m')'; +smb = ncread(ncdata,'smb')'; + +disp(' Interpolating thicknesses'); +md.geometry.thickness=InterpFromGridToMesh(x1,y1,thk,md.mesh.x,md.mesh.y,0); +pos0=find(md.geometry.thickness<=10); +md.geometry.thickness(pos0)=10; + +disp(' Interpolating bedrock topography'); +md.geometry.bed=InterpFromGridToMesh(x1,y1,topg,md.mesh.x,md.mesh.y,0); + +disp(' Constructing surface elevation'); +md.geometry.surface=md.geometry.thickness+md.geometry.bed; + +disp(' Interpolating velocities'); +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; +md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.mesh.z); + +disp(' Interpolating temperatures'); +md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15; %convert to Kelvin + +disp(' Interpolating surface mass balance'); +md.surfaceforcings.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0); +md.surfaceforcings.mass_balance=md.surfaceforcings.mass_balance*md.materials.rho_water/md.materials.rho_ice; + +disp(' Construct basal friction parameters'); +md.friction.coefficient=30*ones(md.mesh.numberofvertices,1); +pos=find(md.mask.vertexonfloatingice); +md.friction.coefficient(pos)=0; %no friction applied on floating ice +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=ones(md.mesh.numberofelements,1); + +disp(' Construct ice rheological properties'); +md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); +md.materials.rheology_B=paterson(md.initialization.temperature); +%Reduce viscosity along the shear margins +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); + +disp(' Set other boundary conditions'); +md=SetMarineIceSheetBC(md,'./Front.exp'); +md.diagnostic.icefront(:,end)=1; %Apply water pressure at the ice front +md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/examples/Jakobshavn/WeakB.exp =================================================================== --- /issm/trunk-jpl/examples/Jakobshavn/WeakB.exp (revision 11676) +++ /issm/trunk-jpl/examples/Jakobshavn/WeakB.exp (revision 11676) @@ -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/trunk-jpl/examples/Jakobshavn/runme.m =================================================================== --- /issm/trunk-jpl/examples/Jakobshavn/runme.m (revision 11676) +++ /issm/trunk-jpl/examples/Jakobshavn/runme.m (revision 11676) @@ -0,0 +1,92 @@ +steps=[1]; + +if any(steps==1) + disp(' Step 1: Mesh creation'); + md=triangle(model,'Domain.exp',2000); + + %Get observed velocity field on mesh nodes + ncdata='Greenland_5km_v1.1.nc'; + if ~exist(ncdata,'file'), + error('Download Greenland5km_v1.1.nc first on http://websrv.cs.umt.edu/isis/index.php/Present_Day_Greenland'); + end + x1 = ncread(ncdata,'x1'); + y1 = ncread(ncdata,'y1'); + velx = ncread(ncdata,'surfvelx'); + vely = 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); + [md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,+1,39,71); + + save JksMesh.mdl md +end +if any(steps==2) + disp(' Step 2: Parameterization'); + md=loadmodel('JksMesh.mdl'); + + md=setmask(md,'',''); + md=parameterize(md,'Jks.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); + + save JksPar.mdl md +end +if any(steps==3) + disp(' Step 3: Control method friction'); + md=loadmodel('JksPar.mdl'); + + md=setflowequation(md,'macayeal','all'); + + %Control general + md.inversion.iscontrol=1; + md.inversion.nsteps=20; + md.inversion.step_threshold=0.99*ones(md.inversion.nsteps,1); + md.inversion.maxiter_per_step=5*ones(md.inversion.nsteps,1); + md.verbose=verbose('solution',true,'control',true); + + %Cost functions + md.inversion.cost_functions=[101*ones(md.inversion.nsteps,1) 103*ones(md.inversion.nsteps,1)]; + md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,2); + md.inversion.cost_functions_coefficients(:,1)=40; + md.inversion.cost_functions_coefficients(:,2)=1; + + %Controls + md.inversion.control_parameters={'FrictionCoefficient'}; + md.inversion.gradient_scaling(1:md.inversion.nsteps)=30; + 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; + + %Go solve + md.cluster=generic('name',oshostname,'np',4); + md.solver=addoptions(md.solver,NoneAnalysisEnum,asmoptions); + md.solver=addoptions(md.solver,DiagnosticVertAnalysisEnum,jacobiasmoptions); + md=solve(md,DiagnosticSolutionEnum); + + save JksControl.mdl md +end +if any(steps==4) + disp(' Plotting') + md=loadmodel('JksControl.mdl'); + + plotmodel(md,'unit#all','km','axis#all','equal',... + 'FontSize#all',12,... + 'data',md.inversion.vel_obs,'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 Index: /issm/trunk-jpl/examples/SquareIceShelf/runme.m =================================================================== --- /issm/trunk-jpl/examples/SquareIceShelf/runme.m (revision 11675) +++ /issm/trunk-jpl/examples/SquareIceShelf/runme.m (revision 11676) @@ -4,3 +4,13 @@ md=parameterize(md,'Square.par'); md=setflowequation(md,'macayeal','all'); -md=solve(md,DiagnosticSolutionEnum); +md.cluster=generic('name',oshostname,'np',4); + +md.initialization.pressure=zeros(md.mesh.numberofvertices,1); +md.transient.isthermal=0; +md.timestepping.final_time=10; +md.timestepping.time_step=.5; +md=solve(md,TransientSolutionEnum); + +md.geometry.surface=md.results.TransientSolution(end).Surface; +md.geometry.bed=md.results.TransientSolution(end).Bed; +md.geometry.thickness=md.results.TransientSolution(end).Thickness;