Index: /issm/workshop/08ISMIP/CheatyIsmipA.par
===================================================================
--- /issm/workshop/08ISMIP/CheatyIsmipA.par	(revision 18038)
+++ /issm/workshop/08ISMIP/CheatyIsmipA.par	(revision 18038)
@@ -0,0 +1,52 @@
+%Parameterization for ISMIP A experiment
+
+%Set the Simulation generic name #md.miscellaneous
+%->
+
+%Geometry
+disp('   Constructing Geometry');
+
+%Define the geometry of the simulation #md.geometry
+%surface is [-x*tan(0.5*pi/180)] #md.mesh
+%->
+md.geometry.surface=-md.mesh.x*tan(0.5*pi/180.);
+%base is [surface-1000+500*sin(x*2*pi/L).*sin(y*2*pi/L)]
+%L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x)
+%->
+L=max(md.mesh.x)-min(md.mesh.x);
+md.geometry.base=md.geometry.surface-1000.0+500.0*sin(md.mesh.x*2.0*pi/L).*sin(md.mesh.y*2.0*pi/L);
+%thickness is the difference between surface and base #md.geometry
+%->
+md.geometry.thickness=md.geometry.surface-md.geometry.base;
+%plot the geometry to check it out
+%->
+plotmodel(md,'data',md.geometry.thickness);
+
+disp('   Defining friction parameters');
+
+%These parameters will not be used but need to be fixed #md.friction
+%one friciton coefficient per node (md.mesh.numberofvertices,1)
+%->
+md.friction.coefficient=200.0*ones(md.mesh.numberofvertices,1);
+%one friciton exponent (p,q) per element
+%->
+md.friction.p=ones(md.mesh.numberofelements,1);
+%->
+md.friction.q=ones(md.mesh.numberofelements,1);
+
+disp('   Construct ice rheological properties');
+
+%The rheology parameters sit in the material section #md.materials
+%B has one value per vertex
+%->
+md.materials.rheology_B=6.8067e7*ones(md.mesh.numberofvertices,1);
+%n has one value per element
+%->
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+
+disp('   Set boundary conditions');
+
+%Set the default boundary conditions for an ice-sheet 
+% #help SetIceSheetBC
+%->
+md=SetIceSheetBC(md);
Index: /issm/workshop/08ISMIP/CheatyIsmipF.par
===================================================================
--- /issm/workshop/08ISMIP/CheatyIsmipF.par	(revision 18038)
+++ /issm/workshop/08ISMIP/CheatyIsmipF.par	(revision 18038)
@@ -0,0 +1,66 @@
+%Parameterization for ISMIP F experiment
+
+%Set the Simulation generic name #md.miscellaneous
+%->
+
+%Geometry
+disp('   Constructing Geometry');
+
+%Define the geometry of the simulation #md.geometry
+%surface is [-x*tan(3.0*pi/180)] #md.mesh
+%->
+md.geometry.surface=md.mesh.x*tan(3.0*pi/180.0);
+%base is [surface-1000+100*exp(-((x-L/2).^2+(y-L/2).^2)/(10000.^2))]
+%L is the size of the side of the square #max(md.mesh.x)-min(md.mesh.x)
+%->
+L=max(md.mesh.x)-min(md.mesh.x);
+%->
+md.geometry.base=md.geometry.surface-1000.0+100.0*exp(-((md.mesh.x-L/2.0).^2.0+(md.mesh.y-L/2.0).^2.0)/(10000.^2.0));
+%thickness is the difference between surface and base #md.geometry
+%->
+md.geometry.thickness=md.geometry.surface-md.geometry.base;
+%plot the geometry to check it out
+%->
+plotmodel(md,'data',md.geometry.thickness);
+
+disp('   Defining friction parameters');
+
+%These parameters will not be used but need to be fixed #md.friction
+%one friciton coefficient per node (md.mesh.numberofvertices,1)
+%conversion form year to seconds with #md.constants.yts
+%->
+md.friction.coefficient=sqrt(md.constants.yts/(1000*2.140373*10^-7))*ones(md.mesh.numberofvertices,1);
+%one friciton exponent (p,q) per element
+%->
+md.friction.p=ones(md.mesh.numberofelements,1);
+%->
+md.friction.q=zeros(md.mesh.numberofelements,1);
+
+disp('   Construct ice rheological properties');
+
+%The rheology parameters sit in the material section #md.materials
+%B has one value per vertex
+%->
+md.materials.rheology_B=(1/(2.140373*10^-7/md.constants.yts))*ones(md.mesh.numberofvertices,1);
+%n has one value per element
+%->
+md.materials.rheology_n=1*ones(md.mesh.numberofelements,1);
+
+disp('   Set boundary conditions');
+
+%Set the default boundary conditions for an ice-sheet 
+% #help SetIceSheetBC
+%->
+md=SetIceSheetBC(md);
+
+disp('   Initializing velocity and pressure');
+
+%initialize the velocity and pressurefields of #md.initialization
+%->
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+%->
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+%->
+md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+%->
+md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
Index: /issm/workshop/08ISMIP/CheatyRunme.m
===================================================================
--- /issm/workshop/08ISMIP/CheatyRunme.m	(revision 18038)
+++ /issm/workshop/08ISMIP/CheatyRunme.m	(revision 18038)
@@ -0,0 +1,232 @@
+%wich steps to perform steps are from 1 to 8 
+%step 7 is specific to ISMIPA
+%step 8 is specific to ISMIPF
+
+steps=[1 2 3 4 5 6 7];
+
+% parameter file to be used, choose between CheatyIsmipA.par or CheatyIsmipF.par
+ParamFile='CheatyIsmipA.par'
+
+name_prefix='ISMIP';
+%Run Steps
+org=organizer('repository','./Models','prefix',name_prefix,'steps',steps);
+
+% {{{ Mesh Generation #1
+if perform(org,'Mesh_Generation')
+
+	%initialize md as a new model #help model
+	%->
+	md=model;
+	% generate a squaremesh #help squaremesh
+	% Side is 80 km long with 20 points
+	%->
+	if(ParamFile=='CheatyIsmipA.par'),
+		md=squaremesh(md,80000,80000,20,20);
+	elseif(ParamFile=='CheatyIsmipF.par'),
+		md=squaremesh(md,100000,100000,30,30);
+  end
+	% plot the given mesh #plotdoc
+	%->
+	plotmodel(md,'data','mesh')
+	% save the given model
+	%->
+	savemodel(org,md);
+end
+% }}}
+
+% {{{ Masks #2
+if perform(org,'SetMask'), 
+
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+	md=loadmodel(org,'Mesh_Generation');
+	% set the mask #help setmask
+	% all MISMIP nodes are grounded
+	%->
+	md=setmask(md,'','');
+	% plot the given mask #md.mask to locate the field
+	%->
+	plotmodel(md,'data',md.mask.groundedice_levelset);
+	% save the given model
+	%->
+	savemodel(org,md);
+end
+% }}}
+
+% {{{ Parameterization #3
+if perform(org,'Parameterization') 
+
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+	md=loadmodel(org,'SetMask');
+	% parametrize the model # help parameterize
+	% you will need to fil-up the parameter file defined by the
+  % ParamFile variable
+	%->
+	md=parameterize(md,ParamFile);
+	% save the given model
+	%->
+	savemodel(org,md);
+end
+% }}}
+
+% {{{ Extrusion #4
+if perform(org,'Extrusion')
+	
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+	md=loadmodel(org,'Parameterization');
+	% vertically extrude the preceding mesh #help extrude
+	% only 5 layers exponent 1
+	%->
+	md=extrude(md,9,1);
+	% plot the 3D geometry #plotdoc
+	%->
+	plotmodel(md,'data',md.geometry.base)
+	% save the given model
+	%->
+	savemodel(org,md);
+end
+% }}}
+
+% {{{ Set the flow computing method #5
+
+if perform(org,'SetFlow')
+
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+	md=loadmodel(org,'Extrusion');
+	% set the approximation for the flow computation #help setflowequation
+	% We will be using the Higher Order Model (HO)
+	%->
+	md=setflowequation(md,'HO','all');
+	% save the given model
+	%->
+	savemodel(org,md);
+end
+% }}}
+
+% {{{ Set Boundary Conditions #6
+if perform(org,'BoundaryCondition')
+
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+	md=loadmodel(org,'SetFlow');
+	% dirichlet boundary condition are known as SPCs
+	% ice frozen to the base, no velocity	#md.stressbalance
+	% SPCs are initialized at NaN one value per vertex
+	%->
+	md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+	%->
+	md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+	%->
+	md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+	% extract the nodenumbers at the base #md.mesh.vertexonbase
+	%->
+	basalnodes=find(md.mesh.vertexonbase);
+  % set the sliding to zero on the bed
+	%->
+	md.stressbalance.spcvx(basalnodes)=0.0;
+	%->
+	md.stressbalance.spcvy(basalnodes)=0.0;
+	% periodic boundaries have to be fixed on the sides
+	% create tabs with the side of the domain
+	% for x
+	% create maxX #help find
+	%->
+	maxX=find(md.mesh.x==max(md.mesh.x));
+	% create minX
+	%->
+	minX=find(md.mesh.x==min(md.mesh.x));
+	% for y, max X and minX should be excluded
+	% create maxY
+	%->
+	maxY=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x));
+	% create minY
+	%->
+	minY=find(md.mesh.y==min(md.mesh.y) & md.mesh.x~=max(md.mesh.x) & md.mesh.x~=min(md.mesh.x));
+	% set the node that should be paired together
+	% #md.stressbalance.vertex_pairing
+	%->
+	md.stressbalance.vertex_pairing=[minX,maxX;minY,maxY];
+	if (ParamFile=='CheatyIsmipF.par')
+		% if we are dealing with IsmipF the solution is in
+		% masstransport
+		md.masstransport.vertex_pairing=md.stressbalance.vertex_pairing;
+  end
+	% save the given model
+	%->
+	savemodel(org,md);
+end
+% }}}
+
+% {{{ Solving #7
+if perform(org,'SolvingISMIPA')
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+	md=loadmodel(org,'BoundaryCondition');
+	% Set cluster #md.cluster
+	% generic parameters #help generic
+	% set only the name and number of process
+	%->
+	md.cluster=generic('name',oshostname(),'np',2);
+	% Set which control message you want to see #help verbose
+	%->
+	md.verbose=verbose('convergence',true);
+	% Solve #help solve
+	% we are solving a StressBalanc
+	%->
+	md=solve(md,StressbalanceSolutionEnum());
+	% save the given model
+	%->
+	savemodel(org,md);
+	% plot the surface velocities #plotdoc
+	%->
+	plotmodel(md,'data',md.results.StressbalanceSolution.Vel)
+end
+% }}}
+
+% {{{ Solving #8
+if perform(org,'SolvingISMIPF')
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+	md=loadmodel(org,'BoundaryCondition');
+	% Set cluster #md.cluster
+	% generic parameters #help generic
+	% set only the name and number of process
+	%->
+	md.cluster=generic('name',oshostname(),'np',2);
+	% Set which control message you want to see #help verbose
+	%->
+	md.verbose=verbose('convergence',true);
+	% set the transient model to ignore the thermal model
+	% #md.transient 
+	%->
+	md.transient.isthermal=0;
+	% define the timestepping scheme
+	% everything here should be provided in years #md.timestepping
+	% give the length of the time_step (4 years)
+	%->
+	md.timestepping.time_step=4;
+	% give final_time (20*4 years time_steps)
+	%->
+	md.timestepping.final_time=4*20;
+	% Solve #help solve
+	% we are solving a TransientSolution
+	%->
+	md=solve(md,TransientSolutionEnum);
+	% save the given model
+	%->
+	savemodel(org,md);
+	% plot the surface velocities #plotdoc
+	%->
+	plotmodel(md,'data',md.results.TransientSolution(20).Vel)
+end
+% }}}
Index: /issm/workshop/08ISMIP/runme.m
===================================================================
--- /issm/workshop/08ISMIP/runme.m	(revision 18038)
+++ /issm/workshop/08ISMIP/runme.m	(revision 18038)
@@ -0,0 +1,226 @@
+%wich steps to perform steps are from 1 to 8 
+%step 7 is specific to ISMIPA
+%step 8 is specific to ISMIPF
+
+steps=[7];
+
+% parameter file to be used, choose between IsmipA.par or IsmipF.par
+ParamFile='IsmipA.par'
+
+name_prefix='ISMIP';
+%Run Steps
+org=organizer('repository','./Models','prefix',name_prefix,'steps',steps);
+
+% {{{ Mesh Generation #1
+if perform(org,'Mesh_Generation')
+
+	%initialize md as a new model #help model
+	%->
+
+	% generate a squaremesh #help squaremesh
+	% Side is 80 km long with 20 points
+	%->
+
+	% plot the given mesh #plotdoc
+	%->
+
+	% save the given model
+	%->
+
+end
+% }}}
+
+% {{{ Masks #2
+if perform(org,'SetMask'), 
+
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+
+	% set the mask #help setmask
+	% all MISMIP nodes are grounded
+	%->
+
+	% plot the given mask #md.mask to locate the field
+	%->
+
+	% save the given model
+	%->
+
+end
+% }}}
+
+% {{{ Parameterization #3
+if perform(org,'Parameterization') 
+
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+
+	% parametrize the model # help parameterize
+	% you will need to fil-up the parameter file (given by the
+  % ParamFile variable)
+	%->
+
+	% save the given model
+	%->
+
+end
+% }}}
+
+% {{{ Extrusion #4
+if perform(org,'Extrusion')
+	
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+
+	% vertically extrude the preceding mesh #help extrude
+	% only 5 layers exponent 1
+	%->
+
+	% plot the 3D geometry #plotdoc
+	%->
+
+	% save the given model
+	%->
+
+end
+% }}}
+
+% {{{ Set the flow computing method #5
+
+if perform(org,'SetFlow')
+
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+
+	% set the approximation for the flow computation #help setflowequation
+	% We will be using the Higher Order Model (HO)
+	%->
+
+	% save the given model
+	%->
+
+end
+% }}}
+
+% {{{ Set Boundary Conditions #6
+if perform(org,'BoundaryCondition')
+
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+
+	% dirichlet boundary condition are known as SPCs
+	% ice frozen to the base, no velocity	#md.stressbalance
+	% SPCs are initialized at NaN one value per vertex
+	%->
+
+	%->
+
+	%->
+
+	% extract the nodenumbers at the base #md.mesh.vertexonbase
+	%->
+
+  % set the sliding to zero on the bed (Vx and Vy)
+	%->
+
+	%->
+
+	% periodic boundaries have to be fixed on the sides
+	% create tabs with the side of the domain
+	% for x
+	% create maxX #help find
+	%->
+
+	% create minX
+	%->
+
+	% for y, max X and minX should be excluded
+	% create maxY
+	%->
+
+	% create minY
+	%->
+
+	% set the node that should be paired together
+	% #md.stressbalance.vertex_pairing
+	%->
+
+	if (ParamFile=='IsmipF.par')
+		% if we are dealing with IsmipF the solution is in masstransport
+		md.masstransport.vertex_pairing=md.stressbalance.vertex_pairing;
+  end
+	% save the given model
+	%->
+
+end
+% }}}
+
+% {{{ Solving #7
+if perform(org,'SolvingISMIPA')
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+
+	% Set cluster #md.cluster
+	% generic parameters #help generic
+	% set only the name and number of process
+	%->
+
+	% Set which control message you want to see #help verbose
+	%->
+
+	% Solve #help solve
+	% we are solving a StressBalanc
+	%->
+
+	% save the given model
+	%->
+
+	% plot the surface velocities #plotdoc
+	%->
+end
+% }}}
+
+% {{{ Solving #8
+if perform(org,'SolvingISMIPF')
+	% load the preceding step #help loadmodel
+	% path is given by the organizer with the name of the given step
+	%->
+
+	% Set cluster #md.cluster
+	% generic parameters #help generic
+	% set only the name and number of process
+	%->
+
+	% Set which control message you want to see #help verbose
+	%->
+
+	% set the transient model to ignore the thermal model
+	% #md.transient 
+	%->
+
+	% define the timestepping scheme
+	% everything here should be provided in years #md.timestepping
+	% give the length of the time_step (4 years)
+	%->
+
+	% give final_time (20*4 years time_steps)
+	%->
+
+	% Solve #help solve
+	% we are solving a TransientSolution
+	%->
+
+	% save the given model
+	%->
+
+	% plot the surface velocities #plotdoc
+	%->
+
+end
+% }}}