%Test Name: EarthSlr Dakota Sampling glaciers.

%mesh earth:
md=model;
md.cluster=generic('name',oshostname(),'np',5);
md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh

%parameterize solidearth solution:
%solidearth loading:  {{{
md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
md.dsl.global_average_thermosteric_sea_level_change=[1;0];
md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
%antarctica
late=sum(md.mesh.lat(md.mesh.elements),2)/3;
longe=sum(md.mesh.long(md.mesh.elements),2)/3;
pos=find(late < -80);
md.solidearth.surfaceload.icethicknesschange(pos)=-100;
%greenland
pos=find(late>70 & late<80 & longe>-60 & longe<-30);
md.solidearth.surfaceload.icethicknesschange(pos)=-100;

%elastic loading from love numbers:
md.solidearth.lovenumbers=lovenumbers('maxdeg',100);

%}}}
%mask:  {{{
mask=gmtmask(md.mesh.lat,md.mesh.long);
icemask=ones(md.mesh.numberofvertices,1);
pos=find(mask==0);
icemask(pos)=-1;
pos=find(sum(mask(md.mesh.elements),2)<3);
icemask(md.mesh.elements(pos,:))=-1;
md.mask.ice_levelset=icemask;
md.mask.ocean_levelset=-icemask;

%make sure that the elements that have loads are fully grounded:
pos=find(md.solidearth.surfaceload.icethicknesschange);
md.mask.ocean_levelset(md.mesh.elements(pos,:))=1;

%make sure wherever there is an ice load, that the mask is set to ice:
pos=find(md.solidearth.surfaceload.icethicknesschange);
md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
% }}}

md.solidearth.settings.ocean_area_scaling=0;

%Geometry for the bed, arbitrary: 
md.geometry.bed=-ones(md.mesh.numberofvertices,1);

%Materials: 
md.materials=materials('hydro');

%Miscellaneous
md.miscellaneous.name='test2006';

%Solution parameters
md.solidearth.settings.reltol=NaN;
md.solidearth.settings.abstol=1e-3;
md.solidearth.settings.computesealevelchange=1;

% max number of iteration reverted back to 10 (i.e., the original default value)
md.solidearth.settings.maxiter=10;

%eustatic + rigid + elastic + rotation run:
md.solidearth.settings.rigid=1;
md.solidearth.settings.elastic=1;
md.solidearth.settings.rotation=1;

%transient settings: 
md.timestepping.start_time=0;
md.timestepping.final_time=10;
md.timestepping.time_step=1;
md.transient.isslr=1;
md.transient.issmb=0;
md.transient.isgia=1;
md.transient.ismasstransport=0;
md.transient.isstressbalance=0;
md.transient.isthermal=0;
dh=md.solidearth.surfaceload.icethicknesschange;
deltathickness=zeros(md.mesh.numberofelements+1,10);
for i=1:10,
	deltathickness(1:end-1,i)=dh*i;
end
deltathickness(end,:)=0:1:9;
md.solidearth.surfaceload.icethicknesschange=deltathickness;

%hack: 
md.geometry.surface=zeros(md.mesh.numberofvertices,1);
md.geometry.thickness=ones(md.mesh.numberofvertices,1);
md.geometry.base=-ones(md.mesh.numberofvertices,1);
md.geometry.bed=md.geometry.base;


%Uncertainty Quantification
%ice sheets {{{
npart=1; nt=1;
partition=-ones(md.mesh.numberofelements,1);
pos=find(late < -80); partition(pos)=0;
pos=find(late>70 & late<80 & longe>-60 & longe<-30); partition(pos)=0;

%variables: 
qmuvar.surfaceload=normal_uncertain('descriptor','scaled_SurfaceloadIceThicknessChange',...
	'mean',1*ones(npart,nt),...
	'stddev',.1*ones(npart,nt),... %10% standard deviation
	'partition',partition,...
	'transient','on',...
	'nsteps',nt);
%}}}

%correlation:
md.qmu.correlation_matrix=[];

%variables final declaration:
md.qmu.variables=struct();
md.qmu.variables.surfaceload=qmuvar.surfaceload;

locations=[1 5 10 15 20];
%responses  % {{{
md.qmu.responses.sealevel1=response_function('descriptor','Outputdefinition1');
md.qmu.responses.sealevel2=response_function('descriptor','Outputdefinition2');
md.qmu.responses.sealevel3=response_function('descriptor','Outputdefinition3');
md.qmu.responses.sealevel4=response_function('descriptor','Outputdefinition4');
md.qmu.responses.sealevel5=response_function('descriptor','Outputdefinition5');

%output definitions: 
for i=1:length(locations),
	ind=locations(i);
	if i==1,
		md.outputdefinition.definitions={nodalvalue('name','SNode','definitionstring','Outputdefinition1', ...
			'model_string','Sealevel','node',ind)}; 
	else
		md.outputdefinition.definitions{end+1}=nodalvalue('name','SNode','definitionstring',['Outputdefinition' num2str(i)], ...
			'model_string','Sealevel','node',ind); 
	end
end 
% }}}

%algorithm:  % {{{
md.qmu.method     =dakota_method('nond_samp');
md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),...
'seed',1234,...
'samples',10,...
'sample_type','random');
md.qmu.output=1; 
%}}}
%parameters % {{{
md.qmu.params.direct=true;
md.qmu.params.interval_type='forward';
md.qmu.params.analysis_driver='matlab';
md.qmu.params.evaluation_scheduling='master';
md.qmu.params.processors_per_evaluation=2;
md.qmu.params.tabular_graphics_data=true; 
md.qmu.isdakota=1;
md.verbose=verbose(0); md.verbose.qmu=1;
% }}}
%qmu statistics %{{{
md.qmu.statistics.nfiles_per_directory=2;
md.qmu.statistics.ndirectories=5;

md.qmu.statistics.method(1).name='Histogram';
md.qmu.statistics.method(1).fields={'Sealevel','BslrIce'};
md.qmu.statistics.method(1).steps=[1:10];
md.qmu.statistics.method(1).nbins=20;

md.qmu.statistics.method(2).name='MeanVariance';
md.qmu.statistics.method(2).fields={'Sealevel','BslrIce'};
md.qmu.statistics.method(2).steps=[1:10];

md.qmu.statistics.method(3).name='SampleSeries';
md.qmu.statistics.method(3).fields={'Sealevel','BslrIce'};
md.qmu.statistics.method(3).steps=[1:10];
md.qmu.statistics.method(3).indices=locations;
%}}}

%run transient dakota solution: 
mds=solve(md,'Transient');

%run without statistics computations:
md.qmu.statistics.method(1).name='None';
md=solve(md,'Transient');

%compare statistics with our own here:
svalues=mds.results.StatisticsSolution(end).SealevelSamples; %all values at locations.

dvalues=zeros(md.qmu.method.params.samples,length(locations));
for i=1:md.qmu.method.params.samples,
	dvalues(i,:)=md.results.dakota.modelresults{i}.TransientSolution(end).Sealevel(locations);
end

samplesnorm=norm(dvalues-svalues,'fro');

%Fields and tolerances to track changes
field_names={'Samples Norm'};
field_tolerances={1e-13};
field_values={samplesnorm};
