1 | %Test Name: SquareSheetShelfDiadSSA3dDakotaSamp
|
---|
2 | md=triangle(model(),'../Exp/Square.exp',150000.);
|
---|
3 | md=setmask(md,'../Exp/SquareShelf.exp','');
|
---|
4 | md=parameterize(md,'../Par/SquareSheetShelf.par');
|
---|
5 | md=setflowequation(md,'SSA','all');
|
---|
6 | md.cluster=generic('name',oshostname(),'np',3);
|
---|
7 | md.materials.rho_ice=10^7; %involved in the mass flux, make it easy
|
---|
8 | md.geometry.thickness(:)=1; %make it easy
|
---|
9 | md.geometry.surface=md.geometry.base+md.geometry.thickness;
|
---|
10 |
|
---|
11 | %constrain all velocities to 1 m/yr, in the y-direction
|
---|
12 | md.stressbalance.spcvx(:)=0;
|
---|
13 | md.stressbalance.spcvy(:)=1;
|
---|
14 | md.stressbalance.spcvz(:)=0;
|
---|
15 |
|
---|
16 | %Dakota options
|
---|
17 |
|
---|
18 | %dakota version
|
---|
19 | version=IssmConfig('_DAKOTA_VERSION_'); version=version(1:3); version=str2num(version);
|
---|
20 |
|
---|
21 | %partitioning
|
---|
22 | npart=20;
|
---|
23 | [partition,md]=partitioner(md,'package','chaco','npart',npart,'weighting','on');
|
---|
24 | partition=partition-1;
|
---|
25 | md.qmu.isdakota=1;
|
---|
26 |
|
---|
27 | md.qmu.variables.drag_coefficient=normal_uncertain('descriptor','scaled_FrictionCoefficient',...
|
---|
28 | 'mean',ones(npart,1),...
|
---|
29 | 'stddev',.01*ones(npart,1),...
|
---|
30 | 'partition',partition);
|
---|
31 |
|
---|
32 | md.qmu.responses.MaxVel=response_function('descriptor','MaxVel');
|
---|
33 | md.qmu.responses.MassFlux1=response_function('descriptor','indexed_MassFlux_1');
|
---|
34 | md.qmu.responses.MassFlux2=response_function('descriptor','indexed_MassFlux_2');
|
---|
35 | md.qmu.responses.MassFlux3=response_function('descriptor','indexed_MassFlux_3');
|
---|
36 | md.qmu.responses.MassFlux4=response_function('descriptor','indexed_MassFlux_4');
|
---|
37 | md.qmu.responses.MassFlux5=response_function('descriptor','indexed_MassFlux_5');
|
---|
38 | md.qmu.responses.MassFlux6=response_function('descriptor','indexed_MassFlux_6');
|
---|
39 | md.qmu.responses.MassFlux7=response_function('descriptor','indexed_MassFlux_7');
|
---|
40 |
|
---|
41 | %mass flux profiles
|
---|
42 | md.qmu.mass_flux_profiles={'../Exp/MassFlux1.exp','../Exp/MassFlux2.exp','../Exp/MassFlux3.exp','../Exp/MassFlux4.exp','../Exp/MassFlux5.exp','../Exp/MassFlux6.exp','../Exp/Square.exp'};
|
---|
43 | md.qmu.mass_flux_profile_directory=pwd;
|
---|
44 |
|
---|
45 |
|
---|
46 | %% nond_sampling study
|
---|
47 |
|
---|
48 | md.qmu.method =dakota_method('nond_samp');
|
---|
49 | md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),...
|
---|
50 | 'seed',1234,...
|
---|
51 | 'samples',20,...
|
---|
52 | 'sample_type','lhs');
|
---|
53 |
|
---|
54 | %% a variety of parameters
|
---|
55 | md.qmu.params.interval_type='forward';
|
---|
56 | md.qmu.params.direct=true;
|
---|
57 | md.qmu.params.tabular_graphics_data=true;
|
---|
58 |
|
---|
59 | if version>=6,
|
---|
60 | md.qmu.params.analysis_driver='matlab';
|
---|
61 | md.qmu.params.evaluation_scheduling='master';
|
---|
62 | md.qmu.params.processors_per_evaluation=2;
|
---|
63 | else
|
---|
64 | md.qmu.params.analysis_driver='stressbalance';
|
---|
65 | md.qmu.params.evaluation_concurrency=1;
|
---|
66 | end
|
---|
67 |
|
---|
68 |
|
---|
69 |
|
---|
70 | md.stressbalance.reltol=10^-5; %tighten for qmu analyses
|
---|
71 |
|
---|
72 | md=solve(md,'Stressbalance','overwrite','y');
|
---|
73 |
|
---|
74 | %Fields and tolerances to track changes
|
---|
75 | md.qmu.results=md.results.dakota;
|
---|
76 |
|
---|
77 | %ok, mass flux of 3 profiles should be -3 Gt/yr -3 Gt/yr and the sum, which is -6 Gt/yr
|
---|
78 | %we recover those mass fluxes through the mean of the response.
|
---|
79 | %also, we recover the max velo, which should be 1m/yr.
|
---|
80 | %we put all that data in the montecarlo field, which we will use to test for success.
|
---|
81 | %also, check that the stddev are 0.
|
---|
82 | md.results.dakota.montecarlo=[];
|
---|
83 | for i=1:8,
|
---|
84 | md.results.dakota.montecarlo=[md.results.dakota.montecarlo md.results.dakota.dresp_out(i).mean];
|
---|
85 | end
|
---|
86 | for i=1:8,
|
---|
87 | md.results.dakota.montecarlo=[md.results.dakota.montecarlo md.results.dakota.dresp_out(i).stddev];
|
---|
88 | end
|
---|
89 | field_names ={'montecarlo'};
|
---|
90 | field_tolerances={1e-11};
|
---|
91 | field_values={...
|
---|
92 | md.results.dakota.montecarlo,...
|
---|
93 | };
|
---|