%Test of a bump

%initialize model with Stokes model
md=model;
md=mesh(md,'DomainOutline.exp',20000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=extrude(md,8,4);
md.pressure=zeros(md.numberofgrids,1);
md.vx=zeros(md.numberofgrids,1);
md.vy=zeros(md.numberofgrids,1);
md.vz=zeros(md.numberofgrids,1);
md=setelementstype(md,'pattyn','all','stokes','all');

%compute solution
md=solve(md,'analysis_type','diagnostic');

%save model stokes
mds=md;
save modelstokes mds

%now run a control method on macayeal
md=model;
md=mesh(md,'DomainOutline.exp',20000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md.pressure=zeros(md.numberofgrids,1);
md.vx=zeros(md.numberofgrids,1);
md.vy=zeros(md.numberofgrids,1);
md.vz=zeros(md.numberofgrids,1);
md=setelementstype(md,'macayeal','all');

%Initialize the CM with average value of vx and vy
md.vx_obs=DepthAverage(mds,mds.results.diagnostic.vx);
md.vy_obs=DepthAverage(mds,mds.results.diagnostic.vy);
md.vel_obs=DepthAverage(mds,mds.results.diagnostic.vel);
md.drag=30*ones(md.numberofgrids,1);

%Compute the CM
md=solve(md,'analysis_type','control')

%%now run a control method on pattyn
%md=model;
%md=mesh(md,'DomainOutline.exp',20000);
%md=geography(md,'','');
%md=parameterize(md,'Square.par');
%md=extrude(md,8,4);
%md.pressure=zeros(md.numberofgrids,1);
%md.vx=zeros(md.numberofgrids,1);
%md.vy=zeros(md.numberofgrids,1);
%md.vz=zeros(md.numberofgrids,1);
%md=setelementstype(md,'pattyn','all');
%
%%Initialize the CM
%md.vx_obs=mds.results.diagnostic.vx;
%md.vy_obs=mds.results.diagnostic.vy;
%md.vel_obs=mds.results.diagnostic.vel;
%
%%Compute the CM
%md=solve(md,'analysis_type','control')
