source: issm/trunk-jpl/test/NightlyRun/test257.m@ 27318

%Test Name: SquareShelfSMBarma
md=triangle(model(),'../Exp/Square.exp',80000.);
md=setmask(md,'all','');
md=parameterize(md,'../Par/SquareShelfConstrained.par');
md=setflowequation(md,'SSA','all');
md.cluster=generic('name',oshostname(),'np',3);
md.transient.requested_outputs={'default','IceVolume','SmbMassBalance'};

ymax = max(md.mesh.y);
xmax = max(md.mesh.x);
%Generate basin IDs for 3 basins
idbasin      = zeros(md.mesh.numberofelements,1);
iid1         = find(md.mesh.y>=2/3*ymax);
iid2         = intersect(find(md.mesh.y<2/3*ymax),find(md.mesh.x>=1/3*xmax));
iid3         = intersect(find(md.mesh.y<2/3*ymax),find(md.mesh.x<1/3*xmax));
for ii=1:md.mesh.numberofelements
    for vertex=1:3
        if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1
            idbasin(ii) = 1;
        end
    end
    if idbasin(ii)==0 %no vertex was found in basin 1
        for vertex=1:3
            if any(iid2==md.mesh.elements(ii,vertex)) %one vertex in basin 2
                idbasin(ii) = 2;
            end
        end
    end
    if idbasin(ii)==0 %no vertex was found in basin 1 and 2
        idbasin(ii) = 3;
    end
end

%SMB parameters
numparams                  = 2;
numbreaks                  = 1;
intercept                  = [0.5,1.0;1.0,0.6;,2.0,3.0]; %intercept values of SMB in basins [m ice eq./yr]
trendlin                   = [0.0,0.0;0.01,0.001;-0.01,0]; %trend values of SMB in basins [m ice eq./yr^2]
polynomialparams           = cat(3,intercept,trendlin);
datebreaks                 = [3;3;3];

md.timestepping.start_time = 0;
md.timestepping.time_step  = 1;
md.timestepping.final_time = 7;
md.smb                     = SMBarma();
md.smb.num_basins          = 3; %number of basins
md.smb.basin_id            = idbasin; %prescribe basin ID number to elements
md.smb.num_params          = numparams; %number of parameters in the polynomial
md.smb.num_breaks          = numbreaks; %number of breakpoints
md.smb.polynomialparams    = polynomialparams;
md.smb.datebreaks          = datebreaks;
md.smb.ar_order            = 4;
md.smb.ma_order            = 1;
md.smb.arma_timestep       = 2.0; %timestep of the ARMA model [yr]
md.smb.arlag_coefs         = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]];
md.smb.malag_coefs         = [1.0;0;0.2];
md.smb.lapserates          = [0.01,0.0;0.01,-0.01;0.0,-0.01];
md.smb.elevationbins       = [100;150;100];

%Stochastic forcing
md.stochasticforcing.isstochasticforcing = 1;
md.stochasticforcing.fields              = [{'SMBarma'}];
md.stochasticforcing.covariance          = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]]; %global covariance among- and between-fields
md.stochasticforcing.randomflag          = 0; %fixed random seeds
md.stochasticforcing.stochastictimestep  = 1.0;

md=solve(md,'Transient');

%Fields and tolerances to track changes
field_names     ={'Vx1','Vy1','Vel1','Thickness1','Volume1','SmbMassBalance1','Vx2','Vy2','Vel2','Thickness2','Volume2','SmbMassBalance2','Vx3','Vy3','Vel3','Thickness3','Volume3','SmbMassBalance3'};
field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,...
                                                1e-13,1e-13,1e-13,1e-13,1e-13,1e-13...
                                                1e-12,1e-12,1e-12,1e-12,1e-12,1e-12};
field_values={...
        (md.results.TransientSolution(1).Vx),...
        (md.results.TransientSolution(1).Vy),...
        (md.results.TransientSolution(1).Vel),...
        (md.results.TransientSolution(1).Thickness),...
        (md.results.TransientSolution(1).IceVolume),...
        (md.results.TransientSolution(1).SmbMassBalance),...
        (md.results.TransientSolution(2).Vx),...
        (md.results.TransientSolution(2).Vy),...
        (md.results.TransientSolution(2).Vel),...
        (md.results.TransientSolution(2).Thickness),...
        (md.results.TransientSolution(2).IceVolume),...
        (md.results.TransientSolution(2).SmbMassBalance),...
        (md.results.TransientSolution(7).Vx),...
        (md.results.TransientSolution(7).Vy),...
        (md.results.TransientSolution(7).Vel),...
        (md.results.TransientSolution(7).Thickness),...
        (md.results.TransientSolution(7).IceVolume),...
        (md.results.TransientSolution(7).SmbMassBalance),...
        };