source: issm/oecreview/Archive/14312-15392/ISSM-14853-14854.diff

../trunk-jpl/test/NightlyRun/test2072.m

@@ +1 @@
+% Benchmark experiments (Figure A2c Ivins and James, 1999, Geophys. J. Int.) 
+md=triangle(model(),'../Exp/RoundFrontEISMINT.exp',200000.);
+md=setmask(md,'','');
+md=parameterize(md,'../Par/GiaBenchmarksCD.par');
+
+%% indicate what you want to compute 
+md.gia.output_rates=1;           % just want "w" solution 
+md.gia.cross_section_shape=1;    % for square-edged x-section 
+
+%% define loading history 
+md.timestepping.start_time=1000.3;  % for t \approx 1 kyr 
+md.timestepping.final_time=2500000; % 2,500 kyr
+% In order to run Ivins1999 benchmarks, for now, we do things manually as follows: 
+% 1. numtimes is fixed to be five 
+% 2. final_time given above is hard-coded in "GiaDeflectionCorex.cpp"
+% 3. evaluation time = start_time 
+% 4. define new load times, according to the chosen evaluation time (care that general load history is not altered). This is really important, because of #1 and also becasue of the fact that futre loading times are not recognized.
+% 5. Ice thickness assoiated with evaluation time is actually the same at final_time
+md.geometry.thickness=[...
+        [md.geometry.thickness*0.0; 0.0],...
+        [md.geometry.thickness/2.0; 0.1],...
+        [md.geometry.thickness; 0.2],...
+        [md.geometry.thickness; 0.25],...
+        [md.geometry.thickness; md.timestepping.start_time],...
+        ];
+
+%% solve for GIA deflection 
+md.cluster=generic('name',oshostname(),'np',3);
+md.verbose=verbose('1111111');
+md=solve(md,GiaSolutionEnum());
+
+%Fields and tolerances to track changes
+field_names     ={'GiaW','GiadWdt'};
+field_tolerances={1e-13,1e-13};
+field_values={...
+        (md.results.GiaSolution.GiaW),...
+        (md.results.GiaSolution.GiadWdt),...
+        };