Index: /issm/trunk-jpl/examples/shakti/moulin.par =================================================================== --- /issm/trunk-jpl/examples/shakti/moulin.par (revision 22370) +++ /issm/trunk-jpl/examples/shakti/moulin.par (revision 22370) @@ -0,0 +1,58 @@ +%Start defining model parameters here + +% Set up bed topography and ice geometry for a tilted 500m thick slab + md.geometry.base = .02*md.mesh.x; + md.geometry.bed = md.geometry.base; + md.geometry.surface = .02*md.mesh.x + 500; + md.geometry.thickness = md.geometry.surface - md.geometry.bed; + +% Define ice sliding velocity (m/yr) + md.initialization.vx = 10^-6*md.constants.yts*ones(md.mesh.numberofvertices,1); + md.initialization.vy = zeros(md.mesh.numberofvertices,1); + md.initialization.pressure=zeros(md.mesh.numberofvertices,1); + +md.initialization.pressure=zeros(md.mesh.numberofvertices,1); + +% Materials + % Ice flow law parameter (note that the standard parameter A=B^(-3)) + md.materials.rheology_B= (5e-25)^(-1/3)*ones(md.mesh.numberofvertices,1); + md.initialization.temperature=(273)*ones(md.mesh.numberofvertices,1); + md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1); + +%Calving +md.calving.calvingrate=zeros(md.mesh.numberofvertices,1); +%md.calving.spclevelset=NaN(md.mesh.numberofvertices,1); + +% Friction law and coefficient + md.friction=frictionsommers(md.friction); + md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1); + + +%Numerical parameters +md.stressbalance.viscosity_overshoot=0.0; +md.masstransport.stabilization=1.; +md.thermal.stabilization=1.; +md.verbose=verbose(0); +md.settings.waitonlock=30; +md.stressbalance.restol=0.05; +md.steadystate.reltol=0.05; +md.stressbalance.reltol=0.05; +md.stressbalance.abstol=NaN; +md.timestepping.time_step=1.; +md.timestepping.final_time=3.; + +%GIA: +md.gia.lithosphere_thickness=100.*ones(md.mesh.numberofvertices,1); % in km +md.gia.mantle_viscosity=1.0*10^21*ones(md.mesh.numberofvertices,1); % in Pa.s +md.materials.lithosphere_shear_modulus=6.7*10^10; % in Pa +md.materials.lithosphere_density=3.32; % in g/cm^-3 +md.materials.mantle_shear_modulus=1.45*10^11; % in Pa +md.materials.mantle_density=3.34; % in g/cm^-3 + +%Boundary conditions: +md=SetIceSheetBC(md); + +%Change name so that no test have the same name +A=dbstack; +if (length(A)>2), md.miscellaneous.name=A(3).file(1:end-2); end + Index: /issm/trunk-jpl/examples/shakti/outline.exp =================================================================== --- /issm/trunk-jpl/examples/shakti/outline.exp (revision 22370) +++ /issm/trunk-jpl/examples/shakti/outline.exp (revision 22370) @@ -0,0 +1,10 @@ +## Name:DomainOutline +## Icon:0 +# Points Count Value +5 1.000000 +# X pos Y pos +0 0 +1000 0 +1000 1000 +0 1000 +0 0 Index: /issm/trunk-jpl/examples/shakti/runme.m =================================================================== --- /issm/trunk-jpl/examples/shakti/runme.m (revision 22370) +++ /issm/trunk-jpl/examples/shakti/runme.m (revision 22370) @@ -0,0 +1,81 @@ +steps=[1:3]; + +if any(steps==1) + disp(' Step 1: Mesh'); + + %Generate unstructured mesh on 1,000 m square with typical element edge length of 20 m + md=triangle(model,'./outline.exp',20); + + save MoulinMesh md +end +if any(steps==2) + disp(' Step 2: Parameterization'); + md=loadmodel('MoulinMesh'); + + md=setmask(md,'',''); + + % Run parameterization script to set up geometry, velocity, material properties, etc. + md=parameterize(md,'moulin.par'); + + % HYDROLOGY SPECIFIC PARAMETERIZATION: + % Change hydrology class to Sommers' SHaKTI model + md.hydrology=hydrologysommers(); + + % Define initial water head such that water pressure is 50% of ice overburden pressure + md.hydrology.head = 0.5*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base; + + % Initial subglacial gap height of 0.01m everywhere + md.hydrology.gap_height = 0.01*ones(md.mesh.numberofelements,1); + + % Typical bed bump bump spacing (2m) + md.hydrology.bump_spacing = 2*ones(md.mesh.numberofelements,1); + + % Typical bed bump height (0.1m) + md.hydrology.bump_height = 0.1*ones(md.mesh.numberofelements,1); + + % Define distributed englacial input to the subglacial system (m/yr) + % Change the value 0.0 to add distributed input + md.hydrology.englacial_input = 0.0*ones(md.mesh.numberofvertices,1); + + % Initial Reynolds number (start at Re=1000 everywhere) + md.hydrology.reynolds= 1000*ones(md.mesh.numberofelements,1); + + % Set up atmospheric pressure Type I boundary condition at left edge of + % domain (outflow, i.e. h=zb at x=xmin) + md.hydrology.spchead = NaN(md.mesh.numberofvertices,1); + pos=find(md.mesh.vertexonboundary & md.mesh.x==max(md.mesh.x)); + md.hydrology.spchead(pos)=md.geometry.base(pos); + + save MoulinParam md; +end +if any(steps==3); + disp(' Step 3: Solve!'); + md=loadmodel('MoulinParam'); + + md.transient=deactivateall(md.transient); + md.transient.ishydrology=1; + + % Specify that you want to run the model on your current computer + % Change the number of processors according to your machine (here np=4) + md.cluster=generic('np',2); + + % Define the time stepping scheme: run for 90 days with a time step of 1 hr + md.timestepping.time_step=3600/md.constants.yts; % Time step (in years) + md.timestepping.final_time=30/365; + + %Add one moulin with steady input at x=500, y=500 + [a,pos] = min(sqrt((md.mesh.x-500).^2+(md.mesh.y-500).^2)); + time=0:md.timestepping.time_step:md.timestepping.final_time; + md.hydrology.moulin_input = zeros(md.mesh.numberofvertices+1,numel(time)); + md.hydrology.moulin_input(end,:)=time; + md.hydrology.moulin_input(pos,:)=4; + + % Specify no-flux Type 2 boundary conditions on all edges (except + % the Type 1 condition set at the outflow above) + md.hydrology.neumannflux=zeros(md.mesh.numberofelements+1,numel(time)); + md.hydrology.neumannflux(end,:)=time; + + md=solve(md,'Transient'); + + save MoulinTransient md +end