md=triangle(model(),'../Exp/Square.exp',100000.);
md=setmask(md,'all','');
md=parameterize(md,'../Par/SquareShelfConstrained.par');
md=setflowequation(md,'SSA','all');
md.cluster=generic('name',oshostname(),'np',3);

md.autodiff.isautodiff=true;
md.verbose.autodiff=true;
md.autodiff.independents={...
	independent('name','Thickness','type','vertex','nods',md.mesh.numberofvertices);
	%independent('name','Thickness','type','vertex','nods',md.mesh.numberofvertices,'fov_forward_indices',(1:md.mesh.numberofvertices)')
	%independent('name','Thickness','type','vertex','nods',md.mesh.numberofvertices,'fos_forward_index',1)
	};

md.autodiff.dependents={...
	dependent('name','MaxVel','type','scalar','fos_reverse_index',1)...
	%dependent('name','IceVolume','type','scalar')...
	%dependent('name','MassFlux','type','scalar','exp','../Exp/MassFlux1.exp','mesh',md.mesh,'index',1)...
	%dependent('name','MassFlux','type','scalar','exp','../Exp/MassFlux2.exp','mesh',md.mesh,'index',2)...
	};

%md.autodiff.driver='fos_forward';
md.autodiff.driver='fos_reverse';

md=solve(md,TransientSolutionEnum());

%plot
nods=md.mesh.numberofvertices;
jac=md.results.TransientSolution(1).AutodiffJacobian;
plotmodel(md,'data',jac(1:nods),'data',jac(nods+1:2*nods),'data',jac(2*nods+1:3*nods),'data',jac(3*nods+1:4*nods));
