source: issm/trunk/test/Validation/EISMINT/MassConservation/runme.m@ 4869

%This test is a test from the EISMINT for Ice shelves 
% Vincent Rommelaere 1996

%The goal is to test the prognostic model
%md=bamg(model,'domain','DomainOutline.exp','hmax',4550);
md=bamg(model,'domain','DomainOutline.exp','hmax',3000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'macayeal','all');

%Evolution of the ice shelf
md.ndt=500;
md.dt=1;
md.artificial_diffusivity=0; %Better result with no artificial diffusivity

%0launch transient solution
%md=solve(md,'analysis_type','transient');
%FOR NOW:
md.cluster=oshostname();
md.np=8;
i=0;
time=0;
md.dummy=struct();
pos=find(md.y>199999.9);
connectivity=full(sparse(md.elements(:),1,1));
while(time<500),
        disp(['step ' num2str(i) '/' num2str(500/md.dt)]);
        i=i+1;
        time=time+md.dt;
        md=solve(md,'analysis_type',PrognosticSolutionEnum);
        thickness   =zeros(md.numberofgrids,1);
        connectivity=sparse(md.results.PrognosticSolution.Thickness.index(:),1,1);
        thickness   =sparse(md.results.PrognosticSolution.Thickness.index(:),1,md.results.PrognosticSolution.Thickness.value(:));
        md.thickness=full(thickness./connectivity);
        md.thickness(pos)=500+100*sin(2*pi*time/200)*ones(size(pos,1),1);
        md.surface=md.bed+md.thickness;
        md.dummy(i).thickness=md.thickness;
        %plotmodel(md,'data',md.thickness,'sectionvalue','CrossLine.exp')
end

%plot results
plotmodel(md,'data',[md.thickness 500+100*sin(2*pi/200*(500-md.y/400))],'sectionvalue','CrossLine.exp')

%Don't forget to add these lines in icetransient2d.m, just before the computation of the thickness to change the thickness on the upper boundary condition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%pos=find(m_p.ys~=0);
%m_p.ys(pos)=500+500/5*sin(2*pi*time/(200*md.yts))*ones(size(pos,1),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%